tnewton.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tnewton.cpp (1761B)
       ---
            1 #include <functional>
            2 #include <armadillo>
            3 #include "qrfunc.h"
            4 
            5 using namespace arma;
            6 
            7 // Newton's method
            8 vec newton(std::function<vec(vec)> f, vec x_0, vec dx, double eps)
            9 {
           10   vec x = x_0;
           11   int n = x.size();
           12   mat A(n,n);
           13   vec y(n);
           14   vec fx(n);
           15   vec fy(n);
           16   vec df(n);
           17   vec Dx(n);
           18 
           19   do {
           20     fx = f(x);
           21     for (int j=0; j<n; ++j) {
           22       x[j] += dx[j];
           23       df = f(x) - fx;
           24 
           25       for (int i=0; i<n; ++i)
           26         A(i,j) = df[i]/dx[j];
           27 
           28       x[j] -= dx[j];
           29     }
           30 
           31     // Perform QR decomposition
           32     QR qr(A);
           33     vec fx_neg = -1.0f*fx;
           34     Dx = qr.backsub(fx_neg);
           35 
           36     double lambda = 2.0f;
           37     do {
           38       lambda /= 2.0f;
           39       y = x + Dx * lambda;
           40       fy = f(y);
           41     } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda > 0.02f);
           42 
           43     x = y;
           44     fx = fy;
           45     std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n';
           46   } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps);
           47   
           48   // Return solution
           49   return x;
           50 }
           51 
           52 // Newton's method - the user supplies his own derivatives
           53 vec newtonJac(std::function<vec(vec)> f, vec x_0, vec dx, double eps, 
           54                   mat (*J)(vec))
           55 {
           56   vec x = x_0;
           57   int n = x.size();
           58   vec y(n);
           59   vec fx(n);
           60   vec fy(n);
           61   vec Dx(n);
           62 
           63   fx = f(x);
           64     
           65   do {
           66     
           67     // Get Jacobian matrix
           68     mat Jx = J(x_0);
           69 
           70     // Perform QR decomposition
           71     QR qr(Jx);
           72     vec fx_neg = -1.0f*fx;
           73     Dx = qr.backsub(fx_neg);
           74 
           75     double lambda = 2.0f;
           76     do {
           77       lambda /= 2.0f;
           78       y = x + Dx * lambda;
           79       fy = f(y);
           80     } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda > 0.02f);
           81 
           82     x = y;
           83     fx = fy;
           84     std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n';
           85   } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps);
           86   
           87   // Return solution
           88   return x;
           89 }