org.apache.commons.math3.linear
Class SymmLQ

java.lang.Object
  extended by org.apache.commons.math3.linear.IterativeLinearSolver
      extended by org.apache.commons.math3.linear.PreconditionedIterativeLinearSolver
          extended by org.apache.commons.math3.linear.SymmLQ

public class SymmLQ
extends PreconditionedIterativeLinearSolver

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 = C · CT that is known to approximate (A - shift · I) in some sense, where systems of the form M · y = x can be solved efficiently. Then SYMMLQ will implicitly solve the system of equations P · (A - shift · I) · PT · xhat = P · b, i.e. Ahat · xhat = bhat, where P = C-1, 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.

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

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 DimensionMismatchException, this class might throw NonSelfAdjointOperatorException if the linear operator or the preconditioner are not symmetric. In this case, the ExceptionContext provides more information

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

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

Since:
3.0
Version:
$Id: SymmLQ.java 1295953 2012-03-01 22:30:26Z erans $

Nested Class Summary
private  class SymmLQ.State
           A simple container holding the non-final variables used in the iterations.
private static class SymmLQ.SymmLQEvent
          The type of all events fired by this implementation of the SYMMLQ method.
 
Field Summary
private static double CBRT_MACH_PREC
          The cubic root of MACH_PREC.
private  boolean check
          true if symmetry of matrix and conditioner must be checked.
private  double delta
          The value of the custom tolerance δ for the default stopping criterion.
private static double MACH_PREC
          The machine precision.
private static String OPERATOR
          Key for the exception context.
private static String THRESHOLD
          Key for the exception context.
private static String VECTOR
          Key for the exception context.
private static String VECTOR1
          Key for the exception context.
private static String VECTOR2
          Key for the exception context.
 
Constructor Summary
SymmLQ(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
private static void checkSymmetry(RealLinearOperator l, RealVector x, RealVector y, RealVector z)
          Performs a symmetry check on the specified linear operator, and throws an exception in case this check fails.
private static void daxpbypz(double a, RealVector x, double b, RealVector y, RealVector z)
          A BLAS-like function, for the operation z ← a · x + b · y + z.
private static void daxpy(double a, RealVector x, RealVector y)
          A clone of the BLAS DAXPY function, which carries out the operation y ← a · x + y.
 boolean getCheck()
          Returns true if symmetry of the matrix, and symmetry as well as positive definiteness of the preconditioner should be checked.
 RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b)
          Returns an estimate of the solution to the linear system A · x = b.
 RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b, boolean goodb, double shift)
          Returns an estimate of the solution to the linear system (A - shift · I) · x = b.
 RealVector solve(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x)
          Returns an estimate of the solution to the linear system A · x = b.
 RealVector solve(RealLinearOperator a, RealVector b)
          Returns an estimate of the solution to the linear system A · x = b.
 RealVector solve(RealLinearOperator a, RealVector b, boolean goodb, double shift)
          Returns the solution to the system (A - shift · I) · x = b.
 RealVector solve(RealLinearOperator a, RealVector b, RealVector x)
          Returns an estimate of the solution to the linear system A · x = b.
 RealVector solveInPlace(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x)
          Returns an estimate of the solution to the linear system A · x = b.
 RealVector solveInPlace(RealLinearOperator a, RealLinearOperator minv, RealVector b, RealVector x, boolean goodb, double shift)
          Returns an estimate of the solution to the linear system (A - shift · I) · x = b.
 RealVector solveInPlace(RealLinearOperator a, RealVector b, RealVector x)
          Returns an estimate of the solution to the linear system A · x = b.
private static void throwNPDLOException(RealLinearOperator l, RealVector v)
          Throws a new NonPositiveDefiniteOperatorException with appropriate context.
 
Methods inherited from class org.apache.commons.math3.linear.PreconditionedIterativeLinearSolver
checkParameters
 
Methods inherited from class org.apache.commons.math3.linear.IterativeLinearSolver
checkParameters, getIterationManager
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Field Detail

CBRT_MACH_PREC

private static final double CBRT_MACH_PREC
The cubic root of MACH_PREC.


MACH_PREC

private static final double MACH_PREC
The machine precision.


OPERATOR

private static final String OPERATOR
Key for the exception context.

See Also:
Constant Field Values

THRESHOLD

private static final String THRESHOLD
Key for the exception context.

See Also:
Constant Field Values

VECTOR

private static final String VECTOR
Key for the exception context.

See Also:
Constant Field Values

VECTOR1

private static final String VECTOR1
Key for the exception context.

See Also:
Constant Field Values

VECTOR2

private static final String VECTOR2
Key for the exception context.

See Also:
Constant Field Values

check

private final boolean check
true if symmetry of matrix and conditioner must be checked.


delta

private final double delta
The value of the custom tolerance δ for the default stopping criterion.

Constructor Detail

SymmLQ

