Implementation of the SYMMLQ iterative linear solver proposed by Paige and Saunders (1975). This implementation is largely based on the FORTRAN code by Pr. Michael A. Saunders, available here.

SYMMLQ is designed to solve the system of linear equations A · x = b where A is an n × n self-adjoint linear operator (defined as a RealLinearOperator), and b is a given vector. The operator A is not required to be positive definite. If A is known to be definite, the method of conjugate gradients might be preferred, since it will require about the same number of iterations as SYMMLQ but slightly less work per iteration.

SYMMLQ is designed to solve the system (A - shift · I) · x = b, where shift is a specified scalar value. If shift and b are suitably chosen, the computed vector x may approximate an (unnormalized) eigenvector of A, as in the methods of inverse iteration and/or Rayleigh-quotient iteration. Again, the linear operator (A - shift · I) need not be positive definite (but must be self-adjoint). The work per iteration is very slightly less if shift = 0.

Preconditioning

Preconditioning may reduce the number of iterations required. The solver may be provided with a positive definite preconditioner M = P^T · P that is known to approximate (A - shift · I)^-1 in some sense, where matrix-vector products of the form M · y = x can be computed efficiently. Then SYMMLQ will implicitly solve the system of equations P · (A - shift · I) · P^T · x_hat = P · b, i.e. A_hat · x_hat = b_hat, where A_hat = P · (A - shift · I) · P^T, b_hat = P · b, and return the solution x = P^T · x_hat. The associated residual is r_hat = b_hat - A_hat · x_hat = P · [b - (A - shift · I) · x] = P · r.

In the case of preconditioning, the IterativeLinearSolverEvents that this solver fires are such that IterativeLinearSolverEvent#getNormOfResidual() returns the norm of the preconditioned, updated residual, ||P · r||, not the norm of the true residual ||r||.

Default stopping criterion

A default stopping criterion is implemented. The iterations stop when || rhat || ≤ δ || Ahat || || xhat ||, where xhat is the current estimate of the solution of the transformed system, rhat the current estimate of the corresponding residual, and δ a user-specified tolerance.

Iteration count

In the present context, an iteration should be understood as one evaluation of the matrix-vector product A · x. The initialization phase therefore counts as one iteration. If the user requires checks on the symmetry of A, this entails one further matrix-vector product in the initial phase. This further product is not accounted for in the iteration count. In other words, the number of iterations required to reach convergence will be identical, whether checks have been required or not.

The present definition of the iteration count differs from that adopted in the original FOTRAN code, where the initialization phase was not taken into account.

Initial guess of the solution

The x parameter in

• solve(RealLinearOperator, RealVector, RealVector),
• solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)},
• solveInPlace(RealLinearOperator, RealVector, RealVector),
• solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector),
• solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, boolean, double),

should not be considered as an initial guess, as it is set to zero in the initial phase. If x₀ is known to be a good approximation to x, one should compute r₀ = b - A · x, solve A · dx = r0, and set x = x₀ + dx.

Exception context

Besides standard DimensionMismatchException, this class might throw NonSelfAdjointOperatorException if the linear operator or the preconditioner are not symmetric. In this case, the ExceptionContext provides more information

• key "operator" points to the offending linear operator, say L,
• key "vector1" points to the first offending vector, say x,
• key "vector2" points to the second offending vector, say y, such that x^T · L · y ≠ y^T · L · x (within a certain accuracy).

NonPositiveDefiniteOperatorException might also be thrown in case the preconditioner is not positive definite. The relevant keys to the ExceptionContext are

• key "operator", which points to the offending linear operator, say L,
• key "vector", which points to the offending vector, say x, such that x^T · L · x < 0.

References

Paige and Saunders (1975)

C. C. Paige and M. A. Saunders, Solution of Sparse Indefinite Systems of Linear Equations, SIAM Journal on Numerical Analysis 12(4): 617-629, 1975