# Find the square root of a matrix

Let $A$ be the matrix $$A = \left(\begin{array}{cc} 41 & 12\\ 12 & 34 \end{array}\right).$$

I want to decompose it into the form of $B^2$.

I tried diagonalization , but can not move one step further.

Any thought on this? Thanks a lot!

ONE STEP FURTHER:

How to find a upper triangular $U$ such that $A = U^T U$?

-
The matrix is symmetric, so it is certainly diagonalizable. Trace and determinant are both positive, so both eigenvalues are positive. So if you can diagonalize, the diagonal form will have a square root, $QAQ^{-1} = D = P^2$, where $Q$ is the change-of-basis matrix. That means that $A = Q^{-1}P^2Q = (Q^{-1}PQ)^2$, so you can let $B=Q^{-1}PQ$. So your idea works; where did you get stuck? – Arturo Magidin Aug 24 '11 at 3:40
@Arturo Magidin I did not figure out how to use the diagonalized form. Your answer is brilliant! Thanks! – BVFanZ Aug 24 '11 at 3:46
en.wikipedia.org/wiki/… – Fixee Aug 24 '11 at 4:02
Re: edit: Cholesky decomposition can be done with an approach similar to Gerry's: write out the expression for the product of a lower triangular matrix with its transpose, equate to your original matrix, and solve the resulting set of equations... – J. M. Aug 24 '11 at 4:03
+1. Nice question! – Pierre-Yves Gaillard Aug 24 '11 at 7:16

This is an expansion of Arturo's comment.

The matrix has eigenvalues $50,25$, and eigenvectors $(4,3),(-3,4)$, so it eigendecomposes to $$A=\begin{pmatrix}4 & -3 \\ 3 & 4\end{pmatrix} \begin{pmatrix}50 & 0 \\ 0 & 25\end{pmatrix} \begin{pmatrix}4 & -3 \\ 3 & 4\end{pmatrix}^{-1}.$$

This is of the form $A=Q\Lambda Q^{-1}$. If this is $B^2$, then there will be a $B$ of the form $Q\Lambda^{1/2} Q^{-1}$ (square this to check this is formally true). A square root of a diagonal matrix is just the square roots of the diagonal entries, so we have

$$B=\begin{pmatrix}4 & -3 \\ 3 & 4\end{pmatrix} \begin{pmatrix}\sqrt{50} & 0 \\ 0 & \sqrt{25}\end{pmatrix} \begin{pmatrix}4 & -3 \\ 3 & 4\end{pmatrix}^{-1}$$

$$=\frac{1}{5}\begin{pmatrix}9+16\sqrt{2} & -12+12\sqrt{2} \\ -12+12\sqrt{2} & 16+9\sqrt{2}\end{pmatrix}.$$

Here we used $\sqrt{50}=5\sqrt{2},\sqrt{25}=5$, and a quick formula for the inverse of a $2\times 2$ matrix:

$$\begin{pmatrix}a&b\\c&d\end{pmatrix}^{-1}=\frac{1}{ad-bc}\begin{pmatrix}d&-b\\-c&a\end{pmatrix}.$$

Keep in mind that matrix square roots are not unique (even up to sign), but this particular method is guaranteed to produce one example of a real matrix square root whenever $A$ has all positive eigenvalues.

Finding an upper triangular $U$ such that $A=U^TU$ is even more straightforward:

$$A=\begin{pmatrix} a&0 \\ b&c \end{pmatrix} \cdot \begin{pmatrix} a&b \\ 0&c \end{pmatrix}$$

This is $a^2=41$ hence $a=\sqrt{41}$, $ab=12$ hence $b=\frac{12}{41}\sqrt{41}$, and $b^2+c^2=34$ hence $c=25\sqrt{\frac{2}{41}}$.

In other words,

$$U=\sqrt{41}\begin{pmatrix}1&\frac{12}{41}\\0&\frac{25}{41}\sqrt{2}\end{pmatrix}.$$

