newton.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 }