IR Basis

class sparse_ir.FiniteTempBasis(statistics, beta, wmax, eps=None, *, max_size=None, kernel=None, sve_result=None)

Intermediate representation (IR) basis for given temperature.

For a continuation kernel from real frequencies, ω ∈ [-ωmax, ωmax], to imaginary time, τ ∈ [0, beta], this class stores the truncated singular value expansion or IR basis:

\[K(\tau, \omega) \approx \sum_{l=0}^{L-1} U_l(\tau) S_l V_l(\omega),\]

where U are the IR basis functions on the imaginary-time axis, stored in u, S are the singular values, stored in s, and V are the IR basis functions on the real-frequency axis, stored in V. The IR basis functions in Matsubara frequency are stored in uhat.

Example

The following example code assumes the spectral function is a single pole at ω = 2.5:

# Compute IR basis for fermions and β = 10, W <= 4.2
import sparse_ir
basis = sparse_ir.FiniteTempBasis(statistics='F', beta=10, wmax=4.2)

# Assume spectrum is a single pole at ω = 2.5, compute G(iw)
# on the first few Matsubara frequencies
gl = basis.s * basis.v(2.5)
giw = gl @ basis.uhat([1, 3, 5, 7])

property accuracy

Accuracy of the basis.

Upper bound to the relative error of reprensenting a propagator with the given number of basis functions (number between 0 and 1).

property beta: Inverse temperature

default_matsubara_sampling_points(*, npoints=None, positive_only=False)

Default sampling points on the imaginary frequency axis

Parameters:

npoints (int) –
Minimum number of sampling points to return.
positive_only (bool) – Only return non-negative frequencies. This is useful if the object to be fitted is symmetric in Matsubura frequency, ghat(w) == ghat(-w).conj(), or, equivalently, real in imaginary time.

default_omega_sampling_points(*, npoints=None)

Return default sampling points in imaginary time.

Parameters:

npoints (int) –

Minimum number of sampling points to return.

default_tau_sampling_points(*, npoints=None)

Default sampling points on the imaginary time axis

Parameters:

npoints (int) –

Minimum number of sampling points to return.

property kernel: Kernel of which this is the singular value expansion

property lambda_: Basis cutoff parameter, Λ == β * wmax, or None if not present

rescale(new_beta)

Return a basis for different temperature.

Uses the same kernel with the same eps, but a different temperature. Note that this implies a different UV cutoff wmax, since lambda_ == beta * wmax stays constant.

property s: ndarray: Vector of singular values of the continuation kernel

property shape: Shape of the basis function set

property significance

Significances of the basis functions

Vector of significance values, one for each basis function. Each value is a number between 0 and 1 which is a a-priori bound on the (relative) error made by discarding the associated coefficient.

property size: Number of basis functions / singular values

property statistics: Quantum statistic (“F” for fermionic, “B” for bosonic)

property u: PiecewiseLegendrePoly

Basis functions on the imaginary time axis.

Set of IR basis functions on the imaginary time (tau) axis, where tau is a real number between zero and beta. To get the l-th basis function at imaginary time tau of some basis basis, use:

ultau = basis.u[l](tau)        # l-th basis function at time tau

Note that u supports vectorization both over l and tau. In particular, omitting the subscript yields a vector with all basis functions, evaluated at that position:

basis.u(tau) == [basis.u[l](tau) for l in range(basis.size)]

Similarly, supplying a vector of tau points yields a matrix A, where A[l,n] corresponds to the l-th basis function evaluated at tau[n]:

tau = [0.5, 1.0]
basis.u(tau) == \
    [[basis.u[l](t) for t in tau] for l in range(basis.size)]

property uhat: PiecewiseLegendreFT

Basis functions on the reduced Matsubara frequency (wn) axis.

Set of IR basis functions reduced Matsubara frequency (wn) axis, where wn is an integer. These are related to u by the following Fourier transform:

