Numerical analytic continuation#
One of the main problems one faces when working with imaginary (Euclidean) time is inferring the
real-time spectra
where again
Using the IR basis expansion, we can “invert” above equation to arrive at:
where
import numpy as np
import matplotlib.pyplot as pl
import sparse_ir
beta = 40.0
wmax = 2.0
basis = sparse_ir.FiniteTempBasis('F', beta, wmax, eps=2e-8)
pl.semilogy(basis.s / basis.s[0], '-+')
pl.title(r'singular values $s_\ell/s_0$ of $K(\tau, \omega)$')
pl.xlabel(r'index $\ell$');

Least squares form#
In order to make meaningful progress, let us reformulate the analytic continutation as a least squares problem:
To simplify and speed up this equation, we want to use the IR basis form of the kernel. For this, remember that for any imaginary-time propagator:
where the error term
allowing us to truncate our analytic continuation problem to (cf. Jarrell and Gubernatis, 1996):
This already is an improvement over many analytical continuation algorithms, as it maximally compresses the observed imaginary-time data
def semicirc_dos(w):
return 2/np.pi * np.sqrt((np.abs(w) < wmax) * (1 - np.square(w/wmax)))
def insulator_dos(w):
return semicirc_dos(8*w/wmax - 4) + semicirc_dos(8*w/wmax + 4)
# For testing, compute exact coefficients g_l for two models
rho1_l = basis.v.overlap(semicirc_dos)
rho2_l = basis.v.overlap(insulator_dos)
g1_l = -basis.s * rho1_l
g2_l = -basis.s * rho2_l
# Put some numerical noise on both of them (30% of basis accuracy)
rng = np.random.RandomState(4711)
noise = 0.3 * basis.s[-1] / basis.s[0]
g1_l_noisy = g1_l + rng.normal(0, noise, basis.size) * np.linalg.norm(g1_l)
g2_l_noisy = g2_l + rng.normal(0, noise, basis.size) * np.linalg.norm(g2_l)
Truncated-SVD regularization#
However, the problem is still ill-posed. As a first attempt to cure it, let us turn to truncated-SVD regularization (Hansen, 1987), by demanding that the spectral function
where
The choice of
# Analytic continuation made (perhaps too) easy
rho1_l_noisy = g1_l_noisy / -basis.s
rho2_l_noisy = g2_l_noisy / -basis.s
w_plot = np.linspace(-wmax, wmax, 1001)
Vmat = basis.v(w_plot).T
Lprime1 = basis.size // 2
def _plot_one(subplot, dos, rho_l, name):
pl.subplot(subplot)
pl.plot(w_plot, dos(w_plot), ":k", label="true")
pl.plot(w_plot, Vmat[:, :Lprime1] @ rho_l[:Lprime1],
label=f"reconstructed ($L' = {Lprime1}$)")
pl.plot(w_plot, Vmat @ rho_l, lw=1,
label=f"reconstructed ($L' = L = {basis.size}$)")
pl.xlabel(r"$\omega$")
pl.title(name)
pl.xlim(-1.02 * wmax, 1.02 * wmax)
pl.ylim(-.1, 1)
_plot_one(121, semicirc_dos, rho1_l_noisy, r"semi-elliptic DOS $\rho(\omega)$")
_plot_one(122, insulator_dos, rho2_l_noisy, r"insulating DOS $\rho(\omega)$")
pl.legend()
pl.gca().set_yticklabels([])
pl.tight_layout(pad=.1, w_pad=.1, h_pad=.1)

Regularization#
Above spectra are, in a sense, the best reconstructions we can achieve without including
any more a priori information about
the spectrum must be non-negative,
, for one orbital, and positive semi-definite, , in general,the spectrum must be a density:
,one may assume that the spectrum may be “sensible”, i.e., not deviate too much from a default model
(MAXENT) or not be too complex in structure (SpM/SOM).
These constraints are often encoded into the least squares problem (1) by restricting the space of valid solutions
All of these constraints act as regularizers.
As a simple example, let us consider Ridge regression in the above problem. We again expand the spectral function in the right singular functions with
where
# Analytic continuation made (perhaps too) easy
alpha = 100 * noise
invsl_reg = -basis.s / (np.square(basis.s) + np.square(alpha))
rho1_l_reg = invsl_reg * g1_l_noisy
rho2_l_reg = invsl_reg * g2_l_noisy
def _plot_one(subplot, dos, rho_l, rho_l_reg, name):
pl.subplot(subplot)
pl.plot(w_plot, dos(w_plot), ":k", label="true")
pl.plot(w_plot, Vmat @ rho_l, lw=1, label=f"t-SVD with $L'=L$")
pl.plot(w_plot, Vmat @ rho_l_reg, label=f"Ridge regression")
pl.xlabel(r"$\omega$")
pl.title(name)
pl.xlim(-1.02 * wmax, 1.02 * wmax)
pl.ylim(-.1, 1)
_plot_one(121, semicirc_dos, rho1_l_noisy, rho1_l_reg, r"semi-elliptic DOS $\rho(\omega)$")
_plot_one(122, insulator_dos, rho2_l_noisy, rho2_l_reg, r"insulating DOS $\rho(\omega)$")
pl.legend()
pl.gca().set_yticklabels([])
pl.tight_layout(pad=.1, w_pad=.1, h_pad=.1)

Real-axis basis#
One problem we are facing when solving the regularized least squares problem (2) is that the regularization might “force” values of
Let us illustrate the decay of
semielliptic (left), where the
decay roughly as , and are thus not compact.discrete set of peaks:
(right), where does not decay at all, signalling the fact that a delta-peak cannot be represented by any finite (or even infinite) expansion in the basis.
dos3 = np.array([-0.6, -0.1, 0.1, 0.6]) * wmax
rho3_l = basis.v(dos3).sum(1)
g3_l = -basis.s * rho3_l
def _plot_one(subplot, g_l, rho_l, title):
pl.subplot(subplot)
n = np.arange(0, g_l.size, 2)
pl.semilogy(n, np.abs(g_l[::2]/g_l[0]), ':+b', label=r'$|G_\ell/G_0|$')
pl.semilogy(n, np.abs(rho_l[::2]/rho_l[0]), ':xr', label=r'$|\rho_\ell/\rho_0|$')
pl.title(title)
pl.xlabel('$\ell$')
pl.ylim(1e-5, 2)
_plot_one(121, g1_l, rho1_l, r'semielliptic $\rho(\omega)$')
pl.legend()
_plot_one(122, g3_l, rho3_l, r'discrete $\rho(\omega)$')
pl.gca().set_yticklabels([])
pl.tight_layout(pad=.1, w_pad=.1, h_pad=.1)

Thus, we need another representation for the real-frequency axis. The simplest one is to choose a
grid
where
It is useful to use some probability distribution as
where
import scipy.stats as sp_stats
f = sp_stats.cauchy(scale=0.1 * np.pi / beta).pdf
pl.plot(w_plot, f(w_plot))
pl.title("Cauchy distribution");

Using this discretization, we finally arrive at a form of the analytic continuation problem suitable for a optimizer:
where
w = np.linspace(-wmax, wmax, 21)
K = -basis.s[:, None] * np.array(
[basis.v.overlap(lambda w: f(w - wi)) for wi in w]).T
Next, we will examine different regularization techniques that build on the concepts in this section…