Fix orthgonality loss in Davidson using Gram-Schmidt reorthogonalisation (#191)

* * Fix orthgonality loss in Davidson using Gram-Schmidt reorthogonalisation in
  case orthogonality loss is detected.

* * Differentiate between block size and number of guesses in Davidson
  procedure. In cases where n_guesses > n_states, this may lead to a decrease
  in the number of matrix applies needed till convergence.

* fix bug in call to estimate_n_guesses

* * Davidson block size as keyword argument 'n_block' and corresponding sanity
  checks.

* adapt tests to changed number of guesses

* move n_block sanity checks to eigsh

* only update entries for the last n_ss_added vectors in the projected matrix

* add more davidson tests

* more explicit n_block test

* add codecov token in CI

---------

Co-authored-by: jonasleitner <[email protected]>

Author: AdrianLDempwolff

.github/workflows/ci.yaml

@@ -89,7 +89,7 @@ jobs:
         if: matrix.documentation
       #
       - name: Upload coverage to codecov
-        # Note: lcov curerntly produces some error and therefore requires the keep-going flag to complete.
+        # Note: lcov currently produces some error and therefore requires the keep-going flag to complete.
         # Since the error might be the result of some gcc/gcov bug, I don't know how to resolve it currently.
         # lcov < 2.0 apparently did hide a lot of errors from the user so the problem might have been around for some time already...
         run: |
@@ -99,6 +99,8 @@ jobs:
           lcov --ignore-errors unused --remove coverage.info '/opt/*' '/Applications/*' '/Library/*' '/usr/*' "${HOME}"'/.cache/*' "${HOME}"'/.local/*' "${PWD}"'/build/*' "${PWD}"'/libadcc/tests/*' --output-file coverage.info
           lcov --list coverage.info
           codecov -X gcov -f coverage.info
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
         if: contains(matrix.os, 'ubuntu')
 
       - name: Upload coverage to coveralls

adcc/solver/SolverStateBase.py

@@ -39,6 +39,7 @@ def __init__(self, matrix):
         self.converged = False            # Flag whether iteration is converged
         self.n_iter = 0                   # Number of iterations
         self.n_applies = 0                # Number of applies
+        self.reortho_triggers = []        # List of reorthogonalisation triggers
         self.timer = Timer()              # Construct a new timer
 
     def describe(self):
@@ -56,6 +57,10 @@ def describe(self):
         text += "| {0:<41s}  {1:>15s} |\n".format(algorithm, conv)
         text += ("| {0:30s} n_iter={1:<3d}  n_applies={2:<5d} |\n"
                  "".format(problem[:30], self.n_iter, self.n_applies))
+        text += ("| n_reortho={0:<7d}  max_overlap_before_reortho={1:<10s}   |\n"
+                 "".format(len(self.reortho_triggers),
+                           "{:<10.4E}".format(max(self.reortho_triggers))
+                           if len(self.reortho_triggers) > 0 else "N/A"))
         text += "+" + 60 * "-" + "+\n"
         text += ("|  #     eigenvalue  res. norm       "
                  "dominant elements       |\n")

adcc/solver/davidson.py

@@ -75,7 +75,7 @@ def default_print(state, identifier, file=sys.stdout):
 
 
 # TODO This function should be merged with eigsh
-def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep,
+def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep, n_block,
                         is_converged, which, callback=None, preconditioner=None,
                         preconditioning_method="Davidson", debug_checks=False,
                         residual_min_norm=None, explicit_symmetrisation=None):
