Add and test iterative solver

[rmsare]

Replace LU_Solver here with

• Gauss Seidel

or

• Overrelaxation

[rmsare]

Work in progress at https://github.com/rmsare/pymcc/tree/update-solver-iterative

[rmsare]

A few issues have come up that make iterative solvers a problematic choice:

• In general, there are zero values on the diagonal of L in the spline equations (for example, the eqns corresponding to the last three rows of L and v)

• If $\lambda = 0$ (no regularization), there may be more zeros than those three

• G-S and SOR cannot solve singular systems of equations (they don't converge)

[rmsare]

There are a couple options:

• Modify the spline eqns. Maybe we can just solve the L[0:p, 0:p], v[0:p] system?

• Regularize. Already doing this.

• Permute the rows of L?

• Another solver? GMRES?

[rmsare]

See https://github.com/rmsare/pymcc/blob/update-solver-iterative/tests/test-solvers.ipynb

[rmsare]

Update: Implemented a conjugate gradient solver and verified it works on a toy example.

It doesn't solve systems produced by the spline functions, though. It looks like it is failing where there are singular points?

[rmsare]

A few notes from today:

• The spline matrices are nonsingular, but they are also never positive definite (last three diagonal entries are always zero)

• Iterative methods require at least diagonally dominant matrices for convergence, so we can't use Jacobi, Gauss-Seidel, or SOR

• The condition number of spline matrices is large (e.g. 14182000944706), so conjugate gradient without preconditioning may not converge quickly

• Common preconditioners for conjugate gradient use the inverse of the diagonal elements (Jacobi) or a factorization (incomplete Cholesky) that requires nonzero diagonal elements. These won't work with the spline matrices

Things to investigate:

• Can we subsample the points for TPS to make it easier to solve?
• What preconditioning is scipy CG using? This is fine compared to my Boost implementation

Nice reference for making TPS more efficient:
https://vision.cornell.edu/se3/wp-content/uploads/2014/09/fulltext4.pdf

[rmsare]

Correction: GMRES solves these systems, CG does not

[rmsare]

Testing GMRES in Boost. One con: it also uses diagonal/Cholesky type preconditioners for faster convergence.

[rmsare]

Modern TPS implementation with sub sampling:

https://elonen.iki.fi/code/tpsdemo/2d-morph.cpp

https://elonen.iki.fi/code/tpsdemo/