Forming Kronecker products for "non-adjacent" vector spaces

Question

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?

Answer 1

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}}
] // Developer`PackedArrayQ
(* True *)
Off@"Packing"

so discussed method can benefit from these optimizations.