Tensor Products in Code

Tensor products become reliable numerical objects only after the order of subsystems and their basis indices is fixed. The same vector can describe different states under different orderings. This page gives an explicit mixed-radix convention and uses it to apply local operators efficiently.

Helpful background. Partial trace supplies the subsystem interpretation; array representations treats single-system storage.

Mixed-radix indexing

For subsystem dimensions $(d_0,\ldots,d_{n-1})$, order basis states as $|i_0\rangle\otimes\cdots\otimes|i_{n-1}\rangle$ and flatten with the last index varying fastest:

$$q=\sum_{a=0}^{n-1}i_a\prod_{b=a+1}^{n-1}d_b.$$

For two subsystems, $q=i_0d_1+i_1$. NumPy’s default C-order reshape matches this convention. Bit-based codes often make site 0 the least significant bit instead; an explicit mapping is required when exchanging data.

A product state has amplitudes $\Psi_{i_0\ldots i_{n-1}}=\prod_a\psi^{(a)}_{i_a}$. Thus a Kronecker product of statevectors reproduces this ordering. Verify the map with unequal dimensions, such as $2\times3$: symmetric qubit examples can conceal swapped axes.

Local operator contractions

An operator on subsystem $a$ is formally $I\otimes\cdots\otimes A\otimes\cdots\otimes I$, but constructing this full matrix is usually unnecessary. Reshape the state, move axis $a$ to the front, flatten the remaining axes, and multiply by $A$:

$$\Psi'_{i_a,\boldsymbol j} =\sum_{k_a}A_{i_ak_a}\Psi_{k_a,\boldsymbol j}.$$

Restore the axes afterward. A two-site operator uses two target axes in a documented order. Swapping their order changes a controlled gate or a nonsymmetric interaction.

import numpy as np

def local_action(psi, op, target, dims):
    tensor = np.asarray(psi).reshape(dims)
    front = np.moveaxis(tensor, target, 0)
    result = op @ front.reshape(dims[target], -1)
    result = result.reshape(front.shape)
    return np.moveaxis(result, 0, target).reshape(-1)

This snippet assumes a square local operator preserving subsystem dimension. For an isometry with different output dimension, update both the shape and the dimension list explicitly.

A complex two-subsystem check

Let $\psi=|0\rangle\otimes(|0\rangle+i|1\rangle)/\sqrt2$, so its flattened vector is $(1,i,0,0)/\sqrt2$. Applying $X$ to the first subsystem gives $(0,0,1,i)/\sqrt2$; applying $X$ to the second gives $(i,1,0,0)/\sqrt2$.

These two expected outputs diagnose subsystem reversal. Test a complex operator such as $Y$ as well: transpose and conjugation errors may survive an $X$-only test. In a small system compare contraction output to the explicit Kronecker matrix.

An unequal-dimension check is even more discriminating. In dimensions $(2,3)$ take $\psi=(1,i,2,0,0,0)/\sqrt6$ and $Y=\begin{pmatrix}0&-i\\i&0\end{pmatrix}$. Acting on subsystem 0 must give $(0,0,0,i,-1,2i)/\sqrt6$. Acting on the qutrit with $\operatorname{diag}(1,2,3)$ instead gives $(1,2i,6,0,0,0)/\sqrt6$; this deliberately nonunitary test must not be renormalized by the kernel.

For ordered targets $(a,b)$, move those axes to the front and flatten their local index as $r=i_a d_b+i_b$. Multiply the $d_ad_b$ rows by the specified two-site matrix, then invert the axis permutation. The targets need not be adjacent. Test both $(a,b)$ and $(b,a)$ with a nonsymmetric gate: the same numerical matrix represents different physical operators under those two interpretations. Reshaping without the required permutation does not relocate an operator.

Cost and physical restrictions

For total dimension $D=\prod_a d_a$, a fixed-size local action costs $O(D)$ up to local dimension factors and needs $O(D)$ state storage. A dense global matrix would require $O(D^2)$ memory. Axis permutations can create copies, so peak memory depends on array layout and library behavior.

A tensor factorization is not automatically a factorization into distinguishable physical particles. Fermionic mode operators need parity signs; applying bare local raising and lowering matrices on a qubit tensor product does not implement the canonical anticommutation relations. Use fermionic basis construction for that case.

Finally, reduced symmetry-sector vectors are not full tensor arrays. Embed them through their basis map before a generic reshape or implement contractions directly in the reduced basis.

Exercise

In dimensions $(2,3,2)$, find the flat index of $(1,2,0)$ and the indices changed by a local $X$ on the final subsystem.

Solution

$q=1(3\cdot2)+2(2)+0=10$. The final bit toggles $q=10$ and $q=11$ while leaving the first two indices unchanged. A bitwise toggle of the most significant bit would instead act on a different subsystem.

References

• QuTiP developers, “Tensor Products and Partial Traces”, QuTiP documentation — subsystem dimensions and operator placement.
• Watrous, J., The Theory of Quantum Information, Cambridge University Press (2018), author’s text — tensor spaces and linear maps.