Computing a partial density matrix from an amplitude vector
It's a matrix multiplication. Not sure how I didn't notice it right away, but the problem is exactly computing $U \cdot U^\dagger$ for a $U$ whose rows are made up of the chunks of amplitude vector.

def density_matrix_of_first_half(amplitude_vector):
  U = new ComplexMatrix(2^n, 2^n)
  for i &lt; 2^n:
    for j &lt; 2^n:
      U[i, j] = amplitude_vector[i + (j &lt;&lt; n)]
  P = U * U.adjoint()  // &lt;-- expensive
  return P


Fancy matrix multiplication does asymptotically better than the naive algorithm

Matrix multiplication can be done in time $O(n^\omega)$, where $n$ is the dimension of the matrix and $\omega$ is known to be below $2.375$. It's conjectured that there may be algorithms where $\omega = 2 + \epsilon$, but no one's managed to find any yet.

In the context of the question with $2^n \times 2^n$ matrices, the complexity is $O(2^{n \omega})$. However, because the asymptotically faster algorithms tend to have large constant factors hidden by the Big-O notation, in practice it's often better to stick with naive multiplication or Strassen multiplication. Also there are probably GPU-optimized matrix multiplication algorithms that would be more relevant given the note at the end of the question.

$U \cdot U^\dagger$ is not easier than $U \cdot V$

General matrix multiplication combines two totally independent matrices $U \cdot V$. Given that our density matrix problem only computes $U \cdot U^\dagger$, where the two inputs are totally dependent, one wonders if this "adjoint-squaring" is easier than general multiplication. Turns out it's not.