Numerical instability of Davidson iterative eigensolver in solving LR-TDDFT excited states

[maki49]

Describe the bug

Now lr_solver=dav_subspace has mainly two problems in solving LR-TDDFT excited states:

1. "Assertion ( psi_norm > 0 ) failed" when doing normalizaiton in cal_grad (also appears in lr_solver=dav)
2. Large nonphysical negative eigenvalues (e.g., ~-1e10)

Both of them already exist before #5199, after which the problem 1 become worse.

To Reproduce

In tests/integrate/291_NO_KP_LR, enlarge lr_nstates (for example, 10), or use more k-points.

Possible Solution (with Refactor)

I'm going to try the following solutions that needs some refactor:

• Theoretically the number of initial-guess subspace basis vectors can be larger than the number of eigenpairs (bands or states) to be solved (expecially when the first several eigenpairs degenerate) , but in the present implementation, only n_bands initial vectors are copied in.
• Pass a precondition function similar to hpsi_func to enable custom preconditioner (which is hard-coded at present), such as sTDA or rid-preconditioner

[Cstandardlib]

• Theoretically the initial guess of subspace basis vectors can be larger than the number of eigenpairs (bands or states) to be solved (expecially when the first several eigenpairs degenerate) , but in the present implementation, only n_bands initial vectors are copied in.

What does it mean to pass larger initial guess of subspace basis than the number of eigenpairs (bands or states) to be solved?

• Is the size of subspace during the main loop n_initial_guess or n_band? For example we can use a subspace of size n_initial_guess when doing Rayleigh-Ritz and use n_band only when checking convergence.

[mohanchen]

Try to solve the issue by yourselves @maki49 @Cstandardlib @haozhihan

[haozhihan]

Will the same numerical instability be triggered for CG and Davidson?

[maki49]

Will the same numerical instability be triggered for CG and Davidson?

Both the two problems occur also in Davidson, but not so severe (or often).

I've tried CG in #5399 , no such problems:

• It canbe slower than direct diagonalization bylapack (with diag_ethr=1e-4), and some eigenvalues differ >0.01 Ry from lapack's result. Is it normal? Maybe the preconditioner can be optimized.

  case	matrix size	t_lapack(s)	t_cg(s)
  1(H2O)	24	45.64	86.12
  2(Si2)	48	31.83	53.14
  3(Si2)	96	60.04	43.04

