Justine Devin / gammalib

Commit 7b177e085f14dba18e72530db0d92eb72a4a2b6a

Authored by knodlseder
1 parent 961ed5d7

add parameter freezing and reduce set lambda_inc=2, lambda_dec=0.5 (this seepds up convergence)

Inline Side-by-side
Showing 1 changed file with 226 additions and 82 deletions
src/opt/GOptimizerLM.cpp
  View file @7b177e0
... ... @@ -20,8 +20,9 @@
20 20 #ifdef HAVE_CONFIG_H
21 21 #include <config.h>
22 22 #endif
23   -#include "GException.hpp"
24 23 #include "GOptimizerLM.hpp"
  24 +#include "GTools.hpp"
  25 +#include "GException.hpp"
25 26  
26 27 /* __ Method name definitions ____________________________________________ */
27 28  
... ... @@ -180,6 +181,53 @@ GModels&amp; GOptimizerLM::operator() (GOptimizerFunction&amp; fct, GModels&amp; m)
180 181 = =
181 182 ==========================================================================*/
182 183  
  184 +/***********************************************************************//**
  185 + * @brief Print optimizer information
  186 + ***************************************************************************/
  187 +std::string GOptimizerLM::print(void) const
  188 +{
  189 + // Initialise result string
  190 + std::string result;
  191 +
  192 + // Append header
  193 + result.append("=== GOptimizerLM ===");
  194 + result.append("\n"+parformat("Optimized function value")+str(m_value));
  195 + result.append("\n"+parformat("Absolute precision")+str(m_eps));
  196 +
  197 + // Append status
  198 + result.append("\n"+parformat("Optimization status"));
  199 + switch (m_status) {
  200 + case G_LM_CONVERGED:
  201 + result.append("converged");
  202 + break;
  203 + case G_LM_STALLED:
  204 + result.append("stalled");
  205 + break;
  206 + case G_LM_SINGULAR:
  207 + result.append("singular curvature matrix encountered");
  208 + break;
  209 + case G_LM_NOT_POSTIVE_DEFINITE:
  210 + result.append("curvature matrix not positive definite");
  211 + break;
  212 + case G_LM_BAD_ERRORS:
  213 + result.append("errors are inaccurate");
  214 + break;
  215 + default:
  216 + result.append("unknown");
  217 + break;
  218 + }
  219 +
  220 + // Append further information
  221 + result.append("\n"+parformat("Number of parameters")+str(m_npars));
  222 + result.append("\n"+parformat("Number of free parameters")+str(m_nfree));
  223 + result.append("\n"+parformat("Number of iterations")+str(m_iter));
  224 + result.append("\n"+parformat("Lambda")+str(m_lambda));
  225 +
  226 + // Return result
  227 + return result;
  228 +}
  229 +
  230 +
