Mathematica cannot find square roots of some matrices?

Question

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?

Answer 1

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} *)

Answer 2

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

{{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}] /. %[[1]]

Answer 3

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) $$