commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From dim...@apache.org
Subject svn commit: r985828 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/optimization/general/ site/xdoc/
Date Mon, 16 Aug 2010 08:45:11 GMT
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.



Mime
View raw message