-
Could you provide a source for this please? That would be much appreciated! – Pubbie Nov 10 '14 at 23:30
@Pubbie What do you want a source for? This answer seems self-contained to me. – anon Nov 11 '14 at 2:00
Writing a paper concerning this. – Pubbie Nov 11 '14 at 4:13
@Pubbie In any paper I'd write I wouldn't bother with a citation for a fact like this. – anon Nov 11 '14 at 15:58
I found it in the book $\textit{Functions of Matrices}:$ Theory and Computation. – Pubbie Nov 12 '14 at 16:21

For the first part of your question, here is a solution that only works for 2-by-2 matrices, but it has the merit that no eigenvalue is needed.

Recall that in the two-dimensional case, there is a magic equation that is useful in many situations. It is $X^2-({\rm tr}X)X+(\det X)I=0$, which arises from the characteristic polynomial of a $2\times2$ matrix $X$. Now, if $X^2=A$, we have $\det X=\pm\sqrt{\det A}=r$ (say). We take the positive value for $r$. Hence $$(\ast):\quad ({\rm tr}X)X=X^2+rI=A+rI$$ and $({\rm tr}X)^2 = {\rm tr}\left(({\rm tr}X)X\right) = {\rm tr}(A+rI) = {\rm tr}A + 2r$. Thus, from $(\ast)$ we obtain $$X = \frac{1}{\sqrt{{\rm tr}A + 2r}}(A+rI)\quad {\rm where}\quad r=\sqrt{\det A}.$$ This method works for all 2-by-2 matrices $A$ when $\det A\ge0$ and ${\rm tr}A + 2\sqrt{\det A}>0$. In particular, it works for positive definite $A$.

For the second part of your question, as the others have pointed out, the decomposition you ask for is a Cholesky decomposition.

-
Dear @user1551: It seems to me your condition ${\rm tr}A + 2\sqrt{\det A}>0$ just means that the eigenvalues are positive. Is this correct? I also think that your formula is exactly the same as Didier's. [Mine is slightly different (but equivalent) because I haven't been smart enough to cover the two cases by a single formula.] – Pierre-Yves Gaillard Aug 24 '11 at 12:44
Yes, the eigenvalues have to be non-negative and at least one of them must be positive, and our formulae are equivalent. I think you are being modest when you said you were not smart enough. Since mathematicians in this forum tend to analyze problems (and generalize the results) from higher perspectives, it is not surprising that you guys do not take a low road as I did. – user1551 Aug 24 '11 at 15:24
I might be modest, ... but less than you! I mention you and Didier in an edit to my answer. I think my road is lower than yours: I try to express things in the high school language of secant and tangent lines. – Pierre-Yves Gaillard Aug 24 '11 at 16:00
No, my road is certainly lower than yours, because my method cannot be generalized to higher dimensions. – user1551 Aug 24 '11 at 16:34
Your method looks very original to me. It's very natural to use Cayley-Hamilton, but I guess most people would have used it for $A$, not for $X$. That is, put $X:=xA+yI$, write $X^2=A$, use CH in the form $A^2=tA-dI$ (with $t=$ trace, $d=$ det) to get rid of the $A^2$ term. You get 2 equations for $x$ and $y$. This can be generalized. – Pierre-Yves Gaillard Aug 24 '11 at 17:47

Still another explicit formula: for every nonnegative real number $\alpha$, $$A^\alpha=\frac{(u^\alpha-v^\alpha)A+(uv^\alpha-vu^\alpha)I}{u-v},$$ where $u$ and $v$ are the two roots of the polynomial $\chi_A(x)=\det(xI-A)$. When $\alpha=1/2$, $$\sqrt{A}=\frac{A+\sqrt{uv}I}{\sqrt{u}+\sqrt{v}}.$$ Note that the coefficients of this formula can be computed directly from the matrix $A$ since $t=\sqrt{uv}$ is simply $t=\sqrt{\det A}$ and $s=\sqrt{u}+\sqrt{v}$ is such that $s^2=u+v+2t$ hence $$s=\sqrt{\text{tr}(A)+2\sqrt{\det(A)}}.$$ In the present case, one can also compute $\chi_A(x)=(x-41)(x-34)-12^2=(x-25)(x-50)$ and use $\sqrt{u}=5\sqrt2$ and $\sqrt{v}=5$.

