Construct Fermionic Kernel Matrix Specified Tau Grid

[cohensbw]

Hi,

I was wondering if you could answer a question for me!

For the analytic continuation problem

$$G(\tau) = \int_{-\infty}^{+\infty} d\omega \ K(\tau,\omega) A(\omega),$$

the Fermionic kernel is given by $K(\tau, \omega) = \frac{e^{-\tau\omega}}{1 + e^{-\beta\omega}}$.
As you know, its convenient to recast this as the matrix equation

$$G(\tau_l) = K(\tau_l, \omega_i) \cdot A(\omega_i).$$

Now assume the imaginary-time grid needs to be evenly spaced in imaginary-time, such that $\tau_l = \Delta\tau \cdot l$ for $l \in [0, L_\tau)$, where $\beta = \Delta\tau \cdot L_\tau$. Can your package help me construct an optimal $K(\tau_l, \omega_i)$ matrix given a energy window $\omega_\min < \omega_i < \omega_\max$ and target precision $\epsilon$? Assume I can evaluate $A(\omega)$ for any frequency grid $\omega_i$ I want. If so, can you show me how to do so?

Thanks for you help,
Ben

[mwallerb]

Hi Ben,

Thanks for reaching out!

I am not exactly sure what exactly it is you are trying to achieve. Are you trying to solve the analytic continuation problem and are you looking for some "optimal" grid in real frequency to help you with that?

If so, let me first note that there is no "optimal" frequency grid in real frequency, that depends on what it is you're trying to achieve. There is the IR real-frequency grid ${\omega_i}$, that you can get as follows:

import sparse_ir, numpy as np
wmax = 10
beta = 30
basis = sparse_ir.FiniteTempBasis('F', beta, wmax)
w = basis.default_omega_sampling_points()

Then, assuming some linear spacing, you can get build the kernel matrix from the SVD of the kernel like this:

tau = np.linspace(0, beta, 6001)
K = basis.u(tau).T * basis.s @ basis.v(w)

But note that this real-frequency grid is in a sense only useful for the imaginary time/frequency side. It tells you that for any spectral function $A(\omega)$, there exists coefficients $g_i$ such that you can approximate the imaginary-frequency Green's function as a linear combination:

$$ G(\tau) \approx \sum_i g_i K(\tau, \omega_i) $$

to the given accuracy target (This is known as the DLR.) But it does not mean that the actual spectral can be approximated as a sum of these poles.

The strategy (see section "Numerical Analytical Continuation" in the tutorial) that works best for analytic continuation is to use the IR basis as a dimensional reduction step before you try to fit the spectral function. To do this, you first create a TauSampling object with your linear discretization as sampling points and fit your Green's function:

tau_sampling = sparse_ir.TauSampling(sampling_points=t)
gl = tau_sampling.fit(gtau)

Then, you have only ~ 100 inputs into your fitting routine. As a second step, you then set up your fitting kernel and your regularization scheme and do the actual fit.

Does that help?

---

PS:
Let me finally note that equally spaced points in imaginary time are dangerous: let me call $\max(|\omega_{\rm min}|, |\omega_{\rm max}|)$ in your notation $\omega_{\rm max}$ to be consistent with this package. Then it is relatively easy to see that you need:

$$ L_\tau \gtrsim \frac12 \beta\omega_{\rm max}\ \log(\epsilon^{-1})$$

points to even be able to "see" spectral features at the boundaries of your energy window. This blows up quickly, in particular in ab initio settings or at low temperatures. This is why we usually use specially crafted "sampling points" in imaginary time (see tutorial).

[cohensbw]

Thank you, I think your comments may end up being helpful! And you are totally correct, evenly spaced imaginary-time points is sub-optimal but I am stuck with them. I am working on trying to analytically continue imaginary-time Green's function data generated by a determinant quantum Monte Carlo (DQMC) simulation, which unavoidably reports the imaginary-time Green's function data on a regular imaginary-time grid. Specifically, I am trying to represent the spectral function as a sum of simple poles in the lower-half complex plane i.e. $A(\omega) = -\frac{1}{\pi}\text{Im}\left[ G(z \rightarrow \omega) \right]$ and

$$G(z) = \sum_p \frac{a_p + \rm{i}b_p}{z - (\epsilon_p - \rm{i} \Gamma_p)},$$

