# Kronecker Square Roots

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[]]]},
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;
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[]]], 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]}
]
]

• 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. Jul 17 at 7:46

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, {4, 4}]] // FullSimplify // Abs // Max
(* 0 *)

• This isn't working for me: a = ArrayReshape[Range, {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. Jul 18 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. Jul 18 at 11:50
• Yeah that works if I use N. Jul 18 at 12:25
• @flinty I've modified to normalise the eigenvector appropriately. Jul 18 at 17:32
• Thanks, any idea how to find other complex solutions? Jul 18 at 18:01