Hybrid Analytic Continuation Algorithm

• Hybrid analytic continuation algorithms are advanced methods that combine techniques like barycentric rational interpolation, sparse modeling, and conformal mapping to reconstruct spectral functions from imaginary-time data.
• They integrate complementary numerical schemes and physical constraints to stabilize the ill-posed inverse problem while preserving analyticity, causality, symmetry, and positivity.
• These methods offer significant computational advantages, achieving speedups up to 100× over traditional techniques in applications such as GW self-energy and multi-orbital spectral analysis.

In recent work, the label hybrid analytic continuation algorithm is used for continuation schemes that blend components such as barycentric rational interpolation with the adaptive Antoulas–Anderson algorithm, sparse modeling with semi-positive definiteness constraints, sparse modeling with Padé approximation, contour deformation with analytic continuation of the screened Coulomb interaction $W$, or conformal maps with constrained Nevanlinna–Pick interpolation (Huang et al., 2024, Motoyama et al., 2024, Motoyama et al., 2021, Duchemin et al., 2019, Bergamaschi et al., 2023). Across these variants, the common objective is to reconstruct real-frequency response functions from imaginary-time or Matsubara data despite the ill-posedness of the inverse problem, the presence of stochastic noise, and the need to preserve analyticity, causality, symmetry, normalization, or positivity when those constraints are physically mandated.

1. Formal problem and analytic structure

Analytic continuation in many-body physics starts from Euclidean or Matsubara correlators and seeks the corresponding retarded real-frequency objects. For fermionic Matsubara Green’s functions sampled at $\omega_n=(2n+1)\pi/\beta$, one is given $\{i\omega_n,G(i\omega_n)\}$ and aims to recover

$$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$

The same inverse problem appears in imaginary-time form through kernels such as $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$, and in matrix-valued settings through

$$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$

with Hermiticity $A(\omega)=A(\omega)^\dagger$ and, in multi-orbital causal formulations, positive semidefiniteness $A(\omega)\succeq 0$ (Huang et al., 2024, Motoyama et al., 2024).

The difficulty is intrinsic. Analytic continuation is an inverse Laplace transform problem, so small perturbations in imaginary-axis data can produce large changes in the recovered spectrum. Quantum Monte Carlo input typically contains stochastic noise and may also contain autocorrelation-related gaps, making continuation numerically delicate (Huang et al., 2024). In a complementary formulation based on analytic function theory, retarded correlators are analytic in the upper half-plane $\mathbb{C}^+$; fermionic correlators can satisfy a Herglotz or Nevanlinna property, while bosonic correlators require different image-domain constraints and branch-cut handling (Bergamaschi et al., 2023). This analytic structure is the basis for essentially all hybrid methods.

2. Hybridization as a design principle

A hybrid method combines components that address different failure modes of continuation. One component may supply stability against noise, another may preserve sharp low-energy structure, another may impose causality or matrix positivity, and another may provide rigorous uncertainty control. This suggests that “hybrid” refers less to a single algorithm than to a recurrent architectural pattern in which complementary numerical and physical priors are coupled within one continuation workflow.

Family	Hybrid components	Typical target
Barycentric rational continuation	AAA + barycentric rational interpolation + optional Prony denoising + pole refinement	Matsubara Green’s functions
Sparse constrained continuation	IR-basis sparse modeling + ADMM + PSD projection	Multi-orbital spectral matrices
Sparse–Padé continuation	SpM regularization + variance-weighted Padé guidance	Noisy QMC spectra
Contour-deformation GW continuation	Exact contour deformation + continuation of $W$	GW self-energies and lifetimes
Conformal interpolation	Cayley or square-root maps + Schur-class interpolation + Wertevorrat bounds	Retarded correlators from Euclidean data
Learning-driven continuation	Neural predictors or evolutionary search + physical constraints or MEM-style refinement	Synthetic and QMC spectral reconstruction

Several papers make the hybrid character explicit. The barycentric method couples adaptive support-point selection with rational interpolation and allows optional Prony denoising and a switch to pole representation for discrete spectra (Huang et al., 2024). The multi-orbital sparse method combines IR-basis $\omega_n=(2n+1)\pi/\beta$0 regularization with semi-positive definiteness constraints enforced either one-shot or self-consistently in ADMM (Motoyama et al., 2024). SpM–Padé augments sparse modeling with a frequency-dependent Padé consistency term weighted by Padé variance (Motoyama et al., 2021). In GW, the hybridization is between contour deformation for the self-energy integral and analytic continuation applied only to $\omega_n=(2n+1)\pi/\beta$1, not to the more structured self-energy $\omega_n=(2n+1)\pi/\beta$2 (Duchemin et al., 2019).