-
The formulae of yours, mine and @Pierre-Yves Gaillard are actually equivalent, but how they are expressed and derived are different. – user1551 Aug 24 '11 at 10:32
@user1551, yes. – Did Aug 24 '11 at 11:01
@didier Does this identity have a name or can you provide a reference? – Steve Apr 12 '12 at 2:23

Let $\lambda$ and $\mu$ be the eigenvalues of your $2$ by $2$ real matrix $A$. (We may have $\lambda=\mu$.) Assume that $\lambda$ and $\mu$ are positive.

If $\lambda\not=\mu$, write the equation of the secant line to the curve $y=\sqrt x$ through the points $(\lambda,\sqrt\lambda)$ and $(\mu,\sqrt\mu)$: $$y=\sqrt\lambda\ \ \frac{x-\mu}{\lambda-\mu}+\sqrt\mu\ \ \frac{x-\lambda}{\mu-\lambda}\quad.$$ The matrix you want is $$\sqrt\lambda\ \ \frac{A-\mu I}{\lambda-\mu}+\sqrt\mu\ \ \frac{A-\lambda I}{\mu-\lambda}\quad,$$ where $I$ is the identity matrix.

If $\lambda=\mu$, write the equation of the tangent line to the curve $y=\sqrt x$ through the point $(\lambda,\sqrt\lambda)$: $$y=\sqrt\lambda+\frac{x-\lambda}{2\sqrt\lambda}\quad.$$ The matrix you want is $$\sqrt\lambda\ I+\frac{A-\lambda I}{2\sqrt\lambda}\quad.$$

Do you see why?

Do you see how to generalize this to $n$ by $n$ matrices?

EDIT 4. This is to just explain why this secant/tangent stuff comes into the picture. Assume to simplify that the eigenvalues $\lambda$ and $\mu$ of your two by two real matrix $A$ are real and distinct. Let $f\in\mathbb R[X]$ be a polynomial, and $s$ the unique polynomial of degree $\le1$ which agrees with $f$ at $\lambda$ and $\mu$. [Graphically, this is a secant line.] Then the characteristic polynomial $$\chi=(X-\lambda)(X-\mu)$$ will divide $f-s$. As $\chi(A)=0$ by the Cayley-Hamilton Theorem, we have $f(A)=s(A)$. But the expression $s(A)$ makes sense whenever $f$ is a (real-valued) function defined at $\lambda$ and $\mu$. Moreover, the map $f\mapsto f(A)$ is compatible with addition and multiplication.

EDIT 1. As noticed by Didier Piau and user1551, there is a cute formula for the "generalized secant line" to the curve $y=\sqrt x$, by which I mean: the secant line if the points are distinct, the tangent line if they coincide. Supposing $\lambda\not=\mu$, the equation of the secant line is $$y=\frac{\sqrt\lambda-\sqrt\mu}{\lambda-\mu}\ \ x+ \frac{\mu\sqrt\lambda-\lambda\sqrt\mu}{\lambda-\mu}= \frac{x+\sqrt{\lambda\mu}}{\sqrt\lambda+\sqrt\mu}\quad,$$ and the miracle is that the last expression makes sense even if $\lambda=\mu$.

EDIT 2. Note that there are other solutions when $\lambda\not=\mu$. Putting $$E:=\frac{A-\lambda I}{\mu-\lambda}\quad,\quad F:=\frac{A-\mu I}{\lambda-\mu}\quad,$$ we get $$E^2=E,\ F^2=F,\ EF=FE=0,\ I=E+F,\ A=\mu E+\lambda F,$$ and thus $$(\pm\sqrt\mu\ E\pm\sqrt\lambda\ F)^2=A$$ for the four choices of signs. [The plus plus choice corresponds to the previous formula.]

