Matrix diagonalization ====================== .. list-table:: Matrix diagonalization methods and options supported by KSSOLV2.0 :widths: 50 25 :header-rows: 1 * - Diagonalization method - eigmethod * - LOBPCG - 'lobpcg' * - Davidson - 'davidson' * - Chebyshev filtering - 'chebyfilt' * - Preconditioned Project Conjugate Gradient - 'ppcg' * - Krylov Schur - 'diagbyeigs' To seek the fixed point of the Kohn-Sham map .. math:: \rho = f_\mu [H(\rho)], we need to evaluate the matrix function :math:`f_\mu[\cdot]`, which is either a step function that jumps from 0 to 1 at :math:`\mu`, or a Fermi-Dirac distribution of the form .. math:: f_\mu [t] = \frac{1}{1+e^{(t-\mu)/k_B T}}, where :math`k_B` is the Boltzman constant and :math:`T` is the temperature. Such an evaluation can be performed by computing the eigenvalues and eigenvectors of the Kohn-Sham hamiltonian :math:`H`. Because the evaluation of :math:`f_{\mu}` does not need to be fully accurate when :math:`\rho` or :math:`X` is away from the ground state density or wavefunctions, and because a good initial guess from the previous SCF iteration is often available, an iterative diagonalization method is more appropriate. KSSOLV2.0 provides several functions/methods for computing the lowest :math:`k` eigenpairs of a fixed Kohn-Sham Hamiltonian :math:`H`, where :math:`k` is specified by the user. It is typically set to either the number of valence electron pairs or a number that is slightly larger (when a finite temperature SCF calculation is performed.) A user can choose a specific solver by setting the :attr:`eigmethod` field of the :attr:`options` class to one of the choices listed in the table above and passing :attr:`options` to the :attr:`updateX()` function. Locally Optimal Block Conjugate Gradient :attr:`lobpcg()` --------------------------------------------------------- LOBPCG is the default eigensolver used in KSSOLV to compute the lowest :math:`k` eigenpairs of a fixed Kohn-Sham Hamiltonian :math:`H`. The algorithm can be viewed as a way to solve the equivalent constrained minimization problem .. math:: \min_{X^*X=I} \mathrm{trace}\left( \frac{1}{2} X^*HX \right). :label: trmin Given a starting guess :math:`X^{(0)}` to the minimizer, LOBPCG iterates using the following updating formula: .. math:: X^{(i+1)} = X^{(i)} C_1 + W^{(i)} C_2 + P^{(i-1)} C_3, :label: update where :math:`C_1`, :math:`C_2`, and :math:`C_3` are coefficient matrices chosen to minimize the trace of :math:`H` within the space spanned by :math:`\{{X^{(i)}, W^{(i)}, P^{(i-1)}}\}` while maintaining the orthonormality constraint :math:`\langle X^{(i+1)}, X^{(i+1)} \rangle = I`. The matrix block :math:`W^{(i)}` is the preconditioned gradient of the Lagrangian :math:`\frac{1}{2} X^* H X - (X^*X - I)\Theta`, where :math:`\Theta = X^* H X`, i.e., .. math:: W^{(i)} = K^{-1} \left( HX^{(i)} - X^{(i)}\Theta^{(i)}\right), where :math:`K` is an appropriately chosen preconditioner. The matrix block :math:`P^{(i)}` defines the :math:`i` th search direction, which satisfies the recurrence formula .. math:: P^{(i)} = W^{(i)} C_2 + P^{(i-1)} C_3, with :math:`P^{(1)}` set to :math:`\emptyset`, and :math:`P^{(2)}` set to :math:`W^{(1)}`. The coefficient matrices :math:`C_1, C_2, C_3` are obtained by solving a projected eigenvalue problem .. math:: \left(Y^{*} H Y\right) C = \left(Y^{*}Y\right) C \Theta, :label: rr where :math:`Y = \left(X^{(i)}, W^{(i)}, P^{(i-1)}\right)`, :math:`C = (C_1^*, C_2^*, C_3^*)^*`, and :math:`\Theta` is a diagonal matrix containing the :math:`k` smallest eigenvalues of the project Hamiltonian. The solution of the projected eigenvalue problem :eq:`rr` and the udpate of the approximate eigenvector in eq:`update` are often collectively called a *Rayleigh Ritz* procedure. The LOBPCG solver can be called as follows .. code-block:: matlab >> [X,lambda] = lobpcg(H, X0, prec, tol, maxit, verbose); where :attr:`H` is a Hamiltonian object, :attr:`X0` is a wavefunction object that stores the initial guess to the desired eigenvectors, :attr:`prec` is a wavefunction object that stores a diagonal preconditioner (in reciprocal space), :attr:`tol` is a user specified convergence tolerance, :attr:`maxit` is the maximum number of iterations allowed, and :attr:`verbose` controls the amount of diagnostic output during the LOBPCG iteration. Krylov Schur :attr:`diagbyeigs()` ------------------------------- KSSOLV2.0 can use MATLAB's sparse eigensolver :attr:`eigs()` which implements the implicitly restarted Lanczos (or Krylov Schur) method. This is not the fastest solver because it does not take full advantage of approximations to several eigenvectors produced from a previous SCF iteration. Also, because the Hamiltonian is multiplied with one wavefunction (coefficients) at a time, there is limited opportunity for parallelization or data movement optimization. The Kyrlov Schur solver can be called as follows .. code-block:: matlab >> [X,lambda] = diagbyeigs(H, neig, tol, maxit); where :attr:`neig` is the number of eigenpairs to be computed, and :attr:`maxit` is the maximum number of implicit restart used in the Krylov Schur algorithm. Chebyshev Filtering :attr:`chebyfilt()` --------------------------------------- In Chebyshev filtering, we use a subspace iteration to compute dominant eigenvalues of :math:`T_d(H;\alpha,\lambda_{\mathrm{ub}})`, where :math:`T_d(t;\alpha,\lambda_{\mathrm{ub}})` is a :math:`d` th degree shifted and scaled Chebyshev polynomial that is bounded by -1 and 1 on the interval :math:`[\alpha,\lambda_{\mathrm{ub}}]`, and :math:`\alpha` is a carefully chosen parameter that is slightly above the largest eigenvalue of :math:`H` that we need to compute (e.g., the largest occupied state), and :math:`\lambda_{\mathrm{ub}}` is an upper bound of the spectrum of :math:`H`. The polynomial is chosen to amplify eigenvalues to the left of :math:`\alpha`. Chebyshev filtering can be called as follows .. code-block:: matlab >> [X,lambda] = chebyfilt(H, X0, degree, lb, ub); where :attr:`degree` is the degree of the polynomial, :attr:`lb` is the parameter :math:`\alpha` described above and :attr:`ub` is the upper bound of the spectrum :math:`\lambda_{\mathrm{ub}}`. Davidson :attr:`davidson()` --------------------------------------- To some extent, the Davidson's algorithm is a generalization of the LOBPCG algorithm. In the *i* th iteration, the LOBPCG algorithm extracts approximations to the desired eigenpairs from a subspace spanned by :math:`\{ W^{(i)}, X^{(i)}, X^{(i-1)} \}`. In a Davidson's algorithm, we can extract approximation from a larger subspace .. math:: \{W^{(i)}, X^{(i)}, X^{(i-1)}, X^{(i-2)},...\} However, as the dimension of the subspace becomes larger, the cost of solving the projected eigenvalue problem as well as maintaining an orthonormal basis of such as subspace can become much higher, and the algorithm can be come unstable due to the difficult in finding an orthonormal basis of the subspace. PPCG :attr:`ppcg()` --------------------------------------- The preconditioned projected conjugate gradient (PPCG) algorithm is a variant of the LOBPCG designed to reduce the cost of solving the projected problem :eq:`rr` when the number of eigenpairs to be computed is relatively large. If :math:`X^{(i)}` is sufficiently close to the solution of :eq:`trmin`, we can partition :math:`X^{(i)}` into several smaller blocks of vectors and apply LOBPCG to each block separately. The orthonormalization of all vectors in :math:`X^{(i)}` and the Rayleigh Ritz procedure is applied periodically to prevent vectors in different blocks from converging to the same set of eigenvectors.