# Simple and efficient construction of tensor product operators

Problem Description

I'm using mathematica to study a problem in quantum mechanics, which is naturally understood in terms of vector spaces. In the simplest example, suppose we have a 16-dimensional vector space built out of the tensor product of four two-dimensional spaces. I need a simple (and ideally, efficient) code to construct operators which act on different subsystems. I think this is doable, but I'm struggling with correctly organizing arrays.

Here is an example. Suppose I have a random unitary matrix $$U_{12}$$ which acts on the first two two-dimensional spaces, but leaves the others alone. So we can write it as $$U_{12}\otimes I_{4}$$, where $$I_{4}$$ is the $$4\times4$$ identity matrix. This is easy to construct in mathematica:

U12 = RandomVariate[CircularUnitaryMatrixDistribution[4];
UBig = KroneckerProduct[U12, IdentityMatrix[4]];


This does what I want it to. What I'm in need of is a code that generates a matrix representation of the analogous operator $$U_{14}$$, which acts non-trivially on subsystems 1 and 4. More generally, it would be nice to have a code which can generate matrix representations of operators that act on subsystems $$i,j$$ when they aren't neighbors (i.e., KroneckerProduct would only work if $$i = j\pm 1$$). Also, my case of interest would involve roughly ten such two-dimensional spaces, so it would be nice if the solution scaled decently.

My attempt at a solution

Here is my thinking as I have tried to find a solution. Nothing has worked yet so it would be great if someone could tell me where I've gone wrong or suggest an alternative. We could start by generating a random unitary matrix that acts on subspaces 1 and 2, and reshaping the array to have the correct index structure:

U12 = RandomVariate[CircularUnitaryMatrixDistribution[4];
UBig = KroneckerProduct[U12, IdentityMatrix[4]];
UBig2 = ArrayReshape[UBig, {2,2,2,2,4,4}]


The thinking here is that I'll want to re-arrange indices such that we end up with U12 acting on subspaces (1,4). This can only happen if we have 2-valued indices to rearrange; the 4-valued indices should be blocks left alone to parameterize the identity for the subspace we aren't touching. Now let's rearrange the indices:

UBig3 = TensorTranspose[UBig2, Cycles[{{2, 5}, {4, 6}}]]


I think this is the correct cycle based on a diagram I've drawn - the dimensions of this tensor are {2,4,2,4,2,2}, where I understand the first and third indices index subspace 1, while indices 5 and 6 index subspace 4, and the 4-valued indices are the identity space (subspaces 2,3). Now I think we want to flatten these indices such that the there is a single "incoming" and "outgoing index for the (1,4) subspace, i.e.,

UBig4 = Flatten[UBig3, {{1,5},{3,6}};


Now we have a 16x16 matrix, which is indeed what I want. Sadly, as you can check, it is not unitary, so something has gone wrong. This operation should preserve unitarity.

This type of computation can be done in two different ways, using TensorProduct or using KroneckerProduct. I find easier to use TensorProduct, but I give here how I'd address this with both methods:

Define the following objects and note how similar the inputs are in the last two lines:

Id2 = IdentityMatrix[2];
U12 = RandomVariate[CircularUnitaryMatrixDistribution[4]];
U12T = Transpose[ArrayReshape[U12, {2, 2, 2, 2}], {1, 3, 2, 4}];
UBigK = KroneckerProduct[U12, Id2, Id2];
UBigT = TensorProduct[U12T, Id2, Id2];


UBigK and UBigT contain the same inofmation, and in fact we can convert each one into the other:

In[]:= Flatten[UBigT, {{1, 3, 5, 7}, {2, 4, 6, 8}}] === UBigK
Out[]= True

In[]:= Transpose[ArrayReshape[UBigK, {2, 2, 2, 2, 2, 2, 2, 2}], {1, 3, 5, 7, 2, 4, 6, 8}] === UBigT
Out[]= True


Let us first discuss UBigT. It is an operator with 8 indices, organized by four consecutive pairs corresponding to operation in your for vector spaces. If we want to exchange the roles of spaces 2 and 4 then we need to change indices 3 and 4 with indices 7 and 8 respectively. Then we flatten again as we did before to get the 16x16 unitary matrix you want:

In[]:= UnitaryMatrixQ[res1 = Flatten[Transpose[UBigT, Cycles[{{3, 7}, {4, 8}}]], {{1, 3, 5, 7}, {2, 4, 6, 8}}]]
Out[]= True


To perform the computation directly on UBigK seems more difficult to me because it requires a nontrivial permutation. A possible way to get it is as follows. Construct a representation of the basis states:

In[]:= states = Ket @@@ Tuples[{0, 1}, 4]
Out[]= {Ket[0, 0, 0, 0], Ket[0, 0, 0, 1], Ket[0, 0, 1, 0], Ket[0, 0, 1, 1],
Ket[0, 1, 0, 0], Ket[0, 1, 0, 1], Ket[0, 1, 1, 0], Ket[0, 1, 1, 1],
Ket[1, 0, 0, 0], Ket[1, 0, 0, 1], Ket[1, 0, 1, 0], Ket[1, 0, 1, 1],
Ket[1, 1, 0, 0], Ket[1, 1, 0, 1], Ket[1, 1, 1, 0], Ket[1, 1, 1, 1]}


Now find the permutation that exchanges slots 2 and 4:

In[]:= perm24 = PermutationList[FindPermutation[states, states /. Ket[a1_, a2_, a3_, a4_] :> Ket[a1, a4, a3, a2]], 16]
Out[]= {1, 5, 3, 7, 2, 6, 4, 8, 9, 13, 11, 15, 10, 14, 12, 16}


Apply to both indices of the matrix:

In[]:= UnitaryMatrixQ[res2 = Part[UBigK, perm24, perm24]]
Out[]= True


Check that the two results coincide:

In[]:= res1 === res2
Out[]= True