partial trace with sparse matrices

Question

Let $\rho_{ABCD}$ be a sparse matrix of 4 systems each in a $d$-dimensional Hilbert space.

For $d<7$ in a reasonable time (few seconds) I able to perform the partial trace $\rho_{AD}$ using the code proposed in http://www3.imperial.ac.uk/people/m.tame/research. I would need an efficient algorithm for calculating $\rho_{AD}$ where $d\geq7$. The algorithm of the site above does not exploit any properties of the matrix and it requires a lot of permutations and rearrangements.

Do you know an efficient algorithm for calculating the partial trace of qudit which uses the fact that the matrix is sparse? It would also be interesting if the algorithm can take advantage of parallel computation.

Thank you very much in advance for your answers.

Regards,

Silvio

This post has been migrated from (A51.SE)

Comments

Strictly speaking, you can only hope for a polynomial scaling with the dimension of the system $d$ (I would not call this efficient in computational complexity terms, unless you set d to be a constant). It would be useful to know how sparse is your density matrix: how does the maximal number-of-coefficients per rows *and* column asymptotically scale with the dimesion of the system $d$?

This post has been migrated from (A51.SE)

---

In general, a partial trace of a sparse matrix does not yield a sparse matrix, so you will get only limited improvement by using the fact that the matrix is sparse (maybe a factor of d/s).

This post has been migrated from (A51.SE)

Answer 1

If you want something specific to Mathematica then I don't know, but in general:

Let $\sigma = \rho_{AD} = \textrm{Tr}_{BC}(\rho)$. $\rho$ is an operator on four subsystems, so it has four inputs and four outputs making it a rank 8 tensor. Let $i,i'$ be the indices corresponding to the input and output of $\rho$ on subsystem $A$, let $j,j'$ correspond to $B$, etc. The components of $\sigma$ are then $\sigma_{ii'll'} = \sum_{jk} \rho_{ii'jjkkll'}$.

If $\rho$ is sparse you only need to sum the entries which are nonzero. If $\rho$ is Hermitian then $\sigma$ will be Hermitian and it is enough to only compute the upper triangle (the lower triangle is then the conjugate of that). I don't believe there are any other optimizations available unless you know of further structure on $\rho$. The sum is trivially parallelizable - each combination $ii'll'$ is completely independent of the others.

This post has been migrated from (A51.SE)

Comments

The efficiency of this method strongly depends on how costly is to list the coefficients of the sparse which are non-zero. It would be good if the OP could clarify this point.

This post has been migrated from (A51.SE)

---

Thank you very much for your answer. Using your suggestion and looking more carefully to the algorithm proposed in Mathematica (link in the question) I was able to parallelize the algorithm and now I am able to arrive until $d\leq9$. Regarding further structures or the asymptotic scaling of the non-zero elements I don't have a so deep analysis because it is not needed for my problem. I do believe that there is not any interesting property that can be used for the partial trace because I apply one unitaries to some states and I see that there are "a lot" of zeros. Thank you very much again.

This post has been migrated from (A51.SE)