decompositions: IncompleteCholesky: Completed tests

Author: Lucas-Haubert

src/module.cpp

@@ -151,37 +151,39 @@ NB_MODULE(nanoeigenpy, m) {
       solvers, "LeastSquareDiagonalPreconditioner");
 #endif
 
+  using Eigen::Lower;
+
   using Eigen::BiCGSTAB;
   using Eigen::ConjugateGradient;
   using Eigen::DiagonalPreconditioner;
   using Eigen::IdentityPreconditioner;
   using Eigen::LeastSquareDiagonalPreconditioner;
   using Eigen::LeastSquaresConjugateGradient;
-  using Eigen::Lower;
   using Eigen::MINRES;
 
-  exposeMINRES<MINRES<Matrix, Lower>>(solvers, "MINRES");
-  exposeBiCGSTAB<BiCGSTAB<Matrix>>(solvers, "BiCGSTAB");
-  exposeConjugateGradient<ConjugateGradient<Matrix, Lower>>(
-      solvers, "ConjugateGradient");
-  exposeLeastSquaresConjugateGradient<LeastSquaresConjugateGradient<Matrix>>(
-      solvers, "LeastSquaresConjugateGradient");
-
-  using IdentityBiCGSTAB = BiCGSTAB<Matrix, IdentityPreconditioner>;
   using IdentityConjugateGradient =
       ConjugateGradient<Matrix, Lower, IdentityPreconditioner>;
   using IdentityLeastSquaresConjugateGradient =
       LeastSquaresConjugateGradient<Matrix, IdentityPreconditioner>;
   using DiagonalLeastSquaresConjugateGradient =
       LeastSquaresConjugateGradient<Matrix, DiagonalPreconditioner<Scalar>>;
+  using IdentityBiCGSTAB = BiCGSTAB<Matrix, IdentityPreconditioner>;
+  using DiagonalMINRES = MINRES<Matrix, Lower, DiagonalPreconditioner<Scalar>>;
 
-  exposeBiCGSTAB<IdentityBiCGSTAB>(solvers, "IdentityBiCGSTAB");
+  exposeConjugateGradient<ConjugateGradient<Matrix, Lower>>(
+      solvers, "ConjugateGradient");
   exposeConjugateGradient<IdentityConjugateGradient>(
       solvers, "IdentityConjugateGradient");
+  exposeLeastSquaresConjugateGradient<LeastSquaresConjugateGradient<Matrix>>(
+      solvers, "LeastSquaresConjugateGradient");
   exposeLeastSquaresConjugateGradient<IdentityLeastSquaresConjugateGradient>(
       solvers, "IdentityLeastSquaresConjugateGradient");
   exposeLeastSquaresConjugateGradient<DiagonalLeastSquaresConjugateGradient>(
       solvers, "DiagonalLeastSquaresConjugateGradient");
+  exposeMINRES<MINRES<Matrix, Lower>>(solvers, "MINRES");
+  exposeMINRES<DiagonalMINRES>(solvers, "DiagonalMINRES");
+  exposeBiCGSTAB<BiCGSTAB<Matrix>>(solvers, "BiCGSTAB");
+  exposeBiCGSTAB<IdentityBiCGSTAB>(solvers, "IdentityBiCGSTAB");
 
   exposeIncompleteLUT<SparseMatrix>(m, "IncompleteLUT");
   exposeIncompleteCholesky<SparseMatrix>(m, "IncompleteCholesky");

tests/test_incomplete_cholesky.py

@@ -2,11 +2,12 @@
 from scipy.sparse import csc_matrix
 import nanoeigenpy
 
+
 dim = 100
-rng = np.random.default_rng(42)
+rng = np.random.default_rng()
 
 A = rng.random((dim, dim))
-A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+A = (A + A.T) * 0.5 + np.diag(5.0 + rng.random(dim))
 A = csc_matrix(A)
 
 ichol = nanoeigenpy.IncompleteCholesky(A)
@@ -35,17 +36,34 @@
 
 ichol.analyzePattern(A)
 ichol.factorize(A)
+ichol.compute(A)
 assert ichol.info() == nanoeigenpy.ComputationInfo.Success
 
-ichol_shift = nanoeigenpy.IncompleteCholesky()
-ichol_shift.setInitialShift(1e-2)
-ichol_shift.compute(A)
-assert ichol_shift.info() == nanoeigenpy.ComputationInfo.Success
-
 L = ichol.matrixL()
-scaling = ichol.scalingS()
+S_diag = ichol.scalingS()
 perm = ichol.permutationP()
+P = perm.toDenseMatrix()
+
 assert isinstance(L, csc_matrix)
-assert isinstance(scaling, np.ndarray)
+assert isinstance(S_diag, np.ndarray)
 assert L.shape == (dim, dim)
