libdlr: Efficient imaginary time calculations using the discrete Lehmann representation

We introduce libdlr, a library implementing the recently introduced discrete Lehmann representation (DLR) of imaginary time Green's functions. The DLR basis consists of a collection of exponentials chosen by the interpolative decomposition to ensure stable and efficient recovery of Green's functions from imaginary time or Matsbuara frequency samples. The library provides subroutines to build the DLR basis and grids, and to carry out various standard operations. The simplicity of the DLR makes it straightforward to incorporate into existing codes as a replacement for less efficient representations of imaginary time Green's functions, and libdlr is intended to facilitate this process. libdlr is written in Fortran, provides a C header interface, and contains a Python module pydlr. We also introduce a stand-alone Julia implementation, Lehmann.jl.


Introduction
Imaginary time Green's functions are a fundamental component of quantum many-body calculations at finite temperature. Since the cost of many algorithms scales with the number of imaginary time degrees of freedom, there has been significant recent interest in developing efficient methods of discretizing the imaginary time domain.
The most common approach used in scientific codes is to discretize imaginary time Green's functions on a uniform grid, or equivalently by a Fourier series. Although it is simple to use, this representation requires O (Λ/ ) degrees of freedom, where is the desired accuracy and Λ = βω max is a dimensionless energy cutoff depending on the inverse temperature β and a real frequency cutoff ω max . This scaling makes calculations at low temperature and high accuracy prohibitively expensive. Representations of Green's functions by orthogonal polynomials require O √ Λ log(1/ ) degrees of freedom, largely addressing the accuracy issue but remaining suboptimal in the low temperature limit [1][2][3]. This has motivated the development of optimized basis sets in which to expand imaginary time Green's functions; namely the intermediate representation  In particular, the DLR makes use of a low rank decomposition of K of the form for smooth functions π l . Substitution of such a decomposition into (4) yields the DLR (3), demonstrating its existence in principle. As we will see, in practice, we do not obtain the DLR expansion of a Green's function in this manner, since the spectral function ρ is typically not known. The DLR basis {K(τ, ω l )} r l=1 is characterized by the DLR frequencies {ω l } r l=1 , which depend only on the high energy cutoff Λ and the desired accuracy , both user-defined parameters. The DLR frequencies are obtained by a two-step procedure. First, the kernel K(τ, ω) is discretized on a high-order, adaptive fine grid (τ f i , ω f j ) M,N i,j=1 ⊂ [0, 1] × [−Λ, Λ], constructed to efficiently resolve small scale features. This yields an accurate discretization A ij = K(τ f i , ω f j ) of the integral operator in (4). Then, a rank revealing column pivoted QR algorithm is applied to A, yielding a minimal collection of columns sufficient to span the full column space of A to accuracy . Since the matrix A provides an accurate discretization of the kernel K(τ, ω), the selected columns, which correspond to selected frequencies {ω l } r l=1 , define an accurate approximate basis {K(τ, ω l )} r l=1 of the column space of the Lehmann representation integral operator. The number r of DLR basis functions is shown as a function of Λ for a few choices of in Figure 1, and an example of the selected DLR frequencies for a given choice of Λ and is shown in the first panel of Figure 2. The number of basis functions, called the DLR rank, is observed to scale as r = O (log (Λ) log (1/ )).
The DLR coefficients g l of a given Green's function G can be recovered directly from samples {G(τ s k )} n k=1 , for some sampling nodes τ s k . This can be done by solving a linear system for k = 1, . . . , n. We consider two possible scenarios: (i) Applying the rank revealing row pivoted QR algorithm to the matrix K(τ f i , ω l ) yields a set of interpolation nodes {τ k } r k=1 , called the DLR imaginary time nodes, such that the coefficients can be recovered from the samples {G(τ k )} r k=1 . Thus, if G(τ ) can be evaluated at arbitrary imaginary time points, then we can take n = r and τ s k = τ k in (6) and solve the resulting small square linear system to obtain an interpolant of G(τ ). An example of the DLR imaginary time nodes is shown in the second panel of Figure 2. (ii) If samples of G are given at a collection of n > r scattered or uniform grid points τ s k , then (6) is an overdetermined system and can be solved by ordinary least squares fitting.
Since the DLR basis functions are given explicitly, a DLR expansion can be analytically transformed to the Matsubara frequency domain. We have with the Matsubara frequencies given by iν n = i(2n + 1)π/β for fermionic Green's functions i(2n)π/β for bosonic Green's functions (8) in unscaled coordinates. In scaled coordinates, we simply set iν n ← iβν n . The DLR expansion in the Matsubara frequency domain is then As in the imaginary time domain, there is a set of Matsubara frequencies {iν n k } r k=1 , the DLR Matsubara frequency nodes, such that the coefficients can be recovered from samples {G(iν n k )} r k=1 by interpolation; or, one can recover the expansion by least squares fitting. An example of the DLR Matsubara frequency nodes for a given choice of Λ and is shown in the third panel of Figure 2.
For further details on the algorithm used to obtain the DLR frequencies, imaginary time nodes, and Matsubara frequencies, as well as a detailed analysis of the accuracy and stability of the method, we again refer the reader to [8].

