4
$\begingroup$

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?

$\endgroup$
18
  • 8
    $\begingroup$ 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}}? $\endgroup$
    – eyorble
    Dec 11, 2022 at 10:32
  • 2
    $\begingroup$ @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. $\endgroup$ Dec 11, 2022 at 15:58
  • 2
    $\begingroup$ @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. $\endgroup$ Dec 11, 2022 at 16:28
  • 1
    $\begingroup$ @userrandrand MatrixPower can compute Sign function on matrices, which is obviously not analytic at zero. $\endgroup$
    – Anixx
    Dec 11, 2022 at 16:38
  • 2
    $\begingroup$ 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 $\endgroup$ Dec 11, 2022 at 17:22

3 Answers 3

9
$\begingroup$

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]

enter image description here

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} *)
$\endgroup$
1
  • 1
    $\begingroup$ (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. $\endgroup$ Dec 11, 2022 at 18:03
7
$\begingroup$

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]]
$\endgroup$
4
  • $\begingroup$ What is this code even supposed to do? $\endgroup$
    – Anixx
    Dec 11, 2022 at 10:43
  • $\begingroup$ 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][[1]] (do not know how to remove intermediate variable) $\endgroup$
    – Anixx
    Dec 11, 2022 at 11:46
  • $\begingroup$ 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. $\endgroup$
    – Anixx
    Dec 11, 2022 at 11:49
  • $\begingroup$ @Anixx For some matrix, the 1/2 power work. MatrixPower[{{2, 3, 0}, {4, 9, 0}, {0, 0, 4}}, 1/2] $\endgroup$
    – cvgmt
    Dec 11, 2022 at 12:16
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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