API reference

adapol.approx_freq_aaa(F, Z, max_n_poles=None, aaa_tol=None, verbose=False)[source]

Approximate frequency data \(F\) sampled at points \(Z\) with a sum of simple poles, by running the AAA algorithm.

Parameters:

F ((N, ...) array_like) – Frequency data to approximate, sampled at points \(Z\).
Z ((N,) array_like) – Sample points in (complex) frequency space, at which \(F\) is sampled.
max_n_poles (int, optional) – Maximum number of poles to use in the approximation.
aaa_tol (float, optional) – Tolerance on the AAA residual, i.e. on the pole step (see the note below).
verbose (bool, optional) – If True, print verbose output during the approximation process.

Returns:

poles ((M,) ndarray) – Poles of the approximating sum of simple poles (M <= max_n_poles).
residues ((M, …) ndarray) – Residues of the approximating sum of simple poles.
error (float) – Maximum absolute error of the approximation at the sample points.

Notes

The approximation is built in two steps:

Pole step: the AAA algorithm is run on the data \((Z, F)\) to determine the pole locations.
Residue step: the residues are then determined by a linear least-squares fit that minimizes the error against the frequency-domain data \(F\) at the sample points \(Z\).

If only max_n_poles is set, AAA runs until at most max_n_poles poles are used.
If only aaa_tol is set, AAA runs until the AAA residual tolerance aaa_tol is reached.
If both are set, AAA stops as soon as either the AAA residual tolerance aaa_tol is reached or max_n_poles poles are used, whichever happens first.

Note

The number of poles produced might be smaller than max_n_poles, for two reasons:

an odd number of poles is always produced, from symmetry considerations, and
the AAA cleanup step might remove poles with small residues.

Note

The tolerance aaa_tol controls the AAA residual, i.e. the maximum absolute deviation from the frequency-domain data \(F\) over the sample points \(Z\) not yet used as AAA support points (pole step); it does not bound the final error, which also depends on the subsequent residue fit.

Examples

Fit a simple function with two poles \(F(Z) = 1/(Z-1) + 0.5/(Z+2)\) sampled on an equispaced imaginary-frequency grid \(Z \in [-10i, 10i]\).

>>> import numpy as np
>>> np.set_printoptions(precision=2, suppress=True)
>>> from adapol.adapol import approx_freq_aaa
>>> # Sample frequency data F at points Z
>>> Z = 1j * np.linspace(-10, 10, 100)
>>> F = 1 / (Z - 1) + 0.5 / (Z + 2)  # Example frequency data with two poles
>>> # Approximate F with a sum of simple poles using AAA
>>> poles, residues, error = approx_freq_aaa(F, Z, aaa_tol=1e-12)
>>> len(poles)
3
>>> physical = np.abs(residues) > 1e-8   # discard the null-residue pole
>>> poles[physical]
array([-2.,  1.])
>>> residues[physical].round(2) + 0.0
array([0.5+0.j, 1. +0.j])
>>> float(error) < 1e-12
True

Note that the fit contains three poles, not two, due to the constrained AAA algorithm. However, the additional pole has a residue that is numerically zero. Since it does not contribute to the fit, its location is not determined by the data: it is fixed only by round-off and therefore varies between platforms and BLAS implementations. For that reason the examples select the poles with a non-negligible residue instead of printing the returned arrays verbatim.

Also tensor valued functions can be fitted, e.g. a 2x2 matrix valued function \(\dim(F(Z)) = 2 \times 2\) with two poles and two matrix residues:

>>> R1 = np.array([[1, 0.1j], [-0.1j, 0]])[None, ...]
>>> R2 = np.array([[0, 0.1], [0.1, 1]])[None, ...]
>>> F = R1 / (Z[:, None, None] - 1) + R2 / (Z[:, None, None] + 2)
>>> poles, residues, error = approx_freq_aaa(F, Z, aaa_tol=1e-12)
>>> weight = np.abs(residues).max(axis=(1, 2))   # residue norm of each pole
>>> poles[weight > 1e-8]
array([-2.,  1.])
>>> residues[weight > 1e-8].round(2) + 0.0  # round-off noise in the zero entries removed
array([[[0. +0.j , 0.1+0.j ],
        [0.1+0.j , 1. +0.j ]],

       [[1. +0.j , 0. +0.1j],
        [0. -0.1j, 0. +0.j ]]])
