Langevin Sampler: Scalable Quantum Tomography

• The paper introduces a Langevin sampler that leverages Burer–Monteiro factorization to reduce dimensionality and enforce the Hermitian PSD structure in quantum tomography.
• It employs a spectral Student–t prior to promote low-rank solutions and directly incorporate physical constraints in the parameterization.
• The unadjusted Langevin algorithm achieves computational efficiency with PAC–Bayesian risk bounds matching state-of-the-art rates and improved mixing over traditional MCMC methods.

A Langevin sampler for quantum tomography is a Bayesian computational approach that leverages the Burer–Monteiro factorization to efficiently estimate quantum states from measurement data, with explicit low-rank structure and scalability guarantees. The method operates by parameterizing a $d \times d$ Hermitian positive semidefinite (PSD) density matrix via a product of a complex matrix and its conjugate transpose, imposing physical constraints directly in the parameter space. This enables the construction of a posterior distribution restricted to matrices of known or bounded rank and, through the use of a spectral Student–$t$ prior, promotes solutions of even lower rank when the true rank is unknown. The posterior is explored via an unadjusted Langevin algorithm (ULA), with rigorous PAC–Bayesian risk bounds that match state-of-the-art rates, and the algorithm achieves substantial computational savings compared to conventional Markov chain Monte Carlo (MCMC) techniques when the target density is low rank (Adel et al., 13 Jan 2026).

1. Parameterization via Burer–Monteiro Factorization

Quantum tomographic inference seeks a density matrix $\rho \in \C^{d \times d}$ satisfying $\rho = \rho^\dagger$, $\rho \succeq 0$, and $\tr(\rho) = 1$. If $\rho$ is known (or assumed) to have rank $r \leq d$, it is parameterized as $\rho = Z Z^\dagger$ with $Z \in \C^{d \times r}$. The unit-trace condition translates to $\|Z\|_F^2 = 1$, so $Z$ resides on the complex hypersphere.

Measurement data is modeled by $A$ observables, each with $S$ possible outcomes. The empirical frequencies $\hat p_{a,s}$ for observable $a$ and outcome $s$ are related to the Born prediction $p_{a,s}(Z) = \tr(\mathcal P_s^a Z Z^\dagger)$, where $\mathcal P_s^a$ denote the POVM elements. A pseudo-likelihood corresponds to a sum-of-squares loss:

$L(Z) = \sum_{a=1}^A \sum_{s=1}^S [\hat p_{a,s} - \tr(\mathcal P_s^a Z Z^\dagger)]^2$

Bayesian inference proceeds with the posterior $$p(Z | \text{data}) \propto \exp(-\lambda L(Z)) p(Z)$$, leading to potential $U(Z) = \lambda L(Z) - \log p(Z)$.

This factorization reduces the ambient parameter space from $d^2$ to $d r$ dimensions and automatically maintains Hermitian PSD structure.

2. Low-Rank–Promoting Spectral Prior

When the rank of $\rho$ is unknown but an upper bound $r$ is available, a spectral Student–$t$ prior is used to promote low-rank solutions:

$$p(Z) \propto \det(\theta^2 I_d + Z Z^\dagger)^{-\alpha},\quad \alpha = \tfrac{2d + r + 2}{2},\;\; \theta > 0$$

The prior decomposes as a product over the singular values $\sigma_i$ of $Z$: $$p(Z) \propto \prod_{i=1}^r (\theta^2 + \sigma_i^2)^{-\alpha}$$, heavily penalizing small $\sigma_i$ and thus favoring low-effective-rank $\rho$. The gradient of $\log p(Z)$ is available in closed form, facilitating efficient implementation:

$$\nabla_Z \log p(Z) = - (2d + r + 2)\, (\theta^2 I_d + Z Z^\dagger)^{-1} Z$$

This prior is a complex generalization of that studied by Dalalyan (2020) for promoting low-rank matrix estimation.

3. Langevin Sampler: Stochastic Dynamics and Discretization

The posterior on $Z$ is sampled via complex Langevin dynamics governed by the SDE:

$$\mathrm{d}Z_t = - \nabla_Z U(Z_t) \, \mathrm{d}t + \sqrt{2}\, \mathrm{d}W_t$$

where $W_t$ is Brownian motion in $\C^{d \times r}$.

Discretization through the unadjusted Langevin algorithm (ULA) with step size $\eta$ yields

$$Z_{k+1} = Z_k - \eta \nabla_Z U(Z_k) + \sqrt{2\eta}\, \Xi_k,\quad \Xi_k \sim \mathcal N(0, I_{dr})$$

