Dual formulation of the maximum entropy method applied to analytic continuation of quantum Monte Carlo data

TL;DR

This work tackles the ill-posed problem of analytic continuation for imaginary-time quantum Monte Carlo data by rethinking maximum entropy methods (MEM) through a dual Newton optimization. By formulating a dual, differentiable KL-regularized least-squares problem that operates in the kernel row-space, the authors retain the full spectral basis while achieving better conditioning and enabling second-order optimization. They provide analytic error bounds linking dual and primal solutions, data perturbations, and the regularization parameter, and demonstrate through LQCD and warm dense matter benchmarks that the dual Newton MEM yields more reliable uncertainty quantification and better resilience to noise than Bryan's algorithm. The results show improved spectral recovery, preserved features such as roton-like dispersions in DSFs, and practical improvements for analyzing authentic PIMC data, underscoring the method's impact for quantum many-body spectroscopy and related applications.

Abstract

Many fields of physics use quantum Monte Carlo techniques, but struggle to estimate dynamic spectra via the analytic continuation of imaginary-time quantum Monte Carlo data. One of the most ubiquitous approaches to analytic continuation is the maximum entropy method (MEM). We supply a dual Newton optimization algorithm to be used within the MEM and provide analytic bounds for the algorithm's error. The MEM is typically used with Bryan's controversial algorithm [Rothkopf, "Bryan's Maximum Entropy Method" Data 5.3 (2020)]. We present new theoretical issues that are not yet in the literature. Our algorithm has all the theoretical benefits of Bryan's algorithm without these theoretical issues. We compare the MEM with Bryan's optimization to the MEM with our dual Newton optimization on test problems from lattice quantum chromodynamics and plasma physics. These comparisons show that in the presence of noise the dual Newton algorithm produces better estimates and error bars; this indicates the limits of Bryan's algorithm's applicability. We use the MEM to investigate authentic quantum Monte Carlo data for the uniform electron gas at warm dense matter conditions and further substantiate the roton-type feature in the dispersion relation.

Figures (7)

• Figure 1: Plots of the parameterized $\rho$-meson spectral function (\ref{['eq:rhomesonSPF']}), the flat Bayesian prior, and the MEM estimates obtained using either Bryan's optimization algorithm (blue) or our dual Newton optimization algorithm (red). Plots keep number of correlator points constant (i.e., $N_\tau=30$) and vary the noise level $\sigma$. From left to right the noise increases as $\sigma^2 = 10^{-4}, 10^{-3}, 10^{-2}$. In the leftmost plot ($\sigma^2 = 10^{-4}$), the estimates and error bands are identical except near $\omega=0$. In the noiseless limit $\sigma^2 = 10^{-4}$ the different optimization algorithms produce similar curves, but when noise is large, MEM with Bryan's algorithm over-regularizes the solution.
• Figure 2: In rows 1 and 2 respectively, we present plots of the solutions proposed by the dual Newton algorithm and Bryan's algorithm at different regularization weights $\alpha$. All inversions use synthetic $\rho$-meson data generated from $N_\tau=30$ and the noise parameter increases from left to right, i.e., $\sigma^2 = 10^{-4}$ (left column), $\sigma^2 = 10^{-3}$ (middle column), $\sigma^2 = 10^{-2}$. (right column). In the bottom row we present the posterior weighting function used to combine the various solutions. The red and blue color bands indicate the $\alpha$ domain $[\alpha_\mathbf{min},\alpha_\mathbf{max}]$ kept by the MEM with dual Newton optimizer and MEM with Bryan's optimizer respectively; see (\ref{['eq:Bayesianweighting']}) for details. In all cases the posterior weighting function associated with the dual Newton method has a sharper peak at smaller $\alpha$ value than Bryan's primal approach. Comparing posterior weighting function across columns, as noise increases the peak of the posterior weighting function for Bryan's primal approach moves rightward. In comparison, the posterior weighting function for the dual Newton changes little.
• Figure 3: We study the theoretical error estimates presented in (\ref{['eq:error-x']}) for the $\rho$-meson problem. We present the upper bound on the error in solid black; $b, b'$ correspond to noisy ($\sigma^2 = 10^{-3}$) and noiseless data. Since $\lVert C^{-1} \rVert$, $\lVert A \rVert$, and $\lVert b -b' \rVert$ are constants, the solid black line shows a $1/\alpha$ scaling. The marked red and blue curves indicate the deviation between the algorithm's solution constructed from noisy data $x$ and the solution constructed from noiseless data $x'$; these two curves are almost identical. Lastly, we plot the deviation between Bryan's solution and the dual Newton's solution from noisy data, indicated by $\lVert x_{Bryan} - x_{dual \, N.}\rVert$, which shows that while the error scaling may be similar, the constructed solutions differ.
• Figure 4: Plots of the estimated dynamic structure factor for the synthetic UEG example with associated posterior weighting function. From left to right, the noise level $\sigma$ increases as $\sigma = 10^{-2}, 10^{-1}, 10^{0}$. The top row contains the plots of the DSF. The completed Mermin (CM) model is the the input spectrum and desired recovery, defined in (\ref{['eq:CMsusceptibility']}) with mean field correction (\ref{['eq:RPAcorrection']}) included. The RPA UEG is the Bayesian prior, defined in (\ref{['eq:idealgas']}) with mean field correction (\ref{['eq:RPAcorrection']}) included. The MEM estimates are in either red or blue. On the bottom row the Bayesian posterior is plotted as a function of $\alpha$. We logarithmically sample alpha in $\alpha \in [10^{0}-10^5]$. The red and blue vertical color bands indicates the $\alpha$ domain is selected for either MEM averaging procedure; see (\ref{['eq:Bayesianweighting']}) for details on how the $\alpha$ domain was selected. Comparing posterior weighting functions, as noise increases the peak of the posterior weighting function for Bryan's primal approach moves rightward. In comparison, the posterior weighting function for the dual Newton changes little.
• Figure 5: Top: From left to right we vary the wavenumber $k$ considered in our DSF, for all $k$ values, the number of correlator points is fixed at $201$ and the relative noise is $\approx 10^{-3}$. The curve labeled "stochastic" is Tobias et al.'s estimate dornheim2018stochasticsamplingalg, not the exact solution. The RPA UEG is the Bayesian prior, defined in (\ref{['eq:idealgas']}), with the mean-field correction (\ref{['eq:RPAcorrection']}) included. Bottom: Plot of the posterior weighting function associated with the above estimates. We collect solutions for $\alpha \in [10^{-3}-10^3]$. The red and blue color bands indicate the $\alpha$ domain kept by the MEM with dual Newton optimizer and MEM with Bryan's optimizer respectively; see (\ref{['eq:Bayesianweighting']}) for details on how the $\alpha$ domain was selected. We find strong overlap between the Bryan and dual Newton approach is readily explained by the small relative noise in the data.