Computing partial traces and reduced density matrices

Taking partial traces for computing reduced density matrices, or related functions, is a ubiquitous procedure in the quantum mechanics of composite systems. In this article, we present a thorough description of this function and analyze the number of elementary operations (ops) needed, under some possible alternative implementations, to compute it on a classical computer. As we notice, it is worthwhile doing some analytical developments in order to avoid making null multiplications and sums, what can considerably reduce the ops. For instance, for a bipartite system $\mathcal{H}_{a}\otimes\mathcal{H}_{b}$ with dimensions $d_{a}=\dim\mathcal{H}_{a}$ and $d_{b}=\dim\mathcal{H}_{b}$ and for $d_{a},d_{b}\gg1$, while a direct use of partial trace definition applied to $\mathcal{H}_{b}$ requires $\mathcal{O}(d_{a}^{6}d_{b}^{6})$ ops, its optimized implementation entails $\mathcal{O}(d_{a}^{2}d_{b})$ ops. In the sequence, we regard the computation of partial traces for general multipartite systems and describe Fortran code provided to implement it numerically. We also consider the calculation of reduced density matrices via Bloch's parametrization with generalized Gell Mann's matrices.