partial trace with sparse matrices

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.

Regards,

Silvio

-
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

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