Physics Stack Exchange is a question and answer site for active researchers, academics and students of physics. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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 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.



share|cite|improve this question
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$? – Juan Bermejo Vega Apr 18 '12 at 7:10
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). – Joe Fitzsimons Apr 18 '12 at 8:20

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.

share|cite|improve this answer
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. – Juan Bermejo Vega Apr 18 '12 at 7:15
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. – Silvio Abruzzo Apr 19 '12 at 8:35

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.