Faster randomized partial trace estimation

TL;DR

The paper tackles the computational bottleneck of obtaining reduced density matrices via partial traces in quantum systems where the total density matrix is exponentially large. It introduces a variance-reduced, matrix-free approach that deflates the top eigenmodes of the matrix function (notably $A=e^{-eta H}$) and employs a Lanczos-based partial-trace approximation to compute $ ext{tr}_{ ext{b}}(f(H))$ efficiently. The main contributions are a deflation-based variance reduction that yields orders-of-magnitude speedups over prior methods, a practical algorithm combining implicit partial-trace estimation with block-Lanczos deflation, and extensive numerical demonstrations on Heisenberg spin systems showing accurate entanglement spectra and ergotropy calculations. This work enables scalable thermodynamic and entanglement analyses of larger quantum many-body subsystems and highlights opportunities for cross-disciplinary collaboration between quantum physics and numerical linear algebra.

Abstract

We develop randomized matrix-free algorithms for estimating partial traces, a generalization of the trace arising in quantum physics and chemistry. Our algorithm improves on the typicality-based approach used in [T. Chen and Y-C. Cheng, \emph{Numerical computation of the equilibrium-reduced density matrix for strongly coupled open quantum systems}, J. Chem. Phys. 157, 064106 (2022)] by deflating important subspaces (e.g. corresponding to the low-energy eigenstates) explicitly. This results in a significant variance reduction, leading to several order-of-magnitude speedups over the previous state of the art. We then apply our algorithm to study the thermodynamics of several Heisenberg spin systems, particularly the entanglement spectrum and ergotropy.

Figures (10)

• Figure 1: Let $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_{d_{\mathrm{t}}}$ be the singular values of $\mathbf{A}$. Then the variance of the partial trace estimator \ref{['eqn:quadratic_partial_trace_estimator']} is bounded above by $2\|\mathbf{A}\|_\mathsf{F}^2 = 2(\sigma_1^2 + \cdots + \sigma_{d_{\mathrm{t}}}^2)$. By deflating the top $k$ singular values, we can reduce the variance to at most $2(\sigma_{k+1}^2 + \cdots + \sigma_{d_{\mathrm{t}}}^2)$ (see \ref{['sec:variance_reduce']}). Here we take $\mathbf{A} = \bm{\rho}_{\mathrm{t}} = \exp(-\beta \mathbf{H}) / \operatorname{tr}(\exp(-\beta \mathbf{H}))$, where $\mathbf{H}$ corresponds to the solvable spin chain with $N=10$ and $h=0.3$ (described in \ref{['sec:solvable_chain']}). We then plot the variance bound $2(\sigma_{k+1}^2 + \cdots + \sigma_{d_{\mathrm{t}}}^2)$ for several values of $k$ and a range of $\beta$. Takeaway: The $k=0$ curve exhibits high variance at low-temperature. Through the use of deflation, the variance can be reduced significantly.
• Figure 1: Using the same example as \ref{['fig:variance_example']}, we plot the approximate variance of the partial trace estimator for $\mathbf{A} = \exp(-\beta \mathbf{H}_{\mathrm{t}}) / \operatorname{tr}(\exp(-\beta \mathbf{H}_{\mathrm{t}}))$ when $\mathbf{Q} = \operatorname{orth}(\exp(-\beta_0 \mathbf{H}_{\mathrm{t}})\mathbf{\Omega})$ is an approximate top subspace. We use $J/\beta_0 = 0.1$ (dotted vertical line) and sample $\mathbf{\Omega}\in\mathbb{R}^{d_{\mathrm{t}} \times k}$ with independent standard normal entries. The light curves are those of \ref{['fig:variance_example']} and correspond to the optimal rank-$k$ subspace. While the approximate top subspace does reduce the variance, the variance does not go to zero in the zero temperature limit $\beta\to\infty$. The right plot shows the normalized singular values of $\exp(-\beta_0 \mathbf{H}_{\mathrm{t}})$, which decay rapidly.
• Figure 1: Kagome-strip chain with $N=20$ sites and periodic boundary conditions. Subsystem (s) is encircled.
• Figure 2: Von Neumann entropy for Kagome-strip chain at varying values of $J_2$, $h$, and $\beta$.
• Figure 3: Entanglement spectrum for Kagome-strip chain at varying values of $J_2$ and $h$ at fixed temperature $J/\beta = 2\times 10^{-2}$.