19
$\begingroup$

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?
$\endgroup$
  • $\begingroup$ BTW does anyone know about factoring exactly? What about different norms? $\endgroup$ – pdmclean Aug 16 '15 at 2:31
  • $\begingroup$ "factoring exactly?" - Pitsianis-Van Loan already covers that. See my first example. $\endgroup$ – J. M. is away Aug 16 '15 at 20:43
21
+100
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Nice! That seems like a fairly general and useful function. $\endgroup$ – chris Aug 15 '15 at 9:46
  • $\begingroup$ Great. Thanks for that! $\endgroup$ – pdmclean Aug 15 '15 at 13:14
  • $\begingroup$ 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 . $\endgroup$ – J. M. is away Jul 29 '16 at 16:07
  • $\begingroup$ @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. $\endgroup$ – march Jul 29 '16 at 19:06
  • 1
    $\begingroup$ @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. $\endgroup$ – march Aug 2 '16 at 17:55

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.