SparseIR

Intermediate representation (IR) for many-body propagators

(polyFT::PiecewiseLegendreFT)(n)

Obtain Fourier transform of polynomial for given frequency index n.

AbstractKernel

Integral kernel K(x, y).

Abstract base type for an integral kernel, i.e. a AbstractFloat binary function $K(x, y)$ used in a Fredhold integral equation of the first kind:

\[ u(x) = ∫ K(x, y) v(y) dy\]

where $x ∈ [x_\mathrm{min}, x_\mathrm{max}]$ and $y ∈ [y_\mathrm{min}, y_\mathrm{max}]$. For its SVE to exist, the kernel must be square-integrable, for its singular values to decay exponentially, it must be smooth.

In general, the kernel is applied to a scaled spectral function $ρ'(y)$ as:

\[ ∫ K(x, y) ρ'(y) dy,\]

where $ρ'(y) = w(y) ρ(y)$.

(kernel::AbstractKernel)(x, y[, x₊, x₋])

Evaluate kernel at point (x, y).

The parameters x₊ and x₋, if given, shall contain the values of x - xₘᵢₙ and xₘₐₓ - x, respectively. This is useful if either difference is to be formed and cancellation expected.

AbstractSVEHints

Discretization hints for singular value expansion of a given kernel.

AbstractSampling

Abstract class for sparse sampling.

Encodes the "basis transformation" of a propagator from the truncated IR basis coefficients G_ir[l] to time/frequency sampled on sparse points G(x[i]) together with its inverse, a least squares fit:

     ________________                   ___________________
    |                |    evaluate     |                   |
    |     Basis      |---------------->|     Value on      |
    |  coefficients  |<----------------|  sampling points  |
    |________________|      fit        |___________________|

CentrosymmSVE <: AbstractSVE

SVE of centrosymmetric kernel in block-diagonal (even/odd) basis.

For a centrosymmetric kernel K, i.e., a kernel satisfying: K(x, y) == K(-x, -y), one can make the following ansatz for the singular functions:

u[l](x) = ured[l](x) + sign[l] * ured[l](-x)
v[l](y) = vred[l](y) + sign[l] * ured[l](-y)

where sign[l] is either +1 or -1. This means that the singular value expansion can be block-diagonalized into an even and an odd part by (anti-)symmetrizing the kernel:

K_even = K(x, y) + K(x, -y)
K_odd  = K(x, y) - K(x, -y)

The lth basis function, restricted to the positive interval, is then the singular function of one of these kernels. If the kernel generates a Chebyshev system [1], then even and odd basis functions alternate.

[1]: A. Karlin, Total Positivity (1968).

Union of several basis functions in the imaginary-time/real-frequency domain domains

Evaluate basis function at position x

Union of several basis functions in the imaginary-frequency domain domains

Evaluate basis function at frequency n

DimensionlessBasis <: AbstractBasis

Intermediate representation (IR) basis in reduced variables.

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

K(x, y) ≈ sum(u[l](x) * s[l] * v[l](y) for l in range(L))

The functions are given in reduced variables, x = 2τ/β - 1 and y = ω/ωmax, which scales both sides to the interval [-1, 1]. The kernel then only depends on a cutoff parameter Λ = β * ωmax.

Examples

The following example code assumes the spectral function is a single pole at x = 0.2. We first compute an IR basis suitable for fermions and β*W ≤ 42. Then we get G(iw) on the first few Matsubara frequencies:

julia> using SparseIR

julia> basis = DimensionlessBasis(fermion, 42);

julia> gl = basis.s .* basis.v(0.2);

julia> giw = transpose(basis.uhat([1, 3, 5, 7])) * gl

Fields

• u::PiecewiseLegendrePolyVector: Set of IR basis functions on the reduced imaginary time (x) axis. These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of points x, you can call the function u(x). To obtain a single basis function, a slice or a subset l, you can use u[l].

• uhat::PiecewiseLegendreFTVector: Set of IR basis functions on the Matsubara frequency (wn) axis.

These objects are stored as a set of Bessel functions.

To obtain the value of all basis functions at a Matsubara frequency or a array of points wn, you can call the function uhat(wn). Note that we expect reduced frequencies, which are simply even/odd numbers for bosonic/fermionic objects. To obtain a single basis function, a slice or a subset l, you can use uhat[l].

• s: Vector of singular values of the continuation kernel

• v::PiecewiseLegendrePolyVector: Set of IR basis functions on the reduced real frequency (y) axis.

These functions are stored as piecewise Legendre polynomials.

To obtain the value of all basis functions at a point or a array of points y, you can call the function v(y). To obtain a single basis function, a slice or a subset l, you can use v[l].

See also FiniteTempBasis for a basis directly in time/frequency.

DimensionlessBasis(statistics, Λ, ε=nothing; kernel=LogisticKernel(Λ), sve_result=compute_sve(kernel; ε))

Construct an IR basis suitable for the given statistics and cutoff Λ.

FiniteTempBasis <: AbstractBasis

Intermediate representation (IR) basis for given temperature.

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

K(τ, ω) ≈ sum(u[l](τ) * s[l] * v[l](ω) for l in 1:L)

This basis is inferred from a reduced form by appropriate scaling of the variables.

Examples

The following example code assumes the spectral function is a single pole at ω = 2.5. We first compute an IR basis suitable for fermions and β = 10, W ≤ 4.2. Then we get G(iw) on the first few Matsubara frequencies:

julia> using SparseIR

julia> basis = FiniteTempBasis(fermion, 42, 4.2);

julia> gl = basis.s .* basis.v(2.5);

julia> giw = transpose(basis.uhat([1, 3, 5, 7])) * gl

Fields

• u::PiecewiseLegendrePolyVector: Set of IR basis functions on the imaginary time (tau) axis. These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of points x, you can call the function u(x). To obtain a single basis function, a slice or a subset l, you can use u[l].

• uhat::PiecewiseLegendreFT: Set of IR basis functions on the Matsubara frequency (wn) axis. These objects are stored as a set of Bessel functions.

  To obtain the value of all basis functions at a Matsubara frequency or a array of points wn, you can call the function uhat(wn). Note that we expect reduced frequencies, which are simply even/odd numbers for bosonic/fermionic objects. To obtain a single basis function, a slice or a subset l, you can use uhat[l].

• s: Vector of singular values of the continuation kernel

• v::PiecewiseLegendrePoly: Set of IR basis functions on the real frequency (w) axis. These functions are stored as piecewise Legendre polynomials.

  To obtain the value of all basis functions at a point or a array of points w, you can call the function v(w). To obtain a single basis function, a slice or a subset l, you can use v[l].

FiniteTempBasis(statistics, β, wmax, ε=nothing; kernel=LogisticKernel(β * wmax), sve_result=compute_sve(kernel; ε))

Construct a finite temperature basis suitable for the given statistics and cutoffs β and wmax.

FiniteTempBasisSet

Class for holding IR bases and sparse-sampling objects.

An object of this class holds IR bases for fermions and bosons and associated sparse-sampling objects.

Fields

• basis_f::FiniteTempBasis: Fermion basis
• basis_b::FiniteTempBasis: Boson basis
• beta::Float64: Inverse temperature
• wmax::Float64: Cut-off frequency
• tau::Vector{Float64}: Sampling points in the imaginary-time domain
• wn_f::Vector{Int}: Sampling fermionic frequencies
• wn_b::Vector{Int}: Sampling bosonic frequencies
• smpltauf::TauSampling: Sparse sampling for tau & fermion
• smpltaub::TauSampling: Sparse sampling for tau & boson
• smplwnf::MatsubaraSampling: Sparse sampling for Matsubara frequency & fermion
• smplwnb::MatsubaraSampling: Sparse sampling for Matsubara frequency & boson
• sve_result::Tuple{PiecewiseLegendrePoly,Vector{Float64},PiecewiseLegendrePoly}: Results of SVE

FiniteTempBasisSet(β, wmax, ε; sve_result=compute_sve(LogisticKernel(β * wmax); ε))

Create basis sets for fermion and boson and associated sampling objects. Fermion and bosonic bases are constructed by SVE of the logistic kernel.

Legendre basis

In the original paper [L. Boehnke et al., PRB 84, 075145 (2011)], they used:

G(\tau) = \sum_{l=0} \sqrt{2l+1} P_l[x(\tau)] G_l/beta,

where P_l[x] is the $l$-th Legendre polynomial.

In this class, the basis functions are defined by

U_l(\tau) \equiv c_l (\sqrt{2l+1}/beta) * P_l[x(\tau)],

where cl are additional l-depenent constant factors. By default, we take cl = 1, which reduces to the original definition.

LogisticKernel <: AbstractKernel

Fermionic/bosonic analytical continuation kernel.

In dimensionless variables $x = 2 τ/β - 1$, $y = β ω/Λ$, the integral kernel is a function on $[-1, 1] × [-1, 1]$:

\[ K(x, y) = \frac{e^{-Λ y (x + 1) / 2}}{1 + e^{-Λ y}}\]

LogisticKernel is a fermionic analytic continuation kernel. Nevertheless, one can model the $τ$ dependence of a bosonic correlation function as follows:

\[ ∫ \frac{e^{-Λ y (x + 1) / 2}}{1 - e^{-Λ y}} ρ(y) dy = ∫ K(x, y) ρ'(y) dy,\]

with

\[ ρ'(y) = w(y) ρ(y),\]

where the weight function is given by

\[ w(y) = \frac{1}{\tanh(Λ y/2)}.\]

LogisticKernelOdd <: AbstractReducedKernel

Fermionic analytical continuation kernel, odd.

In dimensionless variables $x = 2τ/β - 1$, $y = βω/Λ$, the fermionic integral kernel is a function on $[-1, 1] × [-1, 1]$:

\[ K(x, y) = -\frac{\sinh(Λ x y / 2)}{\cosh(Λ y / 2)}\]

Constant term in matsubara-frequency domain

MatsubaraSampling <: AbstractSampling

Sparse sampling in Matsubara frequencies.

Allows the transformation between the IR basis and a set of sampling points in (scaled/unscaled) imaginary frequencies.

MatsubaraSampling(basis[, sampling_points])

Construct a MatsubaraSampling object. If not given, the sampling_points are chosen as the (discrete) extrema of the highest-order basis function in Matsubara. This turns out to be close to optimal with respect to conditioning for this size (within a few percent).

NestedRule{T}

Nested quadrature rule.

PiecewiseLegendreFT <: Function

Fourier transform of a piecewise Legendre polynomial.

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

    p̂(n) == ∫ dx exp(im * π * 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.

PiecewiseLegendrePoly <: Function

Piecewise Legendre polynomial.

Models a function on the interval $[xmin, xmax]$ 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.

PiecewiseLegendrePolyVector{T}

Alias for Vector{PiecewiseLegendrePoly{T}}.

PowerModel

Model from a high-frequency series expansion::

A(iω) == sum(A[n] / (iω)^(n+1) for n in 1:N)

where $iω == i * π/2 * wn$ is a reduced imaginary frequency, i.e., $wn$ is an odd/even number for fermionic/bosonic frequencies.

ReducedKernel

Restriction of centrosymmetric kernel to positive interval.

For a kernel $K$ on $[-1, 1] × [-1, 1]$ that is centrosymmetric, i.e. $K(x, y) = K(-x, -y)$, it is straight-forward to show that the left/right singular vectors can be chosen as either odd or even functions.

Consequentially, they are singular functions of a reduced kernel $K_\mathrm{red}$ on $[0, 1] × [0, 1]$ that is given as either:

\[ K_\mathrm{red}(x, y) = K(x, y) \pm K(x, -y)\]

This kernel is what this class represents. The full singular functions can be reconstructed by (anti-)symmetrically continuing them to the negative axis.

RegularizedBoseKernel <: AbstractKernel

Regularized bosonic analytical continuation kernel.

In dimensionless variables $x = 2 τ/β - 1$, $y = β ω/Λ$, the fermionic integral kernel is a function on $[-1, 1] × [-1, 1]$:

\[ K(x, y) = y \frac{e^{-Λ y (x + 1) / 2}}{e^{-Λ y} - 1}\]

Care has to be taken in evaluating this expression around $y = 0$.

RegularizedBoseKernelOdd <: AbstractReducedKernel

Bosonic analytical continuation kernel, odd.

In dimensionless variables $x = 2 τ / β - 1$, $y = β ω / Λ$, the fermionic integral kernel is a function on $[-1, 1] × [-1, 1]$:

\[ K(x, y) = -y \frac{\sinh(Λ x y / 2)}{\sinh(Λ y / 2)}\]

Rule{T<:AbstractFloat}

Quadrature rule.

Approximation of an integral over [a, b] by a sum over discrete points x with weights w:

\[ ∫ f(x) ω(x) dx ≈ ∑_i f(x_i) w_i\]

where we generally have superexponential convergence for smooth $f(x)$ in the number of quadrature points.

SamplingSVE <: AbstractSVE

SVE to SVD translation by sampling technique [1].

Maps the singular value expansion (SVE) of a kernel kernel onto the singular value decomposition of a matrix A. This is achieved by choosing two sets of Gauss quadrature rules: (x, wx) and (y, wy) and approximating the integrals in the SVE equations by finite sums. This implies that the singular values of the SVE are well-approximated by the singular values of the following matrix:

A[i, j] = √(wx[i]) * K(x[i], y[j]) * √(wy[j])

and the values of the singular functions at the Gauss sampling points can be reconstructed from the singular vectors u and v as follows:

u[l,i] ≈ √(wx[i]) u[l](x[i])
v[l,j] ≈ √(wy[j]) u[l](y[j])

[1] P. Hansen, Discrete Inverse Problems, Ch. 3.1

Sparse pole representation

TauSampling <: AbstractSampling

Sparse sampling in imaginary time.

Allows the transformation between the IR basis and a set of sampling points in (scaled/unscaled) imaginary time.

TauSampling(basis[, sampling_points])

Construct a TauSampling object. If not given, the sampling_points are chosen as the extrema of the highest-order basis function in imaginary time. This turns out to be close to optimal with respect to conditioning for this size (within a few percent).

canonicalize!(u, v)

Canonicalize basis.

Each SVD (u[l], v[l]) pair is unique only up to a global phase, which may differ from implementation to implementation and also platform. We fix that gauge by demanding u[l](1) > 0. This ensures a diffeomorphic connection to the Legendre polynomials as Λ → 0.

choose_accuracy(ε, Twork)

Choose work type and accuracy based on specs and defaults

_compute_unl_inner(poly, wn)

Compute piecewise Legendre to Matsubara transform.

_get_tnl(l, w)

Fourier integral of the l-th Legendre polynomial::

Tₗ(ω) == ∫ dx exp(iωx) Pₗ(x)

_phase_stable(poly, wn)

Phase factor for the piecewise Legendre to Matsubara transform.

Compute the following phase factor in a stable way:

exp.(iπ/2 * wn * cumsum(poly.Δx))

_shift_xmid(knots, Δx)

Return midpoint relative to the nearest integer plus a shift.

Return the midpoints xmid of the segments, as pair (diff, shift), where shift is in (0, 1, -1) and diff is a float such that xmid == shift + diff to floating point accuracy.

_split(poly, x)

Split segment.

Find segment of poly's domain that covers x.

check_domain(kernel, x, y)

Check that (x, y) lies within kernel's domain and return it.

check_reduced_matsubara(n[, ζ])

Checks that n is a reduced Matsubara frequency.

Check that the argument is a reduced Matsubara frequency, which is an integer obtained by scaling the freqency ω[n] as follows:

β / π * ω[n] == 2n + ζ

Note that this means that instead of a fermionic frequency (ζ == 1), we expect an odd integer, while for a bosonic frequency (ζ == 0), we expect an even one. If ζ is omitted, any one is fine.

compute_sve(kernel; 
    ε=nothing, n_sv=typemax(Int), n_gauss=nothing, T=Float64, Twork=nothing,
    sve_strat=iscentrosymmetric(kernel) ? CentrosymmSVE : SamplingSVE,
    svd_strat=nothing)

Perform truncated singular value expansion of a kernel.

Perform a truncated singular value expansion (SVE) of an integral kernel kernel : [xmin, xmax] x [ymin, ymax] -> ℝ:

kernel(x, y) == sum(s[l] * u[l](x) * v[l](y) for l in (1, 2, 3, ...)),

where s[l] are the singular values, which are ordered in non-increasing fashion, u[l](x) are the left singular functions, which form an orthonormal system on [xmin, xmax], and v[l](y) are the right singular functions, which form an orthonormal system on [ymin, ymax].

The SVE is mapped onto the singular value decomposition (SVD) of a matrix by expanding the kernel in piecewise Legendre polynomials (by default by using a collocation).

Arguments

• eps::AbstractFloat: Relative cutoff for the singular values.
• n_sv::Integer: Maximum basis size. If given, only at most the n_sv most

significant singular values and associated singular functions are returned.

• n_gauss::Integer: Order of Legendre polynomials. Defaults to hinted value

by the kernel.

• T: Data type of the result.
• Twork: Working data type. Defaults to a data type with

machine epsilon of at least eps^2, or otherwise most accurate data type available.

• sve_strat: SVE to SVD translation strategy. Defaults to SamplingSVE.
• svd_strat: SVD solver. Defaults to fast (ID/RRQR) based solution

when accuracy goals are moderate, and more accurate Jacobi-based algorithm otherwise.

Return value

Return tuple (u, s, v), where:

• u::PiecewiseLegendrePoly: the left singular functions
• s::Vector: singular values
• v::PiecewiseLegendrePoly: the right singular functions

conv_radius(kernel)

Convergence radius of the Matsubara basis asymptotic model.

For improved relative numerical accuracy, the IR basis functions on the Matsubara axis uhat(basis, n) can be evaluated from an asymptotic expression for abs(n) > conv_radius. If isinf(conv_radius), then the asymptotics are unused (the default).

_default_matsubara_sampling_points(basis; mitigate=true)

Default sampling points on the imaginary frequency axis.

default_omega_sampling_points(basis)

Default sampling points on the real-frequency axis.

default_tau_sampling_points(basis)

Default sampling points on the imaginary time/x axis.

deriv(poly)

Get polynomial for the derivative.

eval_matrix(T, basis, x)

Return evaluation matrix from coefficients to sampling points. T <: AbstractSampling.

evaluate!(buffer::AbstractArray{T,N}, sampling, al; dim=1) where {T,N}

Like evaluate, but write the result to buffer. Please use dim = 1 or N to avoid allocating large temporary arrays internally.

evaluate(sampling, al; dim=1)

Evaluate the basis coefficients al at the sparse sampling points.

findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=_DEFAULT_GRID)

Obtain extrema of fourier-transformed polynomial.

finite_temp_bases(β, wmax, ε, sve_result=compute_sve(LogisticKernel(β * wmax); ε))

Construct FiniteTempBasis objects for fermion and bosons using the same LogisticKernel instance.

fit!(buffer, sampling, al::Array{T,N}; dim=1)

Like fit, but write the result to buffer. Please use dim = 1 or N to avoid allocating large temporary arrays internally. The length of workarry cannot be smaller than the returned value of workarrlengthfit.

fit(sampling, al::AbstractArray{T,N}; dim=1)

Fit basis coefficients from the sparse sampling points Please use dim = 1 or N to avoid allocating large temporary arrays internally.

From IR to SPR

gl: Expansion coefficients in IR

get_symmetrized(kernel, sign)

Construct a symmetrized version of kernel, i.e. kernel(x, y) + sign * kernel(x, -y).

getwmax(basis::FiniteTempBasis)

Real frequency cutoff.

giw(polyFT, wn)

Return model Green's function for reduced frequencies

is_centrosymmetric(kernel)

Return true if kernel(x, y) == kernel(-x, -y) for all values of x and y in range. This allows the kernel to be block-diagonalized, speeding up the singular value expansion by a factor of 4. Defaults to false.

iswellconditioned(basis)

Return true if the sampling is expected to be well-conditioned.

joinrules(rules)

Join multiple Gauss quadratures together.

legder

legendre(n[, T])

Gauss-Legendre quadrature with n points on [-1, 1].

legendre_collocation(rule, n=length(rule.x))

Generate collocation matrix from Gauss-Legendre rule.

legvander(x, deg)

Pseudo-Vandermonde matrix of degree deg.

matop!(buffer, mat, arr::AbstractArray, op, dim)

Apply the operator op to the matrix mat and to the array arr along the first dimension (dim=1) or the last dimension (dim=N).

matop_along_dim!(buffer, mat, arr::AbstractArray, dim::Integer, op)

Apply the operator op to the matrix mat and to the array arr along the dimension dim, writing the result to buffer.

matrices(sve::AbstractSVE)

SVD problems underlying the SVE.

matrix_from_gauss(kernel, gauss_x, gauss_y)

Compute matrix for kernel from Gauss rules.

movedim(arr::AbstractArray, src => dst)

Move arr's dimension at src to dst while keeping the order of the remaining dimensions unchanged.

ngauss(hints)

Gauss-Legendre order to use to guarantee accuracy.

nsvals(hints)

Upper bound for number of singular values.

Upper bound on the number of singular values above the given threshold, i.e. where s[l] ≥ ε * first(s).

overlap(poly::PiecewiseLegendrePoly, f)

Evaluate overlap integral of poly with arbitrary function f.

Given the function f, evaluate the integral::

∫ dx * f(x) * poly(x)

using adaptive Gauss-Legendre quadrature.

piecewise(rule, edges)

Piecewise quadrature with the same quadrature rule, but scaled.

postprocess(sve::AbstractSVE, u, s, v, T=nothing)

Construct the SVE result from the SVD.

quadrature(rule, f)

Approximate f's integral.

reseat(rule, a, b)

Reseat quadrature rule to new domain.

roots(poly)

Find all roots of the piecewise polynomial poly.

scale(rule, factor)

Scale weights by factor.

segments_x(kernel)

Segments for piecewise polynomials on the $x$ axis.

List of segments on the $x$ axis for the associated piecewise polynomial. Should reflect the approximate position of roots of a high-order singular function in $x$.

segments_y(kernel)

Segments for piecewise polynomials on the $y$ axis.

List of segments on the $y$ axis for the associated piecewise polynomial. Should reflect the approximate position of roots of a high-order singular function in $y$.

sve_hints(kernel, ε)

Provide discretisation hints for the SVE routines.

Advises the SVE routines of discretisation parameters suitable in tranforming the (infinite) SVE into an (finite) SVD problem.

See also AbstractSVEHints.

From SPR to IR

g_spr: Expansion coefficients in SPR

truncate(u, s, v[, rtol][, lmax])

Truncate singular value expansion.

Arguments

- `u`, `s`, `v`: Thin singular value expansion
- `rtol` : If given, only singular values satisfying `s[l]/s[0] > rtol` are retained.
- `lmax` : If given, at most the `lmax` most significant singular values are retained.

weight_func(kernel, statistics::Statistics)

Return the weight function for the given statistics.

workarrlength(smpl::AbstractSampling, al; dim=1)

Return length of workarr for fit!.

xrange(kernel)

Return a tuple $(x_\mathrm{min}, x_\mathrm{max})$ delimiting the range of allowed x values.

ypower(kernel)

Power with which the $y$ coordinate scales.

yrange(kernel)

Return a tuple $(y_\mathrm{min}, y_\mathrm{max})$ delimiting the range of allowed y values.

Λ(basis)

Basis cutoff parameter Λ = β * ωmax.