You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2014/02/19 21:31:47 UTC
svn commit: r1569905 - in /commons/proper/math/trunk/src: changes/changes.xml
main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java
Author: luc
Date: Wed Feb 19 20:31:47 2014
New Revision: 1569905
URL: http://svn.apache.org/r1569905
Log:
Make QR in GaussNewton faster and more accurate
Re-factored the code in GaussNewtonOptimizer so that the decomposition
algorithm sees the Jacobian and residuals instead of the normal
equation. This lets the QR algorithm operate directly on the Jacobian
matrix, which is faster and less sensitive to numerical errors. As a
result, one test case that threw a singular matrix exception now passes
with the QR decomposition.
The refactoring also include a speed improvement when computing the
normal matrix for the LU decomposition. Since the normal matrix is
symmetric only half of is computed, which results in a factor of 2 speed
up in computing the normal matrix for problems with many more
measurements than states.
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1569905&r1=1569904&r2=1569905&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Wed Feb 19 20:31:47 2014
@@ -51,6 +51,9 @@ If the output is not quite correct, chec
</properties>
<body>
<release version="3.3" date="TBD" description="TBD">
+ <action dev="luc" type="update" issue="MATH-1099" due-to="Evan Ward">
+ Make QR in GaussNewton faster and more accurate.
+ </action>
<action dev="luc" type="update" issue="MATH-870">
The sparse vector and matrix classes have been un-deprecated. This is a reversal
of a former decision, as we now think we should adopt a generally accepted
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java?rev=1569905&r1=1569904&r2=1569905&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java Wed Feb 19 20:31:47 2014
@@ -21,15 +21,15 @@ import org.apache.commons.math3.exceptio
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
import org.apache.commons.math3.linear.ArrayRealVector;
-import org.apache.commons.math3.linear.BlockRealMatrix;
-import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.LUDecomposition;
+import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.util.Incrementor;
+import org.apache.commons.math3.util.Pair;
/**
* Gauss-Newton least-squares solver. <p/> <p> This class solve a least-square problem by
@@ -46,28 +46,65 @@ public class GaussNewtonOptimizer implem
/** The decomposition algorithm to use to solve the normal equations. */
//TODO move to linear package and expand options?
public static enum Decomposition {
- /** Use {@link LUDecomposition}. */
+ /**
+ * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
+ * using the {@link LUDecomposition}.
+ *
+ * <p> Theoretically this method takes mn<sup>2></sup>/2 operations to compute the
+ * normal matrix and n<sup>3</sup>/3 operations (m > n) to solve the system using
+ * the LU decomposition. </p>
+ */
LU {
@Override
- protected DecompositionSolver getSolver(final RealMatrix matrix) {
- return new LUDecomposition(matrix, SINGULARITY_THRESHOLD).getSolver();
+ protected RealVector solve(final RealMatrix jacobian,
+ final RealVector residuals) {
+ try {
+ final Pair<RealMatrix, RealVector> normalEquation =
+ computeNormalMatrix(jacobian, residuals);
+ final RealMatrix normal = normalEquation.getFirst();
+ final RealVector jTr = normalEquation.getSecond();
+ return new LUDecomposition(normal, SINGULARITY_THRESHOLD)
+ .getSolver()
+ .solve(jTr);
+ } catch (SingularMatrixException e) {
+ throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
+ }
}
},
- /** Use {@link QRDecomposition}. */
+ /**
+ * Solve the linear least squares problem (Jx=r) using the {@link
+ * QRDecomposition}.
+ *
+ * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations
+ * (m > n) and has better numerical accuracy than any method that forms the normal
+ * equations. </p>
+ */
QR {
@Override
- protected DecompositionSolver getSolver(final RealMatrix matrix) {
- return new QRDecomposition(matrix, SINGULARITY_THRESHOLD).getSolver();
+ protected RealVector solve(final RealMatrix jacobian,
+ final RealVector residuals) {
+ try {
+ return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD)
+ .getSolver()
+ .solve(residuals);
+ } catch (SingularMatrixException e) {
+ throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
+ }
}
};
/**
- * Decompose the normal equations.
+ * Solve the linear least squares problem Jx=r.
*
- * @param matrix the normal matrix.
- * @return a solver.
+ * @param jacobian the Jacobian matrix, J. the number of rows >= the number or
+ * columns.
+ * @param residuals the computed residuals, r.
+ * @return the solution x, to the linear least squares problem Jx=r.
+ * @throws ConvergenceException if the matrix properties (e.g. singular) do not
+ * permit a solution.
*/
- protected abstract DecompositionSolver getSolver(RealMatrix matrix);
+ protected abstract RealVector solve(RealMatrix jacobian,
+ RealVector residuals);
}
/**
@@ -132,7 +169,6 @@ public class GaussNewtonOptimizer implem
throw new NullArgumentException();
}
- final int nR = lsp.getObservationSize(); // Number of observed data.
final int nC = lsp.getParameterSize();
final RealVector currentPoint = lsp.getStart();
@@ -160,41 +196,11 @@ public class GaussNewtonOptimizer implem
}
}
- // build the linear problem
- final double[] b = new double[nC];
- final double[][] a = new double[nC][nC];
- for (int i = 0; i < nR; ++i) {
-
- final double[] grad = weightedJacobian.getRow(i);
- final double residual = currentResiduals.getEntry(i);
-
- // compute the normal equation
- //residual is already weighted
- for (int j = 0; j < nC; ++j) {
- b[j] += residual * grad[j];
- }
-
- // build the contribution matrix for measurement i
- for (int k = 0; k < nC; ++k) {
- double[] ak = a[k];
- //Jacobian/gradient is already weighted
- for (int l = 0; l < nC; ++l) {
- ak[l] += grad[k] * grad[l];
- }
- }
- }
-
- try {
- // solve the linearized least squares problem
- RealMatrix mA = new BlockRealMatrix(a);
- DecompositionSolver solver = this.decomposition.getSolver(mA);
- final RealVector dX = solver.solve(new ArrayRealVector(b, false));
- // update the estimated parameters
- for (int i = 0; i < nC; ++i) {
- currentPoint.setEntry(i, currentPoint.getEntry(i) + dX.getEntry(i));
- }
- } catch (SingularMatrixException e) {
- throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
+ // solve the linearized least squares problem
+ final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals);
+ // update the estimated parameters
+ for (int i = 0; i < nC; ++i) {
+ currentPoint.setEntry(i, currentPoint.getEntry(i) + dX.getEntry(i));
}
}
}
@@ -206,4 +212,45 @@ public class GaussNewtonOptimizer implem
'}';
}
+ /**
+ * Compute the normal matrix, J<sup>T</sup>J.
+ *
+ * @param jacobian the m by n jacobian matrix, J. Input.
+ * @param residuals the m by 1 residual vector, r. Input.
+ * @return the n by n normal matrix and the n by 1 J<sup>Tr vector.
+ */
+ private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
+ final RealVector residuals) {
+ //since the normal matrix is symmetric, we only need to compute half of it.
+ final int nR = jacobian.getRowDimension();
+ final int nC = jacobian.getColumnDimension();
+ //allocate space for return values
+ final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
+ final RealVector jTr = new ArrayRealVector(nC);
+ //for each measurement
+ for (int i = 0; i < nR; ++i) {
+ //compute JTr for measurement i
+ for (int j = 0; j < nC; j++) {
+ jTr.setEntry(j, jTr.getEntry(j) +
+ residuals.getEntry(i) * jacobian.getEntry(i, j));
+ }
+
+ // add the the contribution to the normal matrix for measurement i
+ for (int k = 0; k < nC; ++k) {
+ //only compute the upper triangular part
+ for (int l = k; l < nC; ++l) {
+ normal.setEntry(k, l, normal.getEntry(k, l) +
+ jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
+ }
+ }
+ }
+ //copy the upper triangular part to the lower triangular part.
+ for (int i = 0; i < nC; i++) {
+ for (int j = 0; j < i; j++) {
+ normal.setEntry(i, j, normal.getEntry(j, i));
+ }
+ }
+ return new Pair<RealMatrix, RealVector>(normal, jTr);
+ }
+
}
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java?rev=1569905&r1=1569904&r2=1569905&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerWithQRTest.java Wed Feb 19 20:31:47 2014
@@ -49,20 +49,6 @@ public class GaussNewtonOptimizerWithQRT
@Override
@Test
- public void testMoreEstimatedParametersSimple() {
- /*
- * Exception is expected with this optimizer
- */
- try {
- super.testMoreEstimatedParametersSimple();
- fail(optimizer);
- } catch (ConvergenceException e) {
- //expected
- }
- }
-
- @Override
- @Test
public void testMoreEstimatedParametersUnsorted() {
/*
* Exception is expected with this optimizer