7

I'm interested in computing the 'Kronecker Square Root' $\mathbf{A}_{n\times n}$ of a matrix $\mathbf{M}_{n^2\times n^2}=\mathbf{A}\otimes\mathbf{A}$. In case it doesn't exist, I'd like to find the best approximation in some sense.

This question here Nearest Kronecker Product has a good answer that deals with finding two potentially different $\mathbf{A},\mathbf{B}$ such that $\mathbf{M}\approx\mathbf{A}\otimes\mathbf{B}$ but not the square root case.

So far I've just been using NMinimize but is there a better, more efficient way using SingularValueDecomposition?

KrSqrt[mtx_] := With[{dim = Sqrt[Length[mtx[[1]]]]},
  Module[{varsA = Array[a, {dim, dim}]},
   varsA /. Last@NMinimize[

  (** Norm[KroneckerProduct[varsA,varsA]-mtx,"Frobenius"], **)

  (* This Total below is faster than the Fr norm above somehow *)
  Total[(KroneckerProduct[varsA, varsA] - mtx)^2, 2],
  Flatten@varsA]



]]


SeedRandom[1];
m = RandomInteger[20, {4, 4}];
test = KroneckerProduct[m, m];
Round[KrSqrt[test], 10^-6] == m 
(* RETURNS: True *)

Update: You can get slightly better (in "Frobenius" norm sense) approximate square roots on average if $\mathbf{A}$ is allowed to be complex. I'd also like a complex $\mathbf{A}$ option. Here's what I've done based on the same approach as above for Real $\mathbf{A}$:

KrSqrtComplex[mtx_] := 
 Module[{d = Sqrt[Length[mtx[[1]]]], cost, mR, mI},
  cost[m_?(MatrixQ[#, NumericQ] &)] :=
   Total[Abs[KroneckerProduct[m, m] - mtx]^2, 2];
  mR + I mI /. Last@NMinimize[
     cost[mR + I mI], {mR ∈ Matrices[{d, d}, Reals], 
      mI ∈ Matrices[{d, d}, Reals]}
     ]
  ]
flinty
  • 25,147
  • 2
  • 20
  • 86
  • I think that if you dig into what is done there, you will find that you have a symmetric matrix and the solution is given by its principal eigenvector. (Just a guess). This is dependent on M having suitable symmetry, of course. – mikado Jul 17 '21 at 07:46

1 Answers1

2

Here is an efficient implementation that works for some simple test cases.

It is based on Nearest Kronecker Product, but I believe that the intermediate matrix used there is guaranteed to be be symmetric, if the input matrix is a Kronecker square.

kronroot[m_?(MatrixQ[#, NumericQ] &)] := 
 Module[{m1, m2, v1, factor}, 
  m1 = Flatten[Partition[m, Sqrt[Dimensions[m]]], {{2, 1}, {4, 3}}];

(If m has an exact root,m1 will be symmetric.Enforce this)
  m1 = (m1 + Transpose[m1])/2;
  v1 = Normalize[First[Eigenvectors[m1, 1]]];
  factor = Sqrt[v1 . m1 . v1];
  Transpose[Partition[factor*v1, Sqrt[Length[v1]]]]]

To test this

test[m_] := Module[{mk, root},
  mk = KroneckerProduct[m, m];
  root = kronroot[mk];
  mk - KroneckerProduct[root, root]]

test[RandomReal[{-1, 1}, {2, 2}]] // Abs // Max
(* 6.9388910^-17 )


test[ArrayReshape[Range[16], {4, 4}]] // FullSimplify // Abs // Max
(* 0 *)
mikado
  • 16,741
  • 2
  • 20
  • 54
  • This isn't working for me: a = ArrayReshape[Range[16], {4, 4}]; m = KroneckerProduct[a, a]; x = kronroot[m]; KroneckerProduct[x, x] == m - all entries in x are too big by a constant factor of Total[a^2, 2] so I think you need to divide that out. – flinty Jul 18 '21 at 11:40
  • @flinty does it work if you use machine precision numbers e.g. apply N? Could be that I'm assuming the eigenvectors are normalised. – mikado Jul 18 '21 at 11:50
  • Yeah that works if I use N. – flinty Jul 18 '21 at 12:25
  • @flinty I've modified to normalise the eigenvector appropriately. – mikado Jul 18 '21 at 17:32
  • Thanks, any idea how to find other complex solutions? – flinty Jul 18 '21 at 18:01
  • I think that if you might have complex solutions, then the intermediate matrix should be Hermitian. I think you would just need to replace Transpose with ComplexTranspose but it might be more complicated than that. Perhaps you should be considering solutions of the form KroneckerProduct[Conjugate[m],m]? Of course, I don't know why you are trying to solve this equation, so I can't comment on which is more likely to be appropriate. – mikado Jul 18 '21 at 18:05