Stop ConjugateGradientSolver from inverting A on every solve

scala-source/com/barrybecker4/math/linear/ConjugateGradientSolver.scala

@@ -5,7 +5,6 @@ import javax.vecmath.{GMatrix, GVector}
 
 object ConjugateGradientSolver {
   private val DEFAULT_EPS = 0.000001
-  private val DEFAULT_MAX_ITERATIONS = 8
 }
 
 /**
@@ -19,7 +18,8 @@ object ConjugateGradientSolver {
   */
 class ConjugateGradientSolver(val matrix: GMatrix, val b: GVector) {
   private var eps = ConjugateGradientSolver.DEFAULT_EPS
-  private var maxIterations = ConjugateGradientSolver.DEFAULT_MAX_ITERATIONS
+  /** Classic CG may need up to n iterations for an n×n SPD system; allow a small cushion. */
+  private var maxIterations = 2 * matrix.getNumRow
 
   def setEpsilon(e: Double): Unit = { eps = e }
   def setMaxIterations(num: Int): Unit = { maxIterations = num }
@@ -29,7 +29,7 @@ class ConjugateGradientSolver(val matrix: GMatrix, val b: GVector) {
     solve(new GVector(zeros))
   }
 
-  /** Find a solution or return the initial guess if something goes wrong.
+  /** Find a solution or throw if the residual does not fall below [[eps]] within [[maxIterations]].
     * @param initialGuess the initial guess for the solution x, x0
     * @return solution vector
     */
@@ -45,10 +45,7 @@ class ConjugateGradientSolver(val matrix: GMatrix, val b: GVector) {
     val rnew = new GVector(p)
     val pnew = new GVector(p)
     val matrixMultp = new GVector(p)
-    val matrixInverse = new GMatrix(matrix)
-    matrixInverse.invert()
-    var error = 1.0
-    var norm = .0
+    var error = Math.sqrt(r.dot(r))
     var iteration = 0
     while (error > eps && iteration < maxIterations) {
       matrixMultp.mul(matrix, p)
@@ -61,9 +58,7 @@ class ConjugateGradientSolver(val matrix: GMatrix, val b: GVector) {
       p.set(pnew)
       r.set(rnew)
       x.set(xnew)
-      rnew.mul(r, matrixInverse)
-      norm = rnew.dot(r)
-      error = norm * norm
+      error = Math.sqrt(r.dot(r))
       iteration += 1
     }
 

scala-test/com/barrybecker4/math/linear/ConjugateGradientSolverSuite.scala

@@ -52,4 +52,27 @@ class ConjugateGradientSolverSuite extends AnyFunSuite {
     solver.setMaxIterations(3)
     assertThrows[IllegalStateException] { solver.solve }
   }
+
+  /** Diagonal SPD with distinct eigenvalues; classic CG needs up to n steps from x0=0. */
+  test("solveLargerDiagonalSystemWithoutMatrixInverse") {
+    val n = 12
+    val data = Array.fill(n * n)(0.0)
+    val bData = Array.tabulate(n)(i => (i + 1).toDouble)
+    for (i <- 0 until n) data(i * n + i) = (i + 1).toDouble
+    val A = new GMatrix(n, n, data)
+    val b = new GVector(bData)
+    solver = new ConjugateGradientSolver(A, b)
+    // Old default maxIterations=8 cannot finish an n=12 system that needs ~n steps.
+    solver.setMaxIterations(8)
+    assertThrows[IllegalStateException] { solver.solve }
+
+    solver = new ConjugateGradientSolver(A, b)
+    val solution = solver.solve
+    val expected = new GVector(Array.fill(n)(1.0))
+    assertTrue(LinearUtil.appxVectorsEqual(expected, solution, 0.000001), "Unexpected solution:" + solution)
+
+    val ax = new GVector(n)
+    ax.mul(A, solution)
+    assertTrue(LinearUtil.appxVectorsEqual(b, ax, 0.000001), "Residual Ax-b too large for " + solution)
+  }
 }