>>> float(error) < 1e-12
True

adapol.approx_sop_fast(poles, residues, beta, max_n_poles=None, aaa_tol=None, nonlinear_optimization=False, Z=None, verbose=False)[source]

Approximate a sum of simple poles defined by poles and residues with a sum of a (possibly) smaller number of simple poles, by running the AAA algorithm and, optionally, a non-linear optimization step.

Parameters:

poles ((K,) array_like) – Poles \(P_k\) of the original sum of simple poles to approximate.
residues ((K, ...) array_like) – Residues \(R_k\) of the original sum of simple poles to approximate.
beta (float) – Inverse temperature, used to define the L2 norm in imaginary time and the imaginary-frequency grid on which the AAA data is sampled.
max_n_poles (int, optional) – Maximum number of poles to use in the approximation.
aaa_tol (float, optional) – Tolerance on the AAA residual. This is the maximum absolute deviation from the imaginary-frequency-domain data used by the AAA algorithm (the original sum of simple poles evaluated on the grid described below), over the sample points not yet used as AAA support points.
nonlinear_optimization (bool, optional) – If True, run a non-linear optimization step after the AAA approximation, using the AAA poles only as an initial guess.
Z ((N,) array_like, optional) – Custom sample points in (complex) frequency space at which to sample the original sum of simple poles for the AAA algorithm. If not given, a symmetric fermionic Matsubara grid is used by default (see Notes).
verbose (bool, optional) – If True, print verbose output during the approximation process.

Returns:

poles ((M,) ndarray) – Poles of the approximating sum of simple poles.
residues ((M, …) ndarray) – Residues of the approximating sum of simple poles.
error (float) – Normalized imaginary-time \(L^2(\tau)\) norm of the difference between the original and approximating sum of simple poles (see Notes).

Notes

The original sum of simple poles

\[F(Z) = \sum_{k=1}^K \frac{R_k}{Z - P_k}\]

is first evaluated on an equispaced grid \(Z_n\) on the imaginary-frequency axis (see below), to produce the frequency-domain data \(F_n = F(Z_n)\) used by the AAA algorithm.

The approximation is then built in two steps:

Pole step: the AAA algorithm is run on this frequency-domain data \((Z_n, F_n)\) to determine the pole locations.
Residue step: the residues are determined by minimizing the imaginary-time \(L^2(\tau)\) norm of the difference from the original sum of simple poles. By default this is a linear least-squares fit of the residues for the AAA poles; if nonlinear_optimization is True, the pole locations and residues are instead jointly optimized, for details see below.

Note that, unlike approx_freq_aaa, the residues here are fit in imaginary time, not in the frequency domain.

The imaginary-time \(L^2(\tau)\) norm is normalized by the inverse temperature \(\beta\),

\[\lVert f \rVert_{L^2(\tau)} = \left( \frac{1}{\beta} \int_0^\beta |f(\tau)|^2 \, d\tau \right)^{1/2}.\]

If only max_n_poles is set, AAA runs until at most max_n_poles poles are used.
If only aaa_tol is set, AAA runs until the AAA residual tolerance aaa_tol is reached.
If both are set, AAA stops as soon as either the AAA residual tolerance aaa_tol is reached or max_n_poles poles are used, whichever happens first.

Frequency grid: By default, the frequency-domain data used by the AAA algorithm is obtained by evaluating the original sum of simple poles on a symmetric fermionic Matsubara grid \(Z_n = i (2 n + 1) \pi / \beta\) for \(n = -(n_{max}), \ldots, -1, 0, 1, \ldots, n_{max} - 1\), with \(n_{max} = \lfloor 4 \beta \, \omega_{max} / \pi \rfloor + 1\) and \(\omega_{max} = \max_k |poles_k|\) the largest pole magnitude. A custom sample grid may instead be supplied via the optional Z argument, e.g. to use a bosonic Matsubara grid.

Non-linear optimization: The optional non-linear optimization step keeps the number of poles fixed and jointly relocates the pole positions and refits the residues to minimize the imaginary-time L2 norm error between the original and approximating sum of simple poles. It is solved with the L-BFGS-B method (using an analytic gradient) and is run to optimizer convergence, i.e. to reduce the error as much as possible rather than to meet a prescribed error tolerance.

