Include self-symmetric fixed points in symmetrized DLR grids - #20
Include self-symmetric fixed points in symmetrized DLR grids#20Wentzell wants to merge 9 commits into
Conversation
727f3fb to
6f3f2d8
Compare
|
Update 12th Aug: I tested this new PR with CoQui, and this is working for GW in CoQui at least. Although further testing is still required (and ongoing AbInitioQHub/coqui#51 ), i think this PR might fixed the CoQui and GW/GF2 issue already. Weiyi Hi, Nils Thanks for this quick fix for symmetrized grid. I’m going to have a test for CoQui case. Best regards, |
a is already zero from its declaration two lines above. Assisted-by: Claude <noreply@anthropic.com>
build_rf_fine, build_it_fine and geterr_k_it only read their fineparams argument. Source-compatible for callers, but it does change the mangled names, so it is an ABI break for anything linking a prebuilt cppdlr. Assisted-by: Claude <noreply@anthropic.com>
Same quantity, same name as in the eps-cutoff overload: the row where the first mirror pair begins. Assisted-by: Claude <noreply@anthropic.com>
Two commented-out lines and two comments describe an abandoned layout in which the self-paired middle row was placed last; the code places it first. Assisted-by: Claude <noreply@anthropic.com>
Both symcompare examples printed the unsymmetrized rank r on the "Symmetrized DLR" line, understating or overstating the symmetrized rank whenever the two differ. symcompare_it already computed rsym and used it to size its arrays, so only the format argument was wrong; symcompare_if did not compute it at all. Assisted-by: Claude <noreply@anthropic.com>
The symmetrized DLR grid construction needs to select the self-symmetric fixed point of a mirror-symmetric fine grid (omega=0, tau=beta/2) as a single self-paired pivot. Extend pivrgs_sym accordingly: * The eps-cutoff overload accepted only an even number of rows. Allow an odd row count, selecting the middle row (m-1)/2 first as a self-paired pivot and then proceeding in mirror pairs, for an odd epsilon-rank. This mirrors what the fixed-rank overload already does for odd m. The pivot order is built parity-agnostically and the input rows are copied in that order, so the two cannot disagree; results are unchanged for even m. * The fixed-rank overload rejected r = n+1 with n odd. The symmetrized fermionic Matsubara grid needs exactly this: it has no self-symmetric frequency, so it reaches an odd rank n only via n+1 mirror-paired nodes. The remaining checks (r <= min(m, n+1), matching parity of r and m) already cover the meaningful constraints. Assisted-by: Claude <noreply@anthropic.com>
The row ordering was spelled out twice, once for the row copy and once for the pivot vector, each split into an even-m and an odd-m branch. Build the pivot order once and copy the rows through it, as the eps-cutoff overload already does, so the two cannot disagree. Bit-for-bit identical output. Assisted-by: Claude <noreply@anthropic.com>
Each overload spelled out "normalize row j, then orthogonalize every later row against it and update its norm" three times: once for the self-paired middle row and twice per selected mirror pair. Extract it as a lambda. Bit-for-bit identical output. Assisted-by: Claude <noreply@anthropic.com>
The symmetrized construction selected nodes only in exact mirror pairs, so the real-frequency grid had to take a near-degenerate pair +/-delta straddling omega=0, whose basis columns 1/(i*nu_n -/+ delta) are nearly parallel. This ill-conditioned the DLR fit matrices: at eps=1e-10, cond(cf2it) was 2.5e13 and 1.1e14 for lambda=10 and 1000, against ~1e10 unsymmetrized. Harmless on DLR-representable inputs, but workflows applying a DLR-built operator to non-spectral vectors (e.g. the Krylov iterates of an eigensolver) suffer O(1) floating-point error. Instead always include the self-symmetric fixed point as a single self-paired node: omega=0 in real frequency, tau=beta/2 in imaginary time, n=0 in the bosonic Matsubara grid. Conditioning becomes 4.4e10 and 6.8e10 for the same two cases, ranks stay within one of the unsymmetrized rank, and accuracy on spectral inputs is unchanged. Non-symmetrized grids are unchanged bit for bit. Consequences: - The symmetrized DLR rank r is now odd. - The symmetrized fermionic Matsubara grid has no self-symmetric frequency, so it decouples to an even niom = r + 1 of mirror pairs, inverted by least squares (as the symmetrized bosonic case already was). - imtime_ops now rejects an even-rank dlr_rf under SYM, which every symmetrized grid from cppdlr <= 1.3.0 has, so a persisted grid fed back with SYM hits it. Updates the symmetrized DLR grids section of doc/background.rst, which documented the previous pair-only behavior. Still pending: triqs_tprf's DLR Eliashberg tests remain pinned to non-symmetrized meshes (symmetrize=False in test/python/eliashberg/dlr_eliashberg_solver.py), since they were failing with the symmetrized ones. They should be revisited against this change. Assisted-by: Claude <noreply@anthropic.com>
6f3f2d8 to
0c4dc71
Compare
|
I can confirm that with this change, the tprf Py_eliashberg/dlr_eliashberg_solver python test passes also with a symmetric DLR grid, while before it did not. |
|
In addition, I tested PR #20 against triqs/xca higher-order diagrams using native C++ binaries, overriding the hard-coded DLR defaults so that the same xca test suite was compiled once with a non-symmetric grid and once with a symmetric grid. Both configurations passed all 20 CTest executables (82 gtest cases, including 74 assertion-bearing cases) with zero failures. The override was confirmed to be active: the same two-argument DLR construction produced rank 16 with P=48 for the non-symmetric grid and rank 17 with P=51 for the symmetric grid. The OCA and third-order diagram errors were comparable (1.17e-10 versus 1.44e-10), and an additional C++ probe covering XCA diagram orders 2 through 8 across three parameter regimes agreed with the analytic integrals, with a maximum relative error of 3.33e-8. I also compared the current dense and block-sparse Sigma and Green-function evaluators on a common imaginary-time grid through third order; the maximum relative difference between symmetric and non-symmetric results was 6.01e-9. So, as far as these tests show, XCA higher-order integrals work correctly with both symmetric and non-symmetric DLR grids after this fix. (previously, the higher order integration beyond OCA is also failing against symmetric DLR grid in triqs/xca) One thing now I found is xca now directly find the cppdlr dependence through TRIQS core linking version, so unless all further TRIQS build are based on this PR #20, the symmetric grid for higher order integral would not automatically fix. Weiyi |
2 participants
Summary
The symmetrized DLR construction selected nodes only in exact mirror pairs, so the
real-frequency grid was forced onto a near-degenerate pair
±δstraddlingω = 0,whose basis columns
1/(iν_n ∓ δ)are nearly parallel. That ill-conditioned the DLRfit matrices: at
eps=1e-10,cond(cf2it)was2.5e13(λ=10) and1.1e14(λ=1000),against
~1e10unsymmetrized.Instead always include the self-symmetric fixed point as a single self-paired node:
ω = 0in real frequency,τ = β/2in imaginary time (weightless, sincek_it(0.5, ω) = k_it(-0.5, ω)),n = 0in the bosonic Matsubara grid. Conditioningbecomes
4.4e10and6.8e10, ranks stay within one of the unsymmetrized rank,accuracy on spectral inputs is unchanged, and non-symmetrized grids are unchanged bit
for bit.
Addresses #19, which identifies this as the cppdlr-side root cause behind
TRIQS/tprf#50, a diverging GF2 iteration,
and CoQui's self-consistent GW on metals — and asks whether including
ω = 0fixesit. The reproducer from #19 (
tprf_issue, 53b6dd9) is not folded in here and shouldbe re-run before that issue is closed.
pivrgs_symgains the corresponding support: the eps-cutoff overload accepts an oddrow count, selecting the middle row first as a self-paired pivot (and throwing if that
forced pivot is negligible, rather than reporting a misleading rank); the fixed-rank
overload no longer rejects
r = n+1with oddn, which the symmetrized fermionicMatsubara grid needs. Even-
mresults are unchanged. The preparatory commits arebit-for-bit identical refactors, plus a
symcomparefix (it printed the unsymmetrizedrank on the symmetrized line) and the
doc/background.rstupdate.Breaking changes
ris now odd, so code assuming an even rank or purely pairednodes needs updating.
imtime_ops(..., SYM)throws on even-rankdlr_rf— which every symmetrizedgrid from cppdlr ≤ 1.3.0 is, so a persisted grid replayed with
SYMhits it.SYMnowcontains
n = 0, giving a squareniom = r; fermionicSYMhas no self-symmetricMatsubara frequency, so it uses
niom = r + 1mirror pairs inverted by leastsquares, and a fermionic
vals2coefs/coefs2valsround trip recovers samples toO(eps)rather than machine precision.build_rf_fine,build_it_fine,geterr_k_ittakefineparams const &,and the builders gained a defaulted
symmetrizeparameter. Source-compatible, butmangled names change.
Testing
ctest -j 16passes 12/12, with new tests for the fixed points in both the fine andDLR grids, odd-
mand negligible-middle-rowpivrgs_sym, the even-rankSYMthrow,and an HDF5 round trip under
SYMfor both statistics (the fermionic case stores noif2cfLU factors).Follow-up
triqs_tprf's DLR Eliashberg tests are still pinned to
symmetrize=False(
test/python/eliashberg/dlr_eliashberg_solver.py) because they failed withsymmetrized meshes; they should be revisited against this change.