Improved the stopping condition of solve_qp_box_constrained() and...

Improved the stopping condition of solve_qp_box_constrained() and solve_qp_box_constrained_blockdiag().

dlib/optimization/optimization_solve_qp_using_smo.h

@@ -419,8 +419,8 @@ namespace dlib
         matrix<T,NR,NC,MM,L>& alpha,
         const matrix<T,NR,NC,MM,L>& lower,
         const matrix<T,NR,NC,MM,L>& upper,
-        T eps,
-        unsigned long max_iter
+        T eps = 1e-10,
+        unsigned long max_iter = 30000
     )
     {
         // make sure requires clause is not broken
@@ -521,26 +521,35 @@ namespace dlib
 
             v_old = v;
 
-            df = Q*alpha + b;
             // now take a projected gradient step using Nesterov's method.
             v = clamp(alpha - 1.0/lipschitz_bound * df, lower, upper);
             alpha = dlib::clamp((1-gamma)*v + gamma*v_old, lower, upper);
 
+            df = Q*alpha + b;
 
             // check for convergence every 10 iterations
             if (iter%10 == 0)
             {
                 max_df = 0;
+                double absalpha = 0;
+                double thealpha = 0;
                 for (long r = 0; r < Q.nr(); ++r)
                 {
+                    absalpha += std::abs(alpha(r));
                     if (alpha(r) <= lower(r) && df(r) > 0)
                         ;//alpha(r) = lower(r);
                     else if (alpha(r) >= upper(r) && df(r) < 0)
                         ;//alpha(r) = upper(r);
                     else if (std::abs(df(r)) > max_df)
+                    {
                         max_df = std::abs(df(r));
+                        thealpha = std::abs(alpha(r));
+                    }
                 }
-                if (max_df < eps)
+                absalpha /= Q.nr();
+                // Stop when the magnitude of the changes we are making to alpha are eps
+                // smaller than the typical alpha.
+                if (max_df/lipschitz_bound <= eps*std::min(thealpha,absalpha))
                     break;
             }
         }
@@ -613,8 +622,8 @@ namespace dlib
         std::vector<matrix<T,NR,NC,MM,L>>& alphas,
         const std::vector<matrix<T,NR,NC,MM,L>>& lowers,
         const std::vector<matrix<T,NR,NC,MM,L>>& uppers,
-        T eps,
-        unsigned long max_iter
+        T eps = 1e-10,
+        unsigned long max_iter = 30000
     )
     {
         // make sure requires clause is not broken
@@ -834,8 +843,6 @@ namespace dlib
 
             v_old.swap(v);
 
-            //df = Q*alpha + b;
-            compute_df();
 
             // now take a projected gradient step using Nesterov's method.
             for (size_t j = 0; j < alphas.size(); ++j)
@@ -844,11 +851,15 @@ namespace dlib
                 alphas[j] = clamp((1-gamma)*v[j] + gamma*v_old[j], lowers[j], uppers[j]);
             }
 
+            //df = Q*alpha + b;
+            compute_df();
 
             // check for convergence every 10 iterations
             if (iter%10 == 0)
             {
                 max_df = 0;
+                double absalpha = 0;
+                double thealpha = 0;
                 for (size_t r2 = 0; r2 < alphas.size(); ++r2) 
                 {
                     auto& alpha = alphas[r2];
@@ -857,15 +868,22 @@ namespace dlib
                     auto& upper = uppers[r2];
                     for (long r = 0; r < alpha.nr(); ++r)
                     {
+                        absalpha += std::abs(alpha(r));
                         if (alpha(r) <= lower(r) && df_(r) > 0)
                             ;//alpha(r) = lower(r);
                         else if (alpha(r) >= upper(r) && df_(r) < 0)
                             ;//alpha(r) = upper(r);
                         else if (std::abs(df_(r)) > max_df)
+                        {
                             max_df = std::abs(df_(r));
+                            thealpha = std::abs(alpha(r));
+                        }
                     }
                 }
-                if (max_df < eps)
+                absalpha /= num_variables;
+                // Stop when the magnitude of the changes we are making to alpha are eps
+                // smaller than the typical alpha.
+                if (max_df/lipschitz_bound <= eps*std::min(thealpha,absalpha))
                     break;
             }
         }

dlib/optimization/optimization_solve_qp_using_smo_abstract.h

@@ -128,8 +128,8 @@ namespace dlib
         matrix<T,NR,NC,MM,L>& alpha,
         const matrix<T,NR,NC,MM,L>& lower,
         const matrix<T,NR,NC,MM,L>& upper,
