org.apache.commons.math3.linear

## Class SymmLQ

- java.lang.Object
- - org.apache.commons.math3.linear.IterativeLinearSolver
  - - org.apache.commons.math3.linear.PreconditionedIterativeLinearSolver
    - - 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 = 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 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 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

  - 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).

  `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 xT · 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
  Since:
    3.0
  Version:
    $Id: SymmLQ.java 1416643 2012-12-03 19:37:14Z tn $

- - ### Nested Class Summary

    Nested Classes
    Modifier and Type | Class and Description
    --- | ---
    `private static class` | `SymmLQ.State` A simple container holding the non-final variables used in the iterations.

  - ### Field Summary

    Fields
    Modifier and Type | Field and Description
    --- | ---
    `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 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

    Constructors
    Constructor and Description
    `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

    Methods
    Modifier and Type | Method and Description
    --- | ---
    `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 m, RealVector b)` Returns an estimate of the solution to the linear system A · x = b.
    `RealVector` | `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.
    `RealVector` | `solve(RealLinearOperator a, RealLinearOperator m, 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 m, RealVector b, RealVector x)` Returns an estimate of the solution to the linear system A · x = b.
    `RealVector` | `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.
    `RealVector` | `solveInPlace(RealLinearOperator a, RealVector b, RealVector x)` Returns an estimate of the solution to the linear system A · x = b.

    - ### 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

    - #### 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

    - #### 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 m,
                         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
        `m` - the preconditioner, M (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 `m` is not self-adjoint
        `NonPositiveDefiniteOperatorException` - if `m` is not positive definite
        `IllConditionedOperatorException` - if `a` is ill-conditioned
        `NullArgumentException` - if one of the parameters is `null`
        `NonSquareOperatorException` - if `a` or `m` is not square
        `DimensionMismatchException` - if `m` or `b` have dimensions inconsistent with `a`
        `MaxCountExceededException` - at exhaustion of the iteration count, unless a custom `callback` has been set at construction of the `IterationManager`

    - #### solve

          public RealVector solve(RealLinearOperator a,
                         RealLinearOperator m,
                         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
        `m` - the preconditioner, M (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 `m` is not square
        `DimensionMismatchException` - if `m` or `b` have dimensions inconsistent with `a`
        `MaxCountExceededException` - at exhaustion of the iteration count, unless a custom `callback` has been set at construction of the `IterationManager`
        `NonSelfAdjointOperatorException` - if `getCheck()` is `true`, and `a` or `m` is not self-adjoint
        `NonPositiveDefiniteOperatorException` - if `m` is not positive definite
        `IllConditionedOperatorException` - if `a` is ill-conditioned

    - #### solve

          public RealVector solve(RealLinearOperator a,
                         RealLinearOperator m,
                         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
        `m` - the preconditioner, M (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 `m` is not self-adjoint
        `NonPositiveDefiniteOperatorException` - if `m` is not positive definite
        `IllConditionedOperatorException` - if `a` is ill-conditioned
        `NullArgumentException` - if one of the parameters is `null`
        `NonSquareOperatorException` - if `a` or `m` is not square
        `DimensionMismatchException` - if `m`, `b` or `x0` have dimensions inconsistent with `a`
        `MaxCountExceededException` - at exhaustion of the iteration count, unless a custom `callback` has been set at construction of the `IterationManager`

    - #### 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 of the `IterationManager`

    - #### 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 of the `IterationManager`
        `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 of the `IterationManager`

    - #### solveInPlace

          public RealVector solveInPlace(RealLinearOperator a,
                                RealLinearOperator m,
                                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
        `m` - the preconditioner, M (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 `m` is not self-adjoint
        `NonPositiveDefiniteOperatorException` - if `m` is not positive definite
        `IllConditionedOperatorException` - if `a` is ill-conditioned
        `NullArgumentException` - if one of the parameters is `null`
        `NonSquareOperatorException` - if `a` or `m` is not square
        `DimensionMismatchException` - if `m`, `b` or `x0` have dimensions inconsistent with `a`
        `MaxCountExceededException` - at exhaustion of the iteration count, unless a custom `callback` has been set at construction of the `IterationManager`

    - #### solveInPlace

          public RealVector solveInPlace(RealLinearOperator a,
                                RealLinearOperator m,
                                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
        `m` - the preconditioner, M (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 `m` is not square
        `DimensionMismatchException` - if `m`, `b` or `x` have dimensions inconsistent with `a`.
        `MaxCountExceededException` - at exhaustion of the iteration count, unless a custom `callback` has been set at construction of the `IterationManager`
        `NonSelfAdjointOperatorException` - if `getCheck()` is `true`, and `a` or `m` is not self-adjoint
        `NonPositiveDefiniteOperatorException` - if `m` 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 of the `IterationManager`

Copyright (c) 2003-2013 Apache Software Foundation