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
Comments
Post a Comment