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?


migrated from mathematica.stackexchange.com Feb 8 '18 at 18:32

This question came from our site for users of Wolfram Mathematica.

  • $\begingroup$ You should post this question to matheamtics.stackexchange.com. $\endgroup$ – José Antonio Díaz Navas Feb 8 '18 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$ – Sjoerd Smit Feb 8 '18 at 13:39
  • $\begingroup$ @SjoerdSmit The root is requested to be symmetric and Cholesky factors are not (by definition). $\endgroup$ – Henrik Schumacher Feb 8 '18 at 16:42
  • $\begingroup$ Ah, yes. I missed the part about C being symmetric. $\endgroup$ – Sjoerd Smit Feb 8 '18 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]]



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.