Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-406 corrected 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@985589 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Dimitri Pourbaix 2010-08-14 21:59:07 +00:00
parent 6b4085c327
commit 19f58cfa37
3 changed files with 44 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
View File

@ -84,6 +84,12 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
/** Current residuals. */ /** Current residuals. */
protected double[] residuals; protected double[] residuals;
/** Weighted Jacobian */
protected double[][] wjacobian;
/** Weighted residuals */
protected double[] wresiduals;
/** Cost value (square root of the sum of the residuals). */ /** Cost value (square root of the sum of the residuals). */
protected double cost; protected double cost;
@ -189,9 +195,10 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
} }
for (int i = 0; i < rows; i++) { for (int i = 0; i < rows; i++) {
final double[] ji = jacobian[i]; final double[] ji = jacobian[i];
final double factor = -Math.sqrt(residualsWeights[i]); double wi = Math.sqrt(residualsWeights[i]);
for (int j = 0; j < cols; ++j) { for (int j = 0; j < cols; ++j) {
ji[j] *= factor; ji[j] *= -1.0;
wjacobian[i][j] = ji[j]*wi;
} }
} }
} }
@ -219,6 +226,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
for (int i = 0; i < rows; i++) { for (int i = 0; i < rows; i++) {
final double residual = targetValues[i] - objective[i]; final double residual = targetValues[i] - objective[i];
residuals[i] = residual; residuals[i] = residual;
wresiduals[i]= residual*Math.sqrt(residualsWeights[i]);
cost += residualsWeights[i] * residual * residual; cost += residualsWeights[i] * residual * residual;
index += cols; index += cols;
} }
@ -270,7 +278,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
for (int j = i; j < cols; ++j) { for (int j = i; j < cols; ++j) {
double sum = 0; double sum = 0;
for (int k = 0; k < rows; ++k) { for (int k = 0; k < rows; ++k) {
sum += jacobian[k][i] * jacobian[k][j]; sum += wjacobian[k][i] * wjacobian[k][j];
} }
jTj[i][j] = sum; jTj[i][j] = sum;
jTj[j][i] = sum; jTj[j][i] = sum;
@ -342,6 +350,9 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
cols = point.length; cols = point.length;
jacobian = new double[rows][cols]; jacobian = new double[rows][cols];
wjacobian = new double[rows][cols];
wresiduals = new double[rows];
cost = Double.POSITIVE_INFINITY; cost = Double.POSITIVE_INFINITY;
return doOptimize(); return doOptimize();

56
src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
View File

@ -270,7 +270,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
VectorialPointValuePair current = new VectorialPointValuePair(point, objective); VectorialPointValuePair current = new VectorialPointValuePair(point, objective);
while (true) { while (true) {
for (int i=0;i<rows;i++) { for (int i=0;i<rows;i++) {
qtf[i]=residuals[i]; qtf[i]=wresiduals[i];
} }
incrementIterationsCounter(); incrementIterationsCounter();
@ -285,7 +285,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// 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) {
int pk = permutation[k]; int pk = permutation[k];
jacobian[k][pk] = diagR[pk]; wjacobian[k][pk] = diagR[pk];
} }
if (firstIteration) { if (firstIteration) {
@ -318,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] * qtf[i]; sum += wjacobian[i][pj] * qtf[i];
} }
maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost)); maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
} }
@ -387,7 +387,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
double dirJ = lmDir[pj]; double dirJ = lmDir[pj];
work1[j] = 0; work1[j] = 0;
for (int i = 0; i <= j; ++i) { for (int i = 0; i <= j; ++i) {
work1[i] += jacobian[i][pj] * dirJ; work1[i] += wjacobian[i][pj] * dirJ;
} }
} }
double coeff1 = 0; double coeff1 = 0;
@ -514,7 +514,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
int pk = permutation[k]; int pk = permutation[k];
double ypk = lmDir[pk] / diagR[pk]; double ypk = lmDir[pk] / diagR[pk];
for (int i = 0; i < k; ++i) { for (int i = 0; i < k; ++i) {
lmDir[permutation[i]] -= ypk * jacobian[i][pk]; lmDir[permutation[i]] -= ypk * wjacobian[i][pk];
} }
lmDir[pk] = ypk; lmDir[pk] = ypk;
} }
@ -550,7 +550,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
int pj = permutation[j]; int pj = permutation[j];
double sum = 0; double sum = 0;
for (int i = 0; i < j; ++i) { for (int i = 0; i < j; ++i) {
sum += jacobian[i][pj] * work1[permutation[i]]; sum += wjacobian[i][pj] * work1[permutation[i]];
} }
double s = (work1[pj] - sum) / diagR[pj]; double s = (work1[pj] - sum) / diagR[pj];
work1[pj] = s; work1[pj] = s;
@ -565,7 +565,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
int pj = permutation[j]; int pj = permutation[j];
double sum = 0; double sum = 0;
for (int i = 0; i <= j; ++i) { for (int i = 0; i <= j; ++i) {
sum += jacobian[i][pj] * qy[i]; sum += wjacobian[i][pj] * qy[i];
} }
sum /= diag[pj]; sum /= diag[pj];
sum2 += sum * sum; sum2 += sum * sum;
@ -625,7 +625,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
work1[pj] /= work2[j]; work1[pj] /= work2[j];
double tmp = work1[pj]; double tmp = work1[pj];
for (int i = j + 1; i < solvedCols; ++i) { for (int i = j + 1; i < solvedCols; ++i) {
work1[permutation[i]] -= jacobian[i][pj] * tmp; work1[permutation[i]] -= wjacobian[i][pj] * tmp;
} }
} }
sum2 = 0; sum2 = 0;
@ -676,7 +676,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
for (int j = 0; j < solvedCols; ++j) { for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j]; int pj = permutation[j];
for (int i = j + 1; i < solvedCols; ++i) { for (int i = j + 1; i < solvedCols; ++i) {
jacobian[i][pj] = jacobian[j][permutation[i]]; wjacobian[i][pj] = wjacobian[j][permutation[i]];
} }
lmDir[j] = diagR[pj]; lmDir[j] = diagR[pj];
work[j] = qy[j]; work[j] = qy[j];
@ -707,7 +707,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
final double sin; final double sin;
final double cos; final double cos;
double rkk = jacobian[k][pk]; double rkk = wjacobian[k][pk];
if (Math.abs(rkk) < Math.abs(lmDiag[k])) { if (Math.abs(rkk) < Math.abs(lmDiag[k])) {
final double cotan = rkk / lmDiag[k]; final double cotan = rkk / lmDiag[k];
sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); sin = 1.0 / Math.sqrt(1.0 + cotan * cotan);
@ -720,17 +720,17 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// compute the modified diagonal element of R and // compute the modified diagonal element of R and
// the modified element of (Qty,0) // the modified element of (Qty,0)
jacobian[k][pk] = cos * rkk + sin * lmDiag[k]; wjacobian[k][pk] = cos * rkk + sin * lmDiag[k];
final double temp = cos * work[k] + sin * qtbpj; final double temp = cos * work[k] + sin * qtbpj;
qtbpj = -sin * work[k] + cos * qtbpj; qtbpj = -sin * work[k] + cos * qtbpj;
work[k] = temp; work[k] = temp;
// accumulate the tranformation in the row of s // accumulate the tranformation in the row of s
for (int i = k + 1; i < solvedCols; ++i) { for (int i = k + 1; i < solvedCols; ++i) {
double rik = jacobian[i][pk]; double rik = wjacobian[i][pk];
final double temp2 = cos * rik + sin * lmDiag[i]; final double temp2 = cos * rik + sin * lmDiag[i];
lmDiag[i] = -sin * rik + cos * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i];
jacobian[i][pk] = temp2; wjacobian[i][pk] = temp2;
} }
} }
@ -738,8 +738,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
// store the diagonal element of s and restore // store the diagonal element of s and restore
// the corresponding diagonal element of R // the corresponding diagonal element of R
lmDiag[j] = jacobian[j][permutation[j]]; lmDiag[j] = wjacobian[j][permutation[j]];
jacobian[j][permutation[j]] = lmDir[j]; wjacobian[j][permutation[j]] = lmDir[j];
} }
@ -759,7 +759,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
int pj = permutation[j]; int pj = permutation[j];
double sum = 0; double sum = 0;
for (int i = j + 1; i < nSing; ++i) { for (int i = j + 1; i < nSing; ++i) {
sum += jacobian[i][pj] * work[i]; sum += wjacobian[i][pj] * work[i];
} }
work[j] = (work[j] - sum) / lmDiag[j]; work[j] = (work[j] - sum) / lmDiag[j];
} }
@ -800,8 +800,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
for (int k = 0; k < cols; ++k) { for (int k = 0; k < cols; ++k) {
permutation[k] = k; permutation[k] = k;
double norm2 = 0; double norm2 = 0;
for (int i = 0; i < jacobian.length; ++i) { for (int i = 0; i < wjacobian.length; ++i) {
double akk = jacobian[i][k]; double akk = wjacobian[i][k];
norm2 += akk * akk; norm2 += akk * akk;
} }
jacNorm[k] = Math.sqrt(norm2); jacNorm[k] = Math.sqrt(norm2);
@ -815,8 +815,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
double ak2 = Double.NEGATIVE_INFINITY; double ak2 = Double.NEGATIVE_INFINITY;
for (int i = k; i < cols; ++i) { for (int i = k; i < cols; ++i) {
double norm2 = 0; double norm2 = 0;
for (int j = k; j < jacobian.length; ++j) { for (int j = k; j < wjacobian.length; ++j) {
double aki = jacobian[j][permutation[i]]; double aki = wjacobian[j][permutation[i]];
norm2 += aki * aki; norm2 += aki * aki;
} }
if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
@ -837,24 +837,24 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
permutation[k] = pk; permutation[k] = pk;
// choose alpha such that Hk.u = alpha ek // choose alpha such that Hk.u = alpha ek
double akk = jacobian[k][pk]; double akk = wjacobian[k][pk];
double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2); double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2);
double betak = 1.0 / (ak2 - akk * alpha); double betak = 1.0 / (ak2 - akk * alpha);
beta[pk] = betak; beta[pk] = betak;
// transform the current column // transform the current column
diagR[pk] = alpha; diagR[pk] = alpha;
jacobian[k][pk] -= alpha; wjacobian[k][pk] -= alpha;
// transform the remaining columns // transform the remaining columns
for (int dk = cols - 1 - k; dk > 0; --dk) { for (int dk = cols - 1 - k; dk > 0; --dk) {
double gamma = 0; double gamma = 0;
for (int j = k; j < jacobian.length; ++j) { for (int j = k; j < wjacobian.length; ++j) {
gamma += jacobian[j][pk] * jacobian[j][permutation[k + dk]]; gamma += wjacobian[j][pk] * wjacobian[j][permutation[k + dk]];
} }
gamma *= betak; gamma *= betak;
for (int j = k; j < jacobian.length; ++j) { for (int j = k; j < wjacobian.length; ++j) {
jacobian[j][permutation[k + dk]] -= gamma * jacobian[j][pk]; wjacobian[j][permutation[k + dk]] -= gamma * wjacobian[j][pk];
} }
} }
@ -874,11 +874,11 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
int pk = permutation[k]; int pk = permutation[k];
double gamma = 0; double gamma = 0;
for (int i = k; i < rows; ++i) { for (int i = k; i < rows; ++i) {
gamma += jacobian[i][pk] * y[i]; gamma += wjacobian[i][pk] * y[i];
} }
gamma *= beta[pk]; gamma *= beta[pk];
for (int i = k; i < rows; ++i) { for (int i = k; i < rows; ++i) {
y[i] -= gamma * jacobian[i][pk]; y[i] -= gamma * wjacobian[i][pk];
} }
} }
} }

4
src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
View File

@ -472,8 +472,8 @@ extends TestCase {
VectorialPointValuePair optimum = VectorialPointValuePair optimum =
optimizer.optimize(circle, target, weights, new double[] { 0, 0 }); optimizer.optimize(circle, target, weights, new double[] { 0, 0 });
assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1.0e-8); assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1.0e-6);
assertEquals(0.2074999736353867, optimum.getPointRef()[1], 1.0e-8); assertEquals(0.2074999736353867, optimum.getPointRef()[1], 1.0e-6);
assertEquals(0.04268731682389561, optimizer.getRMS(), 1.0e-8); assertEquals(0.04268731682389561, optimizer.getRMS(), 1.0e-8);
} }