Features and layout of the library
We list the main capabilities of libdlr: • Given a choice of Λ and , obtain the DLR frequencies {ω l } r l=1 characterizing the DLR basis • Obtain the DLR imaginary time and Matsubara frequency grids The computational cost of building the DLR basis and grids is typically negligible even for large values of Λ, and subsequent operations involve standard numerical linear algebra procedures with small matrices. libdlr is implemented in Fortran, provides a C header interface, and includes a Python module, pydlr. More information on pydlr is contained in Appendix A. The documentation for the library [12] includes installation intructions, API information, and several examples. The source code is available as a Git repository [13]. The subfolder libdlr/test contains example programs which are run after compilation as tests, and the subfolder libdlr/demo contains additional example programs. An example C program using the C header interface may be found in libdlr/test/ha it.c.
We have made the Julia package Lehmann.jl, which implements functionality similar to libdlr, available in a separate Git repository [14]. More information on Lehmann.jl is contained in Appendix B, and in the documentation [15].

Examples of usage
We describe two examples illustrating the functionality of libdlr. First, we demonstrate recovery of the DLR coefficients of a Green's function both from its values on the DLR imaginary time grid and from noisy data on a uniform grid. Second, we demonstrate the process of solving the Dyson equation self-consistently for the Sachdev-Ye-Kitaev model [16][17][18]. All examples are implemented in Fortran; Appendices A and B contain analogous examples implemented using pydlr and Lehmann.jl, respectively.

Obtaining a DLR expansion from Green's function samples
We consider the Green's function given by the Lehmann representation (1) with spectral density Here, θ is the Heaviside function. In practical calculations, the spectral density is usually not known, so one must recover the DLR from samples of the Green's function itself. We use a Green's function with known spectral density for illustrative purposes, since we can evaluate it with high accuracy by numerical integration of (1) and thereby test our results.
Since supp ρ ⊂ [−1, 1], we can take the frequency support cutoff ω max = 1. Then Λ = βω max ≥ β is a sufficient high energy cutoff. In practice, ω max can typically be estimated on physical grounds, giving an estimate of Λ. Calculations can then be converged with respect to Λ to ensure accuracy. The error tolerance should be chosen based on the desired accuracy, in order to obtain the smallest possible number of basis functions. In this example, we take β = 1000, and fix Λ = 1000.
We note also that libdlr routines work by default with matrix-valued Green's functions G ij , where i and j are typically orbital indices. However, in all examples presented here, we work with scalar-valued Green's functions, and therefore set the parameter n determining the number of orbital indices to 1.

Recovery of DLR from imaginary time grid values
We first consider recovery of the DLR coefficients g l from samples of G(τ ) on the DLR imaginary time nodes τ k . This is the first scenario mentioned in Section 2, and generates an interpolant at the nodes τ k . Figure 3 shows a condensed version of a Fortran code implementing the example using libdlr. In this and all other sample Fortran codes, we do not show variable allocations or steps which do not involve libdlr subroutines, as indicated in the code comments. A complete Fortran code demonstrating this example can be found in the file libdlr/demo/sc it.f90. The file libdlr/demo/sc mf.f90 contains a demonstration of the process of obtaining a DLR from values of a Green's function on the DLR Matsubara frequency nodes.
We first set Λ and , and then build the DLR basis by obtaining the DLR real frequencies ω l , stored in the array dlrrf, using the subroutine dlr it build. This subroutine also produces the r DLR imaginary time nodes τ k , which are stored in the array dlrit.
Next, we assume that the values G(τ k ) have been obtained by some external procedure, and stored in the array g. In this case, we used numerical integration with the known spectral function to obtain the samples. To obtain the DLR coefficients g l , one must solve the linear system (6), with n = r. The subroutine dlr it2cf init initializes this procedure by computing the LU factorization of the system matrix. The LU factors and pivots are stored in the arrays it2cf and it2cfp, respectively. The linear solve is then carried out by the subroutine dlr it2cf, which returns the DLR coefficients g l in the array gc.
We can then evaluate the DLR expansion on output grids in imaginary time and Matsubara frequency. The subroutine eqpts rel generates a uniform grid of imaginary time points in the relative format employed by the library (see Appendix C for an explanation of the relative format) and dlr it eval evaluates the expansion. The subroutine dlr mf eval evaluates the DLR expansion in the Matsubara frequency domain.

