MATH-405 corrected

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@984404 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Dimitri Pourbaix 2010-08-11 13:40:39 +00:00
parent 31fef8fab1
commit 784e4f69ec
3 changed files with 39 additions and 34 deletions

View File

@ -247,12 +247,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
* @return chi-square value * @return chi-square value
*/ */
public double getChiSquare() { public double getChiSquare() {
double chiSquare = 0; return cost*cost;
for (int i = 0; i < rows; ++i) {
final double residual = residuals[i];
chiSquare += residual * residual * residualsWeights[i];
}
return chiSquare;
} }
/** /**

View File

@ -255,6 +255,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
double[] diag = new double[cols]; double[] diag = new double[cols];
double[] oldX = new double[cols]; double[] oldX = new double[cols];
double[] oldRes = new double[rows]; double[] oldRes = new double[rows];
double[] oldObj = new double[rows];
double[] qtf = new double[rows];
double[] work1 = new double[cols]; double[] work1 = new double[cols];
double[] work2 = new double[cols]; double[] work2 = new double[cols];
double[] work3 = new double[cols]; double[] work3 = new double[cols];
@ -267,7 +269,9 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
boolean firstIteration = true; boolean firstIteration = true;
VectorialPointValuePair current = new VectorialPointValuePair(point, objective); VectorialPointValuePair current = new VectorialPointValuePair(point, objective);
while (true) { while (true) {
for (int i=0;i<rows;i++) {
qtf[i]=residuals[i];
}
incrementIterationsCounter(); incrementIterationsCounter();
// compute the Q.R. decomposition of the jacobian matrix // compute the Q.R. decomposition of the jacobian matrix
@ -276,8 +280,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
qrDecomposition(); qrDecomposition();
// compute Qt.res // compute Qt.res
qTy(residuals); qTy(qtf);
// now we don't need Q anymore, // now we don't need Q anymore,
// so let jacobian contain the R matrix with its diagonal elements // so let jacobian contain the R matrix with its diagonal elements
for (int k = 0; k < solvedCols; ++k) { for (int k = 0; k < solvedCols; ++k) {
@ -315,7 +318,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
if (s != 0) { if (s != 0) {
double sum = 0; double sum = 0;
for (int i = 0; i <= j; ++i) { for (int i = 0; i <= j; ++i) {
sum += jacobian[i][pj] * residuals[i]; sum += jacobian[i][pj] * qtf[i];
} }
maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost)); maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
} }
@ -323,6 +326,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
} }
if (maxCosine <= orthoTolerance) { if (maxCosine <= orthoTolerance) {
// convergence has been reached // convergence has been reached
updateResidualsAndCost();
current = new VectorialPointValuePair(point, objective);
return current; return current;
} }
@ -343,9 +348,12 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
double[] tmpVec = residuals; double[] tmpVec = residuals;
residuals = oldRes; residuals = oldRes;
oldRes = tmpVec; oldRes = tmpVec;
tmpVec = objective;
objective = oldObj;
oldObj = tmpVec;
// determine the Levenberg-Marquardt parameter // determine the Levenberg-Marquardt parameter
determineLMParameter(oldRes, delta, diag, work1, work2, work3); determineLMParameter(qtf, delta, diag, work1, work2, work3);
// compute the new point and the norm of the evolution direction // compute the new point and the norm of the evolution direction
double lmNorm = 0; double lmNorm = 0;
@ -357,7 +365,6 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
lmNorm += s * s; lmNorm += s * s;
} }
lmNorm = Math.sqrt(lmNorm); lmNorm = Math.sqrt(lmNorm);
// on the first iteration, adjust the initial step bound. // on the first iteration, adjust the initial step bound.
if (firstIteration) { if (firstIteration) {
delta = Math.min(delta, lmNorm); delta = Math.min(delta, lmNorm);
@ -365,7 +372,6 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// evaluate the function at x + p and calculate its norm // evaluate the function at x + p and calculate its norm
updateResidualsAndCost(); updateResidualsAndCost();
current = new VectorialPointValuePair(point, objective);
// compute the scaled actual reduction // compute the scaled actual reduction
double actRed = -1.0; double actRed = -1.0;
@ -421,6 +427,15 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
xNorm += xK * xK; xNorm += xK * xK;
} }
xNorm = Math.sqrt(xNorm); xNorm = Math.sqrt(xNorm);
current = new VectorialPointValuePair(point, objective);
// tests for convergence.
if (checker != null) {
// we use the vectorial convergence checker
if (checker.converged(getIterations(), previous, current)) {
return current;
}
}
} else { } else {
// failed iteration, reset the previous values // failed iteration, reset the previous values
cost = previousCost; cost = previousCost;
@ -431,24 +446,18 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
tmpVec = residuals; tmpVec = residuals;
residuals = oldRes; residuals = oldRes;
oldRes = tmpVec; oldRes = tmpVec;
tmpVec = objective;
objective = oldObj;
oldObj = tmpVec;
} }
if (checker==null) {
// tests for convergence. if (((Math.abs(actRed) <= costRelativeTolerance) &&
if (checker != null) { (preRed <= costRelativeTolerance) &&
// we use the vectorial convergence checker (ratio <= 2.0)) ||
if (checker.converged(getIterations(), previous, current)) { (delta <= parRelativeTolerance * xNorm)) {
return current; return current;
} }
} else {
// we use the Levenberg-Marquardt specific convergence parameters
if (((Math.abs(actRed) <= costRelativeTolerance) &&
(preRed <= costRelativeTolerance) &&
(ratio <= 2.0)) ||
(delta <= parRelativeTolerance * xNorm)) {
return current;
}
} }
// tests for termination and stringent tolerances // tests for termination and stringent tolerances
// (2.2204e-16 is the machine epsilon for IEEE754) // (2.2204e-16 is the machine epsilon for IEEE754)
if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) { if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {

View File

@ -152,14 +152,14 @@ public class MinpackTest extends TestCase {
minpackTest(new FreudensteinRothFunction(new double[] { 5.0, -20.0 }, minpackTest(new FreudensteinRothFunction(new double[] { 5.0, -20.0 },
12432.833948863, 6.9988751744895, 12432.833948863, 6.9988751744895,
new double[] { new double[] {
11.4121122022341, 11.41300466147456,
-0.8968550851268697 -0.896796038685959
}), false); }), false);
minpackTest(new FreudensteinRothFunction(new double[] { 50.0, -200.0 }, minpackTest(new FreudensteinRothFunction(new double[] { 50.0, -200.0 },
11426454.595762, 6.99887517242903, 11426454.595762, 6.99887517242903,
new double[] { new double[] {
11.412069435091231, 11.412781785788564,
-0.8968582807605691 -0.8968051074920405
}), false); }), false);
} }
@ -325,7 +325,8 @@ public class MinpackTest extends TestCase {
minpackTest(new JennrichSampsonFunction(10, new double[] { 0.3, 0.4 }, minpackTest(new JennrichSampsonFunction(10, new double[] { 0.3, 0.4 },
64.5856498144943, 11.1517793413499, 64.5856498144943, 11.1517793413499,
new double[] { new double[] {
0.2578330049, 0.257829976764542 // 0.2578330049, 0.257829976764542
0.2578199266368004, 0.25782997676455244
}), false); }), false);
} }