3. Rational, Padé, and pole-based hybrids

A major strand of hybrid continuation uses rational approximants. In the barycentric AAA method, the interpolant is

$\omega_n=(2n+1)\pi/\beta$3

with $\omega_n=(2n+1)\pi/\beta$4 and $\omega_n=(2n+1)\pi/\beta$5. The adaptive Antoulas–Anderson algorithm selects support points greedily by maximizing the current residual, computes weights from a linear least-squares problem via the Loewner matrix and an SVD, and stops once the residual over unused samples falls below a tolerance; a default tolerance of $\omega_n=(2n+1)\pi/\beta$6 relative to $\omega_n=(2n+1)\pi/\beta$7 is suggested (Huang et al., 2024). The continued retarded function is then evaluated as $\omega_n=(2n+1)\pi/\beta$8, with $\omega_n=(2n+1)\pi/\beta$9 used in discrete-peak demonstrations. The method can also exploit known symmetries, subtract and restore the Hartree term for self-energies, and switch to a pole representation with BFGS residue optimization when discrete weights are hard to recover. Benchmarks report accurate reconstruction of continuous and discrete spectra, including non-positive-definite spectra, comparable tolerance to intermediate noise relative to MaxEnt, and typical speedups of at least $\{i\omega_n,G(i\omega_n)\}$0 over MaxEnt (Huang et al., 2024).

A second rational line revisits Padé approximation through ensemble averaging. Analytic continuation by averaging Padé approximants varies the number of fitted input points and Padé coefficients independently, accepts only physical continuations, and averages over the accepted set. A similarity-based weighting scheme further suppresses spurious structures by favoring mutually similar spectra among physical least-squares configurations (Schött et al., 2015). In this formulation, the rational approximant is built in least-squares form rather than as a square solve, and the ensemble suppresses random zero–pole artifacts that arise from limited precision and noisy input.

A third hybrid rationale appears earlier in the pipeline: high-precision Matsubara evaluation by Padé decomposition followed by Padé continuation. Replacing slowly convergent Matsubara-frequency summations by Padé-frequency summations substantially improves the precision of the input data used for continuation. In the reported benchmark, $\{i\omega_n,G(i\omega_n)\}$1 Padé frequencies already produced roughly $\{i\omega_n,G(i\omega_n)\}$2 digits for a single Matsubara evaluation, while demanding continuation with $\{i\omega_n,G(i\omega_n)\}$3 required about $\{i\omega_n,G(i\omega_n)\}$4 digits to reproduce exact spectral features; $\{i\omega_n,G(i\omega_n)\}$5 produced perfect agreement in the benchmark spectral reconstruction (Han et al., 2017). Here the hybrid aspect lies in coupling an accurate imaginary-axis computation with a rational real-axis continuation.

A related sparse-pole strategy appears in the continuation of limited noisy Matsubara data. There, the algorithm first interpolates $\{i\omega_n,G(i\omega_n)\}$6 or $\{i\omega_n,G(i\omega_n)\}$7, maps an imaginary-axis interval to the unit circle by a conformal transform, computes Fourier coefficients by FFT, and then applies Prony’s method to a Hankel matrix to recover a small number of poles. The final amplitudes are determined by constrained least squares with non-negativity or imaginary-part constraints (Ying, 2022). This formulation is explicitly targeted at molecule cases with discrete spectra and condensed-matter cases with a quasi-particle prior.

4. Sparse modeling, entropy methods, and convex constrained hybrids

Sparse modeling introduces regularization through the intermediate representation (IR) basis. For multi-orbital continuation, the kernel is factorized by SVD,

$\{i\omega_n,G(i\omega_n)\}$8

and the continuation is posed as an $\{i\omega_n,G(i\omega_n)\}$9-regularized inverse problem in the IR coefficients. The self-consistent variant augments the sparse objective by a semi-positive definiteness penalty,

$$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$0

and solves it by ADMM with Hermiticity enforcement and PSD projection during the iteration (Motoyama et al., 2024). The paper distinguishes a one-shot repair method and a self-consistent method. In the reported two-orbital test, one ADMM sweep costs about $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$1 ms for one-shot PSD repair and about $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$2 ms for self-consistent PSD enforcement; enforcing PSD only every tenth frequency reduces the self-consistent cost to about $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$3 ms per sweep, approximately a $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$4 speedup, with negligible RMSE increase, whereas $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$5 introduces oscillations in the smaller eigenvalue (Motoyama et al., 2024). This is a specifically matrix-valued hybridization: data-driven sparsity is coupled to physics-mandated causality constraints.

SpM–Padé combines sparse modeling and Padé by adding a frequency-dependent quadratic penalty to the sparse objective,

