Class SymmLQ
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 = PT · 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) · PT · xhat = P · b, i.e. Ahat · xhat = bhat, where Ahat = P · (A - shift · I) · PT, bhat = P · b, and return the solution x = PT · xhat. The associated residual is rhat = bhat - Ahat · xhat = P · [b - (A - shift · I) · x] = P · r.
In the case of preconditioning, the IterativeLinearSolverEvent
s 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 countIn 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 x0 is known to be a good approximation to x, one should compute r0 = b - A · x, solve A · dx = r0, and set x = x0 + dx.
Exception context
Besides standard MathIllegalArgumentException
, this class might throw
MathIllegalArgumentException
if the linear operator or the
preconditioner are not symmetric.
- 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 xT · L · y ≠ yT · L · x (within a certain accuracy).
MathIllegalArgumentException
might also be thrown in case the
preconditioner is not positive definite.
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
-
Constructor Summary
ConstructorDescriptionSymmLQ
(int maxIterations, double delta, boolean check) Creates a new instance of this class, with default stopping criterion.SymmLQ
(IterationManager manager, double delta, boolean check) Creates a new instance of this class, with default stopping criterion and custom iteration manager. -
Method Summary
Modifier and TypeMethodDescriptionfinal boolean
Returnstrue
if symmetry of the matrix, and symmetry as well as positive definiteness of the preconditioner should be checked.Returns an estimate of the solution to the linear system A · x = b.solve
(RealLinearOperator a, RealLinearOperator m, RealVector b, boolean goodb, double shift) Returns an estimate of the solution to the linear system (A - shift · I) · x = b.solve
(RealLinearOperator a, RealLinearOperator m, RealVector b, RealVector x) Returns an estimate of the solution to the linear system A · x = b.Returns an estimate of the solution to the linear system A · x = b.solve
(RealLinearOperator a, RealVector b, boolean goodb, double shift) Returns the solution to the system (A - shift · I) · x = b.solve
(RealLinearOperator a, RealVector b, RealVector x) Returns an estimate of the solution to the linear system A · x = b.Returns an estimate of the solution to the linear system A · x = b.solveInPlace
(RealLinearOperator a, RealLinearOperator m, RealVector b, RealVector x, boolean goodb, double shift) Returns an estimate of the solution to the linear system (A - shift · I) · x = b.Returns an estimate of the solution to the linear system A · x = b.Methods inherited from class org.hipparchus.linear.PreconditionedIterativeLinearSolver
checkParameters
Methods inherited from class org.hipparchus.linear.IterativeLinearSolver
checkParameters, getIterationManager
-
Constructor Details
-
SymmLQ
public SymmLQ(int maxIterations, double delta, boolean check) Creates a new instance of this class, with default stopping criterion. Note that settingcheck
totrue
entails an extra matrix-vector product in the initial phase.- Parameters:
maxIterations
- the maximum number of iterationsdelta
- the δ parameter for the default stopping criterioncheck
-true
if self-adjointedness of both matrix and preconditioner should be checked
-
SymmLQ
Creates a new instance of this class, with default stopping criterion and custom iteration manager. Note that settingcheck
totrue
entails an extra matrix-vector product in the initial phase.- Parameters:
manager
- the custom iteration managerdelta
- the δ parameter for the default stopping criterioncheck
-true
if self-adjointedness of both matrix and preconditioner should be checked
-
-
Method Details
-
shouldCheck
public final boolean shouldCheck()Returnstrue
if symmetry of the matrix, and symmetry as well as positive definiteness of the preconditioner should be checked.- Returns:
true
if the tests are to be performed- Since:
- 1.4
-
solve
public RealVector solve(RealLinearOperator a, RealLinearOperator m, RealVector b) throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException, MathIllegalArgumentException Returns an estimate of the solution to the linear system A · x = b.- Overrides:
solve
in classPreconditionedIterativeLinearSolver
- Parameters:
a
- the linear operator A of the systemm
- the preconditioner, M (can benull
)b
- the right-hand side vector- Returns:
- a new vector containing the solution
- Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
orm
is not self-adjointMathIllegalArgumentException
- ifm
is not positive definiteMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-
solve
public RealVector solve(RealLinearOperator a, RealLinearOperator m, RealVector b, boolean goodb, double shift) throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system (A - shift · I) · x = b.If the solution x is expected to contain a large multiple of
b
(as in Rayleigh-quotient iteration), then better precision may be achieved withgoodb
set totrue
; this however requires an extra call to the preconditioner.shift
should be zero if the system A · x = b is to be solved. Otherwise, it could be an approximation to an eigenvalue of A, such as the Rayleigh quotient bT · A · b / (bT · b) corresponding to the vector b. If b is sufficiently like an eigenvector corresponding to an eigenvalue near shift, then the computed x may have very large components. When normalized, x may be closer to an eigenvector than b.- Parameters:
a
- the linear operator A of the systemm
- the preconditioner, M (can benull
)b
- the right-hand side vectorgoodb
- usuallyfalse
, except ifx
is expected to contain a large multiple ofb
shift
- the amount to be subtracted to all diagonal elements of A- Returns:
- a reference to
x
(shallow copy) - Throws:
NullArgumentException
- if one of the parameters isnull
MathIllegalArgumentException
- ifa
orm
is not squareMathIllegalArgumentException
- ifm
orb
have dimensions inconsistent witha
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
orm
is not self-adjointMathIllegalArgumentException
- ifm
is not positive definiteMathIllegalArgumentException
- ifa
is ill-conditioned
-
solve
public RealVector solve(RealLinearOperator a, RealLinearOperator m, RealVector b, RealVector x) throws MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system A · x = b.- Overrides:
solve
in classPreconditionedIterativeLinearSolver
- Parameters:
x
- not meaningful in this implementation; should not be considered as an initial guess (more)a
- the linear operator A of the systemm
- the preconditioner, M (can benull
)b
- the right-hand side vector- Returns:
- a new vector containing the solution
- Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
orm
is not self-adjointMathIllegalArgumentException
- ifm
is not positive definiteMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-
solve
public RealVector solve(RealLinearOperator a, RealVector b) throws MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system A · x = b.- Overrides:
solve
in classPreconditionedIterativeLinearSolver
- Parameters:
a
- the linear operator A of the systemb
- the right-hand side vector- Returns:
- a new vector containing the solution
- Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
is not self-adjointMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-
solve
public RealVector solve(RealLinearOperator a, RealVector b, boolean goodb, double shift) throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException Returns the solution to the system (A - shift · I) · x = b.If the solution x is expected to contain a large multiple of
b
(as in Rayleigh-quotient iteration), then better precision may be achieved withgoodb
set totrue
.shift
should be zero if the system A · x = b is to be solved. Otherwise, it could be an approximation to an eigenvalue of A, such as the Rayleigh quotient bT · A · b / (bT · b) corresponding to the vector b. If b is sufficiently like an eigenvector corresponding to an eigenvalue near shift, then the computed x may have very large components. When normalized, x may be closer to an eigenvector than b.- Parameters:
a
- the linear operator A of the systemb
- the right-hand side vectorgoodb
- usuallyfalse
, except ifx
is expected to contain a large multiple ofb
shift
- the amount to be subtracted to all diagonal elements of A- Returns:
- a reference to
x
- Throws:
NullArgumentException
- if one of the parameters isnull
MathIllegalArgumentException
- ifa
is not squareMathIllegalArgumentException
- ifb
has dimensions inconsistent witha
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
is not self-adjointMathIllegalArgumentException
- ifa
is ill-conditioned
-
solve
public RealVector solve(RealLinearOperator a, RealVector b, RealVector x) throws MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system A · x = b.- Overrides:
solve
in classPreconditionedIterativeLinearSolver
- Parameters:
x
- not meaningful in this implementation; should not be considered as an initial guess (more)a
- the linear operator A of the systemb
- the right-hand side vector- Returns:
- a new vector containing the solution
- Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
is not self-adjointMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-
solveInPlace
public RealVector solveInPlace(RealLinearOperator a, RealLinearOperator m, RealVector b, RealVector x) throws MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system A · x = b. The solution is computed in-place (initial guess is modified).- Specified by:
solveInPlace
in classPreconditionedIterativeLinearSolver
- Parameters:
x
- the vector to be updated with the solution;x
should not be considered as an initial guess (more)a
- the linear operator A of the systemm
- the preconditioner, M (can benull
)b
- the right-hand side vector- Returns:
- a reference to
x0
(shallow copy) updated with the solution - Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
orm
is not self-adjointMathIllegalArgumentException
- ifm
is not positive definiteMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-
solveInPlace
public RealVector solveInPlace(RealLinearOperator a, RealLinearOperator m, RealVector b, RealVector x, boolean goodb, double shift) throws MathIllegalArgumentException, NullArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system (A - shift · I) · x = b. The solution is computed in-place.If the solution x is expected to contain a large multiple of
b
(as in Rayleigh-quotient iteration), then better precision may be achieved withgoodb
set totrue
; this however requires an extra call to the preconditioner.shift
should be zero if the system A · x = b is to be solved. Otherwise, it could be an approximation to an eigenvalue of A, such as the Rayleigh quotient bT · A · b / (bT · b) corresponding to the vector b. If b is sufficiently like an eigenvector corresponding to an eigenvalue near shift, then the computed x may have very large components. When normalized, x may be closer to an eigenvector than b.- Parameters:
a
- the linear operator A of the systemm
- the preconditioner, M (can benull
)b
- the right-hand side vectorx
- the vector to be updated with the solution;x
should not be considered as an initial guess (more)goodb
- usuallyfalse
, except ifx
is expected to contain a large multiple ofb
shift
- the amount to be subtracted to all diagonal elements of A- Returns:
- a reference to
x
(shallow copy). - Throws:
NullArgumentException
- if one of the parameters isnull
MathIllegalArgumentException
- ifa
orm
is not squareMathIllegalArgumentException
- ifm
,b
orx
have dimensions inconsistent witha
.MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
orm
is not self-adjointMathIllegalArgumentException
- ifm
is not positive definiteMathIllegalArgumentException
- ifa
is ill-conditioned
-
solveInPlace
public RealVector solveInPlace(RealLinearOperator a, RealVector b, RealVector x) throws MathIllegalArgumentException, NullArgumentException, MathIllegalArgumentException, MathIllegalStateException Returns an estimate of the solution to the linear system A · x = b. The solution is computed in-place (initial guess is modified).- Overrides:
solveInPlace
in classPreconditionedIterativeLinearSolver
- Parameters:
x
- the vector to be updated with the solution;x
should not be considered as an initial guess (more)a
- the linear operator A of the systemb
- the right-hand side vector- Returns:
- a reference to
x0
(shallow copy) updated with the solution - Throws:
MathIllegalArgumentException
- ifshouldCheck()
istrue
, anda
is not self-adjointMathIllegalArgumentException
- ifa
is ill-conditionedNullArgumentException
- if one of the parameters isnull
MathIllegalStateException
- at exhaustion of the iteration count, unless a customcallback
has been set at construction of theIterationManager
-