Reshaping a matrix for Reverse Kronecker product or Schmidt decomposition

General Usage
,
1

Before I begin, I hope MathJax works on Discourse.
I am trying to do a generalization of Schmidt decomposition on a matrix, which means writing down a matrix H of size m\times n, where m = m_1 \times m_2 and n = n_1\times n_2, as

H = \sum_j \lambda_j H_{1,j}\otimes H_{2,j}\ ,

using the minimum number of terms possible. Here, H_{1,j} and H_{2,j} are matrices of sizes m_1 \times n_1 and m_2 \times n_2 respectively and \otimes is a Kronecker product. What I was trying to do is similar to what is explained in this answer, with a tiny difference at the end.

There are three main steps in this algorithm:

  1. Rearranging elements of H into a matrix of size m_1 n_1 \times m_2 n_2
  2. Doing a Singular Value Decomposition on the result from previous step
  3. Rearranging the outcome of SVD into the form we wanted

The second and third step are straightforward in Julia. For the first step what I initially tried was to reshape the Matrix H into a tensor of rank m_1\times m_2 \times n_1 \times n_2 using reshape and then swap the second and third index using permutedims and then rearrange it back into a matrix using reshape again.

H1 = reshape(H,(m1,m2,n1,n2))
H2 = permutedims(H1,(1,3,2,4))
H3 = reshape(H2,(m1*n1,m2*n2))

However, reshape doesn’t rearrange the elements as I want them to. For example, if I run the procedure above on the example given in the math.stackexchange answer, I get:

\begin{bmatrix} a_{11} & a_{41} & a_{13} & a_{43} \\ a_{21} & a_{51} & a_{23} & a_{53} \\ a_{31} & a_{61} & a_{33} & a_{63} \\ a_{12} & a_{42} & a_{14} & a_{44} \\ a_{22} & a_{52} & a_{24} & a_{54} \\ a_{32} & a_{62} & a_{34} & a_{64} \end{bmatrix}

instead of what is needed:

\begin{bmatrix} a_{11}& a_{21} & a_{12} & a_{22} \\ a_{31}& a_{41} & a_{32} & a_{42} \\ a_{51}& a_{61} & a_{52} & a_{62} \\ a_{13}& a_{23} & a_{14} & a_{24} \\ a_{33}& a_{43} & a_{34} & a_{44} \\ a_{53}& a_{63} & a_{54} & a_{64} \end{bmatrix}

So, my question is: Is there an easy/obvious way of doing this transformation or am I stuck with writing a depth 4 nested for loop?

2

The Kronecker product orders dimensions inversely than how Julia orders tensors. The first dimension in julia corresponds to the last dimension of the Kronecker product. The simplest way to fix you problem would therefore be to swap the order of the two matrices.

In a way, this actually makes sense: it’s like a binary number (or any number really) were the first number is on the right.

Also see https://github.com/JuliaLang/julia/pull/28697 for more discussion. Note that I don’t believe that there is any problem with how things work currently. It’s just a matter of choosing which dimension is “first”.

4

Thanks. I see. So I should do something like this instead:

H1 = reshape(H,(m2,m1,n2,n1))
H2 = permutedims(H1,(1,3,2,4))
H3 = permutedims(reshape(H2,(m2*n2,m1*n1)),(2,1))

As was mentioned in the thread you linked it might be wise for me to write my own tensor product function. Luckily ⊗ is not binded.

5

Try it out! I’ll admit I haven’t looked at exactly what you needed to do, but I figured that the dimension ordering was the problem.