diff --git a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java index 007e13f9d..4f12f6837 100644 --- a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java @@ -50,13 +50,13 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul protected VectorialConvergenceChecker checker; /** - * Jacobian matrix. + * Jacobian matrix of the weighted residuals. *

This matrix is in canonical form just after the calls to * {@link #updateJacobian()}, but may be modified by the solver * in the derived class (the {@link LevenbergMarquardtOptimizer * Levenberg-Marquardt optimizer} does this).

*/ - protected double[][] jacobian; + protected double[][] weightedResidualJacobian; /** Number of columns of the jacobian matrix. */ protected int cols; @@ -81,15 +81,9 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul /** Current objective function value. */ protected double[] objective; - - /** Current residuals. */ - protected double[] residuals; - - /** Weighted Jacobian */ - protected double[][] wjacobian; /** Weighted residuals */ - protected double[] wresiduals; + protected double[] weightedResiduals; /** Cost value (square root of the sum of the residuals). */ protected double cost; @@ -188,17 +182,17 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul */ protected void updateJacobian() throws FunctionEvaluationException { ++jacobianEvaluations; - jacobian = jF.value(point); - if (jacobian.length != rows) { + weightedResidualJacobian = jF.value(point); + if (weightedResidualJacobian.length != rows) { throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, - jacobian.length, rows); + weightedResidualJacobian.length, rows); } for (int i = 0; i < rows; i++) { - final double[] ji = jacobian[i]; + final double[] ji = weightedResidualJacobian[i]; double wi = Math.sqrt(residualsWeights[i]); for (int j = 0; j < cols; ++j) { - ji[j] *= -1.0; - wjacobian[i][j] = ji[j]*wi; + //ji[j] *= -1.0; + weightedResidualJacobian[i][j] = -ji[j]*wi; } } } @@ -225,8 +219,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul int index = 0; for (int i = 0; i < rows; i++) { final double residual = targetValues[i] - objective[i]; - residuals[i] = residual; - wresiduals[i]= residual*Math.sqrt(residualsWeights[i]); + weightedResiduals[i]= residual*Math.sqrt(residualsWeights[i]); cost += residualsWeights[i] * residual * residual; index += cols; } @@ -278,7 +271,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul for (int j = i; j < cols; ++j) { double sum = 0; for (int k = 0; k < rows; ++k) { - sum += wjacobian[k][i] * wjacobian[k][j]; + sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j]; } jTj[i][j] = sum; jTj[j][i] = sum; @@ -343,15 +336,13 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul targetValues = target.clone(); residualsWeights = weights.clone(); this.point = startPoint.clone(); - this.residuals = new double[target.length]; // arrays shared with the other private methods rows = target.length; cols = point.length; - jacobian = new double[rows][cols]; - wjacobian = new double[rows][cols]; - wresiduals = new double[rows]; + weightedResidualJacobian = new double[rows][cols]; + this.weightedResiduals = new double[rows]; cost = Double.POSITIVE_INFINITY; diff --git a/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java index ee8ef904a..43f7d93d8 100644 --- a/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java @@ -81,7 +81,7 @@ public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer { final double[][] a = new double[cols][cols]; for (int i = 0; i < rows; ++i) { - final double[] grad = jacobian[i]; + final double[] grad = weightedResidualJacobian[i]; final double weight = residualsWeights[i]; final double residual = objective[i] - targetValues[i]; diff --git a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java index 286f82c7d..6067de098 100644 --- a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java +++ b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java @@ -270,7 +270,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer { VectorialPointValuePair current = new VectorialPointValuePair(point, objective); while (true) { for (int i=0;i 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2); double betak = 1.0 / (ak2 - akk * alpha); beta[pk] = betak; // transform the current column diagR[pk] = alpha; - wjacobian[k][pk] -= alpha; + weightedResidualJacobian[k][pk] -= alpha; // transform the remaining columns for (int dk = cols - 1 - k; dk > 0; --dk) { double gamma = 0; - for (int j = k; j < wjacobian.length; ++j) { - gamma += wjacobian[j][pk] * wjacobian[j][permutation[k + dk]]; + for (int j = k; j < weightedResidualJacobian.length; ++j) { + gamma += weightedResidualJacobian[j][pk] * weightedResidualJacobian[j][permutation[k + dk]]; } gamma *= betak; - for (int j = k; j < wjacobian.length; ++j) { - wjacobian[j][permutation[k + dk]] -= gamma * wjacobian[j][pk]; + for (int j = k; j < weightedResidualJacobian.length; ++j) { + weightedResidualJacobian[j][permutation[k + dk]] -= gamma * weightedResidualJacobian[j][pk]; } } @@ -874,11 +874,11 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer { int pk = permutation[k]; double gamma = 0; for (int i = k; i < rows; ++i) { - gamma += wjacobian[i][pk] * y[i]; + gamma += weightedResidualJacobian[i][pk] * y[i]; } gamma *= beta[pk]; for (int i = k; i < rows; ++i) { - y[i] -= gamma * wjacobian[i][pk]; + y[i] -= gamma * weightedResidualJacobian[i][pk]; } } } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 37a6616a7..555c599e4 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,8 +52,14 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Bug fixed in Levenberg-Marquardt (handling of weights). + + + Bug fixed in Levenberg-Marquardt (consistency of current). + - Fixed bug in chi-square computation in AbstractLeastSquaresOptimizer. + Bug fixed in chi-square computation in AbstractLeastSquaresOptimizer. Added support for Gaussian curve fitting.