Note

The number of poles produced might be smaller than max_n_poles, for two reasons:

an odd number of poles is always produced, from symmetry considerations, and
the AAA cleanup step might remove poles with small residues.

Note

Setting aaa_tol does not guarantee that the final imaginary-time L2 norm error is below aaa_tol; this is not possible to guarantee within the AAA algorithm alone. However, this function returns the final imaginary-time L2 norm error of approximation (see Returns above).

Examples

Approximate a sum of three simple poles \(s(z) = 1/(z-1) + 0.5/(z+2) + 0.3/(z-0.5)\), specified by its poles and residues, using AAA. Here the input is already minimal, so the three poles are recovered (up to ordering), with the residues re-fit in imaginary time for inverse temperature \(\beta\).

>>> import numpy as np
>>> np.set_printoptions(precision=2, suppress=True)
>>> from adapol.adapol import approx_sop_fast
>>> poles = np.array([1.0, -2.0, 0.5])
>>> residues = np.array([1.0, 0.5, 0.3])
>>> poles, residues, error = approx_sop_fast(poles, residues, beta=20.0, aaa_tol=1e-12)
>>> poles
array([-2. ,  0.5,  1. ])
>>> residues
array([0.5, 0.3, 1. ])
>>> float(error) < 1e-9
True

Unlike approx_freq_aaa, the residues are fit by minimizing the imaginary-time \(L^2(\tau)\) error, so error is an imaginary-time norm rather than a frequency-domain sample error.

Tensor-valued residues are supported as well, e.g. a 2x2 matrix-valued sum of poles with \(\dim(R_k) = 2 \times 2\):

>>> R1 = np.array([[1, 0.1j], [-0.1j, 0]])
>>> R2 = np.array([[0, 0.1], [0.1, 1]])
>>> R3 = np.array([[0.5, 0.0], [0.0, 0.5]])
>>> poles = np.array([1.0, -2.0, 0.3])
>>> residues = np.array([R1, R2, R3])
>>> poles, residues, error = approx_sop_fast(poles, residues, beta=20.0, aaa_tol=1e-12)
>>> poles
array([-2. ,  0.3,  1. ])
>>> residues.round(2) + 0.0  # round-off noise in the zero entries removed
array([[[0. +0.j , 0.1+0.j ],
        [0.1+0.j , 1. +0.j ]],

       [[0.5+0.j , 0. +0.j ],
        [0. +0.j , 0.5+0.j ]],

       [[1. +0.j , 0. +0.1j],
        [0. -0.1j, 0. +0.j ]]])
>>> float(error) < 1e-9
True

adapol.approx_sop_tol(poles, residues, tol, beta, nonlinear_optimization=False, Z=None, verbose=False)[source]

Approximate a sum of simple poles defined by poles and residues with the smallest sum of simple poles whose imaginary-time \(L^2(\tau)\) error is below the tolerance tol.

Parameters:

poles ((K,) array_like) – Poles of the original sum of simple poles to approximate.
residues ((K, ...) array_like) – Residues of the original sum of simple poles to approximate.
tol (float) – Target tolerance on the final imaginary-time \(L^2(\tau)\) norm of the difference between the original and approximating sum of simple poles.
beta (float) – Inverse temperature, used to define the L2 norm in imaginary time and the imaginary-frequency grid on which the AAA data is sampled.
nonlinear_optimization (bool, optional) – If True, the residue step jointly optimizes the pole locations and residues (rather than fitting residues only) to minimize the imaginary-time \(L^2( au)\) error.
Z ((N,) array_like, optional) – Custom sample points in (complex) frequency space at which to sample the original sum of simple poles for the AAA algorithm. If not given, a symmetric fermionic Matsubara grid is used by default (see approx_sop_fast).
verbose (int or bool, optional) – Amount of printed output. 0 (or False) is silent, 1 (or True) prints one line per pass through the AAA and residue fit pipeline, showing the pole count and error of each candidate fit, and 2 additionally prints the indented per step output of the AAA algorithm itself.

Returns:

poles ((M,) ndarray) – Poles of the approximating sum of simple poles.
residues ((M, …) ndarray) – Residues of the approximating sum of simple poles.
error (float) – Normalized imaginary-time \(L^2(\tau)\) norm of the difference between the original and approximating sum of simple poles (with the \(1/\beta\) normalization defined in approx_sop_fast).

Raises:

ValueError – If the target tolerance tol cannot be achieved within the internal maximum number of search steps.

Notes

Unlike approx_sop_fast, where the tolerance only controls the AAA residual in the pole step, here tol is imposed on the final imaginary-time error, i.e. after the residues have been fit. Since AAA only determines pole locations and the residues (and hence the final error) are only known after the residue step, the requested error cannot be reached by a single AAA run. Instead, the AAA + residue pipeline is run repeatedly to search for the minimal number of poles that achieves tol.

Each candidate fit is built in two steps:

Pole step: the original sum of simple poles is evaluated on a symmetric fermionic Matsubara grid (see approx_sop_fast for the grid definition, and the optional Z argument to override it) and AAA is run on this data to determine pole locations.
Residue step: the residues are determined by minimizing the imaginary-time \(L^2(\tau)\) norm of the difference from the original sum of simple poles – a linear least-squares fit of the residues for the AAA poles, or, if nonlinear_optimization is True, a joint optimization of the pole locations and residues.

The pole-count search proceeds in two phases: first the number of AAA poles is increased until the imaginary-time error drops below tol, then a bisection on the number of AAA steps locates the smallest pole count that still meets tol.

Examples

Compress a continuous spectral density into a small sum of simple poles. The density is first discretized as a sum of 200 simple poles on the real axis, which is then compressed to the smallest sum of poles whose imaginary-time \(L^2(\tau)\) error is below tol.

>>> import numpy as np
>>> from adapol.adapol import approx_sop_tol
>>> w = np.linspace(-2, 2, 200)                            # real-frequency grid
>>> dw = w[1] - w[0]
>>> rho = np.sqrt(np.maximum(4 - w**2, 0.0)) / (2 * np.pi)  # semicircle density
>>> poles, residues, error = approx_sop_tol(w, rho * dw, tol=1e-5, beta=20.0)
>>> len(poles) < 30      # 200 input poles compressed to a handful
True
>>> float(error) < 1e-5  # final imaginary-time error is below the tolerance
True

Unlike approx_sop_fast, the number of poles is not prescribed but chosen automatically as the smallest count meeting tol on the final imaginary-time error.

TRIQS interface

adapol.triqs.approx_gf_imfreq_aaa(G_w, max_n_poles=None, aaa_tol=None, verbose=False)[source]

Approximate a Green’s function \(G\) given on an imaginary-frequency mesh with a sum of simple poles, by running the AAA algorithm.

This is the TRIQS front-end to adapol.approx_freq_aaa, taking the frequency data \(F\) and the sample points \(Z\) from the Green’s function and its mesh.

Parameters:

G_w (triqs.gf.Gf) – Green’s function to approximate, defined on an imaginary-frequency mesh (MeshImFreq or MeshDLRImFreq).
max_n_poles (int, optional) – Maximum number of poles to use in the approximation.
aaa_tol (float, optional) – Tolerance on the AAA residual, i.e. on the pole step. It does not bound the final error, which also depends on the subsequent residue fit.
verbose (bool, optional) – If True, print verbose output during the approximation process.

Returns:

poles ((M,) ndarray) – Poles of the approximating sum of simple poles (M <= max_n_poles).
residues ((M, …) ndarray) – Residues of the approximating sum of simple poles.
error (float) – Maximum absolute error of the approximation at the mesh points.

Raises:

ValueError – If neither max_n_poles nor aaa_tol is given.

Notes

If only max_n_poles is set, AAA runs until at most max_n_poles poles are used.
If only aaa_tol is set, AAA runs until the AAA residual tolerance aaa_tol is reached.
If both are set, AAA stops as soon as either the AAA residual tolerance aaa_tol is reached or max_n_poles poles are used, whichever happens first.

See adapol.approx_freq_aaa for details on the algorithm, on why the number of poles produced may be smaller than max_n_poles, and on why aaa_tol does not bound the returned error.