@@ -87,12 +87,15 @@ def davidson_iterations(matrix, state, max_subspace, max_iter, n_ep,
         Matrix to diagonalise
     state
         DavidsonState containing the eigenvector guess
-    max_subspace : int or NoneType, optional
+    max_subspace : int
         Maximal subspace size
-    max_iter : int, optional
+    max_iter : int
         Maximal number of iterations
-    n_ep : int or NoneType, optional
+    n_ep : int
         Number of eigenpairs to be computed
+    n_block : int
+        Davidson block size: the number of vectors that are added to the subspace
+        in each iteration
     is_converged
         Function to test for convergence
     callback : callable, optional
@@ -131,11 +134,11 @@ def callback(state, identifier):
     # The problem size
     n_problem = matrix.shape[1]
 
-    # The block size
-    n_block = len(state.subspace_vectors)
+    # The current subspace size == Number of guesses
+    n_ss_vec = len(state.subspace_vectors)
 
-    # The current subspace size
-    n_ss_vec = n_block
+    # Sanity checks for block size
+    assert n_block >= n_ep and n_block <= n_ss_vec
 
     # The current subspace
     SS = state.subspace_vectors
@@ -157,21 +160,23 @@ def callback(state, identifier):
         Ax = evaluate(matrix @ SS)
         state.n_applies += n_ss_vec
 
+    # Get the worksize view for the first iteration
+    Ass = Ass_cont[:n_ss_vec, :n_ss_vec]
+
+    # Initiall projection of Ax onto the subspace exploiting the hermiticity
+    with state.timer.record("projection"):
+        for i in range(n_ss_vec):
+            for j in range(i, n_ss_vec):
+                Ass[i, j] = SS[i] @ Ax[j]
+                if i != j:
+                    Ass[j, i] = Ass[i, j]
+
     while state.n_iter < max_iter:
         state.n_iter += 1
 
         assert len(SS) >= n_block
         assert len(SS) <= max_subspace
 
-        # Project A onto the subspace, keeping in mind
-        # that the values Ass[:-n_block, :-n_block] are already valid,
-        # since they have been computed in the previous iterations already.
-        with state.timer.record("projection"):
-            Ass = Ass_cont[:n_ss_vec, :n_ss_vec]  # Increase the work view size
-            for i in range(n_block):
-                Ass[:, -n_block + i] = Ax[-n_block + i] @ SS
-            Ass[-n_block:, :] = np.transpose(Ass[:, -n_block:])
-
         # Compute the which(== largest, smallest, ...) eigenpair of Ass
         # and the associated ritz vector as well as residual
         with state.timer.record("rayleigh_ritz"):
@@ -237,7 +242,10 @@ def form_residual(rval, rvec):
                 # Update projection of ADC matrix A onto subspace
                 Ass = Ass_cont[:n_ss_vec, :n_ss_vec]
                 for i in range(n_ss_vec):
-                    Ass[:, i] = Ax[i] @ SS
+                    for j in range(i, n_ss_vec):
+                        Ass[i, j] = SS[i] @ Ax[j]
+                        if i != j:
+                            Ass[j, i] = Ass[i, j]
             # continue to add residuals to space
 
         with state.timer.record("preconditioner"):
@@ -266,12 +274,29 @@ def form_residual(rval, rvec):
             n_ss_added = 0
             for i in range(n_block):
                 pvec = preconds[i]
-                # Project out the components of the current subspace
+                # Project out the components of the current subspace using
+                # conventional Gram-Schmidt (CGS) procedure.
                 # That is form (1 - SS * SS^T) * pvec = pvec + SS * (-SS^T * pvec)
                 coefficients = np.hstack(([1], -(pvec @ SS)))
                 pvec = lincomb(coefficients, [pvec] + SS, evaluate=True)
                 pnorm = np.sqrt(pvec @ pvec)
-                if pnorm > residual_min_norm:
+                if pnorm < residual_min_norm:
+                    continue
+                # Perform reorthogonalisation if loss of orthogonality is
+                # detected; this comes at the expense of computing n_ss_vec
+                # additional scalar products but avoids linear dependence
+                # within the subspace.
+                with state.timer.record("reorthogonalisation"):
+                    ss_overlap = np.array(pvec @ SS)
+                    max_ortho_loss = np.max(np.abs(ss_overlap)) / pnorm
+                    if max_ortho_loss > n_problem * eps:
+                        # Update pvec by instance reorthogonalised against SS
+                        # using a second CGS. Also update pnorm.
+                        coefficients = np.hstack(([1], -ss_overlap))
+                        pvec = lincomb(coefficients, [pvec] + SS, evaluate=True)
+                        pnorm = np.sqrt(pvec @ pvec)
+                        state.reortho_triggers.append(max_ortho_loss)
+                if pnorm >= residual_min_norm:
                     # Extend the subspace
                     SS.append(evaluate(pvec / pnorm))
                     n_ss_added += 1
@@ -284,8 +309,10 @@ def form_residual(rval, rvec):
                 state.subspace_orthogonality = np.max(np.abs(orth))
                 if state.subspace_orthogonality > n_problem * eps:
                     warnings.warn(la.LinAlgWarning(
-                        "Subspace in davidson has lost orthogonality. "
-                        "Expect inaccurate results."
+                        "Subspace in Davidson has lost orthogonality. "
+                        "Max. deviation from orthogonality is {:.4E}. "
+                        "Expect inaccurate results.".format(
+                            state.subspace_orthogonality)
                     ))
 
         if n_ss_added == 0:
@@ -300,12 +327,26 @@ def form_residual(rval, rvec):
                 "be aborted without convergence. Try a different guess."))
             return state
 
+        # Matrix applies for the new vectors
         with state.timer.record("projection"):
             Ax.extend(matrix @ SS[-n_ss_added:])
             state.n_applies += n_ss_added
 
+        # Update the worksize view for the next iteration
+        Ass = Ass_cont[:n_ss_vec, :n_ss_vec]
 
