{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "6cb39adc-f35e-4454-af78-89c5d9a282c1",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (8,5) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"id": "188c4e36",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# Adaptive pole approximation for a continuous spectrum (semicircular density)"
]
},
{
"cell_type": "markdown",
"id": "7480a9dc",
"metadata": {},
"source": [
"In this example we consider the Matsubara Green's function\n",
"\n",
"$$G(\\mathrm{i}\\nu_n) = \\int_{-\\infty}^{\\infty} \\frac{\\rho(\\omega)}{\\mathrm{i}\\nu_n - \\omega}\\, d\\omega$$\n",
"\n",
"for the semicircular spectral density\n",
"\n",
"$$\\rho(\\omega) = \\frac{2}{\\pi}\\sqrt{1-\\omega^2}, \\qquad -1 \\leq \\omega \\leq 1,$$\n",
"\n",
"and $\\rho(\\omega) = 0$ otherwise. Here $\\mathrm{i}\\nu_n = \\mathrm{i}\\,(2n+1)\\pi/\\beta$ are the fermionic Matsubara frequencies at inverse temperature $\\beta$.\n",
"\n",
"For the semicircular density the integral is known in closed form,\n",
"\n",
"$$G(z) = 2\\left(z - \\sqrt{z+1}\\,\\sqrt{z-1}\\right),$$\n",
"\n",
"so we use this analytic formula to generate the Matsubara data, rather than evaluating the integral by adaptive quadrature as in the `semicircle` example (note that this writing of $\\sqrt{z^2-1}$ numerically selects the correct branch in the upper and lower half-planes).\n",
"\n",
"This notebook demonstrates the three functions of the `adapol` public API:\n",
"\n",
"- `approx_freq_aaa` approximates frequency *data* $(Z, F)$ by a sum of simple poles, using the AAA rational approximation algorithm;\n",
"- `approx_sop_fast` approximates an existing sum-of-poles (SOP) expansion by a (hopefully) smaller one, in a single AAA pass, with an optional nonlinear optimization of the pole locations;\n",
"- `approx_sop_tol` finds the *smallest* sum of poles meeting a prescribed error tolerance, by running the AAA (and, optionally, optimization) pipeline repeatedly.\n",
"\n",
"Each section exercises the stopping criteria (a pole budget `max_n_poles` vs. an error tolerance) and, for the two SOP functions, the `nonlinear_optimization` switch. Since the SOP functions take a pole expansion rather than raw data as input, we build such a trial expansion from Matsubara data using a quick and dirty version of the discrete Lehmann representation [[1]](#references)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ea2f2a95",
"metadata": {},
"outputs": [],
"source": [
"import time\n",
"\n",
"import numpy as np\n",
"import scipy.linalg\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from adapol import approx_freq_aaa, approx_sop_fast, approx_sop_tol\n",
"\n",
"np.set_printoptions(precision=6, suppress=True)\n",
"\n",
"beta = 20.0\n",
"\n",
"# Fermionic Matsubara sampling grid Z_n = i(2n+1)*pi/beta for n = -nmax,...,nmax-1.\n",
"nmax = 50\n",
"n_grid = np.arange(-nmax, nmax)\n",
"Z = 1j * (2 * n_grid + 1) * np.pi / beta\n",
"\n",
"\n",
"def G_exact(z):\n",
" \"\"\"Matsubara Green's function for the semicircular density of states.\"\"\"\n",
" return 2.0 * (z - np.sqrt(z + 1.0) * np.sqrt(z - 1.0))\n",
"\n",
"\n",
"def semicircular(w):\n",
" return 2 * np.sqrt(np.maximum(1 - w**2, 0.0)) / np.pi\n",
"\n",
"\n",
"def eval_sop(z, poles, residues):\n",
" \"\"\"Evaluate the sum of simple poles sum_k residues[k] / (z - poles[k]) at points z.\"\"\"\n",
" return (residues[None, :] / (z[:, None] - poles[None, :])).sum(axis=1)\n",
"\n",
"\n",
"G = G_exact(Z)"
]
},
{
"cell_type": "markdown",
"id": "c3ddfa6a",
"metadata": {},
"source": [
"The Matsubara data looks as follows (the real part vanishes by particle-hole symmetry):"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "2af444e4",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(Z.imag, G.real, 'o-', label=r'Re $G(\\mathrm{i}\\nu_n)$')\n",
"plt.plot(Z.imag, G.imag, 'o-', label=r'Im $G(\\mathrm{i}\\nu_n)$')\n",
"plt.xlabel(r'$\\nu_n$')\n",
"plt.ylabel(r'$G(\\mathrm{i}\\nu_n)$')\n",
"plt.legend()\n",
"plt.title('Matsubara data, semicircular density')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "97b37fb1",
"metadata": {},
"source": [
"## Error metric: the $L^2(\\tau)$ norm\n",
"\n",
"Throughout the notebook we measure the quality of a fit by the normalized imaginary time $L^2$ error\n",
"\n",
"$$\\left\\lVert G_{\\mathrm{fit}} - G \\right\\rVert_{L^2(\\tau)} = \\left( \\frac{1}{\\beta}\\int_0^\\beta \\left| G_{\\mathrm{fit}}(\\tau) - G(\\tau) \\right|^2 d\\tau \\right)^{1/2}.$$\n",
"\n",
"A simple pole at $\\omega$ with unit residue has the Matsubara-frequency kernel\n",
"\n",
"$$K(\\mathrm{i}\\nu_n, \\omega) = \\frac{1}{\\mathrm{i}\\nu_n - \\omega},$$\n",
"\n",
"so that a sum of poles reads $G_{\\mathrm{fit}}(\\mathrm{i}\\nu_n) = \\sum_k R_k \\, K(\\mathrm{i}\\nu_n, p_k)$ at the Matsubara frequencies — this is just the sum-of-poles expansion introduced above, evaluated on the imaginary axis. Its Fourier transform is the imaginary time analytic continuation kernel\n",
"\n",
"$$K(\\tau, \\omega) = -\\frac{e^{-\\tau\\omega}}{1 + e^{-\\beta\\omega}},$$\n",
"\n",
"so a sum of poles correspondingly maps to $G_{\\mathrm{fit}}(\\tau) = \\sum_k R_k \\, K(\\tau, p_k)$ in imaginary time.\n",
"\n",
"Below, the exact $G(\\tau) = \\int \\rho(\\omega)\\, K(\\tau,\\omega)\\, d\\omega$ is computed once by Gauss–Chebyshev quadrature — whose weight function $\\sqrt{1-\\omega^2}$ is exactly the semicircle, so a couple hundred nodes give machine precision — and the $\\tau$ integral defining the norm is evaluated by Gauss–Legendre quadrature on $[0, \\beta]$. The kernel is evaluated separately for positive and negative $\\omega$ to avoid overflowing exponentials."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "1ab69afc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"||G||_L2(tau) = 0.204749\n"
]
}
],
"source": [
"def K_tau(tau, omega):\n",
" \"\"\"Imaginary-time kernel K(tau, omega) = -exp(-tau*omega) / (1 + exp(-beta*omega)),\n",
" evaluated in a numerically stable way for both signs of omega.\"\"\"\n",
" K = np.empty((len(tau), len(omega)))\n",
" p = omega > 0\n",
" m = ~p\n",
" t = tau[:, None]\n",
" K[:, p] = -np.exp(-t * omega[p]) / (1 + np.exp(-beta * omega[p]))\n",
" K[:, m] = -np.exp((beta - t) * omega[m]) / (1 + np.exp(beta * omega[m]))\n",
" return K\n",
"\n",
"\n",
"# Gauss-Legendre quadrature on [0, beta], for the tau integral\n",
"x_gl, w_gl = np.polynomial.legendre.leggauss(200)\n",
"tau_q = 0.5 * beta * (x_gl + 1)\n",
"w_q = 0.5 * beta * w_gl\n",
"\n",
"# Exact G(tau) at the quadrature nodes, by Gauss-Chebyshev (2nd kind) quadrature in omega\n",
"M = 200\n",
"j = np.arange(1, M + 1)\n",
"om_cheb = np.cos(j * np.pi / (M + 1))\n",
"w_cheb = np.pi / (M + 1) * np.sin(j * np.pi / (M + 1))**2\n",
"G_tau = (2 / np.pi) * K_tau(tau_q, om_cheb) @ w_cheb\n",
"\n",
"\n",
"def l2tau_error(poles, residues):\n",
" \"\"\"Normalized L2(tau) norm sqrt( (1/beta) int_0^beta |.|^2 dtau )\n",
" of the difference between a sum of simple poles and the exact G.\"\"\"\n",
" fit_tau = K_tau(tau_q, poles) @ residues\n",
" return np.sqrt(np.sum(w_q * np.abs(fit_tau - G_tau)**2) / beta)\n",
"\n",
"\n",
"print(f\"||G||_L2(tau) = {np.sqrt(np.sum(w_q * np.abs(G_tau)**2) / beta):.6f}\")"
]
},
{
"cell_type": "markdown",
"id": "d270c75b",
"metadata": {},
"source": [
"## Fitting frequency data: approx_freq_aaa\n",
"\n",
"`approx_freq_aaa(F, Z, max_n_poles=..., aaa_tol=...)` takes frequency data $F$ sampled at points $Z$ and returns a pole approximation\n",
"\n",
"$$G_{\\mathrm{fit}}(z) = \\sum_{k=1}^{M} \\frac{R_k}{z - p_k}$$\n",
"\n",
"as a tuple `(poles, residues, error)`. It proceeds in two steps: the AAA algorithm determines the pole locations $p_k$ from the data, and the residues $R_k$ are then obtained by a linear least-squares fit against the data. The returned `error` is the **maximum absolute error at the sample points** — a frequency-domain quantity, unlike the $L^2(\\tau)$ metric defined above, so we report both below.\n",
"\n",
"At least one of the two stopping criteria must be given:\n",
"\n",
"- `max_n_poles` bounds the number of poles;\n",
"- `aaa_tol` is the error tolerance at which the AAA iteration stops.\n",
"\n",
"The number of poles produced is always odd (a consequence of the symmetry-constrained AAA variant used here), and can be smaller than `max_n_poles` if AAA converges early, or if the AAA cleanup step removes poles with negligible residues."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "901e983d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"number of poles: 7\n",
"error at samples: 1.39e-07 (max abs, returned by approx_freq_aaa)\n",
"L2(tau) error: 4.45e-08\n",
"poles: [-0.856458 -0.540225 -0.226438 -0. 0.226438 0.540225 0.856458]\n"
]
}
],
"source": [
"poles, residues, error = approx_freq_aaa(G, Z, aaa_tol=1e-6)\n",
"\n",
"print(f\"number of poles: {len(poles)}\")\n",
"print(f\"error at samples: {error:.2e} (max abs, returned by approx_freq_aaa)\")\n",
"print(f\"L2(tau) error: {l2tau_error(poles, residues):.2e}\")\n",
"print(f\"poles: {np.sort(poles)}\")"
]
},
{
"cell_type": "markdown",
"id": "aa917317",
"metadata": {},
"source": [
"We can verify that the maximum error returned by `approx_freq_aaa` matches the error we obtain directly by evaluating the pole expansion:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "6837c491",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"max error from direct evaluation: 1.391e-07\n"
]
}
],
"source": [
"direct_error = np.max(np.abs(eval_sop(Z, poles, residues) - G))\n",
"\n",
"print(f\"max error from direct evaluation: {direct_error:.3e}\")\n"
]
},
{
"cell_type": "markdown",
"id": "6cdd5dc8",
"metadata": {},
"source": [
"### Stopping criterion 1: a pole budget (max_n_poles)\n",
"\n",
"With `max_n_poles` we prescribe a maximum number of poles and take whatever accuracy that budget buys. Sweeping the budget shows the error converging geometrically with the number of poles, down to a roundoff floor of about $10^{-14}$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "38f120a5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"max_n_poles = 2 -> 1 poles, error at samples = 5.56e-01, L2(tau) error = 9.89e-02\n",
"max_n_poles = 4 -> 3 poles, error at samples = 2.31e-02, L2(tau) error = 3.52e-03\n",
"max_n_poles = 6 -> 5 poles, error at samples = 5.49e-05, L2(tau) error = 1.25e-05\n",
"max_n_poles = 8 -> 7 poles, error at samples = 1.39e-07, L2(tau) error = 4.45e-08\n",
"max_n_poles = 10 -> 9 poles, error at samples = 4.04e-11, L2(tau) error = 9.18e-12\n",
"max_n_poles = 12 -> 11 poles, error at samples = 7.01e-15, L2(tau) error = 3.02e-15\n",
"max_n_poles = 14 -> 11 poles, error at samples = 1.29e-14, L2(tau) error = 9.11e-15\n",
"max_n_poles = 16 -> 11 poles, error at samples = 4.12e-14, L2(tau) error = 5.34e-15\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n_poles, err_freq, err_tau = [], [], []\n",
"for max_np in range(2, 18, 2):\n",
" p, r, err = approx_freq_aaa(G, Z, max_n_poles=max_np)\n",
" n_poles.append(len(p))\n",
" err_freq.append(err)\n",
" err_tau.append(l2tau_error(p, r))\n",
" print(f\"max_n_poles = {max_np:2d} -> {len(p):2d} poles, \"\n",
" f\"error at samples = {err:.2e}, L2(tau) error = {err_tau[-1]:.2e}\")\n",
"\n",
"plt.plot(n_poles, err_freq, 'o-', label='max abs error at samples')\n",
"plt.plot(n_poles, err_tau, 's--', label=r'$L^2(\\tau)$ error')\n",
"plt.yscale('log')\n",
"plt.xlabel('Number of poles')\n",
"plt.ylabel('Error')\n",
"plt.legend()\n",
"plt.title('approx_freq_aaa: convergence with pole count')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "12a423c0",
"metadata": {},
"source": [
"Note that the budget is not always saturated: `max_n_poles = 14` and `16` both return 11 poles, since by then AAA has hit the roundoff floor and additional poles would not improve the fit.\n",
"\n",
"### Stopping criterion 2: an AAA tolerance (aaa_tol)\n",
"\n",
"With `aaa_tol`, the AAA iteration instead runs until its internal error estimate drops below the requested tolerance. We saw this above, but will now sweep over values of this parameter. Note that `aaa_tol` controls the *pole step* only: it bounds the AAA interpolation error on the data, not the final $L^2$ error after the residue refit. It is a stopping heuristic rather than a guarantee (though in this case we see that it acts as an upper bound)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "505b264e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"aaa_tol = 1e-03 -> 5 poles, error at samples = 5.49e-05, L2(tau) error = 1.25e-05\n",
"aaa_tol = 1e-06 -> 7 poles, error at samples = 1.39e-07, L2(tau) error = 4.45e-08\n",
"aaa_tol = 1e-09 -> 9 poles, error at samples = 4.04e-11, L2(tau) error = 9.18e-12\n",
"aaa_tol = 1e-12 -> 11 poles, error at samples = 7.01e-15, L2(tau) error = 3.02e-15\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"tols = [1e-3, 1e-6, 1e-9, 1e-12]\n",
"achieved = []\n",
"for tol in tols:\n",
" p, r, err = approx_freq_aaa(G, Z, aaa_tol=tol)\n",
" achieved.append(err)\n",
" print(f\"aaa_tol = {tol:.0e} -> {len(p):2d} poles, error at samples = {err:.2e}, \"\n",
" f\"L2(tau) error = {l2tau_error(p, r):.2e}\")\n",
"\n",
"plt.plot(tols, achieved, 'o-', label='achieved error at samples')\n",
"plt.plot(tols, tols, 'k--', label='achieved = requested')\n",
"plt.xscale('log')\n",
"plt.yscale('log')\n",
"plt.xlabel('aaa_tol')\n",
"plt.ylabel('Error')\n",
"plt.legend()\n",
"plt.title('approx_freq_aaa: requested vs achieved error')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "daaab07e",
"metadata": {},
"source": [
"### Combining both criteria\n",
"\n",
"If both criteria are given, AAA stops at whichever is reached first — useful for requesting a tolerance while capping the cost."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "e5bc44de",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"max_n_poles=20, aaa_tol=1e-8 -> 9 poles, error = 4.04e-11 (tolerance reached first)\n",
"AAA: Warning! Failed to converge after 3 steps. Final residual 7.94E-05 larger than tolerance 1.00E-12.\n",
"max_n_poles=6, aaa_tol=1e-12 -> 5 poles, error = 5.49e-05 (pole budget reached first)\n"
]
}
],
"source": [
"p, r, err = approx_freq_aaa(G, Z, max_n_poles=20, aaa_tol=1e-8)\n",
"print(f\"max_n_poles=20, aaa_tol=1e-8 -> {len(p)} poles, error = {err:.2e} (tolerance reached first)\")\n",
"\n",
"p, r, err = approx_freq_aaa(G, Z, max_n_poles=6, aaa_tol=1e-12)\n",
"print(f\"max_n_poles=6, aaa_tol=1e-12 -> {len(p)} poles, error = {err:.2e} (pole budget reached first)\")"
]
},
{
"cell_type": "markdown",
"id": "a1c77f51",
"metadata": {},
"source": [
"## Compressing a sum of poles: approx_sop_fast\n",
"\n",
"The remaining two functions, `approx_sop_fast` and `approx_sop_tol`, take as input not raw data but an existing *sum-of-poles expansion* — poles $P_k$ and weights $W_k$ with $G(z) \\approx \\sum_{k} W_k / (z - P_k)$. So we first build a trial expansion from the Matsubara data, then compress it.\n",
"\n",
"For the trial expansion, we use a quick and dirty version of the discrete Lehmann representation. On a dense candidate grid $\\omega_j \\in [-1, 1]$, we form the matrix of simple poles $K_{nj} = 1/(\\mathrm{i}\\nu_n - \\omega_j)$ sampled on a dense Matsubara grid, keep a well-conditioned subset of poles via a rank-revealing column-pivoted QR (those with $|R_{jj}| > \\epsilon_{\\mathrm{QR}}|R_{00}|$), and fit their weights to the data for $G$ by least squares."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "e194a2d7",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"selected 21 of 501 candidate poles\n",
"max |Im weights|: 2.9e-07\n",
"max |G - SOP| on the dense grid: 1.9e-14\n",
"L2(tau) error of the input SOP: 2.0e-15\n"
]
}
],
"source": [
"# dense fermionic Matsubara grid and dense candidate pole grid\n",
"n_dense = np.arange(-200, 200)\n",
"Z_dense = 1j * (2 * n_dense + 1) * np.pi / beta\n",
"G_dense = G_exact(Z_dense)\n",
"\n",
"om_cand = np.linspace(-1, 1, 501)\n",
"Kmat = 1.0 / (Z_dense[:, None] - om_cand[None, :])\n",
"\n",
"# rank-revealing pivoted QR selects the linearly independent candidate poles\n",
"qr_tol = 1e-12\n",
"Q, R, piv = scipy.linalg.qr(Kmat, mode='economic', pivoting=True)\n",
"rank = int(np.sum(np.abs(np.diag(R)) > qr_tol * np.abs(R[0, 0])))\n",
"poles_in = np.sort(om_cand[piv[:rank]])\n",
"\n",
"# least-squares fit for the weights\n",
"Ksel = 1.0 / (Z_dense[:, None] - poles_in[None, :])\n",
"weights_in, *_ = np.linalg.lstsq(Ksel, G_dense, rcond=None)\n",
"\n",
"print(f\"selected {rank} of {len(om_cand)} candidate poles\")\n",
"print(f\"max |Im weights|: {np.abs(weights_in.imag).max():.1e}\")\n",
"weights_in = weights_in.real # discard the tiny imaginary parts (see below)\n",
"\n",
"print(f\"max |G - SOP| on the dense grid: {np.abs(eval_sop(Z_dense, poles_in, weights_in) - G_dense).max():.1e}\")\n",
"print(f\"L2(tau) error of the input SOP: {l2tau_error(poles_in, weights_in):.1e}\")"
]
},
{
"cell_type": "markdown",
"id": "564a5448",
"metadata": {},
"source": [
"These 21 poles reproduce $G$ to near machine precision, so the expansion is a faithful proxy whose $L^2(\\tau)$ error, $\\sim 2\\times10^{-15}$, is the accuracy floor for everything below. However, even at this tight tolerance, it is known that a fit using poles *customized to $G$ itself* outperforms the generic DLR-style construction, in the sense that the same accuracy can be achieved with fewer poles [[2]](#references). `approx_sop_fast` and `approx_sop_tol` attempt to compress a given pole approximation.\n",
"\n",
"`approx_sop_fast(poles, residues, beta, ...)` does this in a single pass:\n",
"\n",
"1. **Pole step:** the input expansion is evaluated on a symmetric fermionic Matsubara grid $\\mathrm{i}(2n+1)\\pi/\\beta$, $-n_{\\max} \\le n \\le n_{\\max}-1$, with $n_{\\max} = \\lfloor 4\\beta\\omega_{\\max}/\\pi\\rfloor + 1$ and $\\omega_{\\max} = \\max_k |P_k|$ (here $n_{\\max} = 26$, i.e. 52 points; a custom grid can be supplied via the `Z` argument), and AAA is run on this data to determine the new pole locations.\n",
"2. **Residue step:** the residues are determined by minimizing the $L^2(\\tau)$ norm of the deviation from the input expansion — a *linear* least-squares problem once the poles are fixed. With `nonlinear_optimization=True`, the pole locations and residues are instead optimized *jointly* (L-BFGS-B with an analytic gradient, keeping the number of poles fixed), using the AAA poles as the initial guess.\n",
"\n",
"Note that the residue fit happens in imaginary time — unlike `approx_freq_aaa`, which fits the residues to the frequency data. Accordingly, the returned `error` is the $L^2(\\tau)$ norm of the deviation from the *input expansion*; it is the same norm as our metric, which instead measures against the exact $G$.\n",
"\n",
"The stopping criteria `max_n_poles` and `aaa_tol` work exactly as in `approx_freq_aaa`, with `aaa_tol` again controlling only the AAA pole step."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "338e9b22",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"21 input poles compressed to 9 poles in 4.7 ms\n",
"returned L2(tau) error (vs the input SOP): 9.37e-12\n",
"L2(tau) error vs the exact G: 9.37e-12\n"
]
}
],
"source": [
"t0 = time.perf_counter()\n",
"p, r, err = approx_sop_fast(poles_in, weights_in, beta, aaa_tol=1e-8)\n",
"dt = (time.perf_counter() - t0) * 1e3\n",
"\n",
"print(f\"{len(poles_in)} input poles compressed to {len(p)} poles in {dt:.1f} ms\")\n",
"print(f\"returned L2(tau) error (vs the input SOP): {err:.2e}\")\n",
"print(f\"L2(tau) error vs the exact G: {l2tau_error(p, r):.2e}\")"
]
},
{
"cell_type": "markdown",
"id": "576b59e7",
"metadata": {},
"source": [
"### Convergence, and the effect of nonlinear optimization\n",
"\n",
"We sweep the pole budget with the nonlinear optimization switched off and on. (Timings here and below are rough single-shot measurements.)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "3b575466",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nonlinear_optimization=False max_n_poles= 2 -> 1 poles, error vs exact G = 9.88e-02, t = 4.7 ms\n",
"nonlinear_optimization=False max_n_poles= 4 -> 3 poles, error vs exact G = 3.52e-03, t = 2.0 ms\n",
"nonlinear_optimization=False max_n_poles= 6 -> 5 poles, error vs exact G = 1.24e-05, t = 1.0 ms\n",
"nonlinear_optimization=False max_n_poles= 8 -> 7 poles, error vs exact G = 4.44e-08, t = 7.8 ms\n",
"nonlinear_optimization=False max_n_poles=10 -> 9 poles, error vs exact G = 9.37e-12, t = 3.9 ms\n",
"nonlinear_optimization=False max_n_poles=12 -> 11 poles, error vs exact G = 3.56e-15, t = 4.7 ms\n",
"nonlinear_optimization=False max_n_poles=14 -> 13 poles, error vs exact G = 2.01e-15, t = 8.5 ms\n",
"nonlinear_optimization=False max_n_poles=16 -> 13 poles, error vs exact G = 2.00e-15, t = 2.3 ms\n",
"\n",
"nonlinear_optimization=True max_n_poles= 2 -> 1 poles, error vs exact G = 9.88e-02, t = 1.1 ms\n",
"nonlinear_optimization=True max_n_poles= 4 -> 3 poles, error vs exact G = 2.09e-03, t = 1.7 ms\n",
"nonlinear_optimization=True max_n_poles= 6 -> 5 poles, error vs exact G = 9.25e-06, t = 3.6 ms\n",
"nonlinear_optimization=True max_n_poles= 8 -> 7 poles, error vs exact G = 1.20e-08, t = 11.6 ms\n",
"nonlinear_optimization=True max_n_poles=10 -> 9 poles, error vs exact G = 9.27e-12, t = 5.0 ms\n",
"nonlinear_optimization=True max_n_poles=12 -> 11 poles, error vs exact G = 3.56e-15, t = 6.9 ms\n",
"nonlinear_optimization=True max_n_poles=14 -> 13 poles, error vs exact G = 2.01e-15, t = 7.5 ms\n",
"nonlinear_optimization=True max_n_poles=16 -> 13 poles, error vs exact G = 2.00e-15, t = 4.8 ms\n",
"\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"results = {}\n",
"for nlo in [False, True]:\n",
" ns, errs, times = [], [], []\n",
" for max_np in range(2, 18, 2):\n",
" t0 = time.perf_counter()\n",
" p, r, err = approx_sop_fast(poles_in, weights_in, beta,\n",
" max_n_poles=max_np, nonlinear_optimization=nlo)\n",
" dt = (time.perf_counter() - t0) * 1e3\n",
" ns.append(len(p)); errs.append(err); times.append(dt)\n",
" print(f\"nonlinear_optimization={str(nlo):5s} max_n_poles={max_np:2d} -> \"\n",
" f\"{len(p):2d} poles, error vs exact G = {l2tau_error(p, r):.2e}, \"\n",
" f\"t = {dt:5.1f} ms\")\n",
" results[nlo] = (ns, errs, times)\n",
" print()\n",
"\n",
"plt.plot(results[False][0], results[False][1], 'o-', label='linear residue fit')\n",
"plt.plot(results[True][0], results[True][1], 's--', label='nonlinear optimization')\n",
"plt.yscale('log')\n",
"plt.xlabel('Number of poles')\n",
"plt.ylabel(r'$L^2(\\tau)$ error')\n",
"plt.legend()\n",
"plt.title('approx_sop_fast: convergence of the compression')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "19ea4303",
"metadata": {},
"source": [
"For this smooth, particle-hole-symmetric density the AAA poles are already nearly optimal, so we see there is not much to gain from running the nonlinear optimization step. In general, it is worth experimenting with."
]
},
{
"cell_type": "markdown",
"id": "e6569619",
"metadata": {},
"source": [
"## Compressing to a prescribed accuracy: approx_sop_tol\n",
"\n",
"For `approx_sop_fast`, the tolerance `aaa_tol` is only a stopping heuristic for the pole step. `approx_sop_tol(poles, residues, tol, beta, ...)` instead imposes the tolerance on the **final** $L^2(\\tau)$ error, and returns the *smallest* pole expansion meeting it. Since the final error is only known after the residue fit, the AAA + residue-fit pipeline is run repeatedly: the number of AAA steps is increased until the error drops below `tol`, and a bisection then locates the smallest number of steps that still meets it. If the tolerance cannot be met (for instance, because it lies below the roundoff floor), a `ValueError` is raised.\n",
"\n",
"The search costs several `approx_sop_fast`-like passes, as the timing comparison below shows. With `verbose=True` the progress of the search is printed:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "02e975f2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Adapol: compressing 21 poles, tol 1.00E-06 on the imaginary time L2 error\n",
" beta = 20, 52 Matsubara points, joint pole and residue optimization\n",
"Adapol: [1/2 bracket] growing the AAA step count until the tolerance is met\n",
" 4 steps -> 7 poles, error 1.20E-08 <= tol ok AAA residual 1.36E-07\n",
"Adapol: [2/2 bisect] smallest step count still meeting the tolerance, candidates 1..4\n",
" 2 steps -> 3 poles, error 2.09E-03 > tol fail AAA residual 1.67E-02 candidates 3..4\n",
" 3 steps -> 5 poles, error 9.25E-06 > tol fail AAA residual 7.87E-05 candidates 4..4\n",
"Adapol: done -- 7 poles (4 AAA steps), error 1.20E-08 <= tol 1.00E-06, 3 pipeline passes, 3 AAA runs\n",
"\n",
"approx_sop_tol(tol=1e-6): 7 poles, error = 1.20e-08, t = 14.4 ms\n",
"approx_sop_fast(aaa_tol=1e-6): 7 poles, error = 1.20e-08, t = 7.0 ms\n"
]
}
],
"source": [
"t0 = time.perf_counter()\n",
"p, r, err = approx_sop_tol(poles_in, weights_in, tol=1e-6, beta=beta,\n",
" nonlinear_optimization=True, verbose=True)\n",
"dt_tol = (time.perf_counter() - t0) * 1e3\n",
"print(f\"\\napprox_sop_tol(tol=1e-6): {len(p)} poles, error = {err:.2e}, t = {dt_tol:5.1f} ms\")\n",
"\n",
"# a single fast pass reaching the same accuracy, for timing comparison\n",
"t0 = time.perf_counter()\n",
"p2, r2, err2 = approx_sop_fast(poles_in, weights_in, beta, aaa_tol=1e-6,\n",
" nonlinear_optimization=True)\n",
"dt_fast = (time.perf_counter() - t0) * 1e3\n",
"print(f\"approx_sop_fast(aaa_tol=1e-6): {len(p2)} poles, error = {err2:.2e}, t = {dt_fast:5.1f} ms\")"
]
},
{
"cell_type": "markdown",
"id": "2840dd08",
"metadata": {},
"source": [
"The log shows the search: a first upward sweep in the number of AAA steps until the tolerance is met, followed by re-runs that bisect the step count down to the smallest sufficient one. The result costs a few times more than a single `approx_sop_fast` pass, but comes with an actual guarantee on the final error rather than a heuristic.\n",
"\n",
"### Pole count vs tolerance\n",
"\n",
"Sweeping the tolerance maps out how many poles a given accuracy requires — a staircase, since poles are added in symmetric pairs. In contrast to the heuristic `aaa_tol` of `approx_sop_fast`, the achieved error now always sits at or below the requested tolerance:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "0f0c6167",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"nonlinear_optimization=False tol = 1e-02 -> 3 poles, error = 3.52e-03, t = 3.4 ms\n",
"nonlinear_optimization=False tol = 1e-03 -> 5 poles, error = 1.24e-05, t = 2.7 ms\n",
"nonlinear_optimization=False tol = 1e-04 -> 5 poles, error = 1.24e-05, t = 2.5 ms\n",
"nonlinear_optimization=False tol = 1e-05 -> 7 poles, error = 4.44e-08, t = 2.6 ms\n",
"nonlinear_optimization=False tol = 1e-06 -> 7 poles, error = 4.44e-08, t = 3.2 ms\n",
"nonlinear_optimization=False tol = 1e-07 -> 7 poles, error = 4.44e-08, t = 4.0 ms\n",
"nonlinear_optimization=False tol = 1e-08 -> 9 poles, error = 9.37e-12, t = 4.1 ms\n",
"nonlinear_optimization=False tol = 1e-09 -> 9 poles, error = 9.37e-12, t = 3.9 ms\n",
"nonlinear_optimization=False tol = 1e-10 -> 9 poles, error = 9.37e-12, t = 3.6 ms\n",
"nonlinear_optimization=False tol = 1e-11 -> 9 poles, error = 9.37e-12, t = 4.1 ms\n",
"nonlinear_optimization=False tol = 1e-12 -> 11 poles, error = 3.02e-15, t = 4.7 ms\n",
"\n",
"nonlinear_optimization=True tol = 1e-02 -> 3 poles, error = 2.09e-03, t = 5.3 ms\n",
"nonlinear_optimization=True tol = 1e-03 -> 5 poles, error = 9.25e-06, t = 4.7 ms\n",
"nonlinear_optimization=True tol = 1e-04 -> 5 poles, error = 9.25e-06, t = 5.0 ms\n",
"nonlinear_optimization=True tol = 1e-05 -> 5 poles, error = 9.25e-06, t = 10.5 ms\n",
"nonlinear_optimization=True tol = 1e-06 -> 7 poles, error = 1.20e-08, t = 10.1 ms\n",
"nonlinear_optimization=True tol = 1e-07 -> 7 poles, error = 1.20e-08, t = 12.1 ms\n",
"nonlinear_optimization=True tol = 1e-08 -> 9 poles, error = 9.27e-12, t = 12.6 ms\n",
"nonlinear_optimization=True tol = 1e-09 -> 9 poles, error = 9.27e-12, t = 12.4 ms\n",
"nonlinear_optimization=True tol = 1e-10 -> 9 poles, error = 9.27e-12, t = 12.2 ms\n",
"nonlinear_optimization=True tol = 1e-11 -> 9 poles, error = 9.27e-12, t = 14.2 ms\n",
"nonlinear_optimization=True tol = 1e-12 -> 11 poles, error = 3.96e-15, t = 14.3 ms\n",
"\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"tols = 10.0 ** (-np.arange(2, 13, dtype=float))\n",
"stair = {}\n",
"for nlo in [False, True]:\n",
" ns, errs, times = [], [], []\n",
" for tol in tols:\n",
" t0 = time.perf_counter()\n",
" p, r, err = approx_sop_tol(poles_in, weights_in, tol=tol, beta=beta,\n",
" nonlinear_optimization=nlo)\n",
" dt = (time.perf_counter() - t0) * 1e3\n",
" ns.append(len(p)); errs.append(err); times.append(dt)\n",
" print(f\"nonlinear_optimization={str(nlo):5s} tol = {tol:.0e} -> \"\n",
" f\"{len(p):2d} poles, error = {err:.2e}, t = {dt:5.1f} ms\")\n",
" stair[nlo] = (ns, errs, times)\n",
" print()\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 10))\n",
"ax1.plot(tols, stair[False][0], 'o-', label='linear residue fit')\n",
"ax1.plot(tols, stair[True][0], 's--', label='nonlinear optimization')\n",
"ax1.set_xscale('log')\n",
"ax1.invert_xaxis()\n",
"ax1.set_xlabel('tol')\n",
"ax1.set_ylabel('Number of poles')\n",
"ax1.legend()\n",
"ax1.set_title('poles required vs tolerance')\n",
"\n",
"ax2.plot(tols, stair[False][1], 'o-', label='linear residue fit')\n",
"ax2.plot(tols, stair[True][1], 's--', label='nonlinear optimization')\n",
"ax2.plot(tols, tols, 'k:', label='error = tol')\n",
"ax2.set_xscale('log')\n",
"ax2.set_yscale('log')\n",
"ax2.invert_xaxis()\n",
"ax2.set_xlabel('tol')\n",
"ax2.set_ylabel(r'achieved $L^2(\\tau)$ error')\n",
"ax2.legend()\n",
"ax2.set_title('the tolerance is respected')\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "fa43b3c3",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"We also observe that the nonlinear optimization can sometimes produce a more compact pole approximation at the same tolerance."
]
},
{
"cell_type": "markdown",
"id": "926eaebe",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## Summary\n",
"\n",
"- `approx_freq_aaa` fits raw Matsubara (or other frequency) data directly; its error is measured at the sample points.\n",
"- `approx_sop_fast` is a cheap single-pass compression of an existing pole expansion; `max_n_poles` acts as a budget, `aaa_tol` as a stopping heuristic for the pole search.\n",
"- `approx_sop_tol` returns the smallest expansion meeting a prescribed final $L^2(\\tau)$ accuracy, at the cost of a few fast passes.\n",
"- `nonlinear_optimization=True` refines the AAA pole locations, if possible."
]
},
{
"cell_type": "markdown",
"id": "3b8ae4c4",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## References\n",
"\n",
"[1] J. Kaye, K. Chen, and O. Parcollet, [Discrete Lehmann representation of imaginary time Green's functions](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.235115), Phys. Rev. B **105**, 235115 (2022).\n",
"\n",
"[2] Z. Huang, D. Golež, H. U. R. Strand, J. Kaye, [Automated evaluation of imaginary time strong coupling diagrams by sum-of-exponentials hybridization fitting](https://scipost.org/SciPostPhys.19.5.121), SciPost Phys., 19 (5), 121 (2025)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}