LevenbergMarquardt.php

Go to the documentation of this file.

<?php

// Levenberg-Marquardt in PHP

// http://www.idiom.com/~zilla/Computer/Javanumeric/LM.java

class LevenbergMarquardt {

        function chiSquared($x, $a, $y, $s, $f) {
                $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()


        function solve($x, $a, $y, $s, $vary, $f, $lambda, $termepsilon, $maxiter, $verbose) {
                $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()

}       //      class LevenbergMarquardt