[Bug report] Symmetrized DLR grids develop a near-null ±ω mode and diverge under repeated nonlinear transforms

[weiyiguo9]

Description

With symmetrize=true, the DLR collocation matrices develop a nearly singular mirror-odd mode dominated by the central real-frequency pair ±δ. In a simple self-consistent scalar GF2 iteration, roundoff injected into this mode first appears as rapidly growing, cancelling DLR coefficients and then causes the reconstructed Green's function/self-energy to diverge. The same calculation with symmetrize=false converges to machine precision.

This looks like the cppdlr-side root cause of TRIQS/tprf#50: Eliashberg solver gives incorrect results for symmetrized DLR meshes. It also agrees with the failing reproducer and analysis already present in commit 53b6dd9 on the tprf_issue branch. Moreover, the same issue affects the CoQui self-consistent GW calculation with metallic compounds, which is not affected if using IR basis.

Why GF2, GW, and Eliashberg expose the same cppdlr failure mode

These should not be viewed as three unrelated application-level failures. GF2, self-consistent GW, and the linearized Eliashberg solver all contain the same numerically critical sequence:

1. values on a symmetrized DLR grid are converted to coefficients through the ill-conditioned solve c = K_sym^{-1} v;
2. the coefficients are evaluated in another representation or on another grid;
3. pointwise products are formed;
4. the result is projected back to DLR and the operation is repeated.

The near-null mirror-odd mode associated with the central ±δ pole pair is therefore common to all three calculations. Tiny floating-point components along this mode become large opposite DLR coefficients. They cancel at the original interpolation nodes, but the cancellation is lost after evaluation on another grid and multiplication.

The application-specific products are different,

• GF2: Sigma(tau) ~ G(tau)^2 G(beta-tau);
• GW: P(tau) ~ G(tau)G(-tau) and Sigma(tau) ~ G(tau)W(tau);
• Eliashberg: the kernel applies a G Delta G product,

but they amplify the same cppdlr projection error. In GF2 and self-consistent GW this appears as growth of DLR coefficients followed by divergence of the self-consistency loop. In the Eliashberg solver it appears as loss of floating-point linearity on non-spectral Krylov vectors, producing spurious and non-reproducible eigenvalues.

This interpretation is consistent with TRIQS/tprf#50. It is also supported by two independent controls: the metallic CoQui GW calculation remains stable with the IR basis, and the same GF2 equations remain stable with Lehmann.jl's numerically symmetrized :sym basis. Thus the common discriminator is the current cppdlr symmetrized-grid construction, not the GF2, GW, or Eliashberg physics.

Steps to reproduce

Run the following with TRIQS built against cppdlr 1.3.0 (Julia comparison in attachment):

import numpy as np
from triqs.gfs import (
    Gf, MeshDLRImFreq, MeshDLRImTime,
    make_gf_dlr, make_gf_dlr_imfreq, make_gf_imfreq,
)

BETA, LAMBDA, EPS = 100.0, 1.0e4, 1.0e-13
U, MU, MIXING, NITER = 1.0, 0.15, 0.70, 100


def bethe(z):
    root = np.sqrt(z * z - 1.0)
    root[np.imag(z) * np.imag(root) < 0.0] *= -1.0
    return 2.0 / (z + root)


def run(symmetrize):
    iw_mesh = MeshDLRImFreq(
        BETA, "Fermion", LAMBDA / BETA, EPS, symmetrize
    )
    tau_mesh = MeshDLRImTime(
        BETA, "Fermion", LAMBDA / BETA, EPS, symmetrize
    )
    iw = np.array([complex(point) for point in iw_mesh])
    tau = np.array([float(point) for point in tau_mesh])
    sigma = Gf(mesh=iw_mesh, target_shape=[1, 1])
    sigma.data[:] = 0.0

    for iteration in range(1, NITER + 1):
        g_iw = Gf(mesh=iw_mesh, target_shape=[1, 1])
        g_iw.data[:, 0, 0] = bethe(
            iw + MU - sigma.data[:, 0, 0]
        )
        g = make_gf_dlr(g_iw)
        g_tau = np.array([g(t)[0, 0] for t in tau])
        g_reflected = np.array([g(BETA - t)[0, 0] for t in tau])

        sigma_tau = Gf(mesh=tau_mesh, target_shape=[1, 1])
        sigma_tau.data[:, 0, 0] = U**2 * g_tau**2 * g_reflected
        sigma_new = make_gf_dlr_imfreq(make_gf_dlr(sigma_tau))

        old = make_gf_imfreq(make_gf_dlr(sigma), n_iw=256).data[:, 0, 0]
        new = make_gf_imfreq(make_gf_dlr(sigma_new), n_iw=256).data[:, 0, 0]
        residual = np.max(np.abs(new - old))
        sigma.data[:] = (
            (1.0 - MIXING) * sigma.data + MIXING * sigma_new.data
        )

        dense = make_gf_imfreq(make_gf_dlr(sigma), n_iw=256).data[:, 0, 0]
        coefficient = np.max(np.abs(make_gf_dlr(sigma).data))
        hermiticity = np.max(
            np.abs(dense - np.conjugate(dense[::-1]))
        )
        if iteration in (20, 40, 60, 100):
            print(symmetrize, iteration, residual, coefficient, hermiticity)


