Author: dimpbx
Date: Mon Aug 16 08:45:10 2010
New Revision: 985828
URL: http://svn.apache.org/viewvc?rev=985828&view=rev
Log:
Code simplified in AbstractLeastSquaresOptimizer
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java?rev=985828&r1=985827&r2=985828&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
Mon Aug 16 08:45:10 2010
@@ -50,13 +50,13 @@ public abstract class AbstractLeastSquar
protected VectorialConvergenceChecker checker;
/**
- * Jacobian matrix.
+ * Jacobian matrix of the weighted residuals.
* <p>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).</p>
*/
- protected double[][] jacobian;
+ protected double[][] weightedResidualJacobian;
/** Number of columns of the jacobian matrix. */
protected int cols;
@@ -81,15 +81,9 @@ public abstract class AbstractLeastSquar
/** 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 AbstractLeastSquar
*/
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 AbstractLeastSquar
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 AbstractLeastSquar
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 AbstractLeastSquar
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;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java?rev=985828&r1=985827&r2=985828&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java
Mon Aug 16 08:45:10 2010
@@ -81,7 +81,7 @@ public class GaussNewtonOptimizer extend
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];
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java?rev=985828&r1=985827&r2=985828&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
Mon Aug 16 08:45:10 2010
@@ -270,7 +270,7 @@ public class LevenbergMarquardtOptimizer
VectorialPointValuePair current = new VectorialPointValuePair(point,
objective);
while (true) {
for (int i=0;i<rows;i++) {
- qtf[i]=wresiduals[i];
+ qtf[i]=weightedResiduals[i];
}
incrementIterationsCounter();
@@ -285,7 +285,7 @@ public class LevenbergMarquardtOptimizer
// so let jacobian contain the R matrix with its diagonal elements
for (int k = 0; k < solvedCols; ++k) {
int pk = permutation[k];
- wjacobian[k][pk] = diagR[pk];
+ weightedResidualJacobian[k][pk] = diagR[pk];
}
if (firstIteration) {
@@ -318,7 +318,7 @@ public class LevenbergMarquardtOptimizer
if (s != 0) {
double sum = 0;
for (int i = 0; i <= j; ++i) {
- sum += wjacobian[i][pj] * qtf[i];
+ sum += weightedResidualJacobian[i][pj] * qtf[i];
}
maxCosine = Math.max(maxCosine, Math.abs(sum) / (s *
cost));
}
@@ -345,8 +345,8 @@ public class LevenbergMarquardtOptimizer
oldX[pj] = point[pj];
}
double previousCost = cost;
- double[] tmpVec = residuals;
- residuals = oldRes;
+ double[] tmpVec = weightedResiduals;
+ weightedResiduals = oldRes;
oldRes = tmpVec;
tmpVec = objective;
objective = oldObj;
@@ -387,7 +387,7 @@ public class LevenbergMarquardtOptimizer
double dirJ = lmDir[pj];
work1[j] = 0;
for (int i = 0; i <= j; ++i) {
- work1[i] += wjacobian[i][pj] * dirJ;
+ work1[i] += weightedResidualJacobian[i][pj] * dirJ;
}
}
double coeff1 = 0;
@@ -443,8 +443,8 @@ public class LevenbergMarquardtOptimizer
int pj = permutation[j];
point[pj] = oldX[pj];
}
- tmpVec = residuals;
- residuals = oldRes;
+ tmpVec = weightedResiduals;
+ weightedResiduals = oldRes;
oldRes = tmpVec;
tmpVec = objective;
objective = oldObj;
@@ -514,7 +514,7 @@ public class LevenbergMarquardtOptimizer
int pk = permutation[k];
double ypk = lmDir[pk] / diagR[pk];
for (int i = 0; i < k; ++i) {
- lmDir[permutation[i]] -= ypk * wjacobian[i][pk];
+ lmDir[permutation[i]] -= ypk * weightedResidualJacobian[i][pk];
}
lmDir[pk] = ypk;
}
@@ -550,7 +550,7 @@ public class LevenbergMarquardtOptimizer
int pj = permutation[j];
double sum = 0;
for (int i = 0; i < j; ++i) {
- sum += wjacobian[i][pj] * work1[permutation[i]];
+ sum += weightedResidualJacobian[i][pj] *
work1[permutation[i]];
}
double s = (work1[pj] - sum) / diagR[pj];
work1[pj] = s;
@@ -565,7 +565,7 @@ public class LevenbergMarquardtOptimizer
int pj = permutation[j];
double sum = 0;
for (int i = 0; i <= j; ++i) {
- sum += wjacobian[i][pj] * qy[i];
+ sum += weightedResidualJacobian[i][pj] * qy[i];
}
sum /= diag[pj];
sum2 += sum * sum;
@@ -625,7 +625,7 @@ public class LevenbergMarquardtOptimizer
work1[pj] /= work2[j];
double tmp = work1[pj];
for (int i = j + 1; i < solvedCols; ++i) {
- work1[permutation[i]] -= wjacobian[i][pj] * tmp;
+ work1[permutation[i]] -= weightedResidualJacobian[i][pj] *
tmp;
}
}
sum2 = 0;
@@ -676,7 +676,7 @@ public class LevenbergMarquardtOptimizer
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
for (int i = j + 1; i < solvedCols; ++i) {
- wjacobian[i][pj] = wjacobian[j][permutation[i]];
+ weightedResidualJacobian[i][pj] =
weightedResidualJacobian[j][permutation[i]];
}
lmDir[j] = diagR[pj];
work[j] = qy[j];
@@ -707,7 +707,7 @@ public class LevenbergMarquardtOptimizer
final double sin;
final double cos;
- double rkk = wjacobian[k][pk];
+ double rkk = weightedResidualJacobian[k][pk];
if (Math.abs(rkk) < Math.abs(lmDiag[k])) {
final double cotan = rkk / lmDiag[k];
sin = 1.0 / Math.sqrt(1.0 + cotan * cotan);
@@ -720,17 +720,17 @@ public class LevenbergMarquardtOptimizer
// compute the modified diagonal element of R and
// the modified element of (Qty,0)
- wjacobian[k][pk] = cos * rkk + sin * lmDiag[k];
+ weightedResidualJacobian[k][pk] = cos * rkk + sin *
lmDiag[k];
final double temp = cos * work[k] + sin * qtbpj;
qtbpj = -sin * work[k] + cos * qtbpj;
work[k] = temp;
// accumulate the tranformation in the row of s
for (int i = k + 1; i < solvedCols; ++i) {
- double rik = wjacobian[i][pk];
+ double rik = weightedResidualJacobian[i][pk];
final double temp2 = cos * rik + sin * lmDiag[i];
lmDiag[i] = -sin * rik + cos * lmDiag[i];
- wjacobian[i][pk] = temp2;
+ weightedResidualJacobian[i][pk] = temp2;
}
}
@@ -738,8 +738,8 @@ public class LevenbergMarquardtOptimizer
// store the diagonal element of s and restore
// the corresponding diagonal element of R
- lmDiag[j] = wjacobian[j][permutation[j]];
- wjacobian[j][permutation[j]] = lmDir[j];
+ lmDiag[j] = weightedResidualJacobian[j][permutation[j]];
+ weightedResidualJacobian[j][permutation[j]] = lmDir[j];
}
@@ -759,7 +759,7 @@ public class LevenbergMarquardtOptimizer
int pj = permutation[j];
double sum = 0;
for (int i = j + 1; i < nSing; ++i) {
- sum += wjacobian[i][pj] * work[i];
+ sum += weightedResidualJacobian[i][pj] * work[i];
}
work[j] = (work[j] - sum) / lmDiag[j];
}
@@ -800,8 +800,8 @@ public class LevenbergMarquardtOptimizer
for (int k = 0; k < cols; ++k) {
permutation[k] = k;
double norm2 = 0;
- for (int i = 0; i < wjacobian.length; ++i) {
- double akk = wjacobian[i][k];
+ for (int i = 0; i < weightedResidualJacobian.length; ++i) {
+ double akk = weightedResidualJacobian[i][k];
norm2 += akk * akk;
}
jacNorm[k] = Math.sqrt(norm2);
@@ -815,8 +815,8 @@ public class LevenbergMarquardtOptimizer
double ak2 = Double.NEGATIVE_INFINITY;
for (int i = k; i < cols; ++i) {
double norm2 = 0;
- for (int j = k; j < wjacobian.length; ++j) {
- double aki = wjacobian[j][permutation[i]];
+ for (int j = k; j < weightedResidualJacobian.length; ++j) {
+ double aki = weightedResidualJacobian[j][permutation[i]];
norm2 += aki * aki;
}
if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
@@ -837,24 +837,24 @@ public class LevenbergMarquardtOptimizer
permutation[k] = pk;
// choose alpha such that Hk.u = alpha ek
- double akk = wjacobian[k][pk];
+ double akk = weightedResidualJacobian[k][pk];
double alpha = (akk > 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
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];
}
}
}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=985828&r1=985827&r2=985828&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Aug 16 08:45:10 2010
@@ -52,8 +52,14 @@ The <action> type attribute can be add,u
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="2.2" date="TBD" description="TBD">
+ <action dev="dimpbx" type="fix" issue="MATH-406">
+ Bug fixed in Levenberg-Marquardt (handling of weights).
+ </action>
+ <action dev="dimpbx" type="fix" issue="MATH-405">
+ Bug fixed in Levenberg-Marquardt (consistency of current).
+ </action>
<action dev="dimpbx" type="fix" issue="MATH-377">
- Fixed bug in chi-square computation in AbstractLeastSquaresOptimizer.
+ Bug fixed in chi-square computation in AbstractLeastSquaresOptimizer.
</action>
<action dev="luc" type="add" issue="MATH-400" due-to="J. Lewis Muir">
Added support for Gaussian curve fitting.
