Molassembler  3.0.0
Molecule graph and conformer library
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NewtonRaphson.h
Go to the documentation of this file.
1 
8 #ifndef INCLUDE_TEMPLE_NEWTON_RAPHSON_H
9 #define INCLUDE_TEMPLE_NEWTON_RAPHSON_H
10 
11 #include <Eigen/Core>
12 #include <Eigen/Eigenvalues>
13 #include <unsupported/Eigen/NumericalDiff>
14 #include <iostream>
15 
16 namespace Scine {
17 namespace Molassembler {
18 namespace Temple {
19 
29 template<typename FloatType = double>
30 struct NewtonRaphson {
31  using MatrixType = Eigen::Matrix<FloatType, Eigen::Dynamic, Eigen::Dynamic>;
32  using VectorType = Eigen::Matrix<FloatType, Eigen::Dynamic, 1>;
33 
37  unsigned iterations;
39  FloatType value;
41  VectorType gradient;
42  };
43 
44  template<
45  typename UpdateFunction,
46  typename Checker
47  > OptimizationReturnType minimize(
48  Eigen::Ref<VectorType> parameters,
49  UpdateFunction&& function,
50  Checker&& check
51  ) {
52  const unsigned P = parameters.size();
53  FloatType value;
54  VectorType gradients = Eigen::VectorXd::Zero(P);
55  MatrixType hessian (P, P);
56  hessian.setZero();
57 
58  VectorType numericalGradient(P);
59  {
60  constexpr FloatType h = 1e-4;
61  FloatType a;
62  FloatType b;
63  VectorType diffParameters = parameters;
64  for(unsigned i = 0; i < P; ++i) {
65  diffParameters(i) += h;
66  function(diffParameters, b, gradients, hessian);
67  diffParameters(i) -= 2 * h;
68  function(diffParameters, a, gradients, hessian);
69  diffParameters(i) = parameters(i);
70  numericalGradient(i) = (b - a) / (2 * h);
71  }
72  }
73 
74  function(parameters, value, gradients, hessian);
75 
76  if(numericalGradient.norm() > 1e-4) {
77  if(!numericalGradient.isApprox(gradients, 1e-8)) {
78  std::cout << "MISMATCH: Numerical gradient is " << numericalGradient.transpose() << "\nfunction gradient is : " << gradients.transpose() << "\n";
79  } else {
80  std::cout << "match: norm is " << gradients.norm() << " and diff norm is " << (numericalGradient - gradients).norm() << "\n";
81  }
82  }
83 
84  MatrixType numericalHessian(P, P);
85  {
86  constexpr FloatType h = 1e-4;
87  FloatType a;
88  FloatType b;
89  FloatType c;
90  FloatType d;
91  VectorType diffParameters = parameters;
92  for(unsigned i = 0; i < P; ++i) {
93  // Diagonal second derivative from forward and backward derivative
94  function(diffParameters, b, gradients, hessian);
95  diffParameters(i) += h;
96  function(diffParameters, c, gradients, hessian);
97  diffParameters(i) -= 2 * h;
98  function(diffParameters, a, gradients, hessian);
99  diffParameters(i) = parameters(i);
100  numericalHessian(i, i) = (a - 2 * b + c) / (h * h);
101 
102  // Cross terms from two central derivatives
103  for(unsigned j = i + 1; j < P; ++j) {
104  // (i) + h, (j) + h -> a
105  diffParameters(i) += h;
106  diffParameters(j) += h;
107  function(diffParameters, a, gradients, hessian);
108  // (i) + h, (j) - h - > b
109  diffParameters(j) -= 2 * h;
110  function(diffParameters, b, gradients, hessian);
111  // (i) - h, (j) - h -> d
112  diffParameters(i) -= 2 * h;
113  function(diffParameters, d, gradients, hessian);
114  // (i) - h, (j) + h -> c
115  diffParameters(j) += 2 * h;
116  function(diffParameters, c, gradients, hessian);
117 
118  numericalHessian(i, j) = (a - b - c + d) / (4 * h * h);
119  numericalHessian(j, i) = numericalHessian(i, j);
120  }
121  }
122  }
123 
124  function(parameters, value, gradients, hessian);
125 
126  if(!numericalHessian.isApprox(hessian), 1e-2) {
127  std::cout << "MISMATCH: Hessian is\n" << hessian << "\nNumerical hessian is\n" << numericalHessian << "\n";
128  } else {
129  std::cout << "Hessian matches\n";
130  }
131 
132  unsigned iteration = 0;
133  while(
134  check.shouldContinue(
135  iteration,
136  std::as_const(value),
137  std::as_const(gradients)
138  )
139  ) {
140  // Solve HΔx = g, then apply x_(n+1) = x_n - Δx
141  Eigen::JacobiSVD<MatrixType> svd(
142  hessian,
143  Eigen::ComputeFullV | Eigen::ComputeFullU
144  );
145  svd.setThreshold(svdThreshold);
146  VectorType steps = svd.solve(gradients);
147 
148  // Take a step
149  FloatType rms = std::sqrt(steps.squaredNorm() / P);
150  if (rms > trustRadius) {
151  steps *= trustRadius / rms;
152  }
153  parameters -= steps;
154 
155  function(parameters, value, gradients, hessian);
156  ++iteration;
157  }
158 
159  return {
160  iteration,
161  value,
162  std::move(gradients)
163  };
164  }
165 
167  FloatType svdThreshold = 1.0e-12;
169  FloatType trustRadius = 0.5;
170 };
171 
172 } // namespace Temple
173 } // namespace Molassembler
174 } // namespace Scine
175 
176 #endif
FloatType value
Final function value.
Definition: NewtonRaphson.h:39
unsigned iterations
Number of iterations.
Definition: NewtonRaphson.h:37
Type returned from an optimization.
Definition: NewtonRaphson.h:35
VectorType gradient
Final gradient.
Definition: NewtonRaphson.h:41
FloatType trustRadius
The maximum RMS of a taken step.
Definition: NewtonRaphson.h:169
FloatType svdThreshold
The SVD threshold for the decomposition of the Hessian.
Definition: NewtonRaphson.h:167
A very basic newton-raphson minimizer.
Definition: NewtonRaphson.h:30