-def eigsh(matrix, guesses, n_ep=None, max_subspace=None,
+        # Project Ax onto the subspace, keeping in mind
+        # that the values Ass[:-n_ss_added, :-n_ss_added] are already valid,
+        # since they have been computed in the previous iterations already.
+        with state.timer.record("projection"):
+            for i in range(n_ss_vec - n_ss_added, n_ss_vec):
+                for j in range(i + 1):
+                    Ass[i, j] = SS[i] @ Ax[j]
+                    if i != j:
+                        Ass[j, i] = Ass[i, j]
+
+
+def eigsh(matrix, guesses, n_ep=None, n_block=None, max_subspace=None,
           conv_tol=1e-9, which="SA", max_iter=70,
           callback=None, preconditioner=None,
           preconditioning_method="Davidson", debug_checks=False,
@@ -320,6 +361,9 @@ def eigsh(matrix, guesses, n_ep=None, max_subspace=None,
         Guess vectors (fixes also the Davidson block size)
     n_ep : int or NoneType, optional
         Number of eigenpairs to be computed
+    n_block : int or NoneType, optional
+        The solver block size: the number of vectors that are added to the subspace
+        in each iteration
     max_subspace : int or NoneType, optional
         Maximal subspace size
     conv_tol : float, optional
@@ -364,11 +408,28 @@ def eigsh(matrix, guesses, n_ep=None, max_subspace=None,
     if n_ep is None:
         n_ep = len(guesses)
     elif n_ep > len(guesses):
-        raise ValueError("n_ep cannot exceed the number of guess vectors.")
+        raise ValueError(f"n_ep (= {n_ep}) cannot exceed the number of guess "
+                         f"vectors (= {len(guesses)}).")
+
+    if n_block is None:
+        n_block = n_ep
+    elif n_block < n_ep:
+        raise ValueError(f"n_block (= {n_block}) cannot be smaller than the number "
+                         f"of states requested (= {n_ep}).")
+    elif n_block > len(guesses):
+        raise ValueError(f"n_block (= {n_block}) cannot exceed the number of guess "
+                         f"vectors (= {len(guesses)}).")
+
     if not max_subspace:
         # TODO Arnoldi uses this:
         # max_subspace = max(2 * n_ep + 1, 20)
         max_subspace = max(6 * n_ep, 20, 5 * len(guesses))
+    elif max_subspace < 2 * n_block:
+        raise ValueError(f"max_subspace (= {max_subspace}) needs to be at least "
+                         f"twice as large as n_block (n_block = {n_block}).")
+    elif max_subspace < len(guesses):
+        raise ValueError(f"max_subspace (= {max_subspace}) cannot be smaller than "
+                         f"the number of guess vectors (= {len(guesses)}).")
 
     def convergence_test(state):
         state.residuals_converged = state.residual_norms < conv_tol
@@ -385,7 +446,7 @@ def convergence_test(state):
 
     state = DavidsonState(matrix, guesses)
     davidson_iterations(matrix, state, max_subspace, max_iter,
-                        n_ep=n_ep, is_converged=convergence_test,
+                        n_ep=n_ep, n_block=n_block, is_converged=convergence_test,
                         callback=callback, which=which,
                         preconditioner=preconditioner,
                         preconditioning_method=preconditioning_method,

adcc/tests/backends/backends_crossref_test.py

@@ -65,10 +65,14 @@ def test_adc2_h2o(self, system, case):
             pytest.skip("Veloxchem does not support f-functions. "
                         "Not enough backends available.")
 
-        # fewer states available for fc-fv-cvs
-        n_states = 5
-        if "fc" in case and "fv" in case and "cvs" in case:
-            n_states = 4
+        kwargs = {"n_singlets": 5}
+        # fewer states available for fc-fv-cvs (4) and fv-cvs (5)
+        if "fv" in case and "cvs" in case:
+            kwargs["n_singlets"] = 3
+            kwargs["n_guesses"] = 3
+        elif "cvs" in case:
+            # state 5 and 6 are degenerate -> can't compare the eigenvectors
+            kwargs["n_singlets"] = 4
 
         method = "cvs-adc2" if "cvs" in case else "adc2"
         core_orbitals = system.core_orbitals if "cvs" in case else None
@@ -79,9 +83,9 @@ def test_adc2_h2o(self, system, case):
         for b in backends_test:
             scfres = cached_backend_hf(b, system, conv_tol=1e-10)
             results[b] = adcc.run_adc(
-                scfres, method=method, n_singlets=n_states, conv_tol=1e-9,
+                scfres, method=method, conv_tol=1e-9,
                 core_orbitals=core_orbitals, frozen_core=frozen_core,
-                frozen_virtual=frozen_virtual
+                frozen_virtual=frozen_virtual, **kwargs
             )
             assert results[b].converged
         compare_adc_results(results, 5e-8)

adcc/tests/functionality_test.py

@@ -133,12 +133,16 @@ def test_functionality(self, system: str, case: str, method: str, kind: str,
         n_states = testcases.kinds_to_nstates([kind]).pop()
 
         kwargs = {n_states: 3}
-        # only few states available for h2o sto3g adc0/adc1
-        if system.name == "h2o" and system.basis == "sto-3g" and method.level < 2:
-            if "cvs" in case and "fv" in case:
-                kwargs[n_states] = 1
-            elif "cvs" in case:
-                kwargs[n_states] = 2
+        # only few states available for h2o sto3g
+        if system.name == "h2o" and system.basis == "sto-3g":
+            if method.level < 2:  # adc0/adc1
+                if "cvs" in case and "fv" in case:
+                    kwargs[n_states] = 1
+                elif "cvs" in case:
+                    kwargs[n_states] = 2
+            elif method.level < 4:  # adc2/adc3
+                if "cvs" in case and "fv" in case:  # only 5 states available
+                    kwargs["n_guesses"] = 3
 
         self.base_test(
             system=system, case=case, method=method.name, kind=kind,