EDIT 3. Here is a generalization.

Let $T$ be an $n$ by $n$ complex matrix, and $$p(X)=(X-\lambda_1)^{m(1)}\cdots(X-\lambda_k)^{m(k)}$$ its minimal polynomial (the $\lambda_i$ being distinct and the $m(i)$ positive). Let $A$ be the algebra of those functions $f(z)$ which are holomorphic in a neighborhood of the spectrum $\{ \lambda_1,\dots,\lambda_k \}$ of $T$.

There is a unique $\mathbb C[X]$-algebra morphism from $A$ to $\mathbb C[T]=\mathbb C[X]/(p(X))$. Denote this morphism by $f(z)\mapsto f(T)$. If $f(z)$ is in $A$, then the unique representative of $f(T)$ in $\mathbb C[X]$ of degree less than $\deg p(X)$ is $$\sum_{i=1}^k\ \ \underset{X=\lambda_i}\heartsuit\left( \Big(\ \underset{z=\lambda_i}\heartsuit f(z) \Big)\ \ \frac{(X-\lambda_i)^{m(i)}}{p(X)}\ \right)\ \frac{p(X)}{(X-\lambda_i)^{m(i)}}$$

with $$\underset{u=\lambda_i}\heartsuit\varphi(u):=\sum_{j=0}^{m(i)-1}\frac{\varphi^{(j)}(\lambda_i)}{j!}\ (X-\lambda_i)^j.$$

Moreover, the $\lambda_i$-generalized eigenspace of $T$ is contained in the $f(\lambda_i)$-generalized eigenspace of $f(T)$.

All this follows from the Chinese Remainder Theorem, which says $$\frac{\mathbb C[X]}{(p(X))}=\prod_{i=1}^k\ \ \frac{\mathbb C[X]}{(X-\lambda_i)^{m(i)}}\quad,$$ and from the Taylor Formula.

[There is an Edit 4 above.]

-

Write $$\pmatrix{a&b\cr c&d\cr}^2=\pmatrix{41&12\cr12&34\cr}$$ multiply out the left side, set corresponding entries equal to get four equations in the four unknowns $a,b,c,d$, then see if you can work your way through the algebra to a solution.

-
This method gets ugly in a hurry. Once needs to solve $a^2 + b^2 = 41$, $b(a+d) = 12$, and $b^2 + d^2 = 34$ simultaneously. As far as I can tell, there seems to be no clean algebraic solution. – JavaMan Aug 24 '11 at 5:50
I can work it out on a single sheet of paper (but I have very small handwriting...). – Gerry Myerson Aug 24 '11 at 6:23
I won't ask for all the details, but what is your general process? – JavaMan Aug 24 '11 at 6:33
Subtract third equation from first to get a difference of two squares, and express $a+d$ and $a-d$ in terms of $b$ using the expression for $a+d$ you get from the second equation. This gives you $a$ in terms of $b$ and you can substitute in the first equation to get an equation you can solve for b. It is much easier if you are allowed to guess that the answer can be expressed in integers. – Mark Bennet Aug 24 '11 at 8:11
@DJC, $a^2-d^2=7$, $b/(a-d)=12/7$, $12a=7b+12d$, $b(7b+24d)=144$, $24^2b^4+(24bd)^2=(34)(24)^2b^2$, $24^2b^4+(144-7b^2)^2=(34)(24)^2b^2$ gives a biquadratic for $b$. – Gerry Myerson Aug 25 '11 at 7:01

The second is even simpler than the first question. All you need to to look at the assumed solution, setting unknowns:
$\qquad \qquad \small A = U^t \cdot U = \begin{pmatrix} a&0 \\ b&c \end{pmatrix} \cdot \begin{pmatrix} a&b \\ 0&c \end{pmatrix} \qquad$ and $\small a \cdot a = 41$. Then you can proceed; in Pari/GP it needs something like 5 lines of code...

-