API Documentation pydlr

pydlr

Discrete Lehman Representation (DLR) implementation using Numpy and Scipy.

class pydlr.pydlr.dlr(lamb, eps=1e-15, xi=-1, max_rank=500, nmax=None, verbose=False, python_impl=True)

Discrete Lehmann Representation (DLR) class.

Provides the DRL basis, transforms, and algorithms.

Parameters

• lamb (float) – DLR scale parameter \(\Lambda\).

• eps (float) – Set accuracy of the DLR representation.

• xi (sign, optional,) – Statistical sign \(\xi = \pm 1\) for bosons and fermions respectively.

• max_rank (int, optional) – Maximum rank of the DLR kernel decomposition. Default 500.

• nmax (int, optional) – Maxumum index of the Matsubara frequency grid. Default int(lamb).

• verbose (bool, optional) – Default False.

• python_impl (bool, optional) – Switch between the python and fortran library driver. Default True.

get_dlr_frequencies()

Get real frequency DLR grid \(\omega_x \in [-\Lambda, \Lambda]\)

Returns

w_x – Scaled real frequency points \(\omega_x\) for the DLR representation.

Return type

(n), ndarray

get_tau_over_beta()

Get grid in scaled imaginary time \(\tau/\beta \in [0, 1]\)

Returns

ttau_l – Scaled imaginary time points \(\tilde{\tau}_l = \tau_l / \beta\) for the DLR representation, \(\tilde{\tau}_l \in [0, 1]\).

Return type

(n), ndarray

get_tau(beta)

Get \(\tau\)-grid in imaginary time \(\tau \in [0, \beta]\).

Parameters

beta (float) – Inverse temperature \(\beta\)

Returns

tau_l – Imaginary time points \(\tau_l\) for the DLR representation, \(\tau_l \in [0, \beta]\).

Return type

(n), ndarray

dlr_from_tau(G_laa)

Transform the rank-3 array_like Green’s function G_laa from imaginary time to DLR space.

Parameters

G_laa ((n,m,m), array_like) – Green’s function in imaginary time with \(m \times m\) orbital indices.

Returns

G_xaa – Green’s function i DLR coefficient space with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

tau_from_dlr(G_xaa)

Transform the rank-3 array_like Green’s function G_xaa from DLR space to imaginary time.

Parameters

G_xaa ((n,m,m), array_like) – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

Returns

G_laa – Green’s function in imaginary time with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

tau_from_legendre(G_naa)

Transform the rank-3 array_like Green’s function G_naa from Legendre coefficient space to imaginary time.

Parameters

G_naa ((n,m,m), array_like) – Green’s function in Legendre coefficient space with \(m \times m\) orbital indices.

Returns

G_laa – Green’s function in imaginary time with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

eval_dlr_tau(G_xaa, tau_k, beta)

Evaluate the DLR coefficient Green’s function G_xaa at arbibrary points in imaginary time.

Parameters

• G_xaa ((n,m,m), array_like) – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

• tau_k ((k), array_like) – Imaginary time points \(\tau_k\) where to evaluate the Green’s function.

• beta (float) – Inverse temperature \(\beta\)

Returns

G_kaa – Green’s function at the imaginary time points \(\tau_k\) with \(m \times m\) orbital indices.

Return type

(k,m,m), ndarray

lstsq_dlr_from_tau(tau_i, G_iaa, beta)

Return DLR coefficients by least squares fit to values on arbitrary imaginary time grid.

Parameters

• tau_i ((i), array_like) – Imaginary time points \(\tau_i\) where the Green’s function is sampled.

• G_iaa ((i,m,m), array_like) – Green’s function in imaginary time space \(G(\tau_i)\) with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

Returns

G_xaa – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

Return type

(k,m,m), ndarray

get_matsubara_frequencies(beta)

Get Matsubara frequency grid.

Parameters

beta (float) – Inverse temperature \(\beta\)

Returns

w_q – Matsubara frequency points \(i\omega_q\) for the DLR representation.

Return type

(n), ndarray

dlr_from_matsubara(G_qaa, beta, xi=None)

Transform the rank-3 array_like Green’s function G_qaa from Matsbuara frequency to DLR space.