Recovery of DLR from noisy data
We next consider DLR fitting from noisy data, the second scenario mentioned in Section 2. A condensed sample code is given in Figure 4, and a complete example can be found in the file libdlr/demo/sc it fit.f90. We assume noisy samples G(τ s k ) have been obtained by an external procedure and stored in the array g. The subroutine dlr it fit fits a DLR expansion to the data by solving the overdetermined system (6). The resulting DLR coefficients, stored in the array gc, can be used to evaluate the DLR expansion as in the previous example.  Figure 5a, we plot the Green's function G(τ ) with spectral density (9) and β = 1000, along with noisy data obtained on a uniform grid of n = 2500  points by adding uniform random numbers of magnitude η = 10 −2 to accurate values computed by numerical integration. We see that the DLR fit G DLR to this data with = 10 −2 agrees well with G(τ ). Pointwise errors for this example, and examples with different noise levels η, are given in Figure 5c. In all cases we take = η, and observe that the DLR fitting process is stable; the fitting process does not introduce a significant error above the magnitude of the noise. Pointwise errors for the DLR obtained from numerically exact samples on the DLR imaginary time grid are shown in Figure 5b, and we see that the error is well controlled by .

Solving the Sachdev-Ye-Kitaev model
The Dyson equation for a given self-energy Σ can be written in the Matsubara frequency domain as Here G 0 (iν n ) = (iν n + µ) −1 is the free particle Green's function, with chemical potential µ. This equation can be solved for G by pointwise inversion on the r DLR Matsubara frequency nodes. In typical applications, Σ is a function of G, Σ = Σ [G], which is most easily evaluated in the imaginary time domain. In this case, (10) becomes nonlinear, and must be solved by self-consistent iteration. The standard method is to solve the Dyson equation using (10) for a given iterate of Σ, transform the solution G to the imaginary time domain to evaluate a new iterate of Σ, and then transform Σ back to the Matsubara frequency domain to solve (10) for the next iteration. The standard implementation of this approach represents G and Σ on fine equispaced grids in imaginary time, and uses the fast Fourier transform to move between the imaginary time and Matsubara frequency domains. The process can be carried out significantly more efficiently using the DLR.
As an example, we consider the Dyson equation given by (10)  ! Recovery of DLR from noisy data on uniform grid (assume data 10 ! given at m equispaced grid points t stored in array g) 11 12 call dlr_it_fit(r,n,dlrrf,m,t,g,gc) Figure 4: libdlr Fortran code to obtain DLR from noisy data on a uniform grid.
Here J is a coupling constant. We consider this example with β = 1000, µ = 0, and J = 1, and set Λ = 5000. We note that G is a fermionic Green's function, so we use a fermionic Matsubara frequency grid. Figure 6 gives sample code for a solver which uses a weighted fixed point iteration to handle the nonlinearity: Here, m refers to the iterate, and w is a weighting parameter which can be selected to improve convergence; we use w = 0.3. A complete code demonstrating this example can be found in the file libdlr/demo/syk mf.f90. As in the previous example, the code begins by obtaining the DLR frequencies and imaginary time grid. Here, the DLR Matsubara frequency grid is also built, using the subroutine dlr mf, with a Matsubara frequency cutoff nmax. This cutoff should in principle be taken larger and larger until the selected DLR Matsubara frequencies no longer change, but in practice we find that setting it to Λ is usually sufficient (see the discussion in [8, Sec. IV D]).
The next several lines define physical problem parameters, and parameters for the weighted fixed point iteration. Then, several initialization routines are called which prepare transformations between the imaginary time and Matsubara frequency grid representations of the Green's function, and the DLR coefficients. These transformations are used to convert between imaginary time and Matsubara frequency representations in the self-consistent iteration. The subroutine dlr it2itr init prepares a transformation between a Green's function G(τ ) on the imaginary time grid and its reflection G(β − τ ) on the same grid, which is needed to evaluate the SYK self-energy. if (maxval(abs(gprev-g))<fptol) then ! Check self-consistency  call dlr_it2itr(r,n,it2itr,g,gr) ! Get G(beta-tau) 7 sig = j*j*g*g*gr ! Get Sigma(tau) 8 9 end subroutine sigfun  The free particle Green's function G 0 appearing in (10) is then evaluated on the Matsubara frequency grid, since it appears in the Dyson equation, and on the imaginary time grid, to serve as an initial guess in the iteration. In the weighted fixed point iteration, the self-energy is evaluated on the imaginary time grid, and then transformed to the Matsubara frequency domain, where the Dyson equation is solved. Then the result is transformed back to the imaginary time grid, where the self-consistency of the solver is checked. Figure 7 shows the subroutine used to evaluate the SYK self-energy in imaginary time. Once self-consistency is reached, the solution is returned on the imaginary time grid. It can be expanded in a DLR and evaluated in imaginary time or Matsubara frequency as in the previous example.
Plots of the solution and of the DLR nodes are given in Figure 8, both in imaginary time and Matsubara frequency. For Λ = 5000 and = 10 −10 , there are r = 66 DLR nodes discretizing each domain.
We remark that the Dyson equation (10) may be written in the imaginary time domain as In [8], an efficient method is presented of solving (13) by discretizing the convolutions using the DLR. This approach may be advantageous in certain cases, and can also be implemented using libdlr. A demonstration for the SYK model is given in the file libdlr/demo/syk it.f90. We also note that the convergence properties of nonlinear iteration in the self-consistent solution of the Dyson equation are problem-dependent. For the SYK model (11), which belongs to a particular subclass of approximations [19], we find that the weighted fix point iteration (12) may fail to converge or converge to spurious local fixed points, particularly at extremely low temperatures. Many standard strategies can be used to improve convergence, including the careful selection of an initial guess by adiabatic variation of parameters like β or µ, mixing, and more sophisticated nonlinear iteration procedures. Two other methods have been observed to improve performance in certain difficult cases: the use of explicitly symmetrized DLR grids, and mild oversampling in the Matsubara frequency domain followed by least squares fitting to a DLR expansion. These are topics of our current research, and will be reported on in detail at a later date. Lastly, we have observed that solving (13) directly in the imaginary time domain, as mentioned in the previous paragraph, tends to lead to more robust convergence.

