Levenberg_marquardt_tests.cpp
Go to the documentation of this file.
1 #ifdef MANGO_EIGEN_AVAILABLE // Don't bother doing any testing if Eigen is unavailable.
2 
3 #include <iostream>
4 #include <iomanip>
5 #include "catch.hpp"
7 
8 
9 TEST_CASE("Levenberg_marquardt","[algorithm][Levenberg_marquardt]") {
10  // If N_line_search==1, the normalized lambda grid should be [1.0] for any lambda_step
11  int N_line_search = 1;
12  auto log_step = GENERATE(range(0,5)); // so log_step = 0, 1, 2, 3, or 4.
13  double lambda_step = exp(log_step - 2.0);
14  double lambda_grid;
15  mango::Levenberg_marquardt::compute_lambda_grid(N_line_search, lambda_step, &lambda_grid);
16  CHECK(lambda_grid == Approx(1.0).epsilon(1e-13));
17 }
18 
19 //! An example residual function that is used for unit testing
20 void Levenberg_marquardt_residual_function_1(int* N_parameters, const double* x, int* N_terms, double* f, int* failed_int, mango::Problem* problem, void* user_data) {
21  assert(*N_parameters == 2);
22  assert(*N_terms == 4);
23  for (int j = 0; j < *N_terms; j++) {
24  f[j] = exp(j + x[0] * x[0] - exp(x[1]));
25  }
26  *failed_int = false;
27 }
28 
29 // Create a subclass that handles setup and tear-down:
30 namespace mango {
31  class Levenberg_marquardt_tester : public Least_squares_solver {
32  public:
33  Levenberg_marquardt_tester();
34  ~Levenberg_marquardt_tester();
35  };
36 }
37 
38 mango::Levenberg_marquardt_tester::Levenberg_marquardt_tester() {
39  // The Catch2 macros automatically call the mango::problem() constructor (the version with no arguments).
40  N_parameters = 2;
41  N_terms = 4;
42  best_state_vector = new double[N_parameters];
43  residuals = new double[N_terms]; // We must allocate this variable since the destructor will delete it.
44  //double* base_case_residuals = new double[N_terms];
45  state_vector = new double[N_parameters];
46  targets = new double[N_terms];
47  sigmas = new double[N_terms];
48  best_residual_function = new double[N_terms];
49  //double base_case_objective_function;
50  mpi_partition = new mango::MPI_Partition();
51  residual_function = &Levenberg_marquardt_residual_function_1;
52  function_evaluations = 0;
53  max_function_evaluations = 10000;
54  verbose = 0;
55 
56  // Set the point about which we will compute the derivatives:
57  state_vector[0] = 1.2;
58  state_vector[1] = 0.9;
59 
60  // Initialize targets and sigmas
61  for (int j=0; j<N_terms; j++) {
62  targets[j] = 1.5 + 2 * j;
63  sigmas[j] = 0.8 + 1.3 * j;
64  }
65 
66  centered_differences = true;
67  finite_difference_step_size = 1.0e-7;
68 }
69 
70 mango::Levenberg_marquardt_tester::~Levenberg_marquardt_tester() {
71  delete[] state_vector;
72  delete[] targets;
73  delete[] sigmas;
74  delete[] best_residual_function;
75  //delete[] base_case_residuals;
76  // best_state_vector and residuals will be deleted by destructor.
77 }
78 
79 
80 /////////////////////////////////////////////////////////////////////////////////
81 /////////////////////////////////////////////////////////////////////////////////
82 
83 
84 TEST_CASE_METHOD(mango::Least_squares_solver, "mango::Levenberg_marquardt::process_lambda_grid_results().","[Levenberg_marquardt]") {
85  // First, set up some variables that are needed.
86 
87  // The Catch2 macros automatically call the mango::problem() constructor (the version with no arguments).
88  N_parameters = GENERATE(range(1,6));
89  N_terms = GENERATE(range(1,6));
90  best_state_vector = new double[N_parameters];
91  residuals = new double[N_terms]; // We must allocate this variable since the destructor will delete it.
92  //double* base_case_residuals = new double[N_terms];
93  state_vector = new double[N_parameters];
94  targets = new double[N_terms];
95  sigmas = new double[N_terms];
96  best_residual_function = new double[N_terms];
97  //double base_case_objective_function;
98  mpi_partition = new mango::MPI_Partition();
99  residual_function = &Levenberg_marquardt_residual_function_1;
100  function_evaluations = 0;
101  verbose = 0;
102 
103  // Initialize targets and sigmas
104  for (int j=0; j<N_terms; j++) {
105  //targets[j] = 1.5 + 2 * j;
106  //sigmas[j] = 0.8 + 1.3 * j;
107  targets[j] = GENERATE(take(1,random(-10.0,10.0)));
108  sigmas[j] = GENERATE(take(1,random(0.1,10.0)));
109  }
110 
111  // Set up MPI:
112  mpi_partition->set_N_worker_groups(1); // No reason to scan this, since only proc0_world does any work.
113  mpi_partition->init(MPI_COMM_WORLD);
114 
115  N_line_search = GENERATE(range(1,6));
116  CAPTURE(N_parameters, N_terms, N_line_search);
117 
119  lm.save_lambda_history = false;
120 
121  lm.line_search_succeeded = false;
122  double original_central_lambda = GENERATE(take(1,random(-10.0,10.0)));
123  lm.central_lambda = original_central_lambda;
124  lm.lambda_increase_factor = GENERATE(take(1,random(2.0,30.0)));
125  double original_objective_function = GENERATE(take(1,random(0.1,10.0)));
126  lm.objective_function = original_objective_function;
127 
128  for (int k=0; k<N_line_search; k++) {
129  lm.normalized_lambda_grid[k] = GENERATE(take(1,random(-10.0,10.0)));
130  for (int j=0; j<N_terms; j++) lm.lambda_scan_residuals(j,k) = GENERATE(take(1,random(-10.0,10.0)));
131  for (int j=0; j<N_parameters; j++) lm.lambda_scan_state_vectors(j,k) = GENERATE(take(1,random(-10.0,10.0)));
132  }
133 
134  // Call the subroutine we want to test:
135  if (mpi_partition->get_proc0_worker_groups()) { // Only group leaders should call the subroutine
136  lm.process_lambda_grid_results();
137  }
138 
139  // The subroutine only does computations on proc0_world, so only run the tests on proc0_world:
140  if (lm.proc0_world) {
141 
142  // Verify that function_evaluations was incremented appropriately:
143  CHECK(function_evaluations == N_line_search);
144 
145  // Verify that we succeeded in picking out the best objective function from the set:
146  double best_objective_function = residuals_to_single_objective(lm.lambda_scan_residuals.col(lm.min_objective_function_index).data());
147  double f;
148  for (int j=0; j<N_line_search; j++) {
149  f = residuals_to_single_objective(lm.lambda_scan_residuals.col(j).data());
150  CHECK(f >= best_objective_function);
151  }
152 
153  // Verify that the best-yet objective function in LM is <= the best objective function from the lambda scan:
154  if (lm.line_search_succeeded) {
155  // If we successfully lowered the objective function,
156  // the best-yet objective function in LM should equal the best objective function from the lambda scan:
157  CHECK(lm.objective_function == Approx(best_objective_function));
158  // and lambda should have been decreased:
159  CHECK(lm.central_lambda == Approx(original_central_lambda * lm.normalized_lambda_grid[lm.min_objective_function_index] / lm.lambda_reduction_on_success));
160  } else {
161  // If we failed to lower the objective function,
162  // the best-yet objective function in LM should be smaller than the best objective function from the lambda scan:
163  CHECK(lm.objective_function < best_objective_function);
164  // and lambda should have been increased:
165  CHECK(lm.central_lambda == Approx(original_central_lambda * lm.lambda_increase_factor));
166  }
167  }
168 
169  // Clean up arrays that were allocated:
170  delete[] state_vector;
171  delete[] targets;
172  delete[] sigmas;
173  delete[] best_residual_function;
174 }
175 
176 
177 TEST_CASE_METHOD(mango::Levenberg_marquardt_tester, "mango::Levenberg_marquardt::evaluate_on_lambda_grid() and line_search()",
178  "[Levenberg_marquardt]") {
179  // Validate one piece of the parallelized lambda scan.
180  // We won't actually calculate the true Jacobian, just plug in an ad-hoc matrix for testing.
181  // This test case is a mini-regression test just for one subroutine.
182 
183  // Set up MPI:
184  mpi_partition->set_N_worker_groups(GENERATE(range(1,6)));
185  mpi_partition->init(MPI_COMM_WORLD);
186 
187  N_line_search = 4;
188  //centered_differences = true;
189 
191  lm.check_least_squares_solution = false;
192  lm.save_lambda_history = false;
193 
194  lm.state_vector << 0.6, -0.8;
195 
196  // Here we set central_lambda and normalized_lambda_grid so that if the values for these quantities in
197  // normal operation are changed in Levenberg_marquardt.cpp, it won't break this test.
198  lm.central_lambda = 0.02;
199  lm.normalized_lambda_grid[0] = 0.1;
200  lm.normalized_lambda_grid[1] = 0.3;
201  lm.normalized_lambda_grid[2] = 1.0;
202  lm.normalized_lambda_grid[3] = 3.0;
203 
204  // Not the true Jacobian, just some values I made up.
205  lm.Jacobian <<
206  1.1, 2.5,
207  3.2, 4.0,
208  -1.1, -2.5,
209  -3.2, -4.0;
210 
211  // Not the true residuals, just some values I made up.
212  lm.residuals_extended << 0.7, 0.9, -1.2, -1.9, 0.0, 0.0;
213 
214  //if (mpi_partition->get_proc0_world()) {
215  // std::cout << "Jacobian:" << std::endl << lm.Jacobian << std::endl;
216  // std::cout << "residuals_extended:" << std::endl << lm.residuals_extended << std::endl;
217  //}
218 
219  SECTION("Check mango::Levenberg_marquardt::evaluate_on_lambda_grid()") {
220 
221  // Call the subroutine we want to test:
222  if (mpi_partition->get_proc0_worker_groups()) { // Only group leaders should call the subroutine
223  lm.evaluate_on_lambda_grid();
224  }
225 
226  }
227  SECTION("Verify mango::Levenberg_marquardt::line_search() for a case in which the objective function is reduced already on the 1st lambda_grid") {
228  lm.objective_function = 1.0e200; // Set to any value larger than the min objective function for this lambda scan.
229 
230  // Call the subroutine we want to test:
231  if (mpi_partition->get_proc0_worker_groups()) { // Only group leaders should call the subroutine
232  lm.line_search();
233 
234  CHECK(lm.line_search_succeeded == true);
235  //std::cout << "state_vector after:" << std::endl << std::setprecision(16) << std::scientific << lm.state_vector << std::endl;
236  CHECK(lm.state_vector[0] == Approx(4.5835619960968432e-01));
237  CHECK(lm.state_vector[1] == Approx(-1.0447504305249349e+00));
238  }
239 
240  }
241 
242  if (mpi_partition->get_proc0_world()) {
243  // The reference values later in this subroutine were obtained using these next 2 lines:
244  //std::cout << "lambda_scan_state_vectors:" << std::endl << std::setprecision(16) << std::scientific << lm.lambda_scan_state_vectors << std::endl;
245  //std::cout << "lambda_scan_residuals:" << std::endl << std::setprecision(16) << std::scientific << lm.lambda_scan_residuals << std::endl;
246 
247  Eigen::MatrixXd lambda_scan_state_vectors(2,4);
248  lambda_scan_state_vectors <<
249  6.5966619313657804e-01, 6.2136563284770552e-01, 5.4071043403897379e-01, 4.5835619960968432e-01,
250  -1.1993276056100737e+00, -1.1711321102947181e+00, -1.1107753800090787e+00, -1.0447504305249349e+00;
251 
252  Eigen::MatrixXd lambda_scan_residuals(4,4);
253  lambda_scan_residuals <<
254  1.1431215076343009e+00, 1.0790483326566978e+00, 9.6373805043769645e-01, 8.6789095841964070e-01,
255  3.1073264219230281e+00, 2.9331574746897324e+00, 2.6197116298993373e+00, 2.3591722213560140e+00,
256  8.4465889478040292e+00, 7.9731486634579225e+00, 7.1211145193581968e+00, 6.4128949795174135e+00,
257  2.2960209249278702e+01, 2.1673265127480192e+01, 1.9357196196347253e+01, 1.7432055890638424e+01;
258 
259  for (int j=0; j<N_line_search; j++) {
260  for (int k=0; k<N_parameters; k++) {
261  CHECK(lambda_scan_state_vectors(k,j) == Approx(lm.lambda_scan_state_vectors(k,j)));
262  }
263  for (int k=0; k<N_terms; k++) {
264  CHECK(lambda_scan_residuals(k,j) == Approx(lm.lambda_scan_residuals(k,j)));
265  }
266  }
267  }
268 
269 }
270 
271 
272 TEST_CASE_METHOD(mango::Levenberg_marquardt_tester, "mango::Levenberg_marquardt::solve()",
273  "[Levenberg_marquardt]") {
274  // Take one outer step of Levenberg-Marquardt. This test case is a mini regression test.
275 
276  // Set up MPI:
277  auto N_worker_groups = GENERATE(range(1,6));
278  mpi_partition->set_N_worker_groups(N_worker_groups);
279  CAPTURE(N_worker_groups);
280  mpi_partition->init(MPI_COMM_WORLD);
281 
282  N_line_search = 5;
283 
285  lm.check_least_squares_solution = true;
286  lm.save_lambda_history = false; // Should be false eventually
287  lm.max_outer_iterations = 1; // Do only a single outer iteration, for simplicity.
288  lm.verbose = 0;
289 
290  lm.state_vector << 0.0, 0.0;
291 
292  // Here we set central_lambda, normalized_lambda_grid, and lambda_reduction_on_success so that if the values for these quantities in
293  // normal operation are changed in Levenberg_marquardt.cpp, it won't break this test.
294  lm.central_lambda = 1.0e-2;
295  lm.normalized_lambda_grid[0] = 0.1;
296  lm.normalized_lambda_grid[1] = 0.3;
297  lm.normalized_lambda_grid[2] = 1.0;
298  lm.normalized_lambda_grid[3] = 3.0;
299  lm.normalized_lambda_grid[4] = 10.0;
300  lm.lambda_reduction_on_success = 5.0;
301 
302  // Call the subroutine we want to test:
303  if (mpi_partition->get_proc0_worker_groups()) { // Only group leaders should call the subroutine
304  lm.solve();
305  }
306 
307  if (mpi_partition->get_proc0_world()) {
308  CHECK(lm.Jacobian(0,0) == Approx(0.0));
309  CHECK(lm.Jacobian(1,0) == Approx(0.0));
310  CHECK(lm.Jacobian(2,0) == Approx(0.0));
311  CHECK(lm.Jacobian(3,0) == Approx(0.0));
312  CHECK(lm.Jacobian(0,1) == Approx(-4.5984930134579383e-01));
313  CHECK(lm.Jacobian(1,1) == Approx(-4.7619047620416932e-01));
314  CHECK(lm.Jacobian(2,1) == Approx(-7.9949465544636245e-01));
315  CHECK(lm.Jacobian(3,1) == Approx(-1.5721395948669106e+00));
316 
317  CHECK(lm.state_vector(0) == Approx(0.0));
318  CHECK(lm.state_vector(1) == Approx(-5.3731855765138470e-01));
319 
320  CHECK(lm.lambda_scan_objective_functions(0) == Approx(3.0649073234919824e+00));
321  CHECK(lm.lambda_scan_objective_functions(1) == Approx(3.0650616944702622e+00));
322  CHECK(lm.lambda_scan_objective_functions(2) == Approx(3.0656473368484627e+00));
323  CHECK(lm.lambda_scan_objective_functions(3) == Approx(3.0676889601937392e+00));
324  CHECK(lm.lambda_scan_objective_functions(4) == Approx(3.0784332659074516e+00));
325 
326  CHECK(lm.central_lambda == Approx(2.0e-4));
327 
328  }
329 }
330 
331 
332 #endif // MANGO_EIGEN_AVAILABLE
TEST_CASE
TEST_CASE("does_algorithm_exist(): Verify that return value is true for all the algorithms.","[algorithms]")
Definition: algorithms_tests.cpp:23
mango::Least_squares_solver
Definition: Least_squares_solver.hpp:28
mango::Levenberg_marquardt
Definition: Levenberg_marquardt.hpp:13
TEST_CASE_METHOD
TEST_CASE_METHOD(mango::Solver, "Solver::finite_difference_gradient()","[Solver][finite difference]")
Definition: finite_difference_tests.cpp:40
mango
This C++ namespace contains everything related to MANGO.
Definition: Imfil.hpp:7
mango::MPI_Partition
A class for dividing the set of MPI processes into worker groups.
Definition: mango.hpp:195
Levenberg_marquardt.hpp
mango::Problem
Definition: mango.hpp:442