Parameters

• G_qaa ((n,m,m), array_like) – Green’s function in Matsubara frequency with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

G_xaa – Green’s function i DLR coefficient space with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

matsubara_from_dlr(G_xaa, beta, xi=None)

Transform the rank-3 array_like Green’s function G_xaa from DLR space to Matsbuara frequency.

Parameters

• G_xaa ((n,m,m), array_like) – Green’s function i DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

G_qaa – Green’s function in Matsubara frequency with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

eval_dlr_freq(G_xaa, z, beta, xi=None)

Evaluate the DLR coefficient Green’s function G_xaa at arbibrary points z in frequency space.

Parameters

• G_xaa ((n,m,m), array_like) – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

• z ((k), array_like) – Frequency points \(z_k\) where to evaluate the Green’s function.

• beta (float) – Inverse temperature \(\beta\)

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

G_zaa – Green’s function at the frequency points \(z_k\) with \(m \times m\) orbital indices.

Return type

(k,m,m), ndarray

lstsq_dlr_from_matsubara(w_q, G_qaa, beta)

Return DLR coefficients by least squares fit to values on arbitrary Matsubara frequency grid.

Parameters

• w_q ((q), array_like) – Imaginary frequency points \(i\omega_q\) where the Green’s function is sampled.

• G_qaa ((q,m,m), array_like) – Green’s function in imaginary frequency space \(G(i\omega_q)\) with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

Returns

G_xaa – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

Return type

(k,m,m), ndarray

quadratic_hamiltonian(G_xaa, beta, xi=None)

Get quadratic Hamiltonian contribution of physical Green’s function.

Parameters

• G_xaa ((k,m,m), ndarray) – Green’s function in DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

Returns

H_aa – Quadratic Hamiltonian contribution to the Green’s function.

Return type

(m,m), ndarray

convolution(A_xaa, B_xaa, beta, xi=None)

DLR convolution with \(\mathcal{O}(N^2)\) scaling. Author: Hugo U.R. Strand (2021)

Imaginary time convolution

\[C = A \ast B\]

reformulated in DLR coefficient space. The notation \(C = A \ast B\) is short hand for:

\[C_{ij}(\tau) = \sum_k \int_{0}^\beta d\bar{\tau} A_{ik}(\tau - \bar{\tau}) B_{kj}(\bar{\tau})\]

Parameters

• A_xaa ((n,m,m)) – Green’s function \(A\) in DLR coefficient space with \(m \times m\) orbital indices.

• B_xaa ((n,m,m)) – Green’s function \(B\) in DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

C_xaa – Green’s function \(C\) in DLR coefficient space with \(m \times m\) orbital indices, given by the convolution \(C = A \ast B\).

Return type

(n,m,m), ndarray

convolution_matrix(A_xaa, beta, xi=None)

DLR convolution matrix with \(\mathcal{O}(N^2)\) scaling. Author: Hugo U.R. Strand (2021)

The imaginary time convolution matrix \(M\) is given by

\[M = [A \ast]\]

i.e. the combination of the Green’s function \(A\) and the convolution operator \(\ast\).

Author: Hugo U.R. Strand (2021)

Parameters

• A_xaa ((n,m,m)) – Green’s function \(A\) in DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

M_xaxa – Convolution matrix \([A \ast]\) as a rank-4 tensor in DLR and orbital space.

Return type

(n,m,n,m), ndarray

free_greens_function_dlr(H_aa, beta, S_aa=None, xi=None)

Return the free Green’s function in DLR coefficent space.

The free Green’s function is the solution to the Dyson differential equation

\[(-\partial_\tau - H_{ij}) G_{jk}(\tau) = \delta(\tau)\]

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions, where the Dyson equation takes the form \((-S_{ij} \partial_\tau - H_{ij}) G_{jk}(\tau) = \delta(\tau)\). Default is S_aa = None.

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

G_xaa – Free Green’s function \(G\) in DLR coefficient space with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

Notes

The fastest algorithm for calculation of free Green’s functions is the free_greens_function_tau method.

free_greens_function_tau(H_aa, beta, S_aa=None, xi=None)

Return the free Green’s function in imaginary time.

The free Green’s function is the solution to the Dyson differential equation

\[(-\partial_\tau - H_{ij}) G_{jk}(\tau) = \delta(\tau)\]

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions, where the Dyson equation takes the form \((-S_{ij} \partial_\tau - H_{ij}) G_{jk}(\tau) = \delta(\tau)\). Default is S_aa = None.

• xi (int, optional) – Statistics sign, \(\xi = +1\) for bosons and \(\xi = -1\) for fermions.

Returns

G_laa – Free Green’s function \(G\) in imaginary time with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

free_greens_function_matsubara(H_aa, beta, S_aa=None)

Return the free Green’s function in Matsubara frequency.

The free Green’s function is the solution to the Dyson equation

\[G_{ij}(i \omega_n ) = \left[ i\omega_n - H_{ij} \right]^{-1}\]

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions, where the Dyson equation takes the form \(= G_{jk}(i\omega_n) = (-S_{ij} i\omega_n - H_{ij})^{-1}\). Default is S_aa = None.

Returns

G_qaa – Free Green’s function \(G\) in Matsubara frequency with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

Notes

The fastest algorithm for calculation of free Green’s functions is the free_greens_function_tau method.

dyson_matsubara(H_aa, Sigma_qaa, beta, S_aa=None)

Solve the Dyson equation in Matsubara frequency.

The Dyson equation gives the Green’s function as

\[G_{ij}(i \omega_n ) = \left[ i\omega_n - H_{ij} - \Sigma(i\omega_n) \right]^{-1}\]

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• Sigma_qaa ((n,m,m), ndarray) – Self-energy \(\Sigma\) in Matsubara frequency with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions. Default is S_aa = None.

Returns

G_qaa – Green’s function \(G\) in Matsubara frequency with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

Notes

The Matsubara frequency Dyson solver is the fastest Dyson solver, albeit not as accurate as the DLR solver dyson_dlr.

dyson_dlr(H_aa, Sigma_xaa, beta, S_aa=None, iterative=False, lomem=False, verbose=False, tol=1e-12)

Solve the Dyson equation in DLR coefficient space.

The Dyson equation gives the imaginary time Green’s function \(G_{ij}(\tau)\) as

\[(-\partial_\tau - H_{ij} - \Sigma_{ij} \ast \, ) G_{jk}(\tau) = \delta(\tau)\]

Using the free Green’s function \(g\) this can be rewritten as the integral equation

\[(1 - g \ast \Sigma \ast \, ) G = g\]

which here is solved directly in DLR coefficient space.

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• Sigma_xaa ((n,m,m), ndarray) – Self-energy \(\Sigma\) in DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions. Default is S_aa = None.

• iterative (bool, optional) – Default False

• lomem (bool, optional) – Default False

• tol (float, optional) – Tolerance for iterative GMRES solver.

Returns

G_xaa – Green’s function \(G\) in DLR coefficient space with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

Notes

This DLR space Dyson solver is the most accurate algorithm for solving the Dyson equation.

By default it uses Lapack’s direct solver on matrix problem in combined DLR and orbital space. This is fast and accurate for small problems. For larger problems the \(\mathcal{O}(N^2M^2)\) memory foot-print limits the performance.

Hence, for large problems the solver should be run with iterative=True and lomem=True and will then use GMRES to solve the linear system using an implicit matrix formulation. This formulation gives a memory foot print that is of the same order as storing a single Green’s function. However, the requested GMRES tolerance tol has to be carefully tuned.

dyson_dlr_integrodiff(H_aa, Sigma_xaa, beta, S_aa=None, xi=None)

Solve the Dyson equation in DLR coefficient space.

The Dyson equation gives the imaginary time Green’s function \(G_{ij}(\tau)\) as

\[(-\partial_\tau - H_{ij} - \Sigma_{ij} \ast \, ) G_{jk}(\tau) = \delta(\tau)\]

which here is solved directly in DLR coefficient space.

Parameters

• H_aa ((m,m), array_like) – Single-particle Hamiltonian matrix in \(m \times m\) orbital space.

• Sigma_xaa ((n,m,m), ndarray) – Self-energy \(\Sigma\) in DLR coefficient space with \(m \times m\) orbital indices.

