From: Ian Hincks on
Hello,

I am working on a script to perform quantum process tomography on three qubits. There is a bottleneck in my computation that I don't know how to get around. Here's the setup:

I have an orthonormal basis, called E, for the set of 8x8 matrices, just indexed as E{1}, ..., E{64}.

I need to compute the 4-dimensional array B defined by the following equation:

E{m}*E{j}*E{n} = sum-k( B(m,n,j,k) * E{k} )

Where the first two *'s are matrix multiplication, and the third * is the scalar B{m,n,j,k} multiplied with the matrix E{k}, and sum-k means to sum over all 64 k-indeces.

The naive solution, and also the only one I know, is to exploit the orthonormality of E and use the Hilbert-Schmidt inner product [ <x,y> := trace( (x')*(y) ) where ' is the conjugate transpose] to get:

B(m,n,j,k) = trace(E{k}*(E{m}*E{j}*E{n}))

where the conjugate transpose has been neglected because each E happens to be hermitian. This method involves four nested for loops to get all 64^4 values, and takes a very long time.

Can anyone suggest a better method? I am relatively new to matlab.

FYI: Here is my code:
% Get the E basis
P{1} = [1 0; 0 1];
P{2} = [0 1; 1 0];
P{3} = [0 -i; i, 0];
P{4} = [1 0; 0 -1];

m = 0;
for a = 1:4
for b = 1:4
for c = 1:4
m = m + 1;
% Get basis element, and normalize
E{m} = kron3(P{a}, P{b}, P{c})/sqrt(8);
end
end
end

% Calculate all 64^4 values for B
for m = 1:64
for n = 1:64
for j = 1:64
V = E{m}*E{j}*E{n};
for k = 1:64
B(m,n,j,k) = trace(E{k}*V);
end
end
end
end


And kron3 is just the function kron3(x,y,z) = kron(kron(x,y),z)


Thanks!
Ian
From: Ian Hincks on
Maybe I've bogged this down with too much context. Here is a summary of that which is relevant:

E is a set of 64 8x8 matrices that satisfy the property trace(E{i}*E{j}) = (1 if i==j, 0 else). I need to solve either of the following two (equivalent) problems for B efficiently:

E{m}*E{j}*E{n} = sum(B(m,n,j,k) * E{k}, k=1:64)

OR

B(m,n,j,k) = trace(E{k}*(E{m}*E{j}*E{n}))

Where B is a 64x64x64x64 complex array.

Thanks,
Ian