# Square root of a matrix (or a 2nd order tensor)

Suppose I have two symmetrical matrices $A$ and $C$. Suppose they have the following relation $$A = C^{T} \cdot C = C^{2}$$ and I would like to solve this above equation for $C$ matrix.

Is it correct to write the solution in this form? $$C = A^{1/2}$$

If yes how to define square root of a matrix?

• You should post this question to matheamtics.stackexchange.com. Commented Feb 8, 2018 at 11:58
• You may want to look into trying to get the CholeskyDecomposition of A (assuming it exists), because that would give you a unique matrix C. Commented Feb 8, 2018 at 13:39
• @SjoerdSmit The root is requested to be symmetric and Cholesky factors are not (by definition). Commented Feb 8, 2018 at 16:42
• Ah, yes. I missed the part about C being symmetric. Commented Feb 8, 2018 at 17:01

Nope, it's not correct. If A is positive semi-definite (as it is if it can be written as A = Transpose[B].B), you can do the following:

A = RandomReal[{-1, 1}, {3, 3}];
A = Transpose[A].A;
{λ, U} = Eigensystem[A];
B = Transpose[U].(Sqrt[λ] U);
Max[Abs[B.B - A]]
Max[Abs[Transpose[B].B - A]]


1.55431*10^-15

1.55431*10^-15

Note that B^2 squares each entry and is different from MatrixPower[B,2] (equal B.B).

You can get the "principal" square root using MatrixPower: Using Michael's example:

MatrixPower[{{0,1},{1,1}}, 1/2] //Simplify //TeXForm


$\left( \begin{array}{cc} \frac{\left(-1+\sqrt{5}\right) \sqrt{1+\sqrt{5}}+i \sqrt{-1+\sqrt{5}} \left(1+\sqrt{5}\right)}{2 \sqrt{10}} & \frac{-i \sqrt{-1+\sqrt{5}}+\sqrt{1+\sqrt{5}}}{\sqrt{10}} \\ \frac{-i \sqrt{-1+\sqrt{5}}+\sqrt{1+\sqrt{5}}}{\sqrt{10}} & \frac{i \left(-1+\sqrt{5}\right)^{3/2}+\left(1+\sqrt{5}\right)^{3/2}}{2 \sqrt{10}} \\ \end{array} \right)$

In general, one expects $2^n$ solutions to $C^2 = A$ if $A$ and $C$ are $n \times n$ matrices. So this quickly gets out of hand as $n$ increases. Here is code for find them:

nn = 2;
SeedRandom[4];
aij = RandomInteger[{-2, 2}, nn (nn + 1)/2];
aa = StatisticsLibraryVectorToSymmetricMatrix[Drop[aij, nn], Take[aij, nn], nn];
MatrixForm[aa]


vars = Array[x, nn (nn + 1)/2];
cc = StatisticsLibraryVectorToSymmetricMatrix[Drop[vars, nn], Take[vars, nn], nn];
sols = cc /. Solve[cc.cc == aa, vars];
Map[MatrixForm, Simplify@sols]


Check:

MatrixForm /@ Simplify[#.# & /@ sols]