Problems with CG_METHOD=0 [encountered in many different use cases]

Question

Problems with CG_METHOD=0 [encountered in many different use cases]

beddalumia opened this issue 2 years ago · 2 comments

What's up

As I was performing some initial tests on the incoming changes on fit_overhaul branch, I noticed something is wrong with the Weiss field analytical gradient, as it is currently on master. By wrong I mean that the fitted field is qualitatively different from the original one, so seriously wrong.

Note that this happens only if I request the analytic gradient: switching to CG_METHOD=1 (minimize routine) solves everything, switching only CG_GRAD=1 but remaining with the NR routine gives the correct fit but takes forever to compute, basically freezing for several minutes (✨holy minimize, mighty va10a✨).

Also note that the hybridization gradient appears to work fine, so maybe the actual problem lies within grad_g0and_replica. If that's the case I'll hit the wall very soon, when testing the gradients for the new 'elemental' norm. We'll see...

Evidence

I did not include the real part in the figure for it is basically correct (but not totally: even with the numerical gradient I see some unwanted weight near the origin). Don't know if this could help, but I report it.
No restart file has been used, and here we are at one loop: we do not expect a brilliant fit, but a fair one yes.
I attach down here the input file, so that the problem could be reproduced with the current master head (dc046d8)

Notes

I've already added a warning for the user in commit b6351ff (fit_overhaul), we might want to copy those lines on master without waiting for merge, if this is considered urgent.

Link to the lines