• beta (float) – Inverse temperature \(\beta\)

• S_aa ((m,m), array_like, optional) – Overlap matrix for the generalized case of non-orthogonal basis functions. Default is S_aa = None.

Returns

G_xaa – Green’s function \(G\) in DLR coefficient space with \(m \times m\) orbital indices.

Return type

(n,m,m), ndarray

Notes

This DLR space Dyson solver should not be used in production, use dyson_dlr instead. The differential formulation gives a worse condition number, as compared to the integral equation formulation. The method is kept here for educational purposes.

class pydlr.pydlr.solver_wrapper(solver, verbose=False)

Wrapper for scipy.sparse.linalg iterative linearsystems solvers

kwargs_filter(d)

callback(x)

kernel

Imagniary time analytic continuation kernel routines.

pydlr.kernel.chebyshev_collocation_points_1st_kind(N)

Return the Chebyshev collocation points of the first kind.

The Chebyshev collocation points of the first kind are defined as

\[x_j = \cos \left( \pi \frac{2j + 1}{2N} \right)\]

where \(N\) is the order and \(j = 0, ..., N-1\).

Parameters

N (int) – Order of the collocation point set.

Returns

x_j – Array with the collocation points \(x_j\) sorted in increasing order.

Return type

ndarray

Examples

>>> from pydlr.kernel import chebyshev_collocation_points_1st_kind
>>> x_j = chebyshev_collocation_points_1st_kind(10)
array([-0.98768834, -0.89100652, -0.70710678, -0.4539905 , -0.15643447,
        0.15643447,  0.4539905 ,  0.70710678,  0.89100652,  0.98768834])

pydlr.kernel.chebyshev_barycentric_weights_1st_kind(N)

Return the Chebyshev barycentric interpolation weights of the first kind.

The barycentric interpolation weights are defined as

\[w_j = (-1)^{j} \sin \left( \pi \frac{2j + 1}{2N} \right)\]

where \(N\) is the order and \(j = 0, ..., N-1\).

Parameters

N (int) – Order of the collocation point set.

Returns

w_j – Array with the barycentric weights \(w_j\).

Return type

ndarray

Examples

>>> from pydlr.kernel import chebyshev_barycentric_weights_1st_kind
>>> w_j = chebyshev_barycentric_weights_1st_kind(10)
array([ 0.15643447, -0.4539905 ,  0.70710678, -0.89100652,  0.98768834,
       -0.98768834,  0.89100652, -0.70710678,  0.4539905 , -0.15643447])

pydlr.kernel.barycentric_chebyshev_interpolation(x, x_j, f_j, w_j)

Return the barycentric interpolation of a Chebyshev polynomial on arbitrary points.

The Barycentric interpolation formula has the form

\[f(x) = \left( \sum_{j=0}^{N-1} \frac{f_j w_j}{x - x_j} \right) \Bigg/ \left( \sum_{j=0}^{N-1} \frac{w_j}{x - x_j} \right)\]

Parameters

• x (ndarray) – points \(x\) to evaluate the Chebyshev polynomial at

• x_j (ndarray) – Chebyshev collocation points \(x_j\) of the first kind

• f_j (ndarray) – Values \(f_j\) of the polynomial at the collocation points, \(f_j = f(x_j)\)

• w_j (ndarray) – Chebyshev barycentric interpolation weights \(w_j\)

Returns

f_x – Values \(f(x)\) of the polynomial at the points \(x\)

Return type

ndarray

Examples

>>> import numpy as np
>>> f = lambda x : np.cos(4*np.pi * x)
>>> from pydlr.kernel import chebyshev_collocation_points_1st_kind
>>> from pydlr.kernel import chebyshev_barycentric_weights_1st_kind
>>> x_j = chebyshev_collocation_points_1st_kind(32)
>>> w_j = chebyshev_barycentric_weights_1st_kind(32)
>>> f_j = f(x_j)
>>> from pydlr.kernel import barycentric_chebyshev_interpolation
>>> x = np.linspace(-1, 1, num=1000)
>>> f_x = barycentric_chebyshev_interpolation(x, x_j, f_j, w_j)
>>> np.testing.assert_array_almost_equal(f_x, f(x))
>>> print(np.max(np.abs(f_x - f(x))) < 1e-9)
True

