diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/linear/SymmLQ.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/linear/SymmLQ.java | 1182 |
1 files changed, 1182 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/linear/SymmLQ.java b/src/main/java/org/apache/commons/math3/linear/SymmLQ.java new file mode 100644 index 0000000..bc7be0b --- /dev/null +++ b/src/main/java/org/apache/commons/math3/linear/SymmLQ.java @@ -0,0 +1,1182 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.linear; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.util.ExceptionContext; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.IterationManager; +import org.apache.commons.math3.util.MathUtils; + +/** + * Implementation of the SYMMLQ iterative linear solver proposed by <a href="#PAIG1975">Paige and + * Saunders (1975)</a>. This implementation is largely based on the FORTRAN code by Pr. Michael A. + * Saunders, available <a href="http://www.stanford.edu/group/SOL/software/symmlq/f77/">here</a>. + * + * <p>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 {@link 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. + * + * <p>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 <em>must</em> be self-adjoint). The work per iteration is very slightly + * less if shift = 0. + * + * <h3>Preconditioning</h3> + * + * <p>Preconditioning may reduce the number of iterations required. The solver may be provided with + * a positive definite preconditioner M = P<sup>T</sup> · P that is known to approximate (A - + * shift · I)<sup>-1</sup> 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<sup>T</sup> · x<sub>hat</sub> = P · b, + * i.e. A<sub>hat</sub> · x<sub>hat</sub> = b<sub>hat</sub>, where A<sub>hat</sub> = P + * · (A - shift · I) · P<sup>T</sup>, b<sub>hat</sub> = P · b, and + * return the solution x = P<sup>T</sup> · x<sub>hat</sub>. The associated residual is + * r<sub>hat</sub> = b<sub>hat</sub> - A<sub>hat</sub> · x<sub>hat</sub> = P · [b - (A + * - shift · I) · x] = P · r. + * + * <p>In the case of preconditioning, the {@link IterativeLinearSolverEvent}s that this solver fires + * are such that {@link IterativeLinearSolverEvent#getNormOfResidual()} returns the norm of the + * <em>preconditioned</em>, updated residual, ||P · r||, not the norm of the <em>true</em> + * residual ||r||. + * + * <h3><a id="stopcrit">Default stopping criterion</a></h3> + * + * <p>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. + * + * <h3>Iteration count</h3> + * + * <p>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 <em>not</em> 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. + * + * <p>The present definition of the iteration count differs from that adopted in the original FOTRAN + * code, where the initialization phase was <em>not</em> taken into account. + * + * <h3><a id="initguess">Initial guess of the solution</a></h3> + * + * <p>The {@code x} parameter in + * + * <ul> + * <li>{@link #solve(RealLinearOperator, RealVector, RealVector)}, + * <li>{@link #solve(RealLinearOperator, RealLinearOperator, RealVector, RealVector)}}, + * <li>{@link #solveInPlace(RealLinearOperator, RealVector, RealVector)}, + * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector)}, + * <li>{@link #solveInPlace(RealLinearOperator, RealLinearOperator, RealVector, RealVector, + * boolean, double)}, + * </ul> + * + * should not be considered as an initial guess, as it is set to zero in the initial phase. If + * x<sub>0</sub> is known to be a good approximation to x, one should compute r<sub>0</sub> = b - A + * · x, solve A · dx = r0, and set x = x<sub>0</sub> + dx. + * + * <h3><a id="context">Exception context</a></h3> + * + * <p>Besides standard {@link DimensionMismatchException}, this class might throw {@link + * NonSelfAdjointOperatorException} if the linear operator or the preconditioner are not symmetric. + * In this case, the {@link ExceptionContext} provides more information + * + * <ul> + * <li>key {@code "operator"} points to the offending linear operator, say L, + * <li>key {@code "vector1"} points to the first offending vector, say x, + * <li>key {@code "vector2"} points to the second offending vector, say y, such that x<sup>T</sup> + * · L · y ≠ y<sup>T</sup> · L · x (within a certain accuracy). + * </ul> + * + * <p>{@link NonPositiveDefiniteOperatorException} might also be thrown in case the preconditioner + * is not positive definite. The relevant keys to the {@link ExceptionContext} are + * + * <ul> + * <li>key {@code "operator"}, which points to the offending linear operator, say L, + * <li>key {@code "vector"}, which points to the offending vector, say x, such that x<sup>T</sup> + * · L · x < 0. + * </ul> + * + * <h3>References</h3> + * + * <dl> + * <dt><a id="PAIG1975">Paige and Saunders (1975)</a> + * <dd>C. C. Paige and M. A. Saunders, <a + * href="http://www.stanford.edu/group/SOL/software/symmlq/PS75.pdf"><em> Solution of Sparse + * Indefinite Systems of Linear Equations</em></a>, SIAM Journal on Numerical Analysis 12(4): + * 617-629, 1975 + * </dl> + * + * @since 3.0 + */ +public class SymmLQ extends PreconditionedIterativeLinearSolver { + + /* + * IMPLEMENTATION NOTES + * -------------------- + * The implementation follows as closely as possible the notations of Paige + * and Saunders (1975). Attention must be paid to the fact that some + * quantities which are relevant to iteration k can only be computed in + * iteration (k+1). Therefore, minute attention must be paid to the index of + * each state variable of this algorithm. + * + * 1. Preconditioning + * --------------- + * The Lanczos iterations associated with Ahat and bhat read + * beta[1] = ||P * b|| + * v[1] = P * b / beta[1] + * beta[k+1] * v[k+1] = Ahat * v[k] - alpha[k] * v[k] - beta[k] * v[k-1] + * = P * (A - shift * I) * P' * v[k] - alpha[k] * v[k] + * - beta[k] * v[k-1] + * Multiplying both sides by P', we get + * beta[k+1] * (P' * v)[k+1] = M * (A - shift * I) * (P' * v)[k] + * - alpha[k] * (P' * v)[k] + * - beta[k] * (P' * v[k-1]), + * and + * alpha[k+1] = v[k+1]' * Ahat * v[k+1] + * = v[k+1]' * P * (A - shift * I) * P' * v[k+1] + * = (P' * v)[k+1]' * (A - shift * I) * (P' * v)[k+1]. + * + * In other words, the Lanczos iterations are unchanged, except for the fact + * that we really compute (P' * v) instead of v. It can easily be checked + * that all other formulas are unchanged. It must be noted that P is never + * explicitly used, only matrix-vector products involving are invoked. + * + * 2. Accounting for the shift parameter + * ---------------------------------- + * Is trivial: each time A.operate(x) is invoked, one must subtract shift * x + * to the result. + * + * 3. Accounting for the goodb flag + * ----------------------------- + * When goodb is set to true, the component of xL along b is computed + * separately. From Paige and Saunders (1975), equation (5.9), we have + * wbar[k+1] = s[k] * wbar[k] - c[k] * v[k+1], + * wbar[1] = v[1]. + * Introducing wbar2[k] = wbar[k] - s[1] * ... * s[k-1] * v[1], it can + * easily be verified by induction that wbar2 follows the same recursive + * relation + * wbar2[k+1] = s[k] * wbar2[k] - c[k] * v[k+1], + * wbar2[1] = 0, + * and we then have + * w[k] = c[k] * wbar2[k] + s[k] * v[k+1] + * + s[1] * ... * s[k-1] * c[k] * v[1]. + * Introducing w2[k] = w[k] - s[1] * ... * s[k-1] * c[k] * v[1], we find, + * from (5.10) + * xL[k] = zeta[1] * w[1] + ... + zeta[k] * w[k] + * = zeta[1] * w2[1] + ... + zeta[k] * w2[k] + * + (s[1] * c[2] * zeta[2] + ... + * + s[1] * ... * s[k-1] * c[k] * zeta[k]) * v[1] + * = xL2[k] + bstep[k] * v[1], + * where xL2[k] is defined by + * xL2[0] = 0, + * xL2[k+1] = xL2[k] + zeta[k+1] * w2[k+1], + * and bstep is defined by + * bstep[1] = 0, + * bstep[k] = bstep[k-1] + s[1] * ... * s[k-1] * c[k] * zeta[k]. + * We also have, from (5.11) + * xC[k] = xL[k-1] + zbar[k] * wbar[k] + * = xL2[k-1] + zbar[k] * wbar2[k] + * + (bstep[k-1] + s[1] * ... * s[k-1] * zbar[k]) * v[1]. + */ + + /** + * A simple container holding the non-final variables used in the iterations. Making the current + * state of the solver visible from the outside is necessary, because during the iterations, + * {@code x} does not <em>exactly</em> hold the current estimate of the solution. Indeed, {@code + * x} needs in general to be moved from the LQ point to the CG point. Besides, additional + * upudates must be carried out in case {@code goodb} is set to {@code true}. + * + * <p>In all subsequent comments, the description of the state variables refer to their value + * after a call to {@link #update()}. In these comments, k is the current number of evaluations + * of matrix-vector products. + */ + private static class State { + /** The cubic root of {@link #MACH_PREC}. */ + static final double CBRT_MACH_PREC; + + /** The machine precision. */ + static final double MACH_PREC; + + /** Reference to the linear operator. */ + private final RealLinearOperator a; + + /** Reference to the right-hand side vector. */ + private final RealVector b; + + /** {@code true} if symmetry of matrix and conditioner must be checked. */ + private final boolean check; + + /** The value of the custom tolerance δ for the default stopping criterion. */ + private final double delta; + + /** The value of beta[k+1]. */ + private double beta; + + /** The value of beta[1]. */ + private double beta1; + + /** The value of bstep[k-1]. */ + private double bstep; + + /** The estimate of the norm of P * rC[k]. */ + private double cgnorm; + + /** The value of dbar[k+1] = -beta[k+1] * c[k-1]. */ + private double dbar; + + /** The value of gamma[k] * zeta[k]. Was called {@code rhs1} in the initial code. */ + private double gammaZeta; + + /** The value of gbar[k]. */ + private double gbar; + + /** The value of max(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ + private double gmax; + + /** The value of min(|alpha[1]|, gamma[1], ..., gamma[k-1]). */ + private double gmin; + + /** Copy of the {@code goodb} parameter. */ + private final boolean goodb; + + /** {@code true} if the default convergence criterion is verified. */ + private boolean hasConverged; + + /** The estimate of the norm of P * rL[k-1]. */ + private double lqnorm; + + /** Reference to the preconditioner, M. */ + private final RealLinearOperator m; + + /** The value of (-eps[k+1] * zeta[k-1]). Was called {@code rhs2} in the initial code. */ + private double minusEpsZeta; + + /** The value of M * b. */ + private final RealVector mb; + + /** The value of beta[k]. */ + private double oldb; + + /** The value of beta[k] * M^(-1) * P' * v[k]. */ + private RealVector r1; + + /** The value of beta[k+1] * M^(-1) * P' * v[k+1]. */ + private RealVector r2; + + /** + * The value of the updated, preconditioned residual P * r. This value is given by {@code + * min(}{@link #cgnorm}{@code , }{@link #lqnorm}{@code )}. + */ + private double rnorm; + + /** Copy of the {@code shift} parameter. */ + private final double shift; + + /** The value of s[1] * ... * s[k-1]. */ + private double snprod; + + /** + * An estimate of the square of the norm of A * V[k], based on Paige and Saunders (1975), + * equation (3.3). + */ + private double tnorm; + + /** + * The value of P' * wbar[k] or P' * (wbar[k] - s[1] * ... * s[k-1] * v[1]) if {@code goodb} + * is {@code true}. Was called {@code w} in the initial code. + */ + private RealVector wbar; + + /** + * A reference to the vector to be updated with the solution. Contains the value of xL[k-1] + * if {@code goodb} is {@code false}, (xL[k-1] - bstep[k-1] * v[1]) otherwise. + */ + private final RealVector xL; + + /** The value of beta[k+1] * P' * v[k+1]. */ + private RealVector y; + + /** The value of zeta[1]^2 + ... + zeta[k-1]^2. */ + private double ynorm2; + + /** The value of {@code b == 0} (exact floating-point equality). */ + private boolean bIsNull; + + static { + MACH_PREC = FastMath.ulp(1.); + CBRT_MACH_PREC = FastMath.cbrt(MACH_PREC); + } + + /** + * Creates and inits to k = 1 a new instance of this class. + * + * @param a the linear operator A of the system + * @param m the preconditioner, M (can be {@code null}) + * @param b the right-hand side vector + * @param goodb usually {@code false}, except if {@code x} is expected to contain a large + * multiple of {@code b} + * @param shift the amount to be subtracted to all diagonal elements of A + * @param delta the δ parameter for the default stopping criterion + * @param check {@code true} if self-adjointedness of both matrix and preconditioner should + * be checked + */ + State( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final boolean goodb, + final double shift, + final double delta, + final boolean check) { + this.a = a; + this.m = m; + this.b = b; + this.xL = new ArrayRealVector(b.getDimension()); + this.goodb = goodb; + this.shift = shift; + this.mb = m == null ? b : m.operate(b); + this.hasConverged = false; + this.check = check; + this.delta = delta; + } + + /** + * 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. + * + * @param l the linear operator L + * @param x the candidate vector x + * @param y the candidate vector y = L · x + * @param z the vector z = L · y + * @throws NonSelfAdjointOperatorException when the test fails + */ + private static void checkSymmetry( + final RealLinearOperator l, + final RealVector x, + final RealVector y, + final RealVector z) + throws NonSelfAdjointOperatorException { + final double s = y.dotProduct(y); + final double t = x.dotProduct(z); + final double epsa = (s + MACH_PREC) * CBRT_MACH_PREC; + if (FastMath.abs(s - t) > epsa) { + final NonSelfAdjointOperatorException e; + e = new NonSelfAdjointOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(SymmLQ.OPERATOR, l); + context.setValue(SymmLQ.VECTOR1, x); + context.setValue(SymmLQ.VECTOR2, y); + context.setValue(SymmLQ.THRESHOLD, Double.valueOf(epsa)); + throw e; + } + } + + /** + * Throws a new {@link NonPositiveDefiniteOperatorException} with appropriate context. + * + * @param l the offending linear operator + * @param v the offending vector + * @throws NonPositiveDefiniteOperatorException in any circumstances + */ + private static void throwNPDLOException(final RealLinearOperator l, final RealVector v) + throws NonPositiveDefiniteOperatorException { + final NonPositiveDefiniteOperatorException e; + e = new NonPositiveDefiniteOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(OPERATOR, l); + context.setValue(VECTOR, v); + throw e; + } + + /** + * A clone of the BLAS {@code DAXPY} function, which carries out the operation y ← a + * · x + y. This is for internal use only: no dimension checks are provided. + * + * @param a the scalar by which {@code x} is to be multiplied + * @param x the vector to be added to {@code y} + * @param y the vector to be incremented + */ + private static void daxpy(final double a, final RealVector x, final RealVector y) { + final int n = x.getDimension(); + for (int i = 0; i < n; i++) { + y.setEntry(i, a * x.getEntry(i) + y.getEntry(i)); + } + } + + /** + * A BLAS-like function, for the operation z ← a · x + b · y + z. This is + * for internal use only: no dimension checks are provided. + * + * @param a the scalar by which {@code x} is to be multiplied + * @param x the first vector to be added to {@code z} + * @param b the scalar by which {@code y} is to be multiplied + * @param y the second vector to be added to {@code z} + * @param z the vector to be incremented + */ + private static void daxpbypz( + final double a, + final RealVector x, + final double b, + final RealVector y, + final RealVector z) { + final int n = z.getDimension(); + for (int i = 0; i < n; i++) { + final double zi; + zi = a * x.getEntry(i) + b * y.getEntry(i) + z.getEntry(i); + z.setEntry(i, zi); + } + } + + /** + * Move to the CG point if it seems better. In this version of SYMMLQ, the convergence tests + * involve only cgnorm, so we're unlikely to stop at an LQ point, except if the iteration + * limit interferes. + * + * <p>Additional upudates are also carried out in case {@code goodb} is set to {@code true}. + * + * @param x the vector to be updated with the refined value of xL + */ + void refineSolution(final RealVector x) { + final int n = this.xL.getDimension(); + if (lqnorm < cgnorm) { + if (!goodb) { + x.setSubVector(0, this.xL); + } else { + final double step = bstep / beta1; + for (int i = 0; i < n; i++) { + final double bi = mb.getEntry(i); + final double xi = this.xL.getEntry(i); + x.setEntry(i, xi + step * bi); + } + } + } else { + final double anorm = FastMath.sqrt(tnorm); + final double diag = gbar == 0. ? anorm * MACH_PREC : gbar; + final double zbar = gammaZeta / diag; + final double step = (bstep + snprod * zbar) / beta1; + // ynorm = FastMath.sqrt(ynorm2 + zbar * zbar); + if (!goodb) { + for (int i = 0; i < n; i++) { + final double xi = this.xL.getEntry(i); + final double wi = wbar.getEntry(i); + x.setEntry(i, xi + zbar * wi); + } + } else { + for (int i = 0; i < n; i++) { + final double xi = this.xL.getEntry(i); + final double wi = wbar.getEntry(i); + final double bi = mb.getEntry(i); + x.setEntry(i, xi + zbar * wi + step * bi); + } + } + } + } + + /** + * Performs the initial phase of the SYMMLQ algorithm. On return, the value of the state + * variables of {@code this} object correspond to k = 1. + */ + void init() { + this.xL.set(0.); + /* + * Set up y for the first Lanczos vector. y and beta1 will be zero + * if b = 0. + */ + this.r1 = this.b.copy(); + this.y = this.m == null ? this.b.copy() : this.m.operate(this.r1); + if ((this.m != null) && this.check) { + checkSymmetry(this.m, this.r1, this.y, this.m.operate(this.y)); + } + + this.beta1 = this.r1.dotProduct(this.y); + if (this.beta1 < 0.) { + throwNPDLOException(this.m, this.y); + } + if (this.beta1 == 0.) { + /* If b = 0 exactly, stop with x = 0. */ + this.bIsNull = true; + return; + } + this.bIsNull = false; + this.beta1 = FastMath.sqrt(this.beta1); + /* At this point + * r1 = b, + * y = M * b, + * beta1 = beta[1]. + */ + final RealVector v = this.y.mapMultiply(1. / this.beta1); + this.y = this.a.operate(v); + if (this.check) { + checkSymmetry(this.a, v, this.y, this.a.operate(this.y)); + } + /* + * Set up y for the second Lanczos vector. y and beta will be zero + * or very small if b is an eigenvector. + */ + daxpy(-this.shift, v, this.y); + final double alpha = v.dotProduct(this.y); + daxpy(-alpha / this.beta1, this.r1, this.y); + /* + * At this point + * alpha = alpha[1] + * y = beta[2] * M^(-1) * P' * v[2] + */ + /* Make sure r2 will be orthogonal to the first v. */ + final double vty = v.dotProduct(this.y); + final double vtv = v.dotProduct(v); + daxpy(-vty / vtv, v, this.y); + this.r2 = this.y.copy(); + if (this.m != null) { + this.y = this.m.operate(this.r2); + } + this.oldb = this.beta1; + this.beta = this.r2.dotProduct(this.y); + if (this.beta < 0.) { + throwNPDLOException(this.m, this.y); + } + this.beta = FastMath.sqrt(this.beta); + /* + * At this point + * oldb = beta[1] + * beta = beta[2] + * y = beta[2] * P' * v[2] + * r2 = beta[2] * M^(-1) * P' * v[2] + */ + this.cgnorm = this.beta1; + this.gbar = alpha; + this.dbar = this.beta; + this.gammaZeta = this.beta1; + this.minusEpsZeta = 0.; + this.bstep = 0.; + this.snprod = 1.; + this.tnorm = alpha * alpha + this.beta * this.beta; + this.ynorm2 = 0.; + this.gmax = FastMath.abs(alpha) + MACH_PREC; + this.gmin = this.gmax; + + if (this.goodb) { + this.wbar = new ArrayRealVector(this.a.getRowDimension()); + this.wbar.set(0.); + } else { + this.wbar = v; + } + updateNorms(); + } + + /** + * Performs the next iteration of the algorithm. The iteration count should be incremented + * prior to calling this method. On return, the value of the state variables of {@code this} + * object correspond to the current iteration count {@code k}. + */ + void update() { + final RealVector v = y.mapMultiply(1. / beta); + y = a.operate(v); + daxpbypz(-shift, v, -beta / oldb, r1, y); + final double alpha = v.dotProduct(y); + /* + * At this point + * v = P' * v[k], + * y = (A - shift * I) * P' * v[k] - beta[k] * M^(-1) * P' * v[k-1], + * alpha = v'[k] * P * (A - shift * I) * P' * v[k] + * - beta[k] * v[k]' * P * M^(-1) * P' * v[k-1] + * = v'[k] * P * (A - shift * I) * P' * v[k] + * - beta[k] * v[k]' * v[k-1] + * = alpha[k]. + */ + daxpy(-alpha / beta, r2, y); + /* + * At this point + * y = (A - shift * I) * P' * v[k] - alpha[k] * M^(-1) * P' * v[k] + * - beta[k] * M^(-1) * P' * v[k-1] + * = M^(-1) * P' * (P * (A - shift * I) * P' * v[k] -alpha[k] * v[k] + * - beta[k] * v[k-1]) + * = beta[k+1] * M^(-1) * P' * v[k+1], + * from Paige and Saunders (1975), equation (3.2). + * + * WATCH-IT: the two following lines work only because y is no longer + * updated up to the end of the present iteration, and is + * reinitialized at the beginning of the next iteration. + */ + r1 = r2; + r2 = y; + if (m != null) { + y = m.operate(r2); + } + oldb = beta; + beta = r2.dotProduct(y); + if (beta < 0.) { + throwNPDLOException(m, y); + } + beta = FastMath.sqrt(beta); + /* + * At this point + * r1 = beta[k] * M^(-1) * P' * v[k], + * r2 = beta[k+1] * M^(-1) * P' * v[k+1], + * y = beta[k+1] * P' * v[k+1], + * oldb = beta[k], + * beta = beta[k+1]. + */ + tnorm += alpha * alpha + oldb * oldb + beta * beta; + /* + * Compute the next plane rotation for Q. See Paige and Saunders + * (1975), equation (5.6), with + * gamma = gamma[k-1], + * c = c[k-1], + * s = s[k-1]. + */ + final double gamma = FastMath.sqrt(gbar * gbar + oldb * oldb); + final double c = gbar / gamma; + final double s = oldb / gamma; + /* + * The relations + * gbar[k] = s[k-1] * (-c[k-2] * beta[k]) - c[k-1] * alpha[k] + * = s[k-1] * dbar[k] - c[k-1] * alpha[k], + * delta[k] = c[k-1] * dbar[k] + s[k-1] * alpha[k], + * are not stated in Paige and Saunders (1975), but can be retrieved + * by expanding the (k, k-1) and (k, k) coefficients of the matrix in + * equation (5.5). + */ + final double deltak = c * dbar + s * alpha; + gbar = s * dbar - c * alpha; + final double eps = s * beta; + dbar = -c * beta; + final double zeta = gammaZeta / gamma; + /* + * At this point + * gbar = gbar[k] + * deltak = delta[k] + * eps = eps[k+1] + * dbar = dbar[k+1] + * zeta = zeta[k-1] + */ + final double zetaC = zeta * c; + final double zetaS = zeta * s; + final int n = xL.getDimension(); + for (int i = 0; i < n; i++) { + final double xi = xL.getEntry(i); + final double vi = v.getEntry(i); + final double wi = wbar.getEntry(i); + xL.setEntry(i, xi + wi * zetaC + vi * zetaS); + wbar.setEntry(i, wi * s - vi * c); + } + /* + * At this point + * x = xL[k-1], + * ptwbar = P' wbar[k], + * see Paige and Saunders (1975), equations (5.9) and (5.10). + */ + bstep += snprod * c * zeta; + snprod *= s; + gmax = FastMath.max(gmax, gamma); + gmin = FastMath.min(gmin, gamma); + ynorm2 += zeta * zeta; + gammaZeta = minusEpsZeta - deltak * zeta; + minusEpsZeta = -eps * zeta; + /* + * At this point + * snprod = s[1] * ... * s[k-1], + * gmax = max(|alpha[1]|, gamma[1], ..., gamma[k-1]), + * gmin = min(|alpha[1]|, gamma[1], ..., gamma[k-1]), + * ynorm2 = zeta[1]^2 + ... + zeta[k-1]^2, + * gammaZeta = gamma[k] * zeta[k], + * minusEpsZeta = -eps[k+1] * zeta[k-1]. + * The relation for gammaZeta can be retrieved from Paige and + * Saunders (1975), equation (5.4a), last line of the vector + * gbar[k] * zbar[k] = -eps[k] * zeta[k-2] - delta[k] * zeta[k-1]. + */ + updateNorms(); + } + + /** + * Computes the norms of the residuals, and checks for convergence. Updates {@link #lqnorm} + * and {@link #cgnorm}. + */ + private void updateNorms() { + final double anorm = FastMath.sqrt(tnorm); + final double ynorm = FastMath.sqrt(ynorm2); + final double epsa = anorm * MACH_PREC; + final double epsx = anorm * ynorm * MACH_PREC; + final double epsr = anorm * ynorm * delta; + final double diag = gbar == 0. ? epsa : gbar; + lqnorm = FastMath.sqrt(gammaZeta * gammaZeta + minusEpsZeta * minusEpsZeta); + final double qrnorm = snprod * beta1; + cgnorm = qrnorm * beta / FastMath.abs(diag); + + /* + * Estimate cond(A). In this version we look at the diagonals of L + * in the factorization of the tridiagonal matrix, T = L * Q. + * Sometimes, T[k] can be misleadingly ill-conditioned when T[k+1] + * is not, so we must be careful not to overestimate acond. + */ + final double acond; + if (lqnorm <= cgnorm) { + acond = gmax / gmin; + } else { + acond = gmax / FastMath.min(gmin, FastMath.abs(diag)); + } + if (acond * MACH_PREC >= 0.1) { + throw new IllConditionedOperatorException(acond); + } + if (beta1 <= epsx) { + /* + * x has converged to an eigenvector of A corresponding to the + * eigenvalue shift. + */ + throw new SingularOperatorException(); + } + rnorm = FastMath.min(cgnorm, lqnorm); + hasConverged = (cgnorm <= epsx) || (cgnorm <= epsr); + } + + /** + * Returns {@code true} if the default stopping criterion is fulfilled. + * + * @return {@code true} if convergence of the iterations has occurred + */ + boolean hasConverged() { + return hasConverged; + } + + /** + * Returns {@code true} if the right-hand side vector is zero exactly. + * + * @return the boolean value of {@code b == 0} + */ + boolean bEqualsNullVector() { + return bIsNull; + } + + /** + * Returns {@code true} if {@code beta} is essentially zero. This method is used to check + * for early stop of the iterations. + * + * @return {@code true} if {@code beta < }{@link #MACH_PREC} + */ + boolean betaEqualsZero() { + return beta < MACH_PREC; + } + + /** + * Returns the norm of the updated, preconditioned residual. + * + * @return the norm of the residual, ||P * r|| + */ + double getNormOfResidual() { + return rnorm; + } + } + + /** Key for the exception context. */ + private static final String OPERATOR = "operator"; + + /** Key for the exception context. */ + private static final String THRESHOLD = "threshold"; + + /** Key for the exception context. */ + private static final String VECTOR = "vector"; + + /** Key for the exception context. */ + private static final String VECTOR1 = "vector1"; + + /** Key for the exception context. */ + private static final String VECTOR2 = "vector2"; + + /** {@code true} if symmetry of matrix and conditioner must be checked. */ + private final boolean check; + + /** The value of the custom tolerance δ for the default stopping criterion. */ + private final double delta; + + /** + * Creates a new instance of this class, with <a href="#stopcrit">default stopping + * criterion</a>. Note that setting {@code check} to {@code true} entails an extra matrix-vector + * product in the initial phase. + * + * @param maxIterations the maximum number of iterations + * @param delta the δ parameter for the default stopping criterion + * @param check {@code true} if self-adjointedness of both matrix and preconditioner should be + * checked + */ + public SymmLQ(final int maxIterations, final double delta, final boolean check) { + super(maxIterations); + this.delta = delta; + this.check = check; + } + + /** + * Creates a new instance of this class, with <a href="#stopcrit">default stopping criterion</a> + * and custom iteration manager. Note that setting {@code check} to {@code true} entails an + * extra matrix-vector product in the initial phase. + * + * @param manager the custom iteration manager + * @param delta the δ parameter for the default stopping criterion + * @param check {@code true} if self-adjointedness of both matrix and preconditioner should be + * checked + */ + public SymmLQ(final IterationManager manager, final double delta, final boolean check) { + super(manager); + this.delta = delta; + this.check = check; + } + + /** + * Returns {@code true} if symmetry of the matrix, and symmetry as well as positive definiteness + * of the preconditioner should be checked. + * + * @return {@code true} if the tests are to be performed + */ + public final boolean getCheck() { + return check; + } + + /** + * {@inheritDoc} + * + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * or {@code m} is not self-adjoint + * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive definite + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solve( + final RealLinearOperator a, final RealLinearOperator m, final RealVector b) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + MaxCountExceededException, + NonSelfAdjointOperatorException, + NonPositiveDefiniteOperatorException, + IllConditionedOperatorException { + MathUtils.checkNotNull(a); + final RealVector x = new ArrayRealVector(a.getColumnDimension()); + return solveInPlace(a, m, b, x, false, 0.); + } + + /** + * Returns an estimate of the solution to the linear system (A - shift · I) · x = + * b. + * + * <p>If the solution x is expected to contain a large multiple of {@code b} (as in + * Rayleigh-quotient iteration), then better precision may be achieved with {@code goodb} set to + * {@code true}; this however requires an extra call to the preconditioner. + * + * <p>{@code 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 b<sup>T</sup> + * · A · b / (b<sup>T</sup> · 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. + * + * @param a the linear operator A of the system + * @param m the preconditioner, M (can be {@code null}) + * @param b the right-hand side vector + * @param goodb usually {@code false}, except if {@code x} is expected to contain a large + * multiple of {@code b} + * @param shift the amount to be subtracted to all diagonal elements of A + * @return a reference to {@code x} (shallow copy) + * @throws NullArgumentException if one of the parameters is {@code null} + * @throws NonSquareOperatorException if {@code a} or {@code m} is not square + * @throws DimensionMismatchException if {@code m} or {@code b} have dimensions inconsistent + * with {@code a} + * @throws MaxCountExceededException at exhaustion of the iteration count, unless a custom + * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback} has + * been set at construction of the {@link IterationManager} + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * or {@code m} is not self-adjoint + * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive definite + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + public RealVector solve( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final boolean goodb, + final double shift) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + MaxCountExceededException, + NonSelfAdjointOperatorException, + NonPositiveDefiniteOperatorException, + IllConditionedOperatorException { + MathUtils.checkNotNull(a); + final RealVector x = new ArrayRealVector(a.getColumnDimension()); + return solveInPlace(a, m, b, x, goodb, shift); + } + + /** + * {@inheritDoc} + * + * @param x not meaningful in this implementation; should not be considered as an initial guess + * (<a href="#initguess">more</a>) + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * or {@code m} is not self-adjoint + * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive definite + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solve( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final RealVector x) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + NonPositiveDefiniteOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + MathUtils.checkNotNull(x); + return solveInPlace(a, m, b, x.copy(), false, 0.); + } + + /** + * {@inheritDoc} + * + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * is not self-adjoint + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solve(final RealLinearOperator a, final RealVector b) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + MathUtils.checkNotNull(a); + final RealVector x = new ArrayRealVector(a.getColumnDimension()); + x.set(0.); + return solveInPlace(a, null, b, x, false, 0.); + } + + /** + * Returns the solution to the system (A - shift · I) · x = b. + * + * <p>If the solution x is expected to contain a large multiple of {@code b} (as in + * Rayleigh-quotient iteration), then better precision may be achieved with {@code goodb} set to + * {@code true}. + * + * <p>{@code 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 b<sup>T</sup> + * · A · b / (b<sup>T</sup> · 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. + * + * @param a the linear operator A of the system + * @param b the right-hand side vector + * @param goodb usually {@code false}, except if {@code x} is expected to contain a large + * multiple of {@code b} + * @param shift the amount to be subtracted to all diagonal elements of A + * @return a reference to {@code x} + * @throws NullArgumentException if one of the parameters is {@code null} + * @throws NonSquareOperatorException if {@code a} is not square + * @throws DimensionMismatchException if {@code b} has dimensions inconsistent with {@code a} + * @throws MaxCountExceededException at exhaustion of the iteration count, unless a custom + * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback} has + * been set at construction of the {@link IterationManager} + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * is not self-adjoint + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + public RealVector solve( + final RealLinearOperator a, final RealVector b, final boolean goodb, final double shift) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + MathUtils.checkNotNull(a); + final RealVector x = new ArrayRealVector(a.getColumnDimension()); + return solveInPlace(a, null, b, x, goodb, shift); + } + + /** + * {@inheritDoc} + * + * @param x not meaningful in this implementation; should not be considered as an initial guess + * (<a href="#initguess">more</a>) + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * is not self-adjoint + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solve(final RealLinearOperator a, final RealVector b, final RealVector x) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + MathUtils.checkNotNull(x); + return solveInPlace(a, null, b, x.copy(), false, 0.); + } + + /** + * {@inheritDoc} + * + * @param x the vector to be updated with the solution; {@code x} should not be considered as an + * initial guess (<a href="#initguess">more</a>) + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * or {@code m} is not self-adjoint + * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive definite + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solveInPlace( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final RealVector x) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + NonPositiveDefiniteOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + return solveInPlace(a, m, b, x, false, 0.); + } + + /** + * Returns an estimate of the solution to the linear system (A - shift · I) · x = + * b. The solution is computed in-place. + * + * <p>If the solution x is expected to contain a large multiple of {@code b} (as in + * Rayleigh-quotient iteration), then better precision may be achieved with {@code goodb} set to + * {@code true}; this however requires an extra call to the preconditioner. + * + * <p>{@code 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 b<sup>T</sup> + * · A · b / (b<sup>T</sup> · 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. + * + * @param a the linear operator A of the system + * @param m the preconditioner, M (can be {@code null}) + * @param b the right-hand side vector + * @param x the vector to be updated with the solution; {@code x} should not be considered as an + * initial guess (<a href="#initguess">more</a>) + * @param goodb usually {@code false}, except if {@code x} is expected to contain a large + * multiple of {@code b} + * @param shift the amount to be subtracted to all diagonal elements of A + * @return a reference to {@code x} (shallow copy). + * @throws NullArgumentException if one of the parameters is {@code null} + * @throws NonSquareOperatorException if {@code a} or {@code m} is not square + * @throws DimensionMismatchException if {@code m}, {@code b} or {@code x} have dimensions + * inconsistent with {@code a}. + * @throws MaxCountExceededException at exhaustion of the iteration count, unless a custom + * {@link org.apache.commons.math3.util.Incrementor.MaxCountExceededCallback callback} has + * been set at construction of the {@link IterationManager} + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * or {@code m} is not self-adjoint + * @throws NonPositiveDefiniteOperatorException if {@code m} is not positive definite + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + public RealVector solveInPlace( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final RealVector x, + final boolean goodb, + final double shift) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + NonPositiveDefiniteOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + checkParameters(a, m, b, x); + + final IterationManager manager = getIterationManager(); + /* Initialization counts as an iteration. */ + manager.resetIterationCount(); + manager.incrementIterationCount(); + + final State state; + state = new State(a, m, b, goodb, shift, delta, check); + state.init(); + state.refineSolution(x); + IterativeLinearSolverEvent event; + event = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), x, b, state.getNormOfResidual()); + if (state.bEqualsNullVector()) { + /* If b = 0 exactly, stop with x = 0. */ + manager.fireTerminationEvent(event); + return x; + } + /* Cause termination if beta is essentially zero. */ + final boolean earlyStop; + earlyStop = state.betaEqualsZero() || state.hasConverged(); + manager.fireInitializationEvent(event); + if (!earlyStop) { + do { + manager.incrementIterationCount(); + event = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), x, b, state.getNormOfResidual()); + manager.fireIterationStartedEvent(event); + state.update(); + state.refineSolution(x); + event = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), x, b, state.getNormOfResidual()); + manager.fireIterationPerformedEvent(event); + } while (!state.hasConverged()); + } + event = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), x, b, state.getNormOfResidual()); + manager.fireTerminationEvent(event); + return x; + } + + /** + * {@inheritDoc} + * + * @param x the vector to be updated with the solution; {@code x} should not be considered as an + * initial guess (<a href="#initguess">more</a>) + * @throws NonSelfAdjointOperatorException if {@link #getCheck()} is {@code true}, and {@code a} + * is not self-adjoint + * @throws IllConditionedOperatorException if {@code a} is ill-conditioned + */ + @Override + public RealVector solveInPlace( + final RealLinearOperator a, final RealVector b, final RealVector x) + throws NullArgumentException, + NonSquareOperatorException, + DimensionMismatchException, + NonSelfAdjointOperatorException, + IllConditionedOperatorException, + MaxCountExceededException { + return solveInPlace(a, null, b, x, false, 0.); + } +} |