{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "5308b16e-a9e6-462f-8ee9-c0ef9ba19fe0",
"metadata": {
"editable": true,
"metadata": {
"nbsphinx": "hidden"
},
"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": "5f3fd2e7",
"metadata": {},
"source": [
"# Pole approximation of multi-orbital data with a discrete spectrum\n",
"\n",
"In this example we consider a matrix-valued (multi-orbital) Matsubara function with a discrete spectrum: a sum of $K$ delta functions at random locations $p_k \\in [-1, 1]$, with random rank-one (bath-site-like) residues,\n",
"\n",
"$$F(\\mathrm{i}\\nu_n) = \\sum_{k=1}^{K} \\frac{v_k^{\\phantom\\dagger} v_k^\\dagger}{\\mathrm{i}\\nu_n - p_k}, \\qquad v_k \\in \\mathbb{C}^{n_{\\mathrm{orb}}},$$\n",
"\n",
"with the $p_k$ and $v_k$ drawn from a fixed random seed. Since the exact function is itself a sum of simple poles, the accuracy of any fit can be measured exactly. Note that the goal is *not* to recover the input poles: at a given target accuracy the fitting routines are free to use fewer — and differently placed — poles than the input, and for a compact representation that is exactly what one wants.\n",
"\n",
"We briefly demonstrate the three main `adapol` functions with fixed, representative parameter choices — see the `semicircle` example for a more detailed exploration of the stopping criteria, the `nonlinear_optimization` switch, and the error metric. We then run a small experiment: using `approx_sop_tol`, we measure how many poles are actually required to reach a given accuracy as the number of orbitals grows."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4e6fb011",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"model: 12 poles, 4 orbitals, data shape (100, 4, 4)\n",
"pole locations: [-0.995 -0.967 -0.918 -0.46 0.087 0.213 0.274 0.459 0.627 0.632\n",
" 0.826 0.87 ]\n"
]
}
],
"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=3, suppress=True)\n",
"\n",
"beta = 20.0\n",
"\n",
"# Fermionic Matsubara sampling grid Z_n = i(2n+1)*pi/beta, symmetric about zero\n",
"nmax = 50\n",
"Z = 1j * (2 * np.arange(-nmax, nmax) + 1) * np.pi / beta\n",
"\n",
"\n",
"def make_random_model(n_orb, n_delta, seed_res):\n",
" \"\"\"Random discrete model: n_delta poles uniform on [-1,1] (fixed seed), rank-one\n",
" residues v_k v_k^dag with complex Gaussian v_k, normalized so tr F ~ O(1).\"\"\"\n",
" rng_poles = np.random.default_rng(0)\n",
" rng_res = np.random.default_rng(seed_res)\n",
" poles = np.sort(rng_poles.uniform(-1, 1, n_delta))\n",
" v = (rng_res.standard_normal((n_delta, n_orb))\n",
" + 1j * rng_res.standard_normal((n_delta, n_orb))) / np.sqrt(2 * n_delta)\n",
" residues = np.einsum('ki,kj->kij', v, v.conj())\n",
" return poles, residues\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",
" C = 1.0 / (z[:, None] - poles[None, :])\n",
" return np.einsum('zp,p...->z...', C, residues)\n",
"\n",
"\n",
"n_orb, n_delta = 4, 12\n",
"poles_ex, res_ex = make_random_model(n_orb, n_delta, seed_res=1)\n",
"F = eval_sop(Z, poles_ex, res_ex)\n",
"\n",
"print(f\"model: {n_delta} poles, {n_orb} orbitals, data shape {F.shape}\")\n",
"print(f\"pole locations: {poles_ex}\")"
]
},
{
"cell_type": "markdown",
"id": "a9fa6391",
"metadata": {},
"source": [
"We plot a couple of matrix elements of the data."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b5429b04",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(Z.imag, F[:, 0, 0].imag, 'o-', markersize=3, label=r'Im $F_{00}$')\n",
"plt.plot(Z.imag, F[:, 0, 1].real, 'o-', markersize=3, label=r'Re $F_{01}$')\n",
"plt.plot(Z.imag, F[:, 0, 1].imag, 'o-', markersize=3, label=r'Im $F_{01}$')\n",
"plt.xlabel(r'$\\nu_n$')\n",
"plt.ylabel(r'$F(\\mathrm{i}\\nu_n)$')\n",
"plt.legend()\n",
"plt.title(f'Matsubara data, {n_delta} poles, {n_orb} orbitals')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "b0a1c627",
"metadata": {},
"source": [
"We measure fit quality by the normalized imaginary-time $L^2$ error, with the Frobenius norm over the orbital indices,\n",
"\n",
"$$\\left\\lVert F_{\\mathrm{fit}} - F \\right\\rVert_{L^2(\\tau)} = \\left( \\frac{1}{\\beta}\\int_0^\\beta \\sum_{ij} \\left| \\big(F_{\\mathrm{fit}}\\big)_{ij}(\\tau) - F_{ij}(\\tau) \\right|^2 d\\tau \\right)^{1/2},$$\n",
"\n",
"which is the same norm `approx_sop_fast` and `approx_sop_tol` report. Since the exact function is a known sum of poles, we can evaluate it in imaginary time analytically and compute the error by Gauss–Legendre quadrature (see the `semicircle` example):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d0f3e2d5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"||F||_L2(tau) = 0.3440\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",
"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",
"F_tau = np.einsum('tp,pij->tij', K_tau(tau_q, poles_ex), res_ex)\n",
"\n",
"\n",
"def l2tau_error(poles, residues):\n",
" \"\"\"Normalized Frobenius L2(tau) error of a sum of simple poles vs the exact model.\"\"\"\n",
" fit_tau = np.einsum('tp,p...->t...', K_tau(tau_q, poles), residues)\n",
" d2 = np.abs(fit_tau - F_tau) ** 2\n",
" return np.sqrt(np.sum(w_q * d2.reshape(len(tau_q), -1).sum(axis=1)) / beta)\n",
"\n",
"\n",
"print(f\"||F||_L2(tau) = {np.sqrt(np.einsum('t,tij->', w_q, np.abs(F_tau)**2) / beta):.4f}\")"
]
},
{
"cell_type": "markdown",
"id": "aa309d92",
"metadata": {},
"source": [
"## Fitting the Matsubara data: approx_freq_aaa\n",
"\n",
"`approx_freq_aaa` fits the raw frequency data $(F, Z)$; matrix-valued data of shape `(N, n_orb, n_orb)` is supported directly, and a single set of poles is found for all matrix elements at once. We try a few representative settings. At a moderate tolerance the data — generated by 12 poles — is fit with *fewer* poles: 9 suffice for $\\sim 10^{-6}$ accuracy. A pole budget `max_n_poles` instead caps the cost directly. The fitted residue matrices come out numerically Hermitian."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6d506e8e",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"aaa_tol=0.0001 -> 9 poles, error at samples = 4.2e-06, L2(tau) error = 1.5e-06\n",
"max_n_poles=8 -> 7 poles, error at samples = 5.7e-04, L2(tau) error = 1.3e-04\n",
"aaa_tol=1e-10 -> 13 poles, error at samples = 1.1e-15, L2(tau) error = 2.4e-16\n"
]
}
],
"source": [
"for kwargs in [dict(aaa_tol=1e-4),\n",
" dict(max_n_poles=8),\n",
" dict(aaa_tol=1e-10)]:\n",
" p, r, err = approx_freq_aaa(F, Z, **kwargs)\n",
" opts = ', '.join(f'{k}={v}' for k, v in kwargs.items())\n",
" print(f\"{opts:20s} -> {len(p):2d} poles, error at samples = {err:.1e}, \"\n",
" f\"L2(tau) error = {l2tau_error(p, r):.1e}\")"
]
},
{
"cell_type": "markdown",
"id": "4bb0153a",
"metadata": {},
"source": [
"## Compressing the exact expansion: approx_sop_fast\n",
"\n",
"The two SOP routines take the exact poles and residues as input directly. `approx_sop_fast` makes a single AAA pass (with the residues refit in imaginary time), controlled by `aaa_tol` and/or `max_n_poles`; `nonlinear_optimization=True` additionally refines the pole locations, though here the AAA poles are already essentially optimal and it adds only cost:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "80b985f9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"aaa_tol=0.0001 -> 9 poles, L2(tau) error = 1.5e-06, t = 8.9 ms\n",
"max_n_poles=8 -> 7 poles, L2(tau) error = 1.3e-04, t = 5.3 ms\n",
"max_n_poles=8, nonlinear_optimization=True -> 7 poles, L2(tau) error = 1.3e-04, t = 8.9 ms\n"
]
}
],
"source": [
"for kwargs in [dict(aaa_tol=1e-4),\n",
" dict(max_n_poles=8),\n",
" dict(max_n_poles=8, nonlinear_optimization=True)]:\n",
" t0 = time.perf_counter()\n",
" p, r, err = approx_sop_fast(poles_ex, res_ex, beta, **kwargs)\n",
" dt = (time.perf_counter() - t0) * 1e3\n",
" opts = ', '.join(f'{k}={v}' for k, v in kwargs.items())\n",
" print(f\"{opts:45s} -> {len(p):2d} poles, L2(tau) error = {err:.1e}, t = {dt:5.1f} ms\")"
]
},
{
"cell_type": "markdown",
"id": "2d49780f",
"metadata": {},
"source": [
"## Compression to a prescribed accuracy: approx_sop_tol\n",
"\n",
"`approx_sop_tol` finds the smallest expansion whose final $L^2(\\tau)$ error is below `tol` — so the tolerance directly controls the pole count. Only near roundoff does the required count approach that of the input:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "81f8dd2e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tol = 1e-03 -> 7 poles, L2(tau) error = 1.3e-04\n",
"tol = 1e-06 -> 11 poles, L2(tau) error = 5.1e-09\n",
"tol = 1e-10 -> 13 poles, L2(tau) error = 3.6e-16\n"
]
}
],
"source": [
"for tol in [1e-3, 1e-6, 1e-10]:\n",
" p, r, err = approx_sop_tol(poles_ex, res_ex, tol=tol, beta=beta, verbose=False)\n",
" print(f\"tol = {tol:.0e} -> {len(p):2d} poles, L2(tau) error = {err:.1e}\")"
]
},
{
"cell_type": "markdown",
"id": "4421e5f5",
"metadata": {},
"source": [
"The compressed poles are not simply a subset of the input poles — they are placed wherever the requested accuracy demands. At `tol = 1e-3`, seven poles reproduce the 12-pole function: the cluster of three input poles near $\\omega \\approx -0.95$ is replaced by a *single* pole carrying their combined spectral weight (residue trace $0.50 \\approx 0.19 + 0.18 + 0.13$), the near-degenerate pair at $\\omega \\approx 0.63$ is likewise merged, and the remaining poles shift away from their input locations — yet the function is reproduced to the requested accuracy:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "9c7398c0",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"p7, r7, err7 = approx_sop_tol(poles_ex, res_ex, tol=1e-3, beta=beta, verbose=False)\n",
"idx = np.argsort(p7)\n",
"p7s, tr7 = p7[idx], np.einsum('kii->k', r7[idx]).real\n",
"tr_ex = np.einsum('kii->k', res_ex).real\n",
"\n",
"plt.vlines(poles_ex, 0, tr_ex, color='C0', alpha=0.5)\n",
"plt.plot(poles_ex, tr_ex, 'C0o', label=f'exact ({n_delta} poles)')\n",
"plt.vlines(p7s, 0, tr7, color='C1', alpha=0.8)\n",
"plt.plot(p7s, tr7, 'C1s', markersize=5, label=f'approx_sop_tol(tol=1e-3), {len(p7s)} poles')\n",
"plt.xlabel(r'$\\omega$')\n",
"plt.ylabel(r'tr $R_k$')\n",
"plt.legend()\n",
"plt.title(f'Pole locations and residue traces (L2(tau) error = {err7:.1e})')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7945fcaa",
"metadata": {},
"source": [
"## How does the required number of poles vary with the number of orbitals?\n",
"\n",
"To test this, we fix $K = 40$ delta functions (same random locations throughout) and draw independent random residue vectors for increasing $n_{\\mathrm{orb}}$, then ask `approx_sop_tol` for the minimal pole count at fixed tolerances.\n",
"\n",
"For a single orbital, many of the 40 poles can be merged. As $n_{\\mathrm{orb}}$ grows, the data contains more independent linear combinations of the poles and the required count increases — but it quickly *saturates*, far below 40. The saturation level is the numerical $\\epsilon$-rank of the Matsubara kernel $1/(\\mathrm{i}\\nu_n - \\omega)$ on $\\omega \\in [-1,1]$ (dashed lines, computed with the pivoted-QR construction from the `semicircle` example): a universal pole set of that size can represent *any* function with this spectral support to accuracy $\\epsilon$, no matter how many orbitals. This is the concept underlying the discrete Lehmann representation [[1]](#references): in the limit of a large number of orbitals, the required pole count is set by the spectral support and the target accuracy — not by the number of orbitals (nor by the $n_{\\mathrm{orb}}^2$ matrix elements, nor the 40 input poles)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "a3d70dba",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tol = 1e-04: n_orb = [1, 2, 3, 4, 6, 8, 12, 16, 20, 24] -> poles = [5, 7, 9, 9, 9, 9, 9, 11, 11, 11]\n",
"tol = 1e-06: n_orb = [1, 2, 3, 4, 6, 8, 12, 16, 20, 24] -> poles = [7, 11, 11, 13, 13, 13, 13, 13, 13, 13]\n",
"kernel epsilon-ranks: {0.0001: 11, 1e-06: 13}\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n_delta_sweep = 40\n",
"tols = [1e-4, 1e-6]\n",
"n_orbs = [1, 2, 3, 4, 6, 8, 12, 16, 20, 24]\n",
"\n",
"counts = {}\n",
"for tol in tols:\n",
" counts[tol] = []\n",
" for no in n_orbs:\n",
" poles_s, res_s = make_random_model(no, n_delta_sweep, seed_res=100 + no)\n",
" p, r, err = approx_sop_tol(poles_s, res_s, tol=tol, beta=beta, verbose=False)\n",
" counts[tol].append(len(p))\n",
" print(f\"tol = {tol:.0e}: n_orb = {n_orbs} -> poles = {counts[tol]}\")\n",
"\n",
"# epsilon-rank of the continuum Matsubara kernel on [-1, 1] (cf. the semicircle example)\n",
"Z_dense = 1j * (2 * np.arange(-200, 200) + 1) * np.pi / beta\n",
"om_cand = np.linspace(-1, 1, 501)\n",
"_, R, _ = scipy.linalg.qr(1.0 / (Z_dense[:, None] - om_cand[None, :]),\n",
" mode='economic', pivoting=True)\n",
"diagR = np.abs(np.diag(R))\n",
"ranks = {tol: int(np.sum(diagR > tol * diagR[0])) for tol in tols}\n",
"print(f\"kernel epsilon-ranks: {ranks}\")\n",
"\n",
"for i, tol in enumerate(tols):\n",
" plt.plot(n_orbs, counts[tol], f'C{i}o-', label=f'tol = {tol:.0e}')\n",
" plt.axhline(ranks[tol], color=f'C{i}', linestyle='--', alpha=0.6,\n",
" label=f'kernel rank at {tol:.0e}')\n",
"plt.xlabel('Number of orbitals')\n",
"plt.ylabel('Number of poles required')\n",
"plt.legend()\n",
"plt.title(f'approx_sop_tol: pole count vs orbitals ({n_delta_sweep} input poles)')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "0cad5cb3",
"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)."
]
}
],
"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
}