Conclusion
libdlr facilitates the efficient representation and manipulation of imaginary time Green's functions using the DLR. In this framework, working with imaginary time Green's functions is as simple as with standard discretizations, but involves many fewer degrees of freedom, and user-tuneable accuracy down to the level of machine precision.
As a result, we anticipate the use of the DLR as a basic working tool a variety of equilibrium applications, including continous-time quantum Monte Carlo [20], dynamical mean-field theory [21], self-consistent perturbation theory in both quantum chemistry (GF2) [22][23][24][25] and condensed matter physics [26], as well as Hedin's GW approximation [27,28], including vertex corrections [29,30]. In [31], the DLR was used to discretize imaginary time variables in equilibrium real time contour Green's functions. In nonequilibrium calculations involving two-time Green's functions, it can replace the equispaced imaginary time grids currently in use [32,33], and further improve the efficiency of algorithms making use of low rank compression to reduce the cost of time propagation [34].

C Stable kernel evaluation and relative format for imaginary time points
In libdlr, we work in the dimensionless variables defined at the beginning of Section 2, in which the kernel K(τ, ω) is given by for τ ∈ [0, 1] and ω ∈ (−∞, ∞). When ω ≥ 0, this formula is numerically stable. If ω 0, it can overflow, and we instead use the mathematically equivalent formula This formula, however, leads to another problem: when |ω| is large, catastrophic cancellation in the calculation of 1 − τ for τ near 1 leads to a loss of accuracy in floating point arithmetic.
Let us illustrate the advantage with a concrete example. For simplicity, we assume we are working in three-digit arithmetic. We first consider the evaluation of K(τ, ω) for τ = 0.501e−3 and ω = 1000. No problem arises in this case; we use the formula (14), and obtain the correct value However, if we instead want to calculate K(τ, ω) for τ = 1 − 0.501e−3 = 0.999499 and ω = −1000, the situation is different. In three-digit arithmetic, using the absolute format, we must round to τ = 0.999. Then we find 1 − τ = 0.001 =⇒ (1 − τ )ω = −1, and using (15) directly gives K(τ, ω) = e −1 1 + e −1000 , which is far from the correct value. If we instead use the relative format, τ is stored as τ * = −0.501e−3, and using (16) gives precisely the correct value. In practice, with double precision arithmetic, using the absolute format only leads to a significant loss of accuracy for very large Λ and small . However, to enable calculations to high accuracy in extreme scenarios, all subroutines in libdlr take in imaginary time values in the relative format. Of course, maintaining full relative precision also requires external procedures, such as those used to evaluate Green's functions, to be similarly careful about this issue.
Many users are likely not operating in regimes in which it is important to maintain full relative precision for extreme parameter values. These users can simply ignore this discussion. However, they must still convert imaginary time values to the relative format before using them as inputs to libdlr subroutines (though, of course, imaginary time values which are converted from the absolute format to the relative format are only accurate to the original absolute precision). The subroutine abs2rel performs this conversion. Similarly, the subroutine rel2abs converts imaginary time values in the relative format used by libdlr to the ordinary absolute format.