Author: celestin
Date: Sat Jan 28 15:10:42 2012
New Revision: 1237069
URL: http://svn.apache.org/viewvc?rev=1237069&view=rev
Log:
Added method o.a.c.m.linear.IterativeLinearSolverEvent.getNormOfResidual()
(MATH-735).
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/ConjugateGradient.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/IterativeLinearSolverEvent.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/ConjugateGradientTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/ConjugateGradient.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/ConjugateGradient.java?rev=1237069&r1=1237068&r2=1237069&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/ConjugateGradient.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/ConjugateGradient.java
Sat Jan 28 15:10:42 2012
@@ -99,6 +99,9 @@ public class ConjugateGradient
/** The current estimate of the residual. */
private final RealVector r;
+ /** The current estimate of the norm of the residual. */
+ private final double rnorm;
+
/** The current estimate of the solution. */
private final RealVector x;
@@ -112,12 +115,22 @@ public class ConjugateGradient
* @param x the current estimate of the solution
* @param b the right-hand side vector
* @param r the current estimate of the residual
+ * @param rnorm the norm of the current estimate of the residual
*/
- public ConjugateGradientEvent(final Object source, final int
iterations, final RealVector x, final RealVector b, final RealVector r) {
+ public ConjugateGradientEvent(final Object source, final int
iterations,
+ final RealVector x, final RealVector b, final RealVector r,
+ final double rnorm) {
super(source, iterations);
this.x = RealVector.unmodifiableRealVector(x);
this.b = RealVector.unmodifiableRealVector(b);
this.r = RealVector.unmodifiableRealVector(r);
+ this.rnorm = rnorm;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNormOfResidual() {
+ return rnorm;
}
/** {@inheritDoc} */
@@ -206,7 +219,7 @@ public class ConjugateGradient
final IterationManager manager = getIterationManager();
// Initialization of default stopping criterion
manager.resetIterationCount();
- final double r2max = delta * delta * b.dotProduct(b);
+ final double rmax = delta * b.getNorm();
// Initialization phase counts as one iteration.
manager.incrementIterationCount();
@@ -218,7 +231,7 @@ public class ConjugateGradient
RealVector q = a.operate(p);
final RealVector r = b.combine(1, -1, q);
- double r2 = r.dotProduct(r);
+ double rnorm = r.getNorm();
RealVector z;
if (m == null) {
z = r;
@@ -226,16 +239,16 @@ public class ConjugateGradient
z = null;
}
IterativeLinearSolverEvent evt;
- evt = new ConjugateGradientEvent(this, manager.getIterations(), x, b,
r);
+ evt = new ConjugateGradientEvent(this, manager.getIterations(), x, b,
r, rnorm);
manager.fireInitializationEvent(evt);
- if (r2 <= r2max) {
+ if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
double rhoPrev = 0.;
while (true) {
manager.incrementIterationCount();
- evt = new ConjugateGradientEvent(this, manager.getIterations(), x,
b, r);
+ evt = new ConjugateGradientEvent(this, manager.getIterations(), x,
b, r, rnorm);
manager.fireIterationStartedEvent(evt);
if (m != null) {
z = m.solve(r);
@@ -268,10 +281,10 @@ public class ConjugateGradient
x.combineToSelf(1., alpha, p);
r.combineToSelf(1., -alpha, q);
rhoPrev = rhoNext;
- r2 = r.dotProduct(r);
- evt = new ConjugateGradientEvent(this, manager.getIterations(), x,
b, r);
+ rnorm = r.getNorm();
+ evt = new ConjugateGradientEvent(this, manager.getIterations(), x,
b, r, rnorm);
manager.fireIterationPerformedEvent(evt);
- if (r2 <= r2max) {
+ if (rnorm <= rmax) {
manager.fireTerminationEvent(evt);
return x;
}
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/IterativeLinearSolverEvent.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/IterativeLinearSolverEvent.java?rev=1237069&r1=1237068&r2=1237069&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/IterativeLinearSolverEvent.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/IterativeLinearSolverEvent.java
Sat Jan 28 15:10:42 2012
@@ -54,6 +54,23 @@ public abstract class IterativeLinearSol
public abstract RealVector getRightHandSideVector();
/**
+ * Returns the norm of the residual. The returned value is not required to
+ * be <em>exact</em>. Instead, the norm of the so-called <em>updated</em>
+ * residual (if available) should be returned. For example, the
+ * {@link ConjugateGradient conjugate gradient} method computes a sequence
+ * of residuals, the norm of which is cheap to compute. However, due to
+ * accumulation of round-off errors, this residual might differ from the
+ * true residual after some iterations. See e.g. A. Greenbaum and
+ * Z. Strakos, <em>Predicting the Behavior of Finite Precision Lanzos and
+ * Conjugate Gradient Computations</em>, Technical Report 538, Department
of
+ * Computer Science, New York University, 1991 (available
+ * <a
href="http://www.archive.org/details/predictingbehavi00gree";>here</a>).
+ *
+ * @return an estimate of the norm of the residual
+ */
+ public abstract double getNormOfResidual();
+
+ /**
* Returns the current estimate of the solution to the linear system to be
* solved. This method should return an unmodifiable view, or a deep copy
of
* the actual current solution, in order not to compromise subsequent
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java?rev=1237069&r1=1237068&r2=1237069&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SymmLQ.java
Sat Jan 28 15:10:42 2012
@@ -664,7 +664,7 @@ public class SymmLQ
*/
/** */
- private static final long serialVersionUID = 20120128L;
+ private static final long serialVersionUID = 2012012801L;
/** A reference to the state of this solver. */
private final State state;
@@ -681,6 +681,7 @@ public class SymmLQ
this.state = state;
}
+ /** {@inheritDoc} */
@Override
public int getIterations() {
return getIterationManager().getIterations();
@@ -688,6 +689,12 @@ public class SymmLQ
/** {@inheritDoc} */
@Override
+ public double getNormOfResidual() {
+ return FastMath.min(state.cgnorm, state.lqnorm);
+ }
+
+ /** {@inheritDoc} */
+ @Override
public RealVector getRightHandSideVector() {
return RealVector.unmodifiableRealVector(state.b);
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/ConjugateGradientTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/ConjugateGradientTest.java?rev=1237069&r1=1237068&r2=1237069&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/ConjugateGradientTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/ConjugateGradientTest.java
Sat Jan 28 15:10:42 2012
@@ -20,6 +20,7 @@ import java.util.Arrays;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.MaxCountExceededException;
+import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.util.IterationEvent;
import org.apache.commons.math.util.IterationListener;
import org.junit.Assert;
@@ -175,7 +176,7 @@ public class ConjugateGradientTest {
}
public void iterationPerformed(final IterationEvent e) {
- RealVector v = ((ProvidesResidual)e).getResidual();
+ RealVector v = ((ProvidesResidual) e).getResidual();
r.setSubVector(0, v);
v = ((IterativeLinearSolverEvent) e).getSolution();
x.setSubVector(0, v);
@@ -463,7 +464,6 @@ public class ConjugateGradientTest {
*/
final int[] count = new int[] {0, 0, 0, 0};
final IterationListener listener = new IterationListener() {
-
public void initializationPerformed(final IterationEvent e) {
++count[0];
}
@@ -498,4 +498,97 @@ public class ConjugateGradientTest {
Assert.assertEquals(msg, 1, count[3]);
}
}
+
+ @Test
+ public void testUnpreconditionedNormOfResidual() {
+ final int n = 5;
+ final int maxIterations = 100;
+ final RealLinearOperator a = new HilbertMatrix(n);
+ final IterativeLinearSolver solver;
+ final IterationListener listener = new IterationListener() {
+
+ private void doTestNormOfResidual(final IterationEvent e) {
+ final IterativeLinearSolverEvent evt;
+ evt = (IterativeLinearSolverEvent) e;
+ final RealVector x = evt.getSolution();
+ final RealVector b = evt.getRightHandSideVector();
+ final RealVector r = b.subtract(a.operate(x));
+ final double rnorm = r.getNorm();
+ Assert.assertEquals("iteration performed (residual)",
+ rnorm, evt.getNormOfResidual(),
+ FastMath.max(1E-5 * rnorm, 1E-10));
+ }
+
+ public void initializationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationStarted(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void terminationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+ };
+ solver = new ConjugateGradient(maxIterations, 1E-10, true);
+ solver.getIterationManager().addIterationListener(listener);
+ final RealVector b = new ArrayRealVector(n);
+ for (int j = 0; j < n; j++) {
+ b.set(0.);
+ b.setEntry(j, 1.);
+ solver.solve(a, b);
+ }
+ }
+
+ @Test
+ public void testPreconditionedNormOfResidual() {
+ final int n = 5;
+ final int maxIterations = 100;
+ final RealLinearOperator a = new HilbertMatrix(n);
+ final InvertibleRealLinearOperator m = JacobiPreconditioner.create(a);
+ final PreconditionedIterativeLinearSolver solver;
+ final IterationListener listener = new IterationListener() {
+
+ private void doTestNormOfResidual(final IterationEvent e) {
+ final IterativeLinearSolverEvent evt;
+ evt = (IterativeLinearSolverEvent) e;
+ final RealVector x = evt.getSolution();
+ final RealVector b = evt.getRightHandSideVector();
+ final RealVector r = b.subtract(a.operate(x));
+ final double rnorm = r.getNorm();
+ Assert.assertEquals("iteration performed (residual)",
+ rnorm, evt.getNormOfResidual(),
+ FastMath.max(1E-5 * rnorm, 1E-10));
+ }
+
+ public void initializationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationStarted(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void terminationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+ };
+ solver = new ConjugateGradient(maxIterations, 1E-10, true);
+ solver.getIterationManager().addIterationListener(listener);
+ final RealVector b = new ArrayRealVector(n);
+ for (int j = 0; j < n; j++) {
+ b.set(0.);
+ b.setEntry(j, 1.);
+ solver.solve(a, m, b);
+ }
+ }
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java?rev=1237069&r1=1237068&r2=1237069&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SymmLQTest.java
Sat Jan 28 15:10:42 2012
@@ -624,4 +624,98 @@ public class SymmLQTest {
});
new SymmLQ(100, 1., true).solve(a, m, b);
}
+
+ @Test
+ public void testUnpreconditionedNormOfResidual() {
+ final int n = 5;
+ final int maxIterations = 100;
+ final RealLinearOperator a = new HilbertMatrix(n);
+ final IterativeLinearSolver solver;
+ final IterationListener listener = new IterationListener() {
+
+ private void doTestNormOfResidual(final IterationEvent e) {
+ final IterativeLinearSolverEvent evt;
+ evt = (IterativeLinearSolverEvent) e;
+ final RealVector x = evt.getSolution();
+ final RealVector b = evt.getRightHandSideVector();
+ final RealVector r = b.subtract(a.operate(x));
+ final double rnorm = r.getNorm();
+ Assert.assertEquals("iteration performed (residual)",
+ rnorm, evt.getNormOfResidual(),
+ FastMath.max(1E-5 * rnorm, 1E-10));
+ }
+
+ public void initializationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationStarted(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void terminationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+ };
+ solver = new ConjugateGradient(maxIterations, 1E-10, true);
+ solver.getIterationManager().addIterationListener(listener);
+ final RealVector b = new ArrayRealVector(n);
+ for (int j = 0; j < n; j++) {
+ b.set(0.);
+ b.setEntry(j, 1.);
+ solver.solve(a, b);
+ }
+ }
+
+ @Test
+ public void testPreconditionedNormOfResidual() {
+ final int n = 5;
+ final int maxIterations = 100;
+ final RealLinearOperator a = new HilbertMatrix(n);
+ final InvertibleRealLinearOperator m = JacobiPreconditioner.create(a);
+ final PreconditionedIterativeLinearSolver solver;
+ final IterationListener listener = new IterationListener() {
+
+ private void doTestNormOfResidual(final IterationEvent e) {
+ final IterativeLinearSolverEvent evt;
+ evt = (IterativeLinearSolverEvent) e;
+ final RealVector x = evt.getSolution();
+ final RealVector b = evt.getRightHandSideVector();
+ final RealVector r = b.subtract(a.operate(x));
+ final double rnorm = r.getNorm();
+ Assert.assertEquals("iteration performed (residual)",
+ rnorm, evt.getNormOfResidual(),
+ FastMath.max(1E-5 * rnorm, 1E-10));
+ }
+
+ public void initializationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void iterationStarted(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+
+ public void terminationPerformed(final IterationEvent e) {
+ doTestNormOfResidual(e);
+ }
+ };
+ solver = new ConjugateGradient(maxIterations, 1E-10, true);
+ solver.getIterationManager().addIterationListener(listener);
+ final RealVector b = new ArrayRealVector(n);
+ for (int j = 0; j < n; j++) {
+ b.set(0.);
+ b.setEntry(j, 1.);
+ solver.solve(a, m, b);
+ }
+ }
}
+
