Problem in TCA with AAA compression

[jaemolihm]

During a DFT+DMFT study using XCA solver, I found a case where the combination of (1) TCA (order=3), and (2) compress_hybridization=True breaks XCA. The plot shows the pseudoparticle self-energy at the first evaluation (before entering PPSC loop), and it is very wrong. The problem disappears if I use order = 2 or compress_hybridization=False.

I attached the input data and script to reproduce.

xca_test.zip

[jaemolihm]

It turns out that the problem occurs whenever there are nearly-degenerate poles in the AAA fitting (called Froissart doublets, according to Jason). They are very close to each other with very large residues. For example, see poles 6 and 8:

Per-pole weight norms:
  pole  0: w= +155.851851  ||weight||=6.55e-02
  pole  1: w=  +80.346280  ||weight||=1.17e-01
  pole  2: w=  -58.849499  ||weight||=1.39e-02
  pole  3: w=  +44.370933  ||weight||=2.27e-01
  pole  4: w=  +30.980637  ||weight||=2.32e-01
  pole  5: w=  +14.285597  ||weight||=8.39e-02
  pole  6: w=  -10.490616  ||weight||=2.91e+08
  pole  7: w=  -13.833345  ||weight||=1.70e-02
  pole  8: w=  -10.490616  ||weight||=2.91e+08
  pole  9: w=   +7.103799  ||weight||=2.51e-02
  pole 10: w=   -2.010472  ||weight||=6.04e-03

This is a Claude-suggested change to adapol that merges such degenerate poles. With this change, the TCA+AAA problem is gone.

diff --git a/src/adapol/fit_utils_xca.py b/src/adapol/fit_utils_xca.py
index 39d45a5..348f10c 100644
--- a/src/adapol/fit_utils_xca.py
+++ b/src/adapol/fit_utils_xca.py
@@ -77,6 +77,8 @@ def polefitting(Deltaiw, Z, Deltat,tgrid, Deltat_dense, tgrid_dense, beta,
         pol = np.real(pol)
         weight, _, residue = get_weight_t(pol, tgrid, Deltat,beta)
         pol, weight = aaa_reduce(pol, weight,eps)
+        pol = merge_degenerate_poles(pol)
+        weight, _, residue = get_weight_t(pol, tgrid, Deltat,beta)
 
         def fhere(pole):
             return erroreval_t(pole, tgrid, Deltat,beta, statistics=statistics) 
@@ -110,3 +112,47 @@ def aaa_reduce(pol, R, eps=1e-6):
         Rnorm[i] = np.linalg.norm(R[i])
     nonz_index = Rnorm > eps
     return pol[nonz_index], R[nonz_index]
+
+
+def merge_degenerate_poles(pol, rtol=1e-6):
+    """Merge near-degenerate poles from AAA into single poles.
+
+    The AAA algorithm (via find_pol) can produce exactly degenerate poles
+    from its generalized eigenvalue problem. When two poles coincide, the
+    downstream least-squares weight fitting becomes catastrophically
+    ill-conditioned (two identical columns in the kernel matrix), producing
+    weight matrices with ~1e8 norm that cause blow-up in diagram evaluation.
+
+    Parameters
+    ----------
+    pol : ndarray
+        Pole positions (energy units).
+    rtol : float
+        Relative tolerance for merging. Poles are merged when
+        |pol_i - pol_j| < rtol * (max(|pol|) - min(|pol|) + 1).
+
+    Returns
+    -------
+    pol_merged : ndarray
+        Poles with degenerate groups replaced by their mean.
+    """
+    if len(pol) <= 1:
+        return pol
+
+    scale = np.ptp(np.abs(pol))  # range of |pol|
+    atol = rtol * (scale + 1.0)  # +1 avoids zero scale
+
+    order = np.argsort(pol)
+    pol_sorted = pol[order]
+
+    merged = []
+    i = 0
+    while i < len(pol_sorted):
+        group = [pol_sorted[i]]
+        while i + 1 < len(pol_sorted) and abs(pol_sorted[i + 1] - group[0]) < atol:
+            i += 1
+            group.append(pol_sorted[i])
+        merged.append(np.mean(group))
+        i += 1
+
+    return np.array(merged)

[jasonkaye]

Thanks for point this out @jaemolihm . @Hertz4 Is anything being done in adapol to handle/remove Froissart doublets, as described in the original AAA paper?

[Hertz4]

Not in this version of the code.
I would suggest that for now we could add @jaemolihm 's fix to the current XCA code. For future release of XCA, the pole fitting dependence would be through the standalone adapol package, and we could add the part that removes Froissart doublets properly (as described in the original AAA paper).
@jaemolihm Thanks for catching this!