matrix - Nearest Kronecker Product

Some people (see The ubiquitous Kronecker product by Van Loan) have worked on finding two matrices $\mathbf A$,$\mathbf B$ of specified size whose tensor product $\mathbf A\otimes\mathbf B$ is closest (in a norm) to a given (larger) matrix $\mathbf C$. That is, find $\mathbf A$,$\mathbf B$ which minimizes $\|\mathbf C−\mathbf A\otimes\mathbf B\|$. The algorithm is based on the SVD. There is a MATLAB implementation here. It would be nice to see this algorithm implemented in Mathematica (using SVD, Riffle...)

Has anyone done this?

Possible Extensions are:

• Factorisation - If the error is zero then the algorithm factorises the original matrix.


• Find the nearest Kronecker product over all possible smaller matrices.


• The algorithm uses the Frobenius norm - could other norms be used?

Answer

The Pitsianis-Van Loan algorithm turns out to be surprisingly easy to implement in Mathematica:

nearestKroneckerProductSum[mat_?MatrixQ, dim1_?VectorQ, dim2_?VectorQ,
                           k_Integer?Positive, opts___] /; 
      TrueQ[Dimensions[mat] == dim1 dim2] := Module[{tmp},
      Check[tmp =
            SingularValueDecomposition[Flatten[Partition[mat, dim2], {{2, 1}, {4, 3}}],
                                       k, opts], 
            Return[$Failed], SingularValueDecomposition::take];
      MapThread[#1 MapThread[Composition[Transpose, Partition],
                             {#2, {dim1, dim2}[[All, 1]]}] &,
                {Sqrt[Diagonal[tmp[[2]]]], Transpose[Transpose /@ Delete[tmp, 2]]}]]


nearestKroneckerProduct[mat_?MatrixQ, dim1_?VectorQ, dim2_?VectorQ, opts___] /;
      TrueQ[Dimensions[mat] == dim1 dim2] := 
      First[nearestKroneckerProductSum[mat, dim1, dim2, 1, opts]]

To verify the routine, let's use a manifest Kronecker product as an example:

dim1 = {4, 3}; dim2 = {6, 5};
BlockRandom[SeedRandom[42];
            m1 = RandomReal[{-1, 1}, dim1]; m2 = RandomReal[{-1, 1}, dim2]];
tst = KroneckerProduct[m1, m2];

{t1, t2} = nearestKroneckerProduct[tst, dim1, dim2];
Norm[KroneckerProduct[t1, t2] - tst, "Frobenius"] // Chop
   0

Note, however, that t1 and t2 are not the same as m1 and m2, since the factorization is only unique up to a constant factor.

Now, for an actual random matrix:

BlockRandom[SeedRandom[42]; m3 = RandomReal[{-1, 1}, dim1 dim2]];
{bm, cm} = nearestKroneckerProduct[m3, dim1, dim2];
Norm[KroneckerProduct[bm, cm] - m3, "Frobenius"]
   9.853962754011258

---

As can be ascertained from the implementation given above, the Pitsianis-Van Loan algorithm can in fact solve a more general problem; namely, to express a matrix as a sum of Kronecker products. The underlying algorithm relies on constructing a certain matrix through shuffling the original, and then proceeding to derive a low-rank approximation via SVD, which can then be reshuffled into the Kronecker factors of the terms.

To demonstrate, here is the algorithm applied to m3, and trying to express it as a sum of six Kronecker products:

k = 6;
tx = nearestKroneckerProductSum[m3, dim1, dim2, k];
Norm[m3 - Total[KroneckerProduct @@@ tx], "Frobenius"]
   5.6710996859350775