Stefan Eschbach / gammalib

Commit 9219f873ecb40b6e9baf0a364e3d7eb1accdff9f

Authored by Jürgen Knödlseder
1 parent e7c83368

Allow for function decrease in LM optimizer. 

A new parameter was introduce to specify the acceptable function decrease in an iteration. Acceptable means that the new parameters are kept for continuing the iterations. This behaviour existed already before, but a low threshold of 1e-6 did not really activate it. Now the threshold was raised to 1.0 and this help ctlike to converge quicker. Methods accept_dec() have been implemented to get and set the acceptable threshold.
Inline Side-by-side
Showing 6 changed files with 303 additions and 31 deletions
ChangeLog
  View file @9219f87
1   -2014-10-30 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2014-10-31 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.9.0 released
4 4 ========================
... ... @@ -21,6 +21,7 @@
21 21 Improve CTA response computation for extended models
22 22 Implement CTA stacked analysis
23 23 Optimize CTA binned and stacked analysis
  24 + LM optimizer now allows logL decreases up to 1.0
24 25 Load events only when needed into GCTAObservation
25 26 Implement full automatic typecasting in Python module
26 27 Fill DETX and DETY values in GCTAResponseIrf::mc
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-09-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -26 October 2014
  5 +31 October 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -138,6 +138,13 @@ method and an explicit call to GOptimizerLM::errors() is now required for
138 138 error computation. GOptimizer::errors() has been added as pure virtual
139 139 method for all optimizer classes.
140 140  
  141 +A new parameter was introduce to specify the acceptable function decrease in
  142 +an iteration. Acceptable means that the new parameters are kept for continuing
  143 +the iterations. This behaviour existed already before, but a low threshold of
  144 +1e-6 did not really activate it. Now the threshold was raised to 1.0 and this
  145 +help ctlike to converge quicker. Methods accept_dec() have been implemented
  146 +to get and set the acceptable threshold.
  147 +
141 148  
142 149 12. Sky module
143 150 --------------
... ...
include/GOptimizerLM.hpp
  View file @9219f87
... ... @@ -67,27 +67,29 @@ public:
67 67 virtual std::string classname(void) const;
68 68 virtual void optimize(GOptimizerFunction& fct, GOptimizerPars& pars);
69 69 virtual void errors(GOptimizerFunction& fct, GOptimizerPars& pars);
70   - virtual double value(void) const { return m_value; } //!< @brief Return function value
71   - virtual int status(void) const { return m_status; } //!< @brief Return optimization status
72   - virtual int iter(void) const { return m_iter; } //!< @brief Return number of iterations
  70 + virtual double value(void) const;
  71 + virtual int status(void) const;
  72 + virtual int iter(void) const;
73 73 virtual std::string print(const GChatter& chatter = NORMAL) const;
74 74  
75 75 // Methods
76   - void max_iter(const int& n) { m_max_iter=n; } //!< @brief Set maximum number of iterations
77   - void max_stalls(const int& n) { m_max_stall=n; } //!< @brief Set maximum number of stalls
78   - void max_boundary_hits(const int& n) { m_max_stall=n; } //!< @brief Set maximum number of boundary hits
79   - void lambda_start(const double& val) { m_lambda_start=val; } //!< @brief Set lambda start value
80   - void lambda_inc(const double& val) { m_lambda_inc=val; } //!< @brief Set lambda increment
81   - void lambda_dec(const double& val) { m_lambda_dec=val; } //!< @brief Set lambda decrement
82   - void eps(const double& eps) { m_eps=eps; } //!< @brief Set convergence precisions
83   - int max_iter(void) const { return m_max_iter; } //!< @brief Return maximum number of iterations
84   - int max_stalls(void) const { return m_max_stall; } //!< @brief Return maximum number of stalls
85   - int max_boundary_hits(void) const { return m_max_hit; } //!< @brief Return maximum number of boundary hits
86   - const double& lambda_start(void) const { return m_lambda_start; } //!< @brief Return lambda start value
87   - const double& lambda_inc(void) const { return m_lambda_inc; } //!< @brief Return lambda increment
88   - const double& lambda_dec(void) const { return m_lambda_dec; } //!< @brief Return lambda decrement
89   - const double& lambda(void) const { return m_lambda; } //!< @brief Return lambda value
90   - const double& eps(void) const { return m_eps; } //!< @brief Return convergence precision
  76 + void max_iter(const int& max_iter);
  77 + void max_stalls(const int& max_stalls);
  78 + void max_boundary_hits(const int& max_hit);
  79 + void lambda_start(const double& value);
  80 + void lambda_inc(const double& value);
  81 + void lambda_dec(const double& value);
  82 + void eps(const double& eps);
  83 + void accept_dec(const double& value);
  84 + int max_iter(void) const;
  85 + int max_stalls(void) const;
  86 + int max_boundary_hits(void) const;
  87 + const double& lambda_start(void) const;
  88 + const double& lambda_inc(void) const;
  89 + const double& lambda_dec(void) const;
  90 + const double& lambda(void) const;
  91 + const double& eps(void) const;
  92 + const double& accept_dec(void) const;
91 93  
92 94 protected:
93 95 // Protected methods
... ... @@ -104,6 +106,7 @@ protected:
104 106 double m_lambda_inc; //!< Lambda increase
105 107 double m_lambda_dec; //!< Lambda decrease
106 108 double m_eps; //!< Absolute precision
  109 + double m_accept_dec; //!< Acceptable function decrease
107 110 int m_max_iter; //!< Maximum number of iterations
108 111 int m_max_stall; //!< Maximum number of stalls
109 112 int m_max_hit; //!< Maximum number of successive hits
... ... @@ -134,4 +137,258 @@ std::string GOptimizerLM::classname(void) const
134 137 return ("GOptimizerLM");
135 138 }
136 139  
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Return function value
  143 + *
  144 + * @return Function value.
  145 + ***************************************************************************/
  146 +inline
  147 +double GOptimizerLM::value(void) const
  148 +{
  149 + return (m_value);
  150 +}
  151 +
  152 +
  153 +/***********************************************************************//**
  154 + * @brief Return optimizer status
  155 + *
  156 + * @return Optimizer status.
  157 + ***************************************************************************/
  158 +inline
  159 +int GOptimizerLM::status(void) const
  160 +{
  161 + return (m_status);
  162 +}
  163 +
  164 +
  165 +/***********************************************************************//**
  166 + * @brief Return number of iterations
  167 + *
  168 + * @return Number of iterations.
  169 + ***************************************************************************/
  170 +inline
  171 +int GOptimizerLM::iter(void) const
  172 +{
  173 + return (m_iter);
  174 +}
  175 +
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Set maximum number of iterations
  179 + *
  180 + * @param[in] max_iter Maximum number of iterations.
  181 + ***************************************************************************/
  182 +inline
  183 +void GOptimizerLM::max_iter(const int& max_iter)
  184 +{
  185 + m_max_iter = max_iter;
  186 + return;
  187 +}
  188 +
  189 +
  190 +/***********************************************************************//**
  191 + * @brief Set maximum number of allowed subsequent stalls
  192 + *
  193 + * @param[in] max_stalls Maximum number of allowed subsequent stalls.
  194 + ***************************************************************************/
  195 +inline
  196 +void GOptimizerLM::max_stalls(const int& max_stalls)
  197 +{
  198 + m_max_stall = max_stalls;
  199 + return;
  200 +}
  201 +
  202 +
  203 +/***********************************************************************//**
  204 + * @brief Set maximum number of parameter boundary hits
  205 + *
  206 + * @param[in] max_hit Maximum number of parameter boundary hits.
  207 + ***************************************************************************/
  208 +inline
  209 +void GOptimizerLM::max_boundary_hits(const int& max_hit)
  210 +{
  211 + m_max_hit = max_hit;
  212 + return;
  213 +}
  214 +
  215 +
  216 +/***********************************************************************//**
  217 + * @brief Set lambda starting value
  218 + *
  219 + * @param[in] value Lambda starting value.
  220 + ***************************************************************************/
  221 +inline
  222 +void GOptimizerLM::lambda_start(const double& value)
  223 +{
  224 + m_lambda_start = value;
  225 + return;
  226 +}
  227 +
  228 +
  229 +/***********************************************************************//**
  230 + * @brief Set lambda increment value
  231 + *
  232 + * @param[in] value Lambda increment value.
  233 + ***************************************************************************/
  234 +inline
  235 +void GOptimizerLM::lambda_inc(const double& value)
  236 +{
  237 + m_lambda_inc = value;
  238 + return;
  239 +}
  240 +
  241 +
  242 +/***********************************************************************//**
  243 + * @brief Set lambda decrement value
  244 + *
  245 + * @param[in] value Lambda decrement value.
  246 + ***************************************************************************/
  247 +inline
  248 +void GOptimizerLM::lambda_dec(const double& value)
  249 +{
  250 + m_lambda_dec = value;
  251 + return;
  252 +}
  253 +
  254 +
  255 +/***********************************************************************//**
  256 + * @brief Set requested absolute convergence precision
  257 + *
  258 + * @param[in] eps Requested absolute convergence precision.
  259 + ***************************************************************************/
  260 +inline
  261 +void GOptimizerLM::eps(const double& eps)
  262 +{
  263 + m_eps = eps;
  264 + return;
  265 +}
  266 +
  267 +
  268 +/***********************************************************************//**
  269 + * @brief Set acceptable function decrease
  270 + *
  271 + * @param[in] value Acceptable function decrease.
  272 + *
  273 + * Sets the acceptable function decrease value for which the new solution
  274 + * will be kept and the iterations continue. This strategy provides better
  275 + * convergence in case that a function decrease is encountered.
  276 + ***************************************************************************/
  277 +inline
  278 +void GOptimizerLM::accept_dec(const double& value)
  279 +{
  280 + m_accept_dec = value;
  281 + return;
  282 +}
  283 +
  284 +
  285 +
  286 +
  287 +/***********************************************************************//**
  288 + * @brief Return maximum number of iterations
  289 + *
  290 + * @return Maximum number of iterations.
  291 + ***************************************************************************/
  292 +inline
  293 +int GOptimizerLM::max_iter(void) const
  294 +{
  295 + return (m_max_iter);
  296 +}
  297 +
  298 +
  299 +/***********************************************************************//**
  300 + * @brief Return maximum number of allowed subsequent stalls
  301 + *
  302 + * @return Maximum number of allowed subsequent stalls.
  303 + ***************************************************************************/
  304 +inline
  305 +int GOptimizerLM::max_stalls(void) const
  306 +{
  307 + return (m_max_stall);
  308 +}
  309 +
  310 +
  311 +/***********************************************************************//**
  312 + * @brief Return maximum number of parameter boundary hits
  313 + *
  314 + * @return Maximum number of parameter boundary hits.
  315 + ***************************************************************************/
  316 +inline
  317 +int GOptimizerLM::max_boundary_hits(void) const
  318 +{
  319 + return (m_max_hit);
  320 +}
  321 +
  322 +
  323 +/***********************************************************************//**
  324 + * @brief Return lambda starting value
  325 + *
  326 + * @return Lambda starting value.
  327 + ***************************************************************************/
  328 +inline
  329 +const double& GOptimizerLM::lambda_start(void) const
  330 +{
  331 + return (m_lambda_start);
  332 +}
  333 +
  334 +
  335 +/***********************************************************************//**
  336 + * @brief Return lambda increment value
  337 + *
  338 + * @return Lambda increment value.
  339 + ***************************************************************************/
  340 +inline
  341 +const double& GOptimizerLM::lambda_inc(void) const
  342 +{
  343 + return (m_lambda_inc);
  344 +}
  345 +
  346 +
  347 +/***********************************************************************//**
  348 + * @brief Return lambda decrement value
  349 + *
  350 + * @return Lambda decrement value.
  351 + ***************************************************************************/
  352 +inline
  353 +const double& GOptimizerLM::lambda_dec(void) const
  354 +{
  355 + return (m_lambda_dec);
  356 +}
  357 +
  358 +
  359 +/***********************************************************************//**
  360 + * @brief Return lambda value
  361 + *
  362 + * @return Lambda value.
  363 + ***************************************************************************/
  364 +inline
  365 +const double& GOptimizerLM::lambda(void) const
  366 +{
  367 + return (m_lambda);
  368 +}
  369 +
  370 +
  371 +/***********************************************************************//**
  372 + * @brief Return requested absolute convergence precision
  373 + *
  374 + * @return Requested absolute convergence precision.
  375 + ***************************************************************************/
  376 +inline
  377 +const double& GOptimizerLM::eps(void) const
  378 +{
  379 + return (m_eps);
  380 +}
  381 +
  382 +
  383 +/***********************************************************************//**
  384 + * @brief Return acceptable function decrease
  385 + *
  386 + * @return Acceptable function decrease.
  387 + ***************************************************************************/
  388 +inline
  389 +const double& GOptimizerLM::accept_dec(void) const
  390 +{
  391 + return (m_accept_dec);
  392 +}
  393 +
137 394 #endif /* GOPTIMIZERLM_HPP */
... ...
inst/mwl/test/test_MWL.cpp
  View file @9219f87
... ... @@ -219,13 +219,13 @@ void TestGMWL::test_optimizer(void)
219 219 double fit_results[] = {83.6331, 0,
220 220 22.0145, 0,
221 221 0.00231438040883405, 0.000191759607698647,
222   - -3.12006, 0.128382,
  222 + -3.11815, 0.128382,
223 223 1, 0,
224 224 1, 0,
225 225 83.6331, 0,
226 226 22.0145, 0,
227 227 4.5050091460345e-11, 1.95885803229878e-12,
228   - -2.13765, 0.00578641,
  228 + -2.13705, 0.00578641,
229 229 1000, 0,
230 230 1, 0};
231 231 test_try("Perform LM optimization");
... ...
pyext/GOptimizerLM.i
  View file @9219f87
... ... @@ -55,13 +55,14 @@ public:
55 55 virtual int iter(void) const;
56 56  
57 57 // Methods
58   - void max_iter(const int& n);
59   - void max_stalls(const int& n);
60   - void max_boundary_hits(const int& n);
61   - void lambda_start(const double& val);
62   - void lambda_inc(const double& val);
63   - void lambda_dec(const double& val);
  58 + void max_iter(const int& max_iter);
  59 + void max_stalls(const int& max_stalls);
  60 + void max_boundary_hits(const int& max_hit);
  61 + void lambda_start(const double& value);
  62 + void lambda_inc(const double& value);
  63 + void lambda_dec(const double& value);
64 64 void eps(const double& eps);
  65 + void accept_dec(const double& value);
65 66 int max_iter(void) const;
66 67 int max_stalls(void) const;
67 68 int max_boundary_hits(void) const;
... ... @@ -69,6 +70,7 @@ public:
69 70 const double& lambda_inc(void) const;
70 71 const double& lambda_dec(void) const;
71 72 const double& eps(void) const;
  73 + const double& accept_dec(void);
72 74 };
73 75  
74 76  
... ...
src/opt/GOptimizerLM.cpp
  View file @9219f87
... ... @@ -577,6 +577,8 @@ std::string GOptimizerLM::print(const GChatter&amp; chatter) const
577 577 result.append(gammalib::str(m_value, 3));
578 578 result.append("\n"+gammalib::parformat("Absolute precision"));
579 579 result.append(gammalib::str(m_eps));
  580 + result.append("\n"+gammalib::parformat("Acceptable value decrease"));
  581 + result.append(gammalib::str(std::abs(m_accept_dec)));
580 582  
581 583 // Append status
582 584 result.append("\n"+gammalib::parformat("Optimization status"));
... ... @@ -636,6 +638,7 @@ void GOptimizerLM::init_members(void)
636 638 m_lambda_inc = 10.0;
637 639 m_lambda_dec = 0.1;
638 640 m_eps = 5.0e-3; //!< Changed on 30/10/2014 from 1.0e-6
  641 + m_accept_dec = 1.0; //!< Allow to decrease by 1.0
639 642 m_max_iter = 100; //!< Changed on 30/10/2014 from 1000
640 643 m_max_stall = 10;
641 644 m_max_hit = 3; //!< Maximum successive boundary hits before freeze
... ... @@ -675,6 +678,7 @@ void GOptimizerLM::copy_members(const GOptimizerLM&amp; opt)
675 678 m_lambda_inc = opt.m_lambda_inc;
676 679 m_lambda_dec = opt.m_lambda_dec;
677 680 m_eps = opt.m_eps;
  681 + m_accept_dec = opt.m_accept_dec;
678 682 m_max_iter = opt.m_max_iter;
679 683 m_max_stall = opt.m_max_stall;
680 684 m_step_adjust = opt.m_step_adjust;
... ... @@ -948,9 +952,10 @@ void GOptimizerLM::iteration(GOptimizerFunction&amp; fct, GOptimizerPars&amp; pars)
948 952 }
949 953 }
950 954  
951   - // ... if function worsened slightly then accept new solution and
952   - // increase lambda
953   - else if (m_delta > -1.0e-6) {
  955 + // ... if function worsened by an amount smaller than m_accept_dec
  956 + // then accept new solution and increase lambda for the next
  957 + // iteration
  958 + else if (m_delta >= -std::abs(m_accept_dec)) {
954 959 m_lambda *= m_lambda_inc;
955 960 }
956 961  
... ...