-        T eps,
-        unsigned long max_iter
+        T eps = 1e-10,
+        unsigned long max_iter = 30000
     );
     /*!
         requires
@@ -155,10 +155,10 @@ namespace dlib
             - The solution to the above QP will be stored in #alpha.
             - This function uses a combination of a SMO algorithm along with Nesterov's
               method as the main iteration of the solver.  It starts the algorithm with the
-              given alpha and it works on the problem until the derivative of f(alpha) is
-              smaller than eps for each element of alpha or the alpha value is at a box
-              constraint.  So eps controls how accurate the solution is and smaller values
-              result in better solutions.
+              given alpha and works on the problem until the magnitude of the changes we
+              are making to alpha are eps times smaller than the typical values in alpha.
+              So eps controls how accurate the solution is and smaller values result in
+              better solutions.
             - At most max_iter iterations of optimization will be performed.  
             - returns the number of iterations performed.  If this method fails to
               converge to eps accuracy then the number returned will be max_iter+1.
@@ -176,8 +176,8 @@ namespace dlib
         std::vector<matrix<T,NR,NC,MM,L>>& alphas,
         const std::vector<matrix<T,NR,NC,MM,L>>& lowers,
         const std::vector<matrix<T,NR,NC,MM,L>>& uppers,
-        T eps,
-        unsigned long max_iter
+        T eps = 1e-10,
+        unsigned long max_iter = 30000
     );
     /*!
         requires
@@ -224,10 +224,10 @@ namespace dlib
             - The solution to the above QP will be stored in #alphas.
             - This function uses a combination of a SMO algorithm along with Nesterov's
               method as the main iteration of the solver.  It starts the algorithm with the
-              given alpha and it works on the problem until the derivative of f(alpha) is
-              smaller than eps for each element of alpha or the alpha value is at a box
-              constraint.  So eps controls how accurate the solution is and smaller values
-              result in better solutions.
+              given alpha and works on the problem until the magnitude of the changes we
+              are making to alpha are eps times smaller than the typical values in alpha.
+              So eps controls how accurate the solution is and smaller values result in
+              better solutions.
             - At most max_iter iterations of optimization will be performed.  
             - returns the number of iterations performed.  If this method fails to
               converge to eps accuracy then the number returned will be max_iter+1.

dlib/test/mpc.cpp

@@ -328,8 +328,8 @@ namespace
                 lower = -4;
                 upper = 3;
 
-                solve_qp_box_using_smo(Q,b,alpha,lower, upper, 0.000000001, 500000);
-                solve_qp_box_constrained(Q,b,alpha2,lower, upper, 0.000000001, 50000);
+                solve_qp_box_using_smo(Q,b,alpha,lower, upper, 1e-12, 500000);
+                solve_qp_box_constrained(Q,b,alpha2,lower, upper, 1e-12, 50000);
                 dlog << LINFO << trans(alpha);
                 dlog << LINFO << trans(alpha2);
                 dlog << LINFO << "objective value:  " << 0.5*trans(alpha)*Q*alpha + trans(b)*alpha;

dlib/test/opt_qp_solver.cpp

@@ -590,6 +590,44 @@ namespace
         }
     }
 
+// ----------------------------------------------------------------------------------------
+
+    void test_solve_qp_box_constrained()
+    {
+        dlog << LINFO << "test_solve_qp_box_constrained()";
+        print_spinner();
+        dlib::rand rnd;
+        matrix<double> m = randm(6,6,rnd);
+        m = m*trans(m);
+        matrix<double,0,1> b = randm(6,1,rnd)-0.5;
+
+        matrix<double,0,1> lower, upper, solution;
+        lower = -ones_matrix<double>(6,1)*1e100;
+        upper = ones_matrix<double>(6,1)*1e100;
+
+        solution = zeros_matrix(lower);
+        unsigned long iters = solve_qp_box_constrained(m,b,solution, lower, upper);
+
+        dlog << LINFO << "iters: " << iters;
+
+        matrix<double> true_solution = -pinv(m)*b;
+        DLIB_TEST_MSG(max(abs(solution - true_solution)) < 1e-6, max(abs(solution - true_solution)));
+
+        iters = solve_qp_box_constrained(m,b,solution, lower, upper, 1e-14);
+        dlog << LINFO << "iters: " << iters;
+
+        /*
+        const double obj1 = 0.5*trans(solution)*m*solution + trans(solution)*b;
+        const double obj2 = 0.5*trans(true_solution)*m*true_solution + trans(true_solution)*b;
+        cout << "iters:" << iters << endl;
+        cout << "obj1: "<< obj1 << endl;
+        cout << "obj2: "<< obj2 << endl;
+        cout << "obj1-obj2: "<< obj1-obj2 << endl;
+        */
+
+        DLIB_TEST_MSG(max(abs(solution - true_solution)) < 1e-10, max(abs(solution - true_solution)));
+    }
+
 // ----------------------------------------------------------------------------------------
 
     void test_solve_qp_box_constrained_blockdiag_compact(dlib::rand& rnd, double percent_off_diag_present)
@@ -706,6 +744,7 @@ namespace
         void perform_test(
         )
         {
+            test_solve_qp_box_constrained();
             print_spinner();
             test_solve_qp4_using_smo();
             print_spinner();