Algorithms for optimized maximum entropy and diagnostic tools for analytic continuation

Analytic continuation of numerical data obtained in imaginary time or frequency has become an essential part of many branches of quantum computational physics. It is, however, an ill-conditioned procedure and thus a hard numerical problem. The maximum-entropy approach, based on Bayesian inference, is the most widely used method to tackle that problem. Although the approach is well established and among the most reliable and efficient ones, useful developments of the method and of its implementation are still possible. In addition, while a few free software implementations are available, a well-documented, optimized, general purpose, and user-friendly software dedicated to that specific task is still lacking. Here we analyze all aspects of the implementation that are critical for accuracy and speed and present a highly optimized approach to maximum entropy. Original algorithmic and conceptual contributions include (1) numerical approximations that yield a computational complexity that is almost independent of temperature and spectrum shape (including sharp Drude peaks in broad background, for example) while ensuring quantitative accuracy of the result whenever precision of the data is sufficient, (2) a robust method of choosing the entropy weight $\alpha$ that follows from a simple consistency condition of the approach and the observation that information- and noise-fitting regimes can be identified clearly from the behavior of ${\chi }^2$ with respect to $\alpha$, and (3) several diagnostics to assess the reliability of the result. Benchmarks with test spectral functions of different complexity and an example with an actual physical simulation are presented. Our implementation, which covers most typical cases for fermions, bosons, and response functions, is available as an open source, user-friendly software.

Figure 1

Schematic behavior of $log\left({\chi }^2\right)$ vs $log\left(\alpha \right)$. Three regimes can be identified.

Figure 2

(a) Spectrum obtained with the program at a temperature $T=0.05$, with a constant relative diagonal error of ${10}^{-5}$; (b) ${\chi }^2$ as a function of $\alpha$ for the example in (a); (c) curvature computed from ${log}_{10}{\chi }^2$ as a function of $\gamma {log}_{10}\alpha$, where $\gamma =0.2$, plotted as a function of ${log}_{10}\alpha$; and (d) spectrum as a function of $\alpha$ at sample frequencies corresponding to the three extrema of the spectrum in (a).

Figure 3

Normalized errors and autocorrelations as a function of frequency index at three different values of $\alpha$ for the example given in Fig. 2. (a) Error and (b) autocorrelation at $\alpha =1000{\alpha }^{*}$; (c) error and (d) autocorrelation at $\alpha =10{\alpha }^{*}$; (e) error and (f) autocorrelation at ${\alpha }^{*}$.

Figure 4

(a) Spectrum obtained with the program at a temperature $T=0.002$, with a constant relative diagonal error of ${10}^{-6}$; (b) magnification on the low frequency peak; (c) ${\chi }^2$ as a function of $\alpha$; and (d) spectrum as a function of $\alpha$ at sample frequencies corresponding to the maxima. The value for the $\omega =0$ peak is given on the left axis while the values for the two other peaks are given on the right axis.

Figure 5

(a) Spectrum for single-site DMFT on a square lattice at $U=6$, filling $n=1$ and temperature $T=0.01$; (b) ${\chi }^2$ as a function of $\alpha$ for the above example; (c) spectrum as a function of $\alpha$ at sample frequencies. The relative error in the data is of order $5\times {10}^{-4}$ at most.

Figure 6

Example of frequency grid density $1/\mathrm{\Delta }\omega$ as a function of $\omega$ around the central part of the grid. In this example, the user has provided parameters defining a grid with different constant frequency steps in different intervals in the form of Eq. (26). The program matches smoothly the steps in the different regions using a hyperbolic tangent shape. Then the distance $\mathrm{\Delta }u=1/({\omega }_{i+1}-{\omega }_{\mu })-1/({\omega }_i-{\omega }_{\mu })$ (with ${\omega }_{\mu }$ fixed by the cutoff) is constant outside the central interval $[-5,5]$.

Figure 7

Example of nonuniform Matsubara frequency index grid. The vertical axis is the index of the Matsubara frequency and the horizontal axis numbers the frequencies that are kept.