Nearest Kronecker Product

Question

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?

BTW does anyone know about factoring exactly? What about different norms? — pdmclean, Aug 16 '15 at 02:31
"factoring exactly?" - Pitsianis-Van Loan already covers that. See my first example. — J. M.'s missing motivation, Aug 16 '15 at 20:43

score 27 · Accepted Answer · edited Apr 13 '17 at 12:56

27

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

edited Apr 13 '17 at 12:56

Community

1

answered Aug 15 '15 at 07:14

J. M.'s missing motivation

124,525
11
401
574

Nice! That seems like a fairly general and useful function. – chris Aug 15 '15 at 09:46
Great. Thanks for that! – pdmclean Aug 15 '15 at 13:14
A little confession: the implementation in this answer is actually just a special case of the algorithm; in full generality, the algorithm solves the problem of expressing a matrix as a sum of Kronecker products. In that general case, one just takes multiple singular vectors, from which one can then build the required matrices . – J. M.'s missing motivation Jul 29 '16 at 16:07
@J.M. Would you be willing to post an implementation of the full algorithm if it's not too much work for you? On this post or a new one, either way. That would actually be pretty useful for things that I do. – march Jul 29 '16 at 19:06
@march, if you don't mind waiting for a few days, sure... :) – J. M.'s missing motivation Jul 29 '16 at 19:07
@J.M. That would be great, thanks! No rush, obviously! Do you think you'll post a new self-answered question or will you update this answer? – march Jul 29 '16 at 19:09
@march, have you seen the generalized routine already? – J. M.'s missing motivation Aug 02 '16 at 17:49
1

@J.M. Thank you for the update! I had not seen it. I do theoretical quantum mechanics calculations for part of my living, and systems with multiple degrees of freedom are treated using tensor-products of the different spaces associated with these dof's, and so this kind of decomposition can be useful for me. Thanks again! I'm awarding a bounty as a token of my appreciation. – march Aug 02 '16 at 17:55
@march, you're welcome; it was actually fun to write, and it helped that the paper described the algorithm very clearly. – J. M.'s missing motivation Aug 02 '16 at 20:52
looks like I have a couple of papers to read. :) – rcollyer Aug 02 '16 at 21:11
Very interesting from an image compression standpoint too if m3 is a grayscale image of dimensions {n^2,n^2} and dim1,dim2 are both dimension {n,n}. – flinty Jul 16 '21 at 21:08

Nearest Kronecker Product

1 Answers1

Linked

Related