183 231 /*==========================================================================
184 232 = =
185 233 = Private methods =
... ... @@ -192,14 +240,24 @@ GModels&amp; GOptimizerLM::operator() (GOptimizerFunction&amp; fct, GModels&amp; m)
192 240 void GOptimizerLM::init_members(void)
193 241 {
194 242 // Initialise optimizer parameters
  243 + m_npars = 0;
  244 + m_nfree = 0;
195 245 m_lambda_start = 1.0e-3;
196   - m_lambda_inc = 10000.0;
197   - m_lambda_dec = 0.1;
  246 + //m_lambda_inc = 10000.0; //!< This led to bad convergence in some cases
  247 + m_lambda_inc = 2.0;
  248 + m_lambda_dec = 0.5;
198 249 m_eps = 1.0e-6;
199   - m_max_iter = 100;
  250 + m_max_iter = 1000;
200 251 m_max_stall = 10;
  252 + m_max_hit = 3; //!< Maximum number of successive boundary hits before freeze
201 253 m_step_adjust = true;
202 254  
  255 + // Initialise bookkeeping arrays
  256 + m_hit_boundary.clear();
  257 + m_hit_minimum.clear();
  258 + m_hit_maximum.clear();
  259 + m_par_freeze.clear();
  260 +
203 261 // Initialise optimizer values
204 262 m_lambda = m_lambda_start;
205 263 m_value = 0.0;
... ... @@ -229,6 +287,10 @@ void GOptimizerLM::copy_members(const GOptimizerLM&amp; opt)
229 287 m_max_iter = opt.m_max_iter;
230 288 m_max_stall = opt.m_max_stall;
231 289 m_step_adjust = opt.m_step_adjust;
  290 + m_hit_boundary = opt.m_hit_boundary;
  291 + m_hit_minimum = opt.m_hit_minimum;
  292 + m_hit_maximum = opt.m_hit_maximum;
  293 + m_par_freeze = opt.m_par_freeze;
232 294 m_lambda = opt.m_lambda;
233 295 m_value = opt.m_value;
234 296 m_status = opt.m_status;
... ... @@ -266,19 +328,32 @@ void GOptimizerLM::optimize(GOptimizerFunction* fct, GOptimizerPars* pars)
266 328 if (fct == NULL || pars == NULL)
267 329 continue;
268 330  
269   - // Fall through if there are no parameters to optimize
270   - int npars = pars->npars();
271   - if (npars < 1)
  331 + // Set number of parameters. Fall through if there are no parameters
  332 + // to optimize
  333 + m_npars = pars->npars();
  334 + m_nfree = pars->nfree();
  335 + if (m_nfree < 1)
272 336 continue;
273 337  
274 338 // Initialise optimization parameters
275 339 m_lambda = m_lambda_start;
276   - m_status = 0;
277   -
278   - // Allocate temporary memory for boundary hit bookkeeping
279   - m_hit_boundary = new bool[npars];
280   - for (int i = 0; i < npars; ++i)
281   - m_hit_boundary[i] = false;
  340 + m_status = G_LM_CONVERGED;
  341 +
  342 + // Initialise bookkeeping arrays
  343 + m_hit_boundary.clear();
  344 + m_hit_minimum.clear();
  345 + m_hit_maximum.clear();
  346 + m_par_freeze.clear();
  347 + m_hit_boundary.reserve(m_npars);
  348 + m_hit_minimum.reserve(m_npars);
  349 + m_hit_maximum.reserve(m_npars);
  350 + m_par_freeze.reserve(m_npars);
  351 + for (int i = 0; i < m_npars; ++i) {
  352 + m_hit_boundary.push_back(false);
  353 + m_hit_minimum.push_back(0);
  354 + m_hit_maximum.push_back(0);
  355 + m_par_freeze.push_back(false);
  356 + }
282 357  
283 358 // Initial evaluation
284 359 fct->eval(*pars);
... ... @@ -302,6 +377,9 @@ void GOptimizerLM::optimize(GOptimizerFunction* fct, GOptimizerPars* pars)
302 377 << m_lambda << std::endl;
303 378 #endif
304 379  
  380 + // Initialise iteration flags
  381 + bool check_for_freeze = true;
  382 +
305 383 // Iterative fitting
306 384 for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) {
307 385  
... ... @@ -316,7 +394,10 @@ void GOptimizerLM::optimize(GOptimizerFunction* fct, GOptimizerPars* pars)
316 394 *m_logger << "Iteration " << m_iter << ": ";
317 395 *m_logger << "func=" << m_value << ", ";
318 396 *m_logger << "Lambda=" << m_lambda << ", ";
319   - *m_logger << "delta=" << delta << std::endl;
  397 + *m_logger << "delta=" << delta;
  398 + if (m_lambda > lambda_old)
  399 + *m_logger << " (stalled)";
  400 + *m_logger << std::endl;
320 401 }
321 402 #if G_DEBUG_OPT
322 403 std::cout << "Iteration " << m_iter << ": func="
... ... @@ -330,24 +411,68 @@ void GOptimizerLM::optimize(GOptimizerFunction* fct, GOptimizerPars* pars)
330 411  
331 412 // If function increased while lambda did not increase then stop
332 413 // iterations
333   - if ((m_lambda <= lambda_old) && (delta < 0.0))
334   - break;
  414 + //if ((m_lambda <= lambda_old) && (delta < 0.0))
  415 + // break;
  416 +
  417 + // Stop if convergence was reached. Before stopping, check
  418 + // if some parameters were frozen, and if this was the case,
  419 + // free them now and continue. We do this only once, i.e.
  420 + // the next time a parameter is frozen and convergence is
  421 + // reached we really stop.
  422 + if ((m_lambda <= lambda_old) && (delta < m_eps)) {
  423 +
  424 + // Check for frozen parameters, and if some exist, free
  425 + // them and start over
  426 + if (check_for_freeze) {
  427 +
  428 + // Signal that we won't check frozen parameters
  429 + // anymore
  430 + check_for_freeze = false;
  431 +
  432 + // Free frozen parameters and determine how many
  433 + // have been frozen
  434 + int nfrozen = 0;
  435 + for (int ipar = 0; ipar < m_npars; ++ipar) {
  436 + if (m_par_freeze[ipar]) {
  437 + nfrozen++;
  438 + pars->par(ipar)->free();
  439 + if (m_logger != NULL) {
  440 + *m_logger << " Free parameter \""
  441 + << pars->par(ipar)->name()
  442 + << "\" after convergence was" \
  443 + " reached with frozen" \
  444 + " parameter." << std::endl;
  445 + }
  446 + }
  447 + }
  448 +
  449 + // If there were frozen parameters then start over
  450 + // again (initialise optimizer)
  451 + if (nfrozen > 0) {
  452 + m_lambda = m_lambda_start;
  453 + value_old = m_value;
  454 + lambda_old = m_lambda;
  455 + lambda_inc = 0;
  456 + continue;
  457 + }
  458 + }
335 459  
336   - // Stop if convergence was reached
337   - if ((m_lambda <= lambda_old) && (delta < m_eps))
  460 + // ... otherwise, convergence was reached and we can stop
  461 + // now
338 462 break;
  463 + }
339 464  
340   - // Monitor the number of subsequent increases of lambda and stop if
341   - // the number of increases exceeds threshold
  465 + // Monitor the number of subsequent increases of lambda and
  466 + // stop if the number of increases exceeds threshold
342 467 lambda_inc = (m_lambda > lambda_old) ? lambda_inc + 1 : 0;
343 468 if (lambda_inc > m_max_stall) {
344 469 m_status = G_LM_STALLED;
345 470 break;
346 471 }
347 472  
348   - // Bookkeeping of actual result (we always store the last lambda to
349   - // detect turn arounds in the lambda tendency; however we always keep
350   - // the best function value)
  473 + // Bookkeeping of actual result (we always store the last
  474 + // lambda to detect turn arounds in the lambda tendency; however
  475 + // we always keep the best function value)
351 476 lambda_old = m_lambda;
352 477 if (delta > 0.0)
353 478 value_old = m_value;
... ... @@ -359,9 +484,6 @@ void GOptimizerLM::optimize(GOptimizerFunction* fct, GOptimizerPars* pars)
359 484  
360 485 } while (0); // endwhile: main loop
361 486  
362   - // Free working memory
363   - if (m_hit_boundary != NULL) delete [] m_hit_boundary;
364   -
365 487 // Return
366 488 return;
367 489 }
... ... @@ -388,11 +510,6 @@ void GOptimizerLM::iteration(GOptimizerFunction* fct, GOptimizerPars* pars)
388 510 if (fct == NULL || pars == NULL)
389 511 continue;
390 512  
391   - // Fall through if there are no parameters to optimize
392   - int npars = pars->npars();
393   - if (npars < 1)
394   - continue;
395   -
396 513 // Initialise iteration parameters
397 514 GVector* grad = fct->gradient();
398 515 GSparseMatrix* covar = fct->covar();
... ... @@ -403,12 +520,12 @@ void GOptimizerLM::iteration(GOptimizerFunction* fct, GOptimizerPars* pars)
403 520 GSparseMatrix save_covar = GSparseMatrix(*covar);
404 521  
405 522 // Save parameter values in vector
406   - GVector save_pars(npars);
407   - for (int ipar = 0; ipar < npars; ++ipar)
  523 + GVector save_pars(m_npars);
  524 + for (int ipar = 0; ipar < m_npars; ++ipar)
408 525 save_pars(ipar) = pars->par(ipar)->value();
409 526  
410 527 // Setup matrix and vector for covariance computation
411   - for (int ipar = 0; ipar < npars; ++ipar) {
  528 + for (int ipar = 0; ipar < m_npars; ++ipar) {
412 529 (*covar)(ipar,ipar) *= (1.0 + m_lambda);
413 530 (*grad)(ipar) = -(*grad)(ipar);
414 531 }
... ... @@ -444,36 +561,73 @@ void GOptimizerLM::iteration(GOptimizerFunction* fct, GOptimizerPars* pars)
444 561 double step = step_size(grad, pars);
445 562  
446 563 // Derive new parameter vector
447   - for (int ipar = 0; ipar < npars; ++ipar) {
  564 + for (int ipar = 0; ipar < m_npars; ++ipar) {
448 565  
449   - // Get actual parameter value and limits
450   - double p = pars->par(ipar)->value();
451   - double p_min = pars->par(ipar)->min();
452   - double p_max = pars->par(ipar)->max();
  566 + // Consider only free parameters
  567 + if (pars->par(ipar)->isfree()) {
453 568  
454   - // Compute new parameter value
455   - p += (*grad)(ipar) * step;
  569 + // Get actual parameter value and limits
  570 + double p = pars->par(ipar)->value();
  571 + double p_min = pars->par(ipar)->min();
  572 + double p_max = pars->par(ipar)->max();
456 573  
457   - // Constrain parameter to within the valid range
458   - if (pars->par(ipar)->hasmin() && p < p_min) {
459   - if (m_logger != NULL) {
460   - *m_logger << "... parameter \"" << pars->par(ipar)->name();
461   - *m_logger << "\" hits minimum: ";
462   - *m_logger << p << " < " << p_min << std::endl;
  574 + // Compute new parameter value
  575 + p += (*grad)(ipar) * step;
  576 +
  577 + // Constrain parameter to within the valid range
  578 + if (pars->par(ipar)->hasmin() && p < p_min) {
  579 + if (m_hit_minimum[ipar] >= m_max_hit) {
  580 + if (m_logger != NULL) {
  581 + *m_logger << " Parameter \"" << pars->par(ipar)->name();
  582 + *m_logger << "\" hits minimum " << p_min << " more than ";
  583 + *m_logger << m_max_hit << " times.";
  584 + *m_logger << " Fix parameter at minimum for now." << std::endl;
  585 + }
  586 + m_par_freeze[ipar] = true;
  587 + pars->par(ipar)->fix();
  588 + }
  589 + else {
  590 + if (m_logger != NULL) {
  591 + *m_logger << " Parameter \"" << pars->par(ipar)->name();
  592 + *m_logger << "\" hits minimum: ";
  593 + *m_logger << p << " < " << p_min;
  594 + *m_logger << " (" << m_hit_minimum[ipar]+1 << ")" << std::endl;
  595 + }
  596 + }
  597 + m_hit_minimum[ipar]++;
  598 + p = p_min;
463 599 }
464   - p = p_min;
465   - }
466   - if (pars->par(ipar)->hasmax() && p > p_max) {
467   - if (m_logger != NULL) {
468   - *m_logger << "... parameter \"" << pars->par(ipar)->name();
469   - *m_logger << "\" hits maximum: ";
470   - *m_logger << p << " > " << p_max << std::endl;
  600 + else if (pars->par(ipar)->hasmax() && p > p_max) {
  601 + if (m_hit_maximum[ipar] >= m_max_hit) {
  602 + if (m_logger != NULL) {
  603 + *m_logger << " Parameter \"" << pars->par(ipar)->name();
  604 + *m_logger << "\" hits maximum " << p_max << " more than ";
  605 + *m_logger << m_max_hit << " times.";
  606 + *m_logger << " Fix parameter at maximum for now." << std::endl;
  607 + }
  608 + m_par_freeze[ipar] = true;
  609 + pars->par(ipar)->fix();
  610 + }
  611 + else {
  612 + if (m_logger != NULL) {
  613 + *m_logger << " Parameter \"" << pars->par(ipar)->name();
  614 + *m_logger << "\" hits maximum: ";
  615 + *m_logger << p << " > " << p_max;
  616 + *m_logger << " (" << m_hit_maximum[ipar]+1 << ")" << std::endl;
  617 + }
  618 + }
  619 + m_hit_maximum[ipar]++;
  620 + p = p_max;
471 621 }
472   - p = p_max;
473   - }
  622 + else {
  623 + m_hit_minimum[ipar] = 0;
  624 + m_hit_maximum[ipar] = 0;
  625 + }
  626 +
  627 + // Set new parameter value
  628 + pars->par(ipar)->value(p);
474 629  
475   - // Set new parameter value
476   - pars->par(ipar)->value(p);
  630 + } // endif: Parameter was free
477 631  
478 632 } // endfor: computed new parameter vector
479 633  
... ... @@ -489,7 +643,7 @@ void GOptimizerLM::iteration(GOptimizerFunction* fct, GOptimizerPars* pars)
489 643  
490 644 // Determine how many parameters have changed
491 645 int par_change = 0;
492   - for (int ipar = 0; ipar < npars; ++ipar) {
  646 + for (int ipar = 0; ipar < m_npars; ++ipar) {
493 647 if (pars->par(ipar)->value() != save_pars(ipar))
494 648 par_change++;
495 649 }
... ... @@ -522,7 +676,7 @@ void GOptimizerLM::iteration(GOptimizerFunction* fct, GOptimizerPars* pars)
522 676 m_value = save_value;
523 677 *grad = save_grad;
524 678 *covar = save_covar;
525   - for (int ipar = 0; ipar < npars; ++ipar)
  679 + for (int ipar = 0; ipar < m_npars; ++ipar)
526 680 pars->par(ipar)->value(save_pars(ipar));
527 681 }
528 682  
... ... @@ -578,7 +732,7 @@ double GOptimizerLM::step_size(GVector* grad, GOptimizerPars* pars)
578 732 else if (m_hit_boundary[ipar]) {
579 733 m_hit_boundary[ipar] = false;
580 734 if (m_logger != NULL) {
581   - *m_logger << "... parameter \"";
  735 + *m_logger << " Parameter \"";
582 736 *m_logger << pars->par(ipar)->name();
583 737 *m_logger << "\" does not drive optimization step anymore.";
584 738 *m_logger << std::endl;
... ... @@ -600,7 +754,7 @@ double GOptimizerLM::step_size(GVector* grad, GOptimizerPars* pars)
600 754 else if (m_hit_boundary[ipar]) {
601 755 m_hit_boundary[ipar] = false;
602 756 if (m_logger != NULL) {
603   - *m_logger << "... parameter \"";
  757 + *m_logger << " Parameter \"";
604 758 *m_logger << pars->par(ipar)->name();
605 759 *m_logger << "\" does not drive optimization step anymore.";
606 760 *m_logger << std::endl;
... ... @@ -615,7 +769,7 @@ double GOptimizerLM::step_size(GVector* grad, GOptimizerPars* pars)
615 769 if (ipar_bnd != -1) {
616 770 m_hit_boundary[ipar_bnd] = true;
617 771 if (m_logger != NULL) {
618   - *m_logger << "... parameter \"";
  772 + *m_logger << " Parameter \"";
619 773 *m_logger << pars->par(ipar_bnd)->name();
620 774 *m_logger << "\" drives optimization step (step=";
621 775 *m_logger << step << ")" << std::endl;
... ... @@ -747,20 +901,15 @@ void GOptimizerLM::errors(GOptimizerFunction* fct, GOptimizerPars* pars)
747 901 ==========================================================================*/
748 902  
749 903 /***********************************************************************//**
750   - * @brief Put optimizer in output stream
  904 + * @brief Output operator
751 905 *
752   - * @param[in] os Output stream into which the optimizer will be dumped.
753   - * @param[in] opt Object to be dumped.
  906 + * @param[in] os Output stream.
  907 + * @param[in] opt Optimizer.
754 908 ***************************************************************************/
755 909 std::ostream& operator<< (std::ostream& os, const GOptimizerLM& opt)
756 910 {
757   - // Put optimizer in stream
758   - os << "=== GOptimizerLM ===" << std::endl;
759   - os << " Optimized function value ..: " << opt.m_value << std::endl;
760   - os << " Absolute precision ........: " << opt.m_eps << std::endl;
761   - os << " Optimization status .......: " << opt.m_status << std::endl;
762   - os << " Number of iterations ......: " << opt.m_iter << std::endl;
763   - os << " Lambda ....................: " << opt.m_lambda;
  911 + // Write optimizer in output stream
  912 + os << opt.print();
764 913  
765 914 // Return output stream
766 915 return os;
... ... @@ -768,20 +917,15 @@ std::ostream&amp; operator&lt;&lt; (std::ostream&amp; os, const GOptimizerLM&amp; opt)
768 917  
769 918  
770 919 /***********************************************************************//**
771   - * @brief Write optimizer into logger
  920 + * @brief Log operator
772 921 *
773   - * @param[in] log Logger
774   - * @param[in] opt Optimizer to be writted into logger.
  922 + * @param[in] log Logger.
  923 + * @param[in] opt Optimizer.
775 924 ***************************************************************************/
776 925 GLog& operator<< (GLog& log, const GOptimizerLM& opt)
777 926 {
778 927 // Write optimizer into logger
779   - log << "=== GOptimizerLM ===" << std::endl;
780   - log << " Optimized function value ..: " << opt.m_value << std::endl;
781   - log << " Absolute precision ........: " << opt.m_eps << std::endl;
782   - log << " Optimization status .......: " << opt.m_status << std::endl;
783   - log << " Number of iterations ......: " << opt.m_iter << std::endl;
784   - log << " Lambda ....................: " << opt.m_lambda;
  928 + log << opt.print();
785 929  
786 930 // Return logger
787 931 return log;
... ...