To maintain the trace constraint $\|Z\|_F = 1$, two strategies are used: (1) projected Langevin, normalizing $Z_{k+1}$ after every step, or (2) unconstrained iteration with trace normalization applied only to the final estimator. Empirically, per-step drift in $\|Z\|_F$ is negligible, so the latter is often preferable for computational simplicity.

4. Implementation Workflow

The following high-level pseudocode summarizes the Langevin sampling scheme for quantum tomography:

Input: r, η, λ, θ, N, B, initial Z₀ ∈ ℂ^{d×r}, ‖Z₀‖_F=1 for k = 0 to N−1: G = ∇_Z [ λ L(Z_k) − log p(Z_k) ] Ξ ∼ 𝒩(0, I_{dr}) # i.i.d. complex Gaussian noise Z_{k+1} = Z_k − η G + √(2η) Ξ # Optionally: Z_{k+1} ← Z_{k+1} / ‖Z_{k+1}‖_F end for ρ̄ = (1/(N−B)) ∑_{k=B}^{N−1} Z_k Z_k^† return ρ̂ = ρ̄ / tr(ρ̄)

Key hyperparameter choices:

• λ controls the data-prior trade-off; for complete measurements, λ ≈ m/2 or 3m/8 is recommended.
• θ tunes the rank penalty; small θ strongly penalizes rank, large θ recovers a nearly uniform prior.
• η is chosen empirically for stability, typically $10^{-5}$–$10^{-4}$.
• N (number of iterations) and B (burn-in) are set to ensure convergence and posterior mixing.

5. PAC–Bayesian Risk Guarantees

In the full Pauli measurement setting with $N_{\text{tot}} = 3^n m$ samples and true rank-$r$ density $\rho^0 = Y^0 Y^{0†}$ with $Y^0 \in \C^{d \times r}$, a PAC–Bayesian Frobenius risk bound holds. For any $\epsilon \in (0,1)$ and comparison factor $\bar Y$ of rank $p$ such that

$$\|\bar Y - Y^0\|_F\,\|\bar Y + Y^0\|_F \leq 3^{-3n/2}2^{-n/2} m$$

it holds with probability at least $1-\epsilon$ that \begin{align*} |\hat\rho-\rho^0|_F^2 &\leq \frac{3}{N_{\text{tot}}\bigl(3^{3n/4}2^{(n+6)/4}(r+\sqrt{r}\,|\bar Y|F)+2r/m+1\bigr)} \ &\qquad + \frac{8\cdot3^n}{2^n N{\text{tot}}} \frac{\log(2/\epsilon)+2p(2^{n+1}+r+2) \log(1+|\bar Y|_2/\theta)}{} \end{align*} The analysis leverages exponential moment inequalities for the empirical squared error, KL-divergence bounds for shifted priors, and spectral properties of measurement operators. The leading term (up to log factors) matches the minimax rate $\frac{3^n\,\rank(\rho^0)}{N_{\text{tot}}}$ known from the literature (Mai & Alquier 2017, Mai 2021).

6. Computational Complexity and Empirical Performance

Each Langevin iteration involves:

• Data drift term: $\nabla L(Z)$ costs $O(A d r)$ operations with $A$ measurements (for complete Pauli, $A = 3^n$).
• Prior gradient: solving $(\theta^2 I_d + Z Z^\dagger)^{-1} Z$ is reduced via Sherman–Morrison–Woodbury to $O(d r^2 + r^3)$ for $r \ll d$.
• Noise sampling: $O(d r)$.

Total per-step cost: $O(A d r + d r^2 + r^3)$, dominated by $O(d r^2)$ for low-rank setting and moderate $A$.

Overall, runtime scales with the number of iterations times this per-step cost. Empirical benchmarks indicate:

• Scalability: With small $r$, each Langevin update is vastly faster than in full-rank MCMC schemes.
• Mixing: The chain mixes in fewer than $10^4$ steps, substantially fewer than required by Metropolis–Hastings methods for comparable accuracy.
• Estimation accuracy: The final estimator achieves Frobenius norm errors competitive with, or superior to, existing Bayesian quantum tomography algorithms.

7. Extensions and Practice Considerations

Potential refinements include step-size (η) annealing or adaptive temperature control for improved mixing. Metropolis-adjusted Langevin (MALA) or Riemannian variants may strengthen theoretical convergence guarantees. The method applies directly to process tomography (Choi matrix estimation) and can address incomplete measurement regimes. The reduction from $d^2$ to $d r$ dimensional parameter space, inherent PSD constraint enforcement, and low-rank–favoring prior make the Langevin sampler a scalable and theoretically sound approach for large-scale quantum tomography with explicit risk guarantees (Adel et al., 13 Jan 2026).