# Forming Kronecker products for “non-adjacent” vector spaces

In quantum mechanics, we often want to build operators acting on multiple particles by simply "Kronecker-ing" operators from the single-particle spaces. That is, if I have some matrix $A$ acting on particle one and $B$ on particle two, the appropriate matrix acting on their joint space is just $A\otimes B$. Mathematica gives you an easy way to do this with the function KroneckerProduct.

However, given a matrix $C$ which acts on the join particle space, it is not generically possible to find an $A$ and $B$ such that $C=A\otimes B$. This means that some operators cannot (naively) be built with Mathematica's native KroneckerProduct function: if I have three particles, and an operator $U$ which acts non-trivially on particles one and three but which does not factorize, I do not know how to tell Mathematica to build that matrix. Note that this difficulty does not occur for an operator acting on "adjacent" vector spaces: if $U$ instead acts on spaces two and three, we can write the total matrix as $\mathbb{I}\otimes U$, where $\mathbb{I}$ is the identity on particle 1.

I can think of a few ways around this, but they are all pretty inefficient. Are there any particularly good/clean ways to work this with Mathematica's native matrix functions?

• It depends a lot on the form of the expression for the operator $U$ that you start with. Without an example, it's not likely that you'll get an answer that matches your needs. See e.g. my answer here for spins where the question was about nearest-neighbor interaction but my code can (I think) handle the general case too. Another related answer is here. – Jens Aug 9 '18 at 4:34
• That's a fair point, but the problem is that the matrices I'm dealing with are random, and change with each instance of my code. I'll look at your references though - thanks! – miggle Aug 9 '18 at 5:21
• You can ArrayReshape result of KroneckerProduct, then Flatten it with "shuffled" dimensions. This is equivalent to transposition on corresponding tensor product space. So if you have three spaces $V_1$, $V_2$, $V_3$ with dimensions n1, n2, n3, and u13 matrix with Dimensions: {n1 n3, n1 n3} acting on product space $V_1 \otimes V_3$, then Flatten[ArrayReshape[KroneckerProduct[u13, IdentityMatrix[n2, SparseArray]], {n1, n3, n2, n1, n3, n2}], {{1, 3, 2}, {4, 6, 5}}] will give you corresponding matrix acting on $V_1 \otimes V_2 \otimes V_3$, with identity acting on $V_2$. – jkuczm Aug 10 '18 at 16:04
• Yes, this worked very well - thank you! I was using a similar strategy with TensorProduct but overall this was easier to code/understand. – miggle Aug 10 '18 at 21:15

Any $m \times n$ matrix $u$, with $m = \prod_{i=1}^k m_i$ and $n = \prod_{i=1}^l n_i$ (where $k,l,m_1,\ldots,m_k,n_1,\ldots,n_l \in \mathbb N$), represents operator between product spaces $\bigotimes_{i=1}^k V_i$ and $\bigotimes_{i=1}^l W_i$. Such matrix can be converted to matrix representing corresponding operator on spaces with shuffled constituent spaces: $\bigotimes_{i=1}^k V_{\sigma(i)}$, $\bigotimes_{i=1}^l W_{\delta(i)}$ (where $\sigma$ and $\delta$ are arbitrary permutations of respective index sets).

In Mathematica this can be done using ArrayReshape and Flatten functions:

Flatten[
ArrayReshape[u, {m@1, m@2, ..., m@k, n@1, n@2, ..., n@l}],
{{σ@1, σ@2, ..., σ@k}, {k + δ@1, k + δ@2, ..., k + δ@l}}
]


In case from OP we have three spaces $V_1$, $V_2$, $V_3$ and a matrix u13 representing endomorphism on $V_1 \otimes V_3$. Let's create example matrix u13, not representing a simple tensor, by adding two independent KroneckerProducts:

ClearAll[n, u, v]
n@1 = 2; n@2 = 3; n@3 = 4;
Function[mat, Table[(mat@i = Array[Subsuperscript[mat, Row@{##}, i] &, {n@i, n@i}]) // MatrixForm, {i, 3}]] /@ {u, v}
(u13 = KroneckerProduct[u@1, u@3] + KroneckerProduct[v@1, v@3]) // MatrixForm


Now we would like to create matrix acting on $V_1 \otimes V_2 \otimes V_3$, that acts as u13 on $V_1 \otimes V_3$ "part" and as some arbitrary matrix u@2 on $V_2$. We can do it by first creating matrix acting on $V_1 \otimes V_3 \otimes V_2$ using KroneckerProduct[u13, u@2] and then swap $V_3$ and $V_2$ using ArrayReshape and Flatten. Since in case of our example u13 matrix we know its decomposition to simple tensors, we can easily verify our result by comparing it with sum of appropriate Kronecker products:

Flatten[
ArrayReshape[KroneckerProduct[u13, u@2], {n@1, n@3, n@2, n@1, n@3, n@2}],
{{1, 3, 2}, {4, 6, 5}}
]
% == KroneckerProduct[u@1, u@2, u@3] + KroneckerProduct[v@1, u@2, v@3] // Expand
(* True *)


Finally to get matrix acting as identity on $V_2$ we just need to use IdentityMatrix instead of u@2:

Flatten[
ArrayReshape[
KroneckerProduct[u13, IdentityMatrix@n@2],
{n@1, n@3, n@2, n@1, n@3, n@2}
],
{{1, 3, 2}, {4, 6, 5}}
]


All used functions have specialized versions working with SparseArrays:

Flatten[
ArrayReshape[
KroneckerProduct[u13, IdentityMatrix[n@2, SparseArray]],
{n@1, n@3, n@2, n@1, n@3, n@2}
],
{{1, 3, 2}, {4, 6, 5}}
]
(* SparseArray[...] *)


and packed arrays:

On@"Packing"
Flatten[
ArrayReshape[
KroneckerProduct[
RandomReal[{-1, 1},{n@1 n@3, n@1 n@3}],
IdentityMatrix[n@2, WorkingPrecision -> MachinePrecision]
],
{n@1, n@3, n@2, n@1, n@3, n@2}
],
{{1, 3, 2}, {4, 6, 5}}
] // DeveloperPackedArrayQ
(* True *)
Off@"Packing"
`

so discussed method can benefit from these optimizations.