\[\hat u(n) = \int_0^\beta d\tau \exp(i\pi n \tau/\beta) u(\tau)\]

To get the l-th basis function at some reduced frequency wn of some basis basis, use:

uln = basis.uhat[l](wn)        # l-th basis function at freq wn

uhat supports vectorization both over l and wn. In particular, omitting the subscript yields a vector with all basis functions, evaluated at that position:

basis.uhat(wn) == [basis.uhat[l](wn) for l in range(basis.size)]

Similarly, supplying a vector of wn points yields a matrix A, where A[l,n] corresponds to the l-th basis function evaluated at wn[n]:

wn = [1, 3]
basis.uhat(wn) == \\
    [[basis.uhat[l](wi) for wi in wn] for l in range(basis.size)]

Note

Instead of the value of the Matsubara frequency, these functions expect integers corresponding to the prefactor of pi over beta. For example, the first few positive fermionic frequencies would be specified as [1, 3, 5, 7], and the first bosonic frequencies are [0, 2, 4, 6]. This is also distinct to an index!

property v: PiecewiseLegendrePoly

Basis functions on the real frequency axis.

Set of IR basis functions on the real frequency (omega) axis, where omega is a real number of magnitude less than wmax. To get the l-th basis function at real frequency omega of some basis basis, use:

ulomega = basis.v[l](omega)    # l-th basis function at freq. omega

Note that v supports vectorization both over l and omega. In particular, omitting the subscript yields a vector with all basis functions, evaluated at that position:

basis.v(omega) == [basis.v[l](omega) for l in range(basis.size)]

Similarly, supplying a vector of omega points yields a matrix A, where A[l,n] corresponds to the l-th basis function evaluated at omega[n]:

omega = [0.5, 1.0]
basis.v(omega) == \
    [[basis.v[l](t) for t in omega] for l in range(basis.size)]

property wmax: Real frequency cutoff or None if not present

Piecewise polynomials

class sparse_ir.poly.PiecewiseLegendrePoly(data, knots, dx=None, symm=None)

Piecewise Legendre polynomial.

Models a function on the interval [-1, 1] as a set of segments on the intervals S[i] = [a[i], a[i+1]], where on each interval the function is expanded in scaled Legendre polynomials.

__call__(x): Evaluate polynomial at position x

__getitem__(l): Return part of a set of piecewise polynomials

deriv(n=1): Get polynomial for the n’th derivative

overlap(f, *, rtol=2.3e-16, return_error=False, points=None)

Evaluate overlap integral of this polynomial with function f.

Given the function f, evaluate the integral:

∫ dx * f(x) * self(x)

using piecewise Gauss-Legendre quadrature, where self are the polynomials.

Parameters:

f (callable) – function that is called with a point x and returns f(x) at that position.
points (sequence of floats) – A sequence of break points in the integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities)

Returns:

array-like object with shape (poly_dims, f_dims) poly_dims are the shape of the polynomial and f_dims are those of the function f(x).

roots(alpha=2)

Find all roots of the piecewise polynomial

Assume that between each two knots (pieces) there are at most alpha roots.

value(l, x): Return value for l and x.

class sparse_ir.poly.PiecewiseLegendreFT(poly, freq='even', n_asymp=None, power_model=None)

Fourier transform of a piecewise Legendre polynomial.

For a given frequency index n, the Fourier transform of the Legendre function is defined as:

phat(n) == ∫ dx exp(1j * pi * n * x / (xmax - xmin)) p(x)

The polynomial is continued either periodically (freq='even'), in which case n must be even, or antiperiodically (freq='odd'), in which case n must be odd.

__call__(n): Obtain Fourier transform of polynomial for given frequencies

extrema(*, part=None, grid=None, positive_only=False): Obtain extrema of Fourier-transformed polynomial.

sign_changes(*, part=None, grid=None, positive_only=False): Obtain sign changes of Fourier-transformed polynomial.