-assert scaling.shape == (dim,)
+assert S_diag.shape == (dim,)
+
+L_dense = L.toarray()
+upper_part = np.triu(L_dense, k=1)
+assert np.allclose(upper_part, 0, atol=1e-12)
+
+assert np.all(S_diag > 0)
+
+S = csc_matrix((S_diag, (range(dim), range(dim))), shape=(dim, dim))
+
+PA = P @ A
+PAP = PA @ P.T
+SPAP = S @ PAP
+SPAPS = SPAP @ S
+
+LLT = L @ L.T
+
+diff = SPAPS - LLT
+relative_error = np.linalg.norm(diff.data) / np.linalg.norm(SPAPS.data)
+assert relative_error < 0.5

tests/test_iterative_solvers.py

@@ -6,52 +6,50 @@
 rng = np.random.default_rng()
 MAX_ITER = 8000
 
-_clazzes = [
+_classes = [
     nanoeigenpy.solvers.ConjugateGradient,
     nanoeigenpy.solvers.IdentityConjugateGradient,
     nanoeigenpy.solvers.LeastSquaresConjugateGradient,
     nanoeigenpy.solvers.IdentityLeastSquaresConjugateGradient,
     nanoeigenpy.solvers.DiagonalLeastSquaresConjugateGradient,
     nanoeigenpy.solvers.MINRES,
-    # nanoeigenpy.solvers.BiCGSTAB,
     nanoeigenpy.solvers.IdentityBiCGSTAB,
 ]
 
+_classes_bis = [
+    nanoeigenpy.solvers.DiagonalMINRES,
+    nanoeigenpy.solvers.BiCGSTAB,
+]
+
 
-@pytest.mark.parametrize("cls", _clazzes)
+@pytest.mark.parametrize("cls", _classes)
 def test_solver(cls):
     Q = rng.standard_normal((dim, dim))
     A = 0.5 * (Q.T + Q)
     solver = cls(A)
     solver.setMaxIterations(MAX_ITER)
 
-    # Vector rhs
-
     x = rng.random(dim)
     b = A.dot(x)
     x_est = solver.solve(b)
 
+    assert solver.info() == nanoeigenpy.ComputationInfo.Success
     assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6)
 
-    # Matrix rhs
-
     X = rng.random((dim, 20))
     B = A.dot(X)
     X_est = solver.solve(B)
 
     assert nanoeigenpy.is_approx(B, A.dot(X_est), 1e-6)
 
 
-@pytest.mark.parametrize("cls", _clazzes)
+@pytest.mark.parametrize("cls", _classes)
 def test_solver_with_guess(cls):
     Q = rng.standard_normal((dim, dim))
     A = 0.5 * (Q.T + Q)
     solver = cls(A)
     solver.setMaxIterations(MAX_ITER)
 
-    # With guess
-    # Vector rhs
-
     x = rng.random(dim)
     b = A.dot(x)
     x_est = solver.solveWithGuess(b, x + 0.01)
@@ -60,7 +58,49 @@ def test_solver_with_guess(cls):
     assert nanoeigenpy.is_approx(x, x_est, 1e-6)
     assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6)
 
-    # Matrix rhs
+    X = rng.random((dim, 20))
+    B = A.dot(X)
+    X_est = solver.solveWithGuess(B, X + 0.01)
+
+    assert nanoeigenpy.is_approx(X, X_est, 1e-6)
+    assert nanoeigenpy.is_approx(B, A.dot(X_est), 1e-6)
+
+
+@pytest.mark.parametrize("cls", _classes_bis)
+def test_solver_bis(cls):
+    Q = rng.standard_normal((dim, dim))
+    A = Q.T @ Q + np.eye(dim) * 0.1
+    solver = cls(A)
+    solver.setMaxIterations(MAX_ITER)
+
+    x = rng.random(dim)
+    b = A.dot(x)
+    x_est = solver.solve(b)
+
+    assert solver.info() == nanoeigenpy.ComputationInfo.Success
+    assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6)
+
+    X = rng.random((dim, 20))
+    B = A.dot(X)
+    X_est = solver.solve(B)
+
+    assert nanoeigenpy.is_approx(B, A.dot(X_est), 1e-6)
+
+
+@pytest.mark.parametrize("cls", _classes_bis)
+def test_solver_with_guess_bis(cls):
+    Q = rng.standard_normal((dim, dim))
+    A = Q.T @ Q + np.eye(dim) * 0.1
+    solver = cls(A)
+    solver.setMaxIterations(MAX_ITER)
+
+    x = rng.random(dim)
+    b = A.dot(x)
+    x_est = solver.solveWithGuess(b, x + 0.01)
+
+    assert solver.info() == nanoeigenpy.ComputationInfo.Success
+    assert nanoeigenpy.is_approx(x, x_est, 1e-6)
+    assert nanoeigenpy.is_approx(b, A.dot(x_est), 1e-6)
 
     X = rng.random((dim, 20))
     B = A.dot(X)