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?

  • 1
    $\begingroup$ You should post this question to matheamtics.stackexchange.com. $\endgroup$ Feb 8, 2018 at 11:58
  • $\begingroup$ 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. $\endgroup$ Feb 8, 2018 at 13:39
  • $\begingroup$ @SjoerdSmit The root is requested to be symmetric and Cholesky factors are not (by definition). $\endgroup$ Feb 8, 2018 at 16:42
  • $\begingroup$ Ah, yes. I missed the part about C being symmetric. $\endgroup$ Feb 8, 2018 at 17:01

3 Answers 3


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]]



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;
aij = RandomInteger[{-2, 2}, nn (nn + 1)/2];
aa = Statistics`Library`VectorToSymmetricMatrix[Drop[aij, nn], Take[aij, nn], nn];

Mathematica graphics

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

Mathematica graphics


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

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.