$$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$6

Padé contributes a low-energy guide where its variance is small, while sparse modeling maintains robustness to noise and enforces positivity and sum rules (Motoyama et al., 2021). The reported outcome is low-variance and low-bias continuation at almost the same computational cost as sparse modeling alone.

Maximum-entropy methods have also been reformulated in hybrid terms. The dual Newton formulation of MEM recasts the entropy-regularized inverse problem into a smooth strongly convex dual optimization in $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$7 variables, solves the finite-temperature normalization $$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$8 explicitly, and retains all singular vectors rather than truncating to a Bryan subspace (Chuna et al., 3 Jan 2025). The paper argues that this preserves the theoretical benefits of Bryan’s MEM while avoiding theoretical issues associated with hard singular-vector truncation, and reports better estimates and error bars under noise on test problems from lattice QCD and plasma physics (Chuna et al., 3 Jan 2025).

5. Conformal, interpolation-theoretic, and domain-specific hybrids

Another line of hybridization proceeds through complex analysis rather than direct inverse optimization. In the conformal-map approach, the upper half-plane is mapped to the unit disk by a Cayley transform,

$$G^R(\omega)=\lim_{\eta\to 0^+}G(\omega+i\eta),\qquad A(\omega)=-\frac{1}{\pi}\operatorname{Im}G^R(\omega).$$9

while bosonic correlators can additionally require a square-root map

$$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$0

After mapping Euclidean data to Schur-class interpolation data in the disk, the full family of admissible interpolants is characterized by constrained Nevanlinna–Pick theory, and the value set at an interior point is a disk $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$1 with explicitly computable center and radius (Bergamaschi et al., 2023). Mapping these disks back to $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$2 gives rigorous bounds for the smeared spectral function

$$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$3

The unsmeared limit $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$4 remains ill-posed: on the boundary of the disk, the Wertevorrat fills $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$5, and the bounds become infinite (Bergamaschi et al., 2023). This hybrid structure couples conformal geometry, constrained interpolation, and rigorous uncertainty quantification.

In GW theory, the hybrid object is the self-energy evaluation itself. The contour-deformation formula splits $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$6 into an integral along the imaginary axis and explicit residue terms from poles of $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$7. Instead of continuing the full self-energy $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$8, the method continues matrix elements of the screened Coulomb interaction $$K(\tau,\omega)=e^{-\tau\omega}/(1\pm e^{-\beta\omega})$$9, which are much smoother than $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$0, and inserts those values into the exact contour-deformation decomposition (Duchemin et al., 2019). For frontier quasiparticles, about $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$1 imaginary-axis points are sufficient for meV-level accuracy in difficult cases, while for deeper valence and core states the method augments the imaginary-axis references by a coarse grid of points parallel to the real axis with $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$2 eV and height $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$3 (Duchemin et al., 2019). The paper reports that continuing $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$4 is far more robust than continuing $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$5, especially because $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$6 inherits a dense pole structure from both $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$7 and $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$8.

6. Learning-based and evolutionary hybrids

Machine-learning approaches bring a different form of hybridization: a learned prior is combined with physical constraints and, in some cases, with subsequent MEM-style refinement. A supervised ANN continuation framework uses imaginary-time inputs projected to low-dimensional representations—three principal components for an oscillator problem and sixty-four Legendre coefficients for a fermionic Green’s-function problem—and reconstructs spectra on a $$G_{ab}(i\omega_n)=\int_{-\infty}^{\infty}\frac{A_{ab}(\omega)}{i\omega_n-\omega}\,d\omega,$$9-point frequency grid through an MLP with batch normalization, ReLU activations, and a softmax output layer (Fournier et al., 2018). The softmax enforces positivity and unit normalization. In the reported comparison, the ANN reached the same level of accuracy as MaxEnt for low-noise input, performed significantly better at higher noise, and processed $A(\omega)=A(\omega)^\dagger$0 pairs in about $A(\omega)=A(\omega)^\dagger$1 s including library load, versus about $A(\omega)=A(\omega)^\dagger$2 min for MaxEnt, which is described as an almost three orders of magnitude speedup (Fournier et al., 2018). The hybrid pathways described for this framework use the ANN prediction as a MaxEnt default model or as a denoising and projection step before conventional continuation (Fournier et al., 2018).