for symmetric_dlr in (False, True):
    run(symmetric_dlr)

Observed at iteration 100:

grid	rank	residual	max DLR coefficient	Hermiticity error
non-symmetrized	88	5.66e-16	2.38e-2	9.95e-15
symmetrized	88	4.38e3	2.34e14	2.53e3

For the symmetrized run, the coefficient already reaches 9.13e2 at iteration 20 and 3.02e8 at iteration 40, while the non-symmetrized run remains stable.

The failure is not caused by particle-hole asymmetry: repeating the test with MU=0 also diverges.

Cross-implementation check: Lehmann.jl

I ported the same GF2 equations and parameters to Julia using Lehmann.jl 0.2.7. The like-for-like comparison is Lehmann's symmetry=:sym, which uses a numerically symmetrized grid/basis without imposing physical particle-hole symmetry. Lehmann's :ph mode is a different, physically constrained kernel and was not used here (the test has MU=0.15).

After 100 iterations:

Lehmann.jl grid	rank	residual	max DLR coefficient	Hermiticity error
:none	89	1.10e-15	1.32e-1	1.30e-13
:sym	94	2.41e-15	1.30e-1	8.89e-16

Neither Lehmann.jl run shows coefficient growth or divergence. The early fixed-point residuals agree with the TRIQS/cppdlr calculation through iteration 15 to the displayed precision, which checks that the GF2 equation, the G(beta-tau) sign, and the Bethe square-root branch were ported consistently. This cross-check indicates that the instability is specific to cppdlr's current symmetric-grid construction/conditioning, rather than an unavoidable property of a symmetrized DLR representation.
see attached file for Julia reproduce

Conditioning evidence

For the parameters above, the coefficient-to-imaginary-time matrix has:

grid	cond(K_tau)	smallest singular value
non-symmetrized	5.02e13	4.04e-13
symmetrized	1.25e16	1.62e-15

The smallest symmetric singular vector is mirror-odd and dominated by the near-zero pole pair ±0.00130680955. Large opposite coefficients cancel at the DLR nodes but lose cancellation after evaluation in another representation and nonlinear multiplication.

The effect is not restricted to the tight GF2 parameters. At the smaller parameters used by the Eliashberg reproducer (beta=2, lambda=10, eps=1e-8), I obtain cond(K_tau)=8.75e9 for the symmetric grid versus 2.16e8 for the non-symmetric grid.

A possibly contributing source-level issue is that the tolerance-based pivrgs_sym checks the residual norm of the first member of a mirror pair, but normalizes the second member without an equivalent small-norm check. This can retain a nearly dependent second direction. The resulting matrices are then inverted with an unregularized LU solve.

Expected behavior

Symmetrizing the numerical node set should not introduce a nearly singular mode that destroys repeated DLR transforms or changes a converged physical result.

Actual behavior

The symmetric basis can hide enormous mirror-odd coefficients behind cancellation at its nodes. Repeated transforms and nonlinear operations amplify this direction until residuals, Hermiticity, and the physical solution diverge.

Proposed direction

The fixed-point construction described in the tprf_issue branch looks appropriate: include each symmetry fixed point as a single node (omega=0, tau=beta/2, and bosonic n=0). This makes the symmetric real-frequency rank odd; fermionic Matsubara nodes can then use r+1 mirror-paired samples with a least-squares solve.

Versions

• cppdlr 1.3.0, git hash 5f9adf91d7620ce145ead8e4e906151493038af3
• TRIQS 4.0.0, git hash 86c3b4d6842a706b040dc5a364671097d8ce2473
• Python 3.14.6
• macOS 15.7.3, arm64

gf2_sym_vs_nonsym_lehmann.zip

Additional Information

This bug detection and comparison is also suggested by Xiansheng Cai

[jasonkaye]

Thanks @weiyiguo9, this is is very interesting. Do we have concrete evidence that explicitly including omega = 0 in the DLR frequency selection fixes this issue? Has someone prototyped this fix?