Average spectrum method for analytic continuation: Efficient blocked-mode sampling and dependence on the discretization grid

The average spectrum method is a promising approach for the analytic continuation of imaginary time or frequency data to the real axis. It determines the analytic continuation of noisy data from a functional average over all admissible spectral functions, weighted by how well they fit the data. Its main advantage is the apparent lack of adjustable parameters and smoothness constraints, using instead the information on the statistical noise in the data. Its main disadvantage is the enormous computational cost of performing the functional integral. Here we introduce an efficient implementation, based on the singular value decomposition of the integral kernel, eliminating this problem. It allows us to analyze the behavior of the average spectrum method in detail. We find that the discretization of the real-frequency grid, on which the spectral function is represented, biases the results. The distribution of the grid points plays the role of a default model while the number of grid points acts as a regularization parameter. We give a quantitative explanation for this behavior, point out the crucial role of the default model and provide a practical method for choosing it, making the average spectrum method a reliable and efficient technique for analytic continuation.

Figure 1

Schematic contour plot of the Gaussian probability density $exp(-{\chi }^2\left(\mathbf{f}\right)/2)$ in the plane of two values $f_1$ and $f_2$. The unphysical region $\mathbf{f}<0$ is shaded in gray. In components sampling the moves $f_i\to f_i^{\prime }$ are proposed parallel to the coordinate axes, resulting in narrow Gaussians of widths that are of the order of $1/\max \left({\mathrm{d}}_{\mathrm{i}}\right)$. In modes sampling, moves $e_i\to e_i^{\prime }$ are proposed along the principal axes of the multivariate Gaussian, so that the moves in directions corresponding to small singular values can take large steps. Note that for ill-conditioned problems the singular values $d_i$ vary over many orders of magnitude.

Figure 2

Example of the hierarchy of grid partitionings used in blocked-mode sampling. At the highest level (top), the grid on which $\mathbf{f}$ is represented forms a single block. Sampling on this block is modes sampling. At the level below, the grid is split into two blocks. If going to a lower level we split the blocks in half, there would always be a block boundary at the center of the grid. To avoid this, we shift the intervals at every other level by half their width. At the lowest level (bottom), the blocks are the individual intervals $f_n$. Sampling on these blocks is components sampling.

Figure 3

Efficiency of different sampling methods comparing the changes in $\mathbf{f}$ between Monte Carlo steps, defined as ${\mathrm{\Delta }}_i=\parallel {\mathbf{f}}^{(i+1)}-{\mathbf{f}}^{\left(i\right)}\parallel /\parallel {\mathbf{f}}^{\left(i\right)}\parallel$. For a grid with small cutoff (top), modes and blocked-mode sampling are of comparable efficiency, while the steps taken when sampling components are exceedingly small. When the cutoff is large enough so that the tail of the spectral function is represented on the grid (bottom) the steps taken when sampling modes also become extremely small. Only blocked-mode sampling always efficiently samples the space of models $\mathbf{f}$.

Figure 4

Optical conductivity $\sigma \left(\omega \right)$ obtained by analytic continuation to uniform grids ${\omega }_n=(n-1/2)\mathrm{\Delta }\omega$ for $n\in \{1,...,N\}$ with fixed grid spacing $\mathrm{\Delta }\omega =0.25$ and increasing number of grid points $N=32,64,128$, and 256, corresponding to a cutoff ${\omega }_{\text{max}}\approx 8,16,32$, and 64. For comparison, the dashed line shows the optical conductivity from which the imaginary-frequency data for the analytic continuation was calculated. Even though all functions are essentially zero for $\omega \gtrsim 8$, the result depends very strongly on the length of the grid: the result of the analytic continuation develops spurious peaks that get sharper with increasing cutoff.

Figure 5

Grid mapping from uniform grid on the interval [0,1] to a nonuniform grid on the half infinite interval $[0,\infty )$. Grid points are indicated as dots, limits of intervals by bars.

Figure 6

Optical conductivity $\sigma \left(\omega \right)$ obtained by analytic continuation to different grids of $N=64$ points. The uniform grid has a spacing $\mathrm{\Delta }\omega =0.125$, corresponding to a cutoff ${\omega }_{\max }\approx 8$. The width parameter of the Gaussian grid was chosen $\alpha =4$, for the exponential $\beta =3$, and for the Lorentzian $\gamma =2.5$. The dashed line shows the exact result. Removing the cutoff by going to a nonuniform grid improves the result significantly, and the results depend much less on the chosen width parameter than on the cutoff. Still, the average spectra obtained for different grid densities differ by more than their error bars. This grid dependence becomes somewhat stronger for larger noise in the data [(top) ${\sigma }_{\mathrm{\Pi }}=0.001$ and (bottom) 0.01].

Figure 7

Simulating one finite grid on another. For the same problem as in the upper panel of Fig. 6, we compare the result for an exponential grid with $\beta =3$ and $N=128$ points with the simulation on the same exponential but with $\tilde{N}=64$ (top) and vice versa (bottom). The results agree within error bars. In addition, runs on the $\tilde{N}$ grid are shown, where the noise in the data is scaled by $N/\tilde{N}$. We do not plot the large symbols distinguishing the curves as they would obscure the near perfect agreement.

Figure 8

Dependence of the average spectrum on the number $N$ of grid points for the same problem as in Fig. 7 on a Lorentzian with $\gamma =2.5$. For larger $N,$ the average spectra get worse. The reason becomes clear from the histograms in the bottom panel, showing the contribution of spectra with a given fit $\chi$ to ${\tilde{\mathbf{f}}}_{\mathrm{ASM}}$: with increasing $N$ the histogram moves to the right, i.e., worse fits.

Figure 9

Determining a reasonable default model by fitting the grid density $\rho \left(x\right)$ (Gaussian of standard deviation $\alpha$, exponential with rate $\beta$, or Lorentzian of width $\gamma$, uniform with cutoff ${\omega }_c)$ to the data (dotted line) for the same problem as in Fig. 8. The inset shows $\sum {\left(\mathrm{\Pi }\left({\omega }_m\right)-{\mathrm{\Pi }}_{\text{def}}\left({\omega }_m\right)\right)}^2$ as a function of the grid parameter, highlighting the importance of choosing a reasonable default model.

Figure 10

Dependence of the average spectrum on the number $N$ of grid points for the same problem as in Fig. 8 on an optimized Gaussian grid with $\alpha =2.4$. The average spectra are largely independent of $N$ and the histograms show a consistently good fit.