1 #ifdef MANGO_EIGEN_AVAILABLE // Don't bother doing any testing if Eigen is unavailable.
9 TEST_CASE(
"Levenberg_marquardt",
"[algorithm][Levenberg_marquardt]") {
11 int N_line_search = 1;
12 auto log_step = GENERATE(range(0,5));
13 double lambda_step = exp(log_step - 2.0);
15 mango::Levenberg_marquardt::compute_lambda_grid(N_line_search, lambda_step, &lambda_grid);
16 CHECK(lambda_grid == Approx(1.0).epsilon(1e-13));
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]));
31 class Levenberg_marquardt_tester :
public Least_squares_solver {
33 Levenberg_marquardt_tester();
34 ~Levenberg_marquardt_tester();
38 mango::Levenberg_marquardt_tester::Levenberg_marquardt_tester() {
42 best_state_vector =
new double[N_parameters];
43 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];
51 residual_function = &Levenberg_marquardt_residual_function_1;
52 function_evaluations = 0;
53 max_function_evaluations = 10000;
57 state_vector[0] = 1.2;
58 state_vector[1] = 0.9;
61 for (
int j=0; j<N_terms; j++) {
62 targets[j] = 1.5 + 2 * j;
63 sigmas[j] = 0.8 + 1.3 * j;
66 centered_differences =
true;
67 finite_difference_step_size = 1.0e-7;
70 mango::Levenberg_marquardt_tester::~Levenberg_marquardt_tester() {
71 delete[] state_vector;
74 delete[] best_residual_function;
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];
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];
99 residual_function = &Levenberg_marquardt_residual_function_1;
100 function_evaluations = 0;
104 for (
int j=0; j<N_terms; j++) {
107 targets[j] = GENERATE(take(1,random(-10.0,10.0)));
108 sigmas[j] = GENERATE(take(1,random(0.1,10.0)));
112 mpi_partition->set_N_worker_groups(1);
113 mpi_partition->init(MPI_COMM_WORLD);
115 N_line_search = GENERATE(range(1,6));
116 CAPTURE(N_parameters, N_terms, N_line_search);
119 lm.save_lambda_history =
false;
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;
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)));
135 if (mpi_partition->get_proc0_worker_groups()) {
136 lm.process_lambda_grid_results();
140 if (lm.proc0_world) {
143 CHECK(function_evaluations == N_line_search);
146 double best_objective_function = residuals_to_single_objective(lm.lambda_scan_residuals.col(lm.min_objective_function_index).data());
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);
154 if (lm.line_search_succeeded) {
157 CHECK(lm.objective_function == Approx(best_objective_function));
159 CHECK(lm.central_lambda == Approx(original_central_lambda * lm.normalized_lambda_grid[lm.min_objective_function_index] / lm.lambda_reduction_on_success));
163 CHECK(lm.objective_function < best_objective_function);
165 CHECK(lm.central_lambda == Approx(original_central_lambda * lm.lambda_increase_factor));
170 delete[] state_vector;
173 delete[] best_residual_function;
177 TEST_CASE_METHOD(mango::Levenberg_marquardt_tester,
"mango::Levenberg_marquardt::evaluate_on_lambda_grid() and line_search()",
178 "[Levenberg_marquardt]") {
184 mpi_partition->set_N_worker_groups(GENERATE(range(1,6)));
185 mpi_partition->init(MPI_COMM_WORLD);
191 lm.check_least_squares_solution =
false;
192 lm.save_lambda_history =
false;
194 lm.state_vector << 0.6, -0.8;
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;
212 lm.residuals_extended << 0.7, 0.9, -1.2, -1.9, 0.0, 0.0;
219 SECTION(
"Check mango::Levenberg_marquardt::evaluate_on_lambda_grid()") {
222 if (mpi_partition->get_proc0_worker_groups()) {
223 lm.evaluate_on_lambda_grid();
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;
231 if (mpi_partition->get_proc0_worker_groups()) {
234 CHECK(lm.line_search_succeeded ==
true);
236 CHECK(lm.state_vector[0] == Approx(4.5835619960968432e-01));
237 CHECK(lm.state_vector[1] == Approx(-1.0447504305249349e+00));
242 if (mpi_partition->get_proc0_world()) {
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;
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;
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)));
263 for (
int k=0; k<N_terms; k++) {
264 CHECK(lambda_scan_residuals(k,j) == Approx(lm.lambda_scan_residuals(k,j)));
272 TEST_CASE_METHOD(mango::Levenberg_marquardt_tester,
"mango::Levenberg_marquardt::solve()",
273 "[Levenberg_marquardt]") {
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);
285 lm.check_least_squares_solution =
true;
286 lm.save_lambda_history =
false;
287 lm.max_outer_iterations = 1;
290 lm.state_vector << 0.0, 0.0;
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;
303 if (mpi_partition->get_proc0_worker_groups()) {
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));
317 CHECK(lm.state_vector(0) == Approx(0.0));
318 CHECK(lm.state_vector(1) == Approx(-5.3731855765138470e-01));
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));
326 CHECK(lm.central_lambda == Approx(2.0e-4));
332 #endif // MANGO_EIGEN_AVAILABLE