# Mathematica cannot find square roots of some matrices?

We know that

$$\left( \begin{array}{cccc} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ \end{array} \right)^2=\left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array} \right)$$

Why then Mathematica cannot find the square root of the later matrix? How to make it able to do it?

• Likely because there are infinitely many matrices that, when squared, result in the latter matrix. How is Mathematica supposed to know you want that one instead of, e.g. {{0,0,0,1/17},{17,0,13,0},{0,0,0,0},{0,0,17,0}}? Dec 11, 2022 at 10:32
• @user293787 MatrixPower[...,1/2] uses power series and in the possible issues sections for MatrixFunction one sees that it will not work if the scalar function is not analytical at the eigenvalues. For a nilpotent matrix all of the eigenvalues are 0 and Sqrt is not analytic at 0. Dec 11, 2022 at 15:58
• @user293787 sorry I retract what I said, actually the documentation is not clear on this as it can still compute MatrixPower[...,1/2] for a diagonal matrix with a 0 eigenvalue. I suppose that applies if the matrix is not diagonalizable and it needs to compute the power series. In any case power series will not work because one would have to apply Sqrt[1+(#-1)]& and for a Nilpotent matrix N, N-1 does not have a matrix norm smaller than 1. Dec 11, 2022 at 16:28
• @userrandrand MatrixPower can compute Sign function on matrices, which is obviously not analytic at zero. Dec 11, 2022 at 16:38
• Ohh quite nice so it respects the algebra for complex numbers on a matrix representation. In any case I think Mathematica is just applying the formula in section 7 of this wikipedia page en.wikipedia.org/wiki/Jordan_normal_form#Matrix_functions Dec 11, 2022 at 17:22

A brute force approach is possible:

A = Array[a, {4, 4}];

rhs = DiagonalMatrix[{1, 1}, 2];

eqns = A . A == rhs // LogicalExpand;

solns = Solve[Reduce[eqns]];

Length[solns]
(* 3 *)

MatrixForm /@ Simplify[A /. solns] Verify that these do indeed solve the equation

A . A == rhs /. solns // Simplify
(* {True, True, True} *)


The first solution is a function of 4 parameters. Clearly, the second solution is simply the special case of the first with a[4,4]->0. I wonder if the third is a special limiting case of the first, but I haven't tried too hard to find it.

Some clues are available as to why Mathematica fails with more direct approaches. For example

MatrixFunction[Sqrt, rhs, Method -> "Jordan"]
(* MatrixFunction::fnanc: The function Sqrt[#1]& is not analytic at 0. *)
(* MatrixFunction[Sqrt, {{0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0,
0}, {0, 0, 0, 0}}, Method -> "Jordan"] *)


In fact, MatrixFunction documentation states

MatrixFunction does not return a result when the scalar function or any of its initial derivatives are not defined at matrix eigenvalues

In this case, we have

Eigenvalues[rhs]
(* {0, 0, 0, 0} *)

• (sorry for so many edits). I think Mathematica just uses the formula on this wikipedia page en.wikipedia.org/wiki/Jordan_normal_form#Matrix_functions. In particular, that would explain why Mathematica can find a square root of {{1,0},{0,0}} even though it has a zero eigenvalue. Dec 11, 2022 at 18:03

Edit

• Since there are many 4*4 matrix satisfies such matrix equation m^2==m1,we use FindInstance to find one of such m.
m = Array[a, {4, 4}];
m1 = SparseArray[Band[{1, 3}] -> 1, {4, 4}];
m1//MatrixForm
sol = FindInstance[MatrixPower[m, 2] == m1, Flatten@m, Integers]
m /. sol[]


{{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}}

• If we remove the Integers in above code, we get

{{0, 0, -1, -1}, {-1, 0, 1, 1}, {0, 0, 0, 0}, {0, 0, -1, 0}}

It gave another example which satisfies the original matrix equation.

Original

FindInstance[
MatrixPower[Array[a, {4, 4}], 2] ==
SparseArray[Band[{1, 3}] -> 1, {4, 4}],
Flatten@Array[a, {4, 4}], Integers]
Array[a, {4, 4}] /. %[]

• What is this code even supposed to do? Dec 11, 2022 at 10:43
• OK, modified the code to this: MSqrt[A_] := Array[a, {Length[A], Length[A]}] /. FindInstance[ MatrixPower[Array[a, {Length[A], Length[A]}] m, 2] == A, Flatten@m, Integers][] (do not know how to remove intermediate variable) Dec 11, 2022 at 11:46
• With Integers the code fails on matrices which have no integer roots, and without it it gives undesirable answer in the example from the OP. Dec 11, 2022 at 11:49
• @Anixx For some matrix, the 1/2 power work. MatrixPower[{{2, 3, 0}, {4, 9, 0}, {0, 0, 4}}, 1/2] Dec 11, 2022 at 12:16

Here's an optimization approach where I minimize the square error between $$A^2$$ and $$B$$, regularized with the 1-norm to encourage sparse solutions.:

$$\underset{A_{ij}}{\mathrm{argmin}}\left(\left\Vert \mathrm{vec}(A^2-B)\right\Vert_2+\lambda\cdot \left\Vert \mathrm{vec}(A)\right\Vert_1\right)$$

objective[X_?(MatrixQ[#, NumericQ]&), Y_?(MatrixQ[#, NumericQ]&)] := Norm[Flatten[X.X-Y], 2]

regularizer[X_?(MatrixQ[#, NumericQ]&)] := Norm[Flatten@X,1];

b = {{0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}};
A = Array[a, {4, 4}];
result = With[{λ = 0.1},
NMinimize[objective[A, b] + λ*regularizer[A], Flatten[A],
Method -> "DifferentialEvolution"]
] // Last;

tolerance = 1*^-7;
(a1 = Chop[A /. result, tolerance]) // MatrixForm
Chop[b1 = a1 . a1, tolerance] // MatrixForm


$$\left( \begin{array}{cccc} 0 & 0 & 0 & -1.42473 \\ -0.701888 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & -0.701888 & 0 \\ \end{array} \right)$$