Suppose you're given the list of $2^{2n}$ complex coefficients defining the state of a $2n$ qubit register. You want to compute the $2^n \times 2^n$ density matrix of the first $n$ qubits.
Is it possible to compute the result in $O(4^n)$ time? That would be optimal, since the output has size $\Theta(4^n)$.
It's easy to hit $O(8^n)$. The naive algorithm that simply computes the density matrix for each possible value of the unrelated qubits and then adds those together achieves that:
def density_matrix_of_first_half(amplitude_vector):
v = amplitude_vector
M = new complex[2^n, 2^n]
for r < 2^n:
offset = r * 2^n
for i < 2^n:
for j < 2^n:
M[i, j] += (v[i + offset] * v[j + offset]).conjugate()
return M
But I was hoping it was possible to do better. Are there any papers about this problem?
(The reason I care so much more about $4^n$ vs $8^n$ is because I'm actually computing things with a GPU. My textures can have $2^{2n} = 4^{n}$ pixels, so I get the $4^n$ inner-loop in the application of a single shader. But $8^n$ is far over the max texture size, so the $8^n$ algorithm would require a pipeline of $2^n$ shaders. Or at least might force me to drop the qubit limit another 30%.)