The Feature Learning Network develops this further by separating spectral-feature extraction from Matsubara-to-feature inference. FL-net uses two encoders and one decoder: one encoder maps $A(\omega)=A(\omega)^\dagger$3 to latent features $A(\omega)=A(\omega)^\dagger$4, another maps Matsubara data $A(\omega)=A(\omega)^\dagger$5 to the same latent space, and the decoder reconstructs $A(\omega)=A(\omega)^\dagger$6 from $A(\omega)=A(\omega)^\dagger$7 with a softmax output (Zhao et al., 2024). In a single-Gaussian dataset, the learned latent coordinates were shown to be locally equivalent to the physical parameters $A(\omega)=A(\omega)^\dagger$8 through a full-rank Jacobian, and the model achieved an improvement of at least $A(\omega)=A(\omega)^\dagger$9 over MEM and previous neural-network approaches on the tested synthetic datasets (Zhao et al., 2024). The paper also derives a robustness measure from the Jacobian $A(\omega)\succeq 0$0, with singular-value decomposition $A(\omega)\succeq 0$1 and sensitivity

$A(\omega)\succeq 0$2

and shows that increasing hidden dimensionality lowers the loss while decreasing robustness (Zhao et al., 2024).

Differential Evolution for Analytic Continuation introduces a parameter-free evolutionary search in which the spectral weights, crossover probability $A(\omega)\succeq 0$3, and differential weight $A(\omega)\succeq 0$4 are all embedded in the genome and updated self-adaptively (Nichols et al., 2022). Mutation and crossover operate frequency-wise, the fitness is the usual $A(\omega)\succeq 0$5 misfit, and the stopping criterion is $A(\omega)\succeq 0$6, with $A(\omega)\succeq 0$7 for simulated data and $A(\omega)\succeq 0$8 for superfluid helium (Nichols et al., 2022). In the reported CPU-time comparison at the large-noise level $A(\omega)\succeq 0$9, DEAC reduced CPU hours by as much as $\mathbb{C}^+$0 relative to MEM and $\mathbb{C}^+$1 relative to FESOM, while reproducing the phonon–roton spectrum in bulk $\mathbb{C}^+$2He, including a maxon near $\mathbb{C}^+$3 at $\mathbb{C}^+$4 meV and a roton near $\mathbb{C}^+$5 at $\mathbb{C}^+$6 meV (Nichols et al., 2022). The accompanying hybrid guidance describes DEAC combined with MaxEnt, SAC, and physics-informed parameterizations (Nichols et al., 2022).

7. Applications, performance envelopes, and persistent limitations

Hybrid continuation methods are now applied across scalar and matrix Green’s functions, optical conductivity, anomalous Nambu components, self-energies, GW quasiparticle energies and lifetimes, lattice-QCD correlators, hadronic vacuum polarization, plasma spectra, and dynamic structure factors (Huang et al., 2024, Duchemin et al., 2019, Bergamaschi et al., 2023, Chuna et al., 3 Jan 2025, Nichols et al., 2022). Performance claims are method-dependent rather than universal. The barycentric AAA method is reported to resolve discrete peak positions and weights far from the Fermi level and to handle non-positive-definite spectra without special constraints (Huang et al., 2024). By contrast, the multi-orbital sparse method explicitly requires $\mathbb{C}^+$7 to preserve causality in matrix spectra (Motoyama et al., 2024). This corrects a common misconception: positivity or semidefiniteness is indispensable in some continuation problems and irrelevant or even inappropriate in others.

The limitations are equally method-specific. In the barycentric method, too few Matsubara points can produce oscillatory continuations; the reported practical guidance states that accuracy improves with the number of Matsubara points up to about $\mathbb{C}^+$8, gains saturate beyond that, and with fewer than $\mathbb{C}^+$9 points AAA may oscillate (Huang et al., 2024). Sharp band edges remain difficult for both BarRat and MaxEnt, and both can show oscillations and edge or bandwidth bias in gapped systems (Huang et al., 2024). In the multi-orbital sparse method, aggressive skipping of PSD enforcement eventually degrades the smaller eigenvalue and introduces oscillations (Motoyama et al., 2024). In GW, imaginary-axis-only continuation of $W$0 can still yield errors approaching $W$1 eV for states far from the gap, while adding about $W$2 coarse first-quadrant reference points reduces errors below $W$3 eV within $W$4 gaps (Duchemin et al., 2019). In FL-net, increasing hidden dimension past the practical optimum can reduce robustness even as training loss decreases (Zhao et al., 2024).

No single hybrid construction removes the fundamental ill-posedness of unsmeared continuation. The conformal-interpolation program makes this explicit by showing that rigorous uncertainty disks remain finite only for interior points $W$5 with $W$6, whereas the boundary limit $W$7 restores the ill-posed problem (Bergamaschi et al., 2023). This suggests that the enduring role of hybridization is not to circumvent ill-posedness, but to redistribute it: away from the most unstable representation and into a combination of rational approximation, variational regularization, analytic structure, learned priors, or uncertainty geometry that is better matched to the data and to the physics.