public SymmLQ(int maxIterations,
              double delta,
              boolean check)
Creates a new instance of this class, with default stopping criterion. Note that setting check to true entails an extra matrix-vector product in the initial phase.

Parameters:
maxIterations - the maximum number of iterations
delta - the δ parameter for the default stopping criterion
check - true if self-adjointedness of both matrix and preconditioner should be checked

SymmLQ

public SymmLQ(IterationManager manager,
              double delta,
              boolean check)
Creates a new instance of this class, with default stopping criterion and custom iteration manager. Note that setting check to true entails an extra matrix-vector product in the initial phase.

Parameters:
manager - the custom iteration manager
delta - the δ parameter for the default stopping criterion
check - true if self-adjointedness of both matrix and preconditioner should be checked
Method Detail

checkSymmetry

private static void checkSymmetry(RealLinearOperator l,
                                  RealVector x,
                                  RealVector y,
                                  RealVector z)
                           throws NonSelfAdjointOperatorException
Performs a symmetry check on the specified linear operator, and throws an exception in case this check fails. Given a linear operator L, and a vector x, this method checks that x' · L · y = y' · L · x (within a given accuracy), where y = L · x.

Parameters:
l - the linear operator L
x - the candidate vector x
y - the candidate vector y = L · x
z - the vector z = L · y
Throws:
NonSelfAdjointOperatorException - when the test fails

daxpbypz

private static void daxpbypz(double a,
                             RealVector x,
                             double b,
                             RealVector y,
                             RealVector z)
A BLAS-like function, for the operation z ← a · x + b · y + z. This is for internal use only: no dimension checks are provided.

Parameters:
a - the scalar by which x is to be multiplied
x - the first vector to be added to z
b - the scalar by which y is to be multiplied
y - the second vector to be added to z
z - the vector to be incremented

daxpy

private static void daxpy(double a,
                          RealVector x,
                          RealVector y)
A clone of the BLAS DAXPY function, which carries out the operation y ← a · x + y. This is for internal use only: no dimension checks are provided.

Parameters:
a - the scalar by which x is to be multiplied
x - the vector to be added to y
y - the vector to be incremented

throwNPDLOException

private static void throwNPDLOException(RealLinearOperator l,
                                        RealVector v)
                                 throws NonPositiveDefiniteOperatorException
Throws a new NonPositiveDefiniteOperatorException with appropriate context.

Parameters:
l - the offending linear operator
v - the offending vector
Throws:
NonPositiveDefiniteOperatorException - in any circumstances

getCheck

public final boolean getCheck()
Returns true 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

solve

public RealVector solve(RealLinearOperator a,
                        RealLinearOperator minv,
                        RealVector b)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        MaxCountExceededException,
                        NonSelfAdjointOperatorException,
                        NonPositiveDefiniteOperatorException,
                        IllConditionedOperatorException
Returns an estimate of the solution to the linear system A · x = b.

Overrides:
solve in class PreconditionedIterativeLinearSolver
Parameters:
a - the linear operator A of the system
minv - the inverse of the preconditioner, M-1 (can be null)
b - the right-hand side vector
Returns:
a new vector containing the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a or minv is not self-adjoint
NonPositiveDefiniteOperatorException - if minv is not positive definite
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a or minv is not square
DimensionMismatchException - if minv or b have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction

solve

public RealVector solve(RealLinearOperator a,
                        RealLinearOperator minv,
                        RealVector b,
                        boolean goodb,
                        double shift)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        MaxCountExceededException,
                        NonSelfAdjointOperatorException,
                        NonPositiveDefiniteOperatorException,
                        IllConditionedOperatorException
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 with goodb set to true; 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 system
minv - the inverse of the preconditioner, M-1 (can be null)
b - the right-hand side vector
goodb - usually false, except if x is expected to contain a large multiple of b
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 is null
NonSquareOperatorException - if a or minv is not square
DimensionMismatchException - if minv or b have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction
NonSelfAdjointOperatorException - if getCheck() is true, and a or minv is not self-adjoint
NonPositiveDefiniteOperatorException - if minv is not positive definite
IllConditionedOperatorException - if a is ill-conditioned

solve

public RealVector solve(RealLinearOperator a,
                        RealLinearOperator minv,
                        RealVector b,
                        RealVector x)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        NonSelfAdjointOperatorException,
                        NonPositiveDefiniteOperatorException,
                        IllConditionedOperatorException,
                        MaxCountExceededException
Returns an estimate of the solution to the linear system A · x = b.

Overrides:
solve in class PreconditionedIterativeLinearSolver
Parameters:
x - not meaningful in this implementation; should not be considered as an initial guess (more)
a - the linear operator A of the system
minv - the inverse of the preconditioner, M-1 (can be null)
b - the right-hand side vector
Returns:
a new vector containing the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a or minv is not self-adjoint
NonPositiveDefiniteOperatorException - if minv is not positive definite
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a or minv is not square
DimensionMismatchException - if minv, b or x0 have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction

solve

public RealVector solve(RealLinearOperator a,
                        RealVector b)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        NonSelfAdjointOperatorException,
                        IllConditionedOperatorException,
                        MaxCountExceededException
Returns an estimate of the solution to the linear system A · x = b.

Overrides:
solve in class PreconditionedIterativeLinearSolver
Parameters:
a - the linear operator A of the system
b - the right-hand side vector
Returns:
a new vector containing the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a is not self-adjoint
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a is not square
DimensionMismatchException - if b has dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction

solve

public RealVector solve(RealLinearOperator a,
                        RealVector b,
                        boolean goodb,
                        double shift)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        NonSelfAdjointOperatorException,
                        IllConditionedOperatorException,
                        MaxCountExceededException
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 with goodb set to true.

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 system
b - the right-hand side vector
goodb - usually false, except if x is expected to contain a large multiple of b
shift - the amount to be subtracted to all diagonal elements of A
Returns:
a reference to x
Throws:
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a is not square
DimensionMismatchException - if b has dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction
NonSelfAdjointOperatorException - if getCheck() is true, and a is not self-adjoint
IllConditionedOperatorException - if a is ill-conditioned

solve

public RealVector solve(RealLinearOperator a,
                        RealVector b,
                        RealVector x)
                 throws NullArgumentException,
                        NonSquareOperatorException,
                        DimensionMismatchException,
                        NonSelfAdjointOperatorException,
                        IllConditionedOperatorException,
                        MaxCountExceededException
Returns an estimate of the solution to the linear system A · x = b.

Overrides:
solve in class PreconditionedIterativeLinearSolver
Parameters:
x - not meaningful in this implementation; should not be considered as an initial guess (more)
a - the linear operator A of the system
b - the right-hand side vector
Returns:
a new vector containing the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a is not self-adjoint
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a is not square
DimensionMismatchException - if b or x0 have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction

solveInPlace

public RealVector solveInPlace(RealLinearOperator a,
                               RealLinearOperator minv,
                               RealVector b,
                               RealVector x)
                        throws NullArgumentException,
                               NonSquareOperatorException,
                               DimensionMismatchException,
                               NonSelfAdjointOperatorException,
                               NonPositiveDefiniteOperatorException,
                               IllConditionedOperatorException,
                               MaxCountExceededException
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 class PreconditionedIterativeLinearSolver
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 system
minv - the inverse of the preconditioner, M-1 (can be null)
b - the right-hand side vector
Returns:
a reference to x0 (shallow copy) updated with the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a or minv is not self-adjoint
NonPositiveDefiniteOperatorException - if minv is not positive definite
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a or minv is not square
DimensionMismatchException - if minv, b or x0 have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction.

solveInPlace

public RealVector solveInPlace(RealLinearOperator a,
                               RealLinearOperator minv,
                               RealVector b,
                               RealVector x,
                               boolean goodb,
                               double shift)
                        throws NullArgumentException,
                               NonSquareOperatorException,
                               DimensionMismatchException,
                               NonSelfAdjointOperatorException,
                               NonPositiveDefiniteOperatorException,
                               IllConditionedOperatorException,
                               MaxCountExceededException
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 with goodb set to true; 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 system
minv - the inverse of the preconditioner, M-1 (can be null)
b - the right-hand side vector
x - the vector to be updated with the solution; x should not be considered as an initial guess (more)
goodb - usually false, except if x is expected to contain a large multiple of b
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 is null
NonSquareOperatorException - if a or minv is not square
DimensionMismatchException - if minv, b or x have dimensions inconsistent with a.
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction
NonSelfAdjointOperatorException - if getCheck() is true, and a or minv is not self-adjoint
NonPositiveDefiniteOperatorException - if minv is not positive definite
IllConditionedOperatorException - if a is ill-conditioned

solveInPlace

public RealVector solveInPlace(RealLinearOperator a,
                               RealVector b,
                               RealVector x)
                        throws NullArgumentException,
                               NonSquareOperatorException,
                               DimensionMismatchException,
                               NonSelfAdjointOperatorException,
                               IllConditionedOperatorException,
                               MaxCountExceededException
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 class PreconditionedIterativeLinearSolver
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 system
b - the right-hand side vector
Returns:
a reference to x0 (shallow copy) updated with the solution
Throws:
NonSelfAdjointOperatorException - if getCheck() is true, and a is not self-adjoint
IllConditionedOperatorException - if a is ill-conditioned
NullArgumentException - if one of the parameters is null
NonSquareOperatorException - if a is not square
DimensionMismatchException - if b or x0 have dimensions inconsistent with a
MaxCountExceededException - at exhaustion of the iteration count, unless a custom callback has been set at construction


Copyright (c) 2003-2013 Apache Software Foundation