7
$\begingroup$

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]}
     ]
  ]
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – mikado
    Jul 17 at 7:46
1
$\begingroup$

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 *)
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – flinty
    Jul 18 at 11:40
  • $\begingroup$ @flinty does it work if you use machine precision numbers e.g. apply N? Could be that I'm assuming the eigenvectors are normalised. $\endgroup$
    – mikado
    Jul 18 at 11:50
  • $\begingroup$ Yeah that works if I use N. $\endgroup$
    – flinty
    Jul 18 at 12:25
  • $\begingroup$ @flinty I've modified to normalise the eigenvector appropriately. $\endgroup$
    – mikado
    Jul 18 at 17:32
  • $\begingroup$ Thanks, any idea how to find other complex solutions? $\endgroup$
    – flinty
    Jul 18 at 18:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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