• Since my HamiltLR does not inherit Hamilt, DiagoCG cannot benefit from diagH_subspace (which depends on Hamilt) as the subspace_func but can only use need_subspace=false now. (What's the improvement of subspace diagonalization in CG? )

  • Some refactor sugesstions are given in Feature: apply CG eigensolver to LR-TDDFT #5399.

Here's an example of full matrix which faces problem 1 in dav_subspace, for someone who tries to debug:

2.54191 0 -0.05226 0 0 -0.074615 0 0.0470716 0 0 -0.0211977 0 0.0202894 0 0.0840934 0 0 0.0502221 0 0 0 0.283082 0 0 
0 2.7113 0 0 -0.0073687 0 -0.0560734 0 -0.0167916 0 0 -0.0431683 0 -0.0422384 0 0 -0.0133889 0 0 0 0 0 0 0 
-0.05226 0 3.14399 0 0 -0.0365811 0 -0.0300141 0 0 -0.0532967 0 -0.0340724 0 0.0294075 0 0 0.0231556 0 0 0 -0.0462096 0 0 
0 0 0 3.12617 0 0 0 0 0 0 0 0 0 0 0 0.0111669 0 0 0.154337 0 -0.0471049 0 0 0.0586978 
0 -0.0073687 0 0 3.24875 0 0.0390991 0 -0.0376704 0 0 -0.00822859 0 0.017133 0 0 0.0551303 0 0 0 0 0 0 0 
-0.074615 0 -0.0365811 0 0 3.74872 0 0.0279213 0 0 -0.0630207 0 0.0274256 0 0.0670718 0 0 -0.0175326 0 0 0 0.0766089 0 0 
0 -0.0560734 0 0 0.0390991 0 1.1967 0 -0.0721625 0 0 -0.0703711 0 0.00916651 0 0 0.0280813 0 0 0 0 0 0 0 
0.0470716 0 -0.0300141 0 0 0.0279213 0 1.61779 0 0 -0.0290441 0 0.0655451 0 0.204046 0 0 0.0755009 0 0 0 0.272002 0 0 
0 -0.0167916 0 0 -0.0376704 0 -0.0721625 0 1.85009 0 0 -0.0574844 0 0.00280023 0 0 -0.00323539 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 1.82894 0 0 0 0 0 0 0 0 0 0.0520155 0 0 -0.0340305 0 
-0.0211977 0 -0.0532967 0 0 -0.0630207 0 -0.0290441 0 0 2.03901 0 -0.00669861 0 -0.0714484 0 0 -0.0249481 0 0 0 -0.123083 0 0 
0 -0.0431683 0 0 -0.00822859 0 -0.0703711 0 -0.0574844 0 0 2.58974 0 0.079856 0 0 0.0750716 0 0 0 0 0 0 0 
0.0202894 0 -0.0340724 0 0 0.0274256 0 0.0655451 0 0 -0.00669861 0 0.957622 0 0.0270841 0 0 -0.0797534 0 0 0 0.0622309 0 0 
0 -0.0422384 0 0 0.017133 0 0.00916651 0 0.00280023 0 0 0.079856 0 1.16322 0 0 0.0403892 0 0 0 0 0 0 0 
0.0840934 0 0.0294075 0 0 0.0670718 0 0.204046 0 0 -0.0714484 0 0.0270841 0 1.76232 0 0 0.0181566 0 0 0 0.328707 0 0 
0 0 0 0.0111669 0 0 0 0 0 0 0 0 0 0 0 1.52936 0 0 -0.0206693 0 0.0856325 0 0 0.00342213 
0 -0.0133889 0 0 0.0551303 0 0.0280813 0 -0.00323539 0 0 0.0750716 0 0.0403892 0 0 1.7002 0 0 0 0 0 0 0 
0.0502221 0 0.0231556 0 0 -0.0175326 0 0.0755009 0 0 -0.0249481 0 -0.0797534 0 0.0181566 0 0 2.24998 0 0 0 0.0973974 0 0 
0 0 0 -0.0125944 0 0 0 0 0 0 0 0 0 0 0 0.00654741 0 0 2.17657 0 0.0597456 0 0 0.185599 
0 0 0 0 0 0 0 0 0 0.0090608 0 0 0 0 0 0 0 0 0 2.59722 0 0 -0.00665076 0 
0 0 0 0.00111033 0 0 0 0 0 0 0 0 0 0 0 0.00610819 0 0 0.0597456 0 3.05256 0 0 0.0617511 
0.11615 0 0.00200569 0 0 0.0860229 0 0.229048 0 0 -0.085706 0 0.0894476 0 0.249182 0 0 0.0841823 0 0 0 3.52626 0 0 
0 0 0 0 0 0 0 0 0 0.00334601 0 0 0 0 0 0 0 0 0 -0.00665076 0 0 3.00593 0 
0 0 0 0.0681118 0 0 0 0 0 0 0 0 0 0 0 -0.009793 0 0 0.185599 0 0.0617511 0 0 3.97099 

[haozhihan]

#4068

[haozhihan]

#563

[Cstandardlib]

• It canbe slower than direct diagonalization bylapack (with diag_ethr=1e-4), and some eigenvalues differ >0.01 Ry from lapack's result. Is it normal? Maybe the preconditioner can be optimized.

Iterative solvers are usually intended for relatively large problems and many of them use a LAPACK syevd/hegvd routine for small projected eigenproblem inside.

For example the function diagH_subspace takes an initial guess block psi (I'll call it space V here) to perform a first so-called "subspace diag"
-- actually a Rayleigh-Ritz procedure, solution of a projected eigenvalue problem to extract approximations to the eigenpairs from somehow constructed subspace, here the initual guess.

We can see from the code that what diagH_subspace does is:

void DiagoIterAssist<T, Device>::diagH_subspace(const hamilt::Hamilt<T, Device>* const pHamilt, // hamiltonian operator carrier
                                                const psi::Psi<T, Device>& psi,     // [in] wavefunction
                                                psi::Psi<T, Device>& evc,           // [out] wavefunction
                                                Real* en,                           // [out] eigenvalues
                                                int n_band // [in] number of bands to be calculated, also number of rows
                                                           // of evc, if set to 0, n_band = nstart, default 0
)
{
    // two case:
    // 1. pw base: nstart = n_band, psi(nbands * npwx)
    // 2. lcao_in_pw base: nstart >= n_band, psi(NLOCAL * npwx)
...

    { // code block to calculate hcc and scc
        setmem_complex_op()(ctx, temp, 0, nstart * dmax);

        T* hphi = temp;
        // do hPsi for all bands
        psi::Range all_bands_range(1, psi.get_current_k(), 0, nstart - 1);
        hpsi_info hpsi_in(&psi, all_bands_range, hphi);
        pHamilt->ops->hPsi(hpsi_in);

        gemm_op<T, Device>()('C', 'N', psi.get_pointer(),hphi,  hcc);

        T* sphi = temp;
        // do sPsi for all bands
        pHamilt->sPsi(psi.get_pointer(), sphi, dmax, dmin, nstart);

        gemm_op<T, Device>()('C', 'N', psi, sphi, scc);
    }

    // after generation of H and S matrix, diag them
    DiagoIterAssist::diagH_LAPACK(nstart, n_band, hcc, scc, nstart, en, vcc);

    { // code block to calculate evc
        gemm_op<T, Device>()('N', 'N', psi.get_pointer(), vcc, evc);
    }
}

I extract only the core operation of this function.
It consists of 3 parts:

1. First construct the Gram Matrix $hcc=V^THV$ and $scc=V^TSV$
2. Second Call the diagH_LAPACK (wrapper of zhegvd inside) to solve the dense projected eigenproblem $(V^THV) y =(V^TSV) y \Lambda$ where $\Lambda$ is en , and $y$ is vcc in the code
3. Last update eigenvector approximation block evc by doing evc = psi * vcc to extract a new set of approximations $V_{new} = Vy$ from subspace V(i.e. psi input, or initial guess space)

This technique is usually applied to a subspace to get eigenpair approximations in all kinds of block iterative methods including Davidson,
and in my opinion it is nothing different than what we do in the main loop of Davidson, and I don't know why it's an independent part named as diagH_subspace. Maybe someone can help?

[maki49]

@Cstandardlib Thanks. So can I use diagH_subspace_init here, or use this instead of diagH_subspace as subspace_func in the CG call in hsolver_pw?

[Cstandardlib]

@Cstandardlib Thanks. So can I use diagH_subspace_init here, or use this instead of diagH_subspace as subspace_func in the CG call in hsolver_pw?

What's the difference between diagH_subspace_init and diagH_subspace?
Typically Rayleigh-Ritz procedure is used on iterative subspaces like Krylov subspace generated during iteration steps and classical CG method does not have this subspace diag procedure. Why do we introduce this part?

[maki49]

@Cstandardlib I refer to the code from here:

abacus-develop/source/module_hsolver/hsolver_pw.cpp

Lines 375 to 403 in e5768ec

	auto subspace_func = [hm, ngk_pointer](const ct::Tensor& psi_in, ct::Tensor& psi_out) {
	// psi_in should be a 2D tensor:
	// psi_in.shape() = [nbands, nbasis]
	const auto ndim = psi_in.shape().ndim();
	REQUIRES_OK(ndim == 2, "dims of psi_in should be less than or equal to 2");
	// Convert a Tensor object to a psi::Psi object
	auto psi_in_wrapper = psi::Psi<T, Device>(psi_in.data<T>(),
	1,
	psi_in.shape().dim_size(0),
	psi_in.shape().dim_size(1),
	ngk_pointer);
	auto psi_out_wrapper = psi::Psi<T, Device>(psi_out.data<T>(),
	1,
	psi_out.shape().dim_size(0),
	psi_out.shape().dim_size(1),
	ngk_pointer);
	auto eigen = ct::Tensor(ct::DataTypeToEnum<Real>::value,
	ct::DeviceType::CpuDevice,
	ct::TensorShape({psi_in.shape().dim_size(0)}));
	DiagoIterAssist<T, Device>::diagH_subspace(hm, psi_in_wrapper, psi_out_wrapper, eigen.data<Real>());
	};
	DiagoCG<T, Device> cg(this->basis_type,
	this->calculation_type,
	this->need_subspace,
	subspace_func,
	this->diag_thr,
	this->diag_iter_max,
	this->nproc_in_pool);

[Cstandardlib]

@maki49 Seems that this first-round subspace diagonalization is used to improve the initial guess passed to the iterative solvers.

I found changes to these two functions in #4047.

diagH_subspace_init is referred only in hsolver_lcaopw.cpp:

        /// solve eigenvector and eigenvalue for H(k)
        hsolver::DiagoIterAssist<T>::diagH_subspace_init(
            pHamilt,                 // interface to hamilt
            transform.get_pointer(), // transform matrix between lcao and pw
            transform.get_nbands(),
            transform.get_nbasis(),
            psi,                                  // psi in pw basis
            eigenvalues.data() + ik * pes->ekb.nc // eigenvalues
#ifdef __EXX
            ,
            add_exx_to_subspace_hamilt,
            set_exxlip_lcaowfc
#endif
        );

and it takes much more parameters as well as supporting more devices than diagH_subspace

void DiagoIterAssist<T, Device>::diagH_subspace_init(hamilt::Hamilt<T, Device>* pHamilt,
    const T* psi,
    int psi_nr,
    int psi_nc,
    psi::Psi<T, Device>& evc,
    Real* en,
    const std::function<void(T*, const int)>& add_to_hcc,
    const std::function<void(const T* const, const int, const int)>& export_vcc)
{
...
if (base_device::get_device_type(ctx) == base_device::GpuDevice)
...

In which field do these parameters affect this process? What's the motivation to replace old diagH_subspace with diagH_subspace_init?

[maki49]

Oh no... both diagH_subspace and diagH_subspace_init depends on Hamilt* and Psi ... I see.
So I will not use any of them for a while. Let's focus on how to improve the Davidson Algorithm.

[maki49]

Thanks for the help of @haozhihan @Cstandardlib. Now there seems no bug in the usage of dav_subspace solving LR-TDDFT any more.

However, the hard-coded formula here is unreasonable for LR-TDDFT and slows down the convergence, which needs further refactor:

abacus-develop/source/module_hsolver/diago_dav_subspace.cpp

Lines 342 to 344 in 7c46674

	// pre[i] = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]);
	double x = std::abs(this->precondition[i] - (*eigenvalue_iter)[m]);
	pre[i] = 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0)));