inputHM2D.conf

 WMIXING=0.500000000                           !Mixing bath parameter
 TS=2.500000000E-01                            !hopping parameter
 NX=1                                          !Number of cluster sites in x direction
 NY=1                                          !Number of cluster sites in y direction
 NKX=30                                        !Number of kx point for BZ integration
 NKY=30                                        !Number of ky point for BZ integration
 NLAT=1                                        !Number of cluster sites
 NORB=1                                        !Number of impurity orbitals (max 5).
 NBATH=7                                       !Number of bath sites:(normal=>Nbath per orb)(hybrid=>Nbath total)(replica=>Nbath=Nreplica)
 NSPIN=1                                       !Number of spin degeneracy (max 2)
 ULOC=1.000000000,0.d0,0.d0,0.d0,0.d0          !Values of the local interaction per orbital (max 5)
 UST=0.d0                                      !Value of the inter-orbital interaction term
 JH=0.d0                                       !Hunds coupling
 JX=0.d0                                       !S-E coupling
 JP=0.d0                                       !P-H coupling
 BETA=300.000000000                            !Inverse temperature, at T=0 is used as a IR cut-off.
 XMU=0.d0                                      !Chemical potential. If HFMODE=T, xmu=0 indicates half-filling condition.
 NLOOP=1                                       !Max number of DMFT iterations.
 DMFT_ERROR=1.000000000E-05                    !Error threshold for DMFT convergence
 SB_FIELD=1.000000000E-01                      !Value of a symmetry breaking field for magnetic solutions.
 GF_FLAG=T                                     !flag to evaluate GFs and related quantities.
 DM_FLAG=F                                     !flag to evaluate the cluster density matrix \rho_IMP = Tr_BATH(\rho))
 ED_TWIN=F                                     !flag to reduce (T) or not (F,default) the number of visited sector using twin symmetry.
 ED_SECTORS=F                                  !flag to reduce sector scan for the spectrum to specific sectors +/- ed_sectors_shift.
 ED_SECTORS_SHIFT=1                            !shift to ed_sectors
 ED_SPARSE_H=T                                 !flag to select  storage of sparse matrix H (mem--, cpu++) if TRUE, or direct on-the-fly H*v product (mem++, cpu--) if FALSE
 ED_GF_SYMMETRIC=F                             !flag to assume Gij = Gji
 ED_PRINT_SIGMA=T                              !flag to print impurity Self-energies
 ED_PRINT_G=T                                  !flag to print impurity Greens function
 ED_PRINT_G0=T                                 !flag to print non-interacting impurity Greens function
 ED_VERBOSE=5                                  !Verbosity level: 0=almost nothing --> 5:all. Really: all
 NSUCCESS=1                                    !Number of successive iterations below threshold for convergence
 LMATS=5000                                    !Number of Matsubara frequencies.
 LREAL=5000                                    !Number of real-axis frequencies.
 LTAU=1024                                     !Number of imaginary time points.
 LFIT=1000                                     !Number of Matsubara frequencies used in the \Chi2 fit.
 NREAD=0.d0                                    !Objective density for fixed density calculations.
 NERR=1.000000000E-04                          !Error threshold for fixed density calculations.
 NDELTA=1.000000000E-01                        !Initial step for fixed density calculations.
 NCOEFF=1.000000000                            !multiplier for the initial ndelta read from a file (ndelta-->ndelta*ncoeff).
 WINI=-5.000000000                             !Smallest real-axis frequency
 WFIN=5.000000000                              !Largest real-axis frequency
 HFMODE=T                                      !Flag to set the Hartree form of the interaction (n-1/2). see xmu.
 EPS=1.000000000E-02                           !Broadening on the real-axis.
 CUTOFF=1.000000000E-09                        !Spectrum cut-off, used to determine the number states to be retained.
 GS_THRESHOLD=1.000000000E-09                  !Energy threshold for ground state degeneracy loop up
 HWBAND=2.000000000                            !half-bandwidth for the bath initialization: flat in -hwband:hwband
 LANC_METHOD=arpack                            !select the lanczos method to be used in the determination of the spectrum. ARPACK (default), LANCZOS (T=0 only), DVDSON (no MPI)
 LANC_NSTATES_SECTOR=2                         !Initial number of states per sector to be determined.
 LANC_NSTATES_TOTAL=1                          !Initial number of total states to be determined.
 LANC_NSTATES_STEP=2                           !Number of states added to the spectrum at each step.
 LANC_NCV_FACTOR=10                            !Set the size of the block used in Lanczos-Arpack by multiplying the required Neigen (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
 LANC_NCV_ADD=0                                !Adds up to the size of the block to prevent it to become too small (Ncv=lanc_ncv_factor*Neigen+lanc_ncv_add)
 LANC_NITER=512                                !Number of Lanczos iteration in spectrum determination.
 LANC_NGFITER=200                              !Number of Lanczos iteration in GF determination. Number of momenta.
 LANC_TOLERANCE=1.000000000E-12                !Tolerance for the Lanczos iterations as used in Arpack and plain lanczos.
 LANC_DIM_THRESHOLD=1024                       !Min dimension threshold to use Lanczos determination of the spectrum rather than Lapack based exact diagonalization.
 CG_METHOD=0                                   !Conjugate-Gradient method: 0=NR, 1=minimize.
 CG_GRAD=0                                     !Gradient evaluation method: 0=analytic (default), 1=numeric.
 CG_FTOL=1.000000000E-05                       !Conjugate-Gradient tolerance.
 CG_STOP=0                                     !Conjugate-Gradient stopping condition: 0-3, 0=C1.AND.C2, 1=C1, 2=C2 with C1=|F_n-1 -F_n|<tol*(1+F_n), C2=||x_n-1 -x_n||<tol*(1+||x_n||).
 CG_NITER=500                                  !Max. number of Conjugate-Gradient iterations.
 CG_WEIGHT=1                                   !Conjugate-Gradient weight form: 1=1.0, 2=1/n , 3=1/w_n.
 CG_SCHEME=weiss                               !Conjugate-Gradient fit scheme: delta or weiss.
 CG_POW=2                                      !Fit power for the calculation of the Chi distance function as 1/L*|G0 - G0and|**cg_pow
 CG_MINIMIZE_VER=F                             !Flag to pick old/.false. (Krauth) or new/.true. (Lichtenstein) version of the minimize CG routine
 CG_MINIMIZE_HH=1.000000000E-04                !Unknown parameter used in the CG minimize procedure.
 HFILE=hamiltonian                             !File where to retrieve/store the bath parameters.
 IMPHFILE=inputHLOC.in                         !File read the input local H.
 LOGFILE=6                                     !LOG unit.

Answer 1 · 2022-06-17T05:53:21.000Z

NEWS (relative to fit_overhaul branch)

Commits 1bab32a and 3240d40 have shown that the issue is wider (hence the title change). Here I report some evidence and try to wrap a brief recap.

RECAP

With numerical gradients everything works, with both "Frobenius" and "Elemental" definition of $\chi^2$ distance.

	Deltaˆ	Weiss
Elemental	🟡	🟢
Frobenius	🟡	🟢

ˆThe Delta fits are tagged "yellow" for they have worse quality wrt the Weiss ones, if comparing with NR-CG results with DMFT_ED code (NR-CG freezes with CDMFT code, consistently. I don't know why). Yet minimize-CG results are all on the same level accross the two codes, and very similar to the "yellow" ones (so we can take them as "fairly good"). More info on the freezing in 3240d40 commit message; plots of fitted functions are reported below. Note that instead minimize and NR give the very same results for the Weiss field, hence tagged green.

With analytical gradients we see some problems if we a) leave unchanged Frobenius gradients (wrt master branch) and b) port analytic gradients for elemental $\chi^2$ from current version in LIB_DMFT_ED. The situation is:

	Delta	Weiss
Elemental	🔴	🔴
Frobenius	🟡	🔴

If we swap sign in the elemental implementation of $\nabla\chi^2$ (thus diverging from LIB_DMFT_ED!) we get:

	Delta	Weiss
Elemental	🟡	🟡
Frobenius	🟡	🔴

If we further notice that commit 38bb300 had fixed the sign of \grad\chi^2(\Delta), but left untouched the sign of \grad\chi^2(g_0), so that applying the missing fix, we get:

	Delta	Weiss
Elemental	🟡	🟡
Frobenius	🟡	🟡

DETAILS

All plots with same input other than CG options, and at one loop.
Solid Line: FG
Dash-Dot: Fit

What I mean with "🔴" is

Delta	Weiss

What I mean with "🟡" is

Delta	Weiss

What I mean with "🟢" is

Deltaˆ	Weiss

ˆThis plot is the only one generated with DMFT_ED code, for this quality appears to be unreachable with minimize-CG, and NR-CG is de facto unavailable within CDMFT_ED.

Answer 2 · 2022-06-28T15:05:09.000Z

Relevant update in SciFortran: QcmPlab/SciFortran@c742471

Effects on this issue to be tested (it could probably fix the freezing with CG_METHOD=0 and CG_GRAD=1).

edit: it does not.