Kronecker Square Roots

Question

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]}
     ]
  ]

Answer 1

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.93889*10^-17 *)

test[ArrayReshape[Range[16], {4, 4}]] // FullSimplify // Abs // Max
(* 0 *)