where it is assumed that $\Gamma_p > 0$ so that the pole is in the lower-half plane. If $b_p = 0$ then that term in the sum corresponds to a lorentzian, which has very fat tails. I need to then calculate $G(\tau)$ given this definition for the spectral function, and these relatively fat Lorentzian tails make it sort of trick to do well. I am currently doing Gauss-Hermite quadrature by doing a variable transform to express the integral in terms of the quantiles of a lorentzian given by assuming $b_p = 0$. However, this requires evaluating the integral for each pole $p$ separately instead of all together. I was hoping to figure out how to intelligently construct a $K(\tau_l,\omega_i)$ matrix to address this problem. Being able to express the forward process of $A(\omega) \mapsto G(\tau)$ as a simple equation like

$$G(\tau_l) = K(\tau_l, \omega_i) \cdot A(\omega_i).$$

is also important because I need to be able to back-propagate through the $A(\omega) \mapsto G(\tau)$ process. Anyways, thanks for responding to my questions, and I will try playing around with your package to see if I can get something that can replace what I am currently doing.

[mwallerb]

Ah, that makes sense!

If you're stuck with equispaced points, do have a look at work coming out of Emanuel Gull's (@egull) group on applying pole fitting algorithms like ESPRIT to analytic continuation: basically, you assume your Green's function is generated by a weakly stationary $AR(p)$ process, and you're trying to figure out its parameters, which you then relate to $p$ poles in the lower complex plane. It builds on half a century of signal theory and is quite robust. The cost is unfortunately a bit bad: $O(L^3)$ if you do it naively or $O(L^2 p)$ if you're smart about how to do the SVD.

For the reverse problem, i.e., computing $G(\tau)$ given a set of $a_p, b_p, \epsilon_p, \Gamma_p$ there is a simple trick:

1. construct the appropriate IR basis, and a Matsubara sampling object, which gives you a set of sampling points { $i\omega_n$ } and a fitting matrix.
2. Evaluate $\hat G(z = iw_n)$ by your equation (optionally speed it up using FMM if you have lots of poles) (cost $O(N+p)$, where $N$ is the size of the IR basis)
3. Fit the IR coefficients $g_l = \hat F \hat G$ using the .fit() routine. (cost $O(N^2)$ )
4. Profit: You can now evaluate the IR coefficients anywhere in time or frequency. You can also autograd through the whole procedure, since everything is linear.

By using these 4 steps, you can even turn this into an analytic continuation algorithm: you first fit the IR coefficients of your QMC result in time, expand it into the sampling frequencies, and then your $\chi^2$ term is the difference between that and the pole form, with the covariance given by the fitting matrix:

$$ r_n({\bf a, z}) := G(i\omega_n) - \sum_p \frac{a_p}{i\omega_n - z_p} $$

$$ \chi^2({\bf a, z}) = -\frac12 \bf r^\dagger F^\dagger F r $$

[cohensbw]

Sorry for not responding sooner, but I finally got around to trying to follow your advice, but I wanted to make sure I understood one detail correctly. Consider evaluating $G(z=i\omega_n)$ for some $i\omega_n$ value given that $G(z) = \frac{a+ib}{z-(\epsilon-i\Gamma)}$. Then I would want to transform $G(i\omega_n) \mapsto G(\tau)$. This would require doing calculating a fit using the fit() function then? So each time I were to change $(a, b, \epsilon, \Gamma)$ I would need to refit the IR coefficients? There is not static matrix $A(\tau_l, i\omega_n)$ I could construct once at the "beginning" to evaluate $G(\tau_l) = A(\tau_l, i\omega_n) G(i\omega_n)$ that I could then reuse, even as $(a, b, \epsilon, \Gamma)$ are changed? Thank you again for humoring! I have also been spending some time reading your papers to try and get a better grasp of how the IR decomposition encodes the information, as I am new to trying to use this type of representation.

[mwallerb]

You can also extract the matrices mediating the transformation, and then combine them.

basis = ...
s_tau = sparse_ir.MatsubaraSampling(basis)
s_iw = sparse_ir.TauSampling(basis)
A = np.linalg.lstsq(np.transpose(s_mat.matrix), np.transpose(s_tau.matrix))[0].T