pydlr.kernel.fermi_function(E, beta)

Evaluate the Fermi distribution function at energy \(E\) and inverse temperature \(eta\).

The Fermi distribution function \(f_\beta(E)\) has the form

\[f_\beta(E) = \frac{1}{1 + e^{\beta E}}\]

the evaluation is stabilized using separate formulas for \(E \lessgtr 0\).

Parameters

• E (ndarray) – Energies

• beta (float) – Inverse temperature

Returns

f – The Fermi distribution function evaluated at \(E\)

Return type

ndarray

pydlr.kernel.kernel(tau, omega)

Evaluate the imaginary time and real frequency analytical continuation kernel \(K(\tau, \omega)\).

The Fermionic analytical continuation kernel has the form

\[K(\tau, \omega) = \frac{e^{-\omega \tau}}{1 + e^{\omega}}\]

in normalized imaginary time \(\tau \in [0, 1]\) and frequency \(\omega\).

The evaluation is stabilized using separate formulas for \(E \lessgtr 0\).

Parameters

• tau (ndarray) – Points in imaginary time \(\tau_j\)

• omega (ndarray) – Points in real frequency \(\omega_k\)

Returns

kernel – Kernel \(K_{jk}\) evaluated on the \(\tau_j\) and \(\omega_k\) grids, \(K_{jk} = K(\tau_j, \omega_k)\)

Return type

ndarray

pydlr.kernel.gridparams(lamb, order=24, lambda_scale=1.0)

Empirical discretization parameters of \(K(\tau, \omega)\) for given \(\Lambda\)

Parameters

• lamb (float) – Cutoff parameter \(\Lambda\)

• order (int, optional) – Order of the polynomial expansion

Returns

• order (int) – polynomial order

• npt (int) – number of panel refinements in imaginary time

• npo (int) – number of panel refinements in frequency

• nt (int) – total number of imaginary time points

• no (int) – total number of real frenquency points

pydlr.kernel.kernel_discretization(lamb, error_est=False)

Return kernel discretization correct to machine prescision for given \(\Lambda\)

Parameters

• lamb (float) – Cut-off parameter \(\Lambda\)

• error_est (bool, optional) – Perform convergence check of discretization (NB, slow)

Returns

• kmat (ndarray) – Discretization of the analytical continuation kernel \(K_{jk}\)

• t (ndarray) – Discretized imaginary time grid \(\tau_j\)

• w (ndarray) – Discretized real frequency grid \(\omega_k\)

class pydlr.kernel.KernelInterpolativeDecoposition(lamb, eps=1e-15, xi=-1, max_rank=500, nmax=None, verbose=False)

Interpolative decomposition class for the imaginary-time analytical continuation kernel.

This is the SciPy based driver class for the Discrete Lehmann Representation (DLR).

Parameters

• lamb (float) – DLR scale parameter \(\Lambda\).

• eps (float) – Set accuracy of the DLR representation.

• xi (sign, optional,) – Statistical sign \(\xi = \pm 1\) for bosons and fermions respectively.

• max_rank (int, optional) – Maximum rank of the DLR kernel decomposition. Default 500.

• nmax (int, optional) – Maxumum index of the Matsubara frequency grid. Default int(lamb).

• verbose (bool, optional) – Default False.

utils

Utilities functions used in testing and demonstration.

pydlr.utils.analytic_bethe_G_tau(tau, beta)

Get Bethe graph imaginary time Green’s function.

Obtained by semi-analytic evaluation using adaptive integration of

\[G(\tau) = -\frac{2}{\pi} \int_{-1}^{1} K\left(\frac{\tau}{\beta}, \beta \omega\right) \sqrt{1 - \omega^2} \, d\omega\]

Parameters

• tau ((n), ndarray) – Imaginary time \(\tau\) points to evaulate the Bethe Green’s function on.

• beta (float) – Inverse temperature \(\beta\)

Returns

G_iaa – Imaginary time Green’s function \(G(\tau)\) evaluated on the given imaginary times.

Return type

(n, 1, 1), ndarray