ILIAS  Release_4_0_x_branch Revision 61816
 All Data Structures Namespaces Files Functions Variables Groups Pages
Public Member Functions
LevenbergMarquardt Class Reference
+ Collaboration diagram for LevenbergMarquardt:

Public Member Functions

 chiSquared ($x, $a, $y, $s, $f)
 Calculate the current sum-squared-error.
 solve ($x, $a, $y, $s, $vary, $f, $lambda, $termepsilon, $maxiter, $verbose)
 Minimize E = sum {(y[k] - f(x[k],a)) / s[k]}^2 The individual errors are optionally scaled by s[k].

Detailed Description

Definition at line 7 of file LevenbergMarquardt.php.

Member Function Documentation

LevenbergMarquardt::chiSquared (   $x,
  $a,
  $y,
  $s,
  $f 
)

Calculate the current sum-squared-error.

Chi-squared is the distribution of squared Gaussian errors, thus the name.

Parameters
double[][]$x
double[]$a
double[]$y,
double[]$s,
object$f

Definition at line 21 of file LevenbergMarquardt.php.

References $d, $f, $x, and $y.

Referenced by solve().

{
$npts = count($y);
$sum = 0.0;
for ($i = 0; $i < $npts; ++$i) {
$d = $y[$i] - $f->val($x[$i], $a);
$d = $d / $s[$i];
$sum = $sum + ($d*$d);
}
return $sum;
} // function chiSquared()

+ Here is the caller graph for this function:

LevenbergMarquardt::solve (   $x,
  $a,
  $y,
  $s,
  $vary,
  $f,
  $lambda,
  $termepsilon,
  $maxiter,
  $verbose 
)

Minimize E = sum {(y[k] - f(x[k],a)) / s[k]}^2 The individual errors are optionally scaled by s[k].

Note that LMfunc implements the value and gradient of f(x,a), NOT the value and gradient of E with respect to a!

Parameters
xarray of domain points, each may be multidimensional
ycorresponding array of values
athe parameters/state of the model
varyfalse to indicate the corresponding a[k] is to be held fixed
s2sigma^2 for point i
lambdablend between steepest descent (lambda high) and jump to bottom of quadratic (lambda zero). Start with 0.001.
termepsilontermination accuracy (0.01)
maxiterstop and return after this many iterations if not done
verboseset to zero (no prints), 1, 2
Returns
the new lambda for future iterations. Can use this and maxiter to interleave the LM descent with some other task, setting maxiter to something small.

Definition at line 57 of file LevenbergMarquardt.php.

References $f, $verbose, $x, $y, and chiSquared().

{
$npts = count($y);
$nparm = count($a);
if ($verbose > 0) {
print("solve x[".count($x)."][".count($x[0])."]");
print(" a[".count($a)."]");
println(" y[".count(length)."]");
}
$e0 = $this->chiSquared($x, $a, $y, $s, $f);
//double lambda = 0.001;
$done = false;
// g = gradient, H = hessian, d = step to minimum
// H d = -g, solve for d
$H = array();
$g = array();
//double[] d = new double[nparm];
$oos2 = array();
for($i = 0; $i < $npts; ++$i) {
$oos2[$i] = 1./($s[$i]*$s[$i]);
}
$iter = 0;
$term = 0; // termination count test
do {
++$iter;
// hessian approximation
for( $r = 0; $r < $nparm; ++$r) {
for( $c = 0; $c < $nparm; ++$c) {
for( $i = 0; $i < $npts; ++$i) {
if ($i == 0) $H[$r][$c] = 0.;
$xi = $x[$i];
$H[$r][$c] += ($oos2[$i] * $f->grad($xi, $a, $r) * $f->grad($xi, $a, $c));
} //npts
} //c
} //r
// boost diagonal towards gradient descent
for( $r = 0; $r < $nparm; ++$r)
$H[$r][$r] *= (1. + $lambda);
// gradient
for( $r = 0; $r < $nparm; ++$r) {
for( $i = 0; $i < $npts; ++$i) {
if ($i == 0) $g[$r] = 0.;
$xi = $x[$i];
$g[$r] += ($oos2[$i] * ($y[$i]-$f->val($xi,$a)) * $f->grad($xi, $a, $r));
}
} //npts
// scale (for consistency with NR, not necessary)
if ($false) {
for( $r = 0; $r < $nparm; ++$r) {
$g[$r] = -0.5 * $g[$r];
for( $c = 0; $c < $nparm; ++$c) {
$H[$r][$c] *= 0.5;
}
}
}
// solve H d = -g, evaluate error at new location
//double[] d = DoubleMatrix.solve(H, g);
// double[] d = (new Matrix(H)).lu().solve(new Matrix(g, nparm)).getRowPackedCopy();
//double[] na = DoubleVector.add(a, d);
// double[] na = (new Matrix(a, nparm)).plus(new Matrix(d, nparm)).getRowPackedCopy();
// double e1 = chiSquared(x, na, y, s, f);
// if (verbose > 0) {
// System.out.println("\n\niteration "+iter+" lambda = "+lambda);
// System.out.print("a = ");
// (new Matrix(a, nparm)).print(10, 2);
// if (verbose > 1) {
// System.out.print("H = ");
// (new Matrix(H)).print(10, 2);
// System.out.print("g = ");
// (new Matrix(g, nparm)).print(10, 2);
// System.out.print("d = ");
// (new Matrix(d, nparm)).print(10, 2);
// }
// System.out.print("e0 = " + e0 + ": ");
// System.out.print("moved from ");
// (new Matrix(a, nparm)).print(10, 2);
// System.out.print("e1 = " + e1 + ": ");
// if (e1 < e0) {
// System.out.print("to ");
// (new Matrix(na, nparm)).print(10, 2);
// } else {
// System.out.println("move rejected");
// }
// }
// termination test (slightly different than NR)
// if (Math.abs(e1-e0) > termepsilon) {
// term = 0;
// } else {
// term++;
// if (term == 4) {
// System.out.println("terminating after " + iter + " iterations");
// done = true;
// }
// }
// if (iter >= maxiter) done = true;
// in the C++ version, found that changing this to e1 >= e0
// was not a good idea. See comment there.
//
// if (e1 > e0 || Double.isNaN(e1)) { // new location worse than before
// lambda *= 10.;
// } else { // new location better, accept new parameters
// lambda *= 0.1;
// e0 = e1;
// // simply assigning a = na will not get results copied back to caller
// for( int i = 0; i < nparm; i++ ) {
// if (vary[i]) a[i] = na[i];
// }
// }
} while(!$done);
return $lambda;
} // function solve()

+ Here is the call graph for this function:

The documentation for this class was generated from the following file: