Added Powell's BOBYQA algorithm for optimization without derivatives.

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403214

dlib/optimization.h

@@ -4,6 +4,7 @@
 #define DLIB_SVM_HEADER
 
 #include "optimization/optimization.h"
+#include "optimization/optimization_bobyqa.h"
 
 #endif // DLIB_OPTIMIZATIOn_HEADER
 

dlib/optimization/optimization_bobyqa_abstract.h

+// Copyright (C) 2009  Davis E. King (davisking@users.sourceforge.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#undef DLIB_OPTIMIZATIOn_BOBYQA_ABSTRACT_H__
+#ifdef DLIB_OPTIMIZATIOn_BOBYQA_ABSTRACT_H__
+
+#include "dlib/matrix.h"
+
+// ----------------------------------------------------------------------------------------
+
+/*
+    This file defines the dlib interface to the BOBYQA software developed by M.J.D Powell.
+    BOBYQA is a method for optimizing a function in the absence of derivative information.  
+    Powell described it as a method that seeks the least value of a function of many 
+    variables, by applying a trust region method that forms quadratic models by 
+    interpolation.  There is usually some freedom in the interpolation conditions, 
+    which is taken up by minimizing the Frobenius norm of the change to the second 
+    derivative of the model, beginning with the zero matrix. The values of the variables 
+    are constrained by upper and lower bounds.  
+
+
+    The following paper, published in 2009 by Powell, describes the
+    detailed working of the BOBYQA algorithm.  
+
+        The BOBYQA algorithm for bound constrained optimization 
+        without derivatives by M.J.D. Powell
+*/
+
+// ----------------------------------------------------------------------------------------
+
+namespace dlib
+{
+    class bobyqa_failure : public error;
+    /*!
+        This is the exception class used by the functions defined in this file.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename funct,
+        typename T, 
+        typename U
+        >
+    void find_min_bobyqa (
+        const funct& f,
+        T& x,
+        long npt,
+        const U& x_lower,
+        const U& x_upper,
+        const double rho_begin,
+        const double rho_end,
+        const long max_f_evals
+    );
+    /*!
+        requires
+            - f(x) must be a valid expression that evaluates to a double
+            - is_col_vector(x) == true
+            - is_col_vector(x_lower) == true
+            - is_col_vector(x_upper) == true
+            - x.size() == x_lower.size() == x_upper.size()
+            - x.size() > 1
+            - x.size() + 2 <= npt <= (x.size()+1)*(x.size()+2)/2
+            - 0 < rho_end < rho_begin
+            - min(x_upper - x_lower) > 2*rho_begin
+              (i.e. the lower and upper bounds on each x element must be larger than 2*rho_begin)
+            - min(x - x_lower) >= 0 && min(x_upper - x) >= 0
+              (i.e. the given x should be within the bounds defined by x_lower and x_upper)
+            - max_f_evals > 1
+        ensures
+            - Performs a constrained minimization of the function f() starting from 
+              the initial point x.  
+            - The BOBYQA algorithm uses a number of interpolating points to perform its
+              work.  The npt argument controls how many points get used.  Typically,
+              a good value to use is 2*x.size()+1.
+            - #x == the value of x (within the bounds defined by x_lower and x_upper) that 
+              was found to minimize f().  More precisely:
+                - min(#x - x_lower) >= 0 && min(x_upper - #x) >= 0
+            - rho_begin and rho_end are used as the initial and final values of a trust 
+              region radius.  Typically, rho_begin should be about one tenth of the greatest 
+              expected change to a variable, while rho_end should indicate the accuracy that 
+              is required in the final values of the variables. 
+        throws
+            - bobyqa_failure
+                This exception is thrown if the algorithm is unable to make progress towards
+                solving the problem.  This may occur because the algorithm detects excessive
+                numerical errors or because max_f_evals of f() have occurred without reaching
+                convergence.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename funct,
+        typename T, 
+        typename U
+        >
+    void find_max_bobyqa (
+        const funct& f,
+        T& x,
+        long npt,
+        const U& x_lower,
+        const U& x_upper,
+        const double rho_begin,
+        const double rho_end,
+        const long max_f_evals
+    );
+    /*!
+        This function is identical to the find_min_bobyqa() routine defined above
+        except that it negates the f() function before performing optimization.  
+        Thus this function will attempt to find the maximizer of f() rather than 
+        the minimizer.
+    !*/
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_OPTIMIZATIOn_BOBYQA_ABSTRACT_H__
+

dlib/test/optimization.cpp

@@ -169,6 +169,20 @@ namespace
         }
 
 
+        if (p.size() < 100)
+        {
+            total_count = 0;
+            x = p;
+            find_min_bobyqa(wrap_function(apq<T>), x, 2*x.size()+1,
+                            uniform_matrix<double>(x.size(),1,-1e100),
+                            uniform_matrix<double>(x.size(),1,1e100),
+                            (max(abs(x))+1)/10,
+                            1e-6,
+                            10000);
+            DLIB_TEST_MSG(dlib::equal(x,opt, 1e-5),opt-x);
+            dlog << LINFO << "find_min_bobyqa(): got apq in " << total_count;
+        }
+
         total_count = 0;
         x = p;
         find_min(lbfgs_search_strategy(10), 
@@ -312,6 +326,19 @@ namespace
                                                &powell, x, minf, 1e-10);
         DLIB_TEST_MSG(dlib::equal(x,opt, 1e-1),opt-x);
         dlog << LINFO << "find_min() cg: got powell/noder2 in " << total_count;
+
+
+        total_count = 0;
+        x = p;
+        find_min_bobyqa(&powell, x, 2*x.size()+1,
+                        uniform_matrix<double>(x.size(),1,-1e100),
+                        uniform_matrix<double>(x.size(),1,1e100),
+                        (max(abs(x))+1)/10,
+                        1e-7,
+                        10000);
+        DLIB_TEST_MSG(dlib::equal(x,opt, 1e-3),opt-x);
+        dlog << LINFO << "find_min_bobyqa(): got powell in " << total_count;
+
     }
 
 
@@ -419,6 +446,18 @@ namespace
         DLIB_TEST_MSG(dlib::equal(x,opt, 1e-5),opt-x);
         dlog << LINFO << "find_min() cg: got simple/noder2 in " << total_count;
 
+
+        total_count = 0;
+        x = p;
+        find_min_bobyqa(&simple, x, 2*x.size()+1,
+                        uniform_matrix<double>(x.size(),1,-1e100),
+                        uniform_matrix<double>(x.size(),1,1e100),
+                        (max(abs(x))+1)/10,
+                        1e-6,
+                        10000);
+        DLIB_TEST_MSG(dlib::equal(x,opt, 1e-5),opt-x);
+        dlog << LINFO << "find_min_bobyqa(): got simple in " << total_count;
+
     }
 
 
@@ -514,6 +553,21 @@ namespace
                                                &rosen, x, minf);
         DLIB_TEST_MSG(dlib::equal(x,opt, 1e-4),opt-x);
         dlog << LINFO << "find_min() cg: got rosen/noder2 in " << total_count;
+
+
+        if (max(abs(p)) < 1000)
+        {
+            total_count = 0;
+            x = p;
+            find_min_bobyqa(&rosen, x, 2*x.size()+1,
+                            uniform_matrix<double>(x.size(),1,-1e100),
+                            uniform_matrix<double>(x.size(),1,1e100),
+                            (max(abs(x))+1)/10,
+                            1e-6,
+                            10000);
+            DLIB_TEST_MSG(dlib::equal(x,opt, 1e-5),opt-x);
+            dlog << LINFO << "find_min_bobyqa(): got rosen in " << total_count;
+        }
     }
 
 
@@ -561,6 +615,18 @@ namespace
             objective_delta_stop_strategy(eps), &neg_rosen, x, maxf);
         DLIB_TEST_MSG(dlib::equal(x,opt, 1e-7),opt-x);
         dlog << LINFO << "find_max() cg: got neg_rosen/noder2 in " << total_count;
+
+
+        total_count = 0;
+        x = p;
+        find_max_bobyqa(&neg_rosen, x, 2*x.size()+1,
+                        uniform_matrix<double>(x.size(),1,-1e100),
+                        uniform_matrix<double>(x.size(),1,1e100),
+                        (max(abs(x))+1)/10,
+                        1e-6,
+                        10000);
+        DLIB_TEST_MSG(dlib::equal(x,opt, 1e-5),opt-x);
+        dlog << LINFO << "find_max_bobyqa(): got neg_rosen in " << total_count;
     }
 
 // ----------------------------------------------------------------------------------------