From 2a3cf8732b2c9e16ebf4d79df82b203cebb8e810 Mon Sep 17 00:00:00 2001 From: Taylor Braun-Jones Date: Thu, 25 Feb 2016 12:47:49 -0500 Subject: [PATCH 1/5] Use FindThreads module to properly find pthreads for cross-builds --- dlib/CMakeLists.txt | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/dlib/CMakeLists.txt b/dlib/CMakeLists.txt index 5762ffd29..6d70f4ab0 100644 --- a/dlib/CMakeLists.txt +++ b/dlib/CMakeLists.txt @@ -132,6 +132,13 @@ if (NOT TARGET dlib) stack_trace.cpp ) + set(dlib_needed_libraries) + if(UNIX) + set(CMAKE_THREAD_PREFER_PTHREAD ON) + find_package(Threads REQUIRED) + set(dlib_needed_libraries ${dlib_needed_libraries} ${CMAKE_THREAD_LIBS_INIT}) + endif() + # we want to link to the right stuff depending on our platform. if (WIN32 AND NOT CYGWIN) ############################################################################### if (DLIB_NO_GUI_SUPPORT) @@ -140,9 +147,6 @@ if (NOT TARGET dlib) set (dlib_needed_libraries ws2_32 winmm comctl32 gdi32 imm32) endif() elseif(APPLE) ############################################################################ - find_library(pthreadlib pthread) - set (dlib_needed_libraries ${pthreadlib}) - if (NOT DLIB_NO_GUI_SUPPORT) find_package(X11 QUIET) if (X11_FOUND) @@ -178,9 +182,6 @@ if (NOT TARGET dlib) mark_as_advanced(pthreadlib xlib xlib_path x11_path) else () ################################################################################## - find_library(pthreadlib pthread) - set (dlib_needed_libraries ${pthreadlib}) - # link to the nsl library if it exists. this is something you need sometimes find_library(nsllib nsl) if (nsllib) From 378431eeeae500e48e363c34ca701db1d34c58ae Mon Sep 17 00:00:00 2001 From: Taylor Braun-Jones Date: Thu, 25 Feb 2016 12:51:13 -0500 Subject: [PATCH 2/5] Support building dlib from a top-level CMakeLists file This follows convention which simplifies using dlib as a CMake ExternalProject. --- CMakeLists.txt | 2 ++ dlib/CMakeLists.txt | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) create mode 100644 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 000000000..3f9c43e68 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,2 @@ +cmake_minimum_required(VERSION 2.8.4) +add_subdirectory(dlib) diff --git a/dlib/CMakeLists.txt b/dlib/CMakeLists.txt index 6d70f4ab0..6fe72b09e 100644 --- a/dlib/CMakeLists.txt +++ b/dlib/CMakeLists.txt @@ -432,7 +432,7 @@ if (NOT TARGET dlib) ARCHIVE DESTINATION lib) endif() - install(DIRECTORY ${CMAKE_SOURCE_DIR}/ DESTINATION include/dlib + install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/ DESTINATION include/dlib FILES_MATCHING PATTERN "*.h" REGEX "${CMAKE_CURRENT_BINARY_DIR}" EXCLUDE) From 9e9fc7409282dffc084136c3b575af9527e54596 Mon Sep 17 00:00:00 2001 From: Davis King Date: Sat, 5 Mar 2016 10:40:23 -0500 Subject: [PATCH 3/5] Added solve_qp_box_constrained() --- .../optimization_solve_qp_using_smo.h | 147 ++++++++++++++++++ ...optimization_solve_qp_using_smo_abstract.h | 49 ++++++ dlib/test/mpc.cpp | 22 +++ 3 files changed, 218 insertions(+) diff --git a/dlib/optimization/optimization_solve_qp_using_smo.h b/dlib/optimization/optimization_solve_qp_using_smo.h index fd15c5291..8a37bb6f4 100644 --- a/dlib/optimization/optimization_solve_qp_using_smo.h +++ b/dlib/optimization/optimization_solve_qp_using_smo.h @@ -404,6 +404,153 @@ namespace dlib return iter+1; } +// ---------------------------------------------------------------------------------------- + + template < + typename EXP1, + typename EXP2, + typename T, long NR, long NC, typename MM, typename L + > + unsigned long solve_qp_box_constrained ( + const matrix_exp& _Q, + const matrix_exp& _b, + matrix& alpha, + matrix& lower, + matrix& upper, + T eps, + unsigned long max_iter + ) + { + const_temp_matrix Q(_Q); + const_temp_matrix b(_b); + + // make sure requires clause is not broken + DLIB_ASSERT(Q.nr() == Q.nc() && + alpha.size() == lower.size() && + alpha.size() == upper.size() && + is_col_vector(b) && + is_col_vector(alpha) && + is_col_vector(lower) && + is_col_vector(upper) && + b.size() == alpha.size() && + b.size() == Q.nr() && + alpha.size() > 0 && + 0 <= min(alpha-lower) && + 0 <= max(upper-alpha) && + eps > 0 && + max_iter > 0, + "\t unsigned long solve_qp_box_constrained()" + << "\n\t Invalid arguments were given to this function" + << "\n\t Q.nr(): " << Q.nr() + << "\n\t Q.nc(): " << Q.nc() + << "\n\t is_col_vector(b): " << is_col_vector(b) + << "\n\t is_col_vector(alpha): " << is_col_vector(alpha) + << "\n\t is_col_vector(lower): " << is_col_vector(lower) + << "\n\t is_col_vector(upper): " << is_col_vector(upper) + << "\n\t b.size(): " << b.size() + << "\n\t alpha.size(): " << alpha.size() + << "\n\t lower.size(): " << lower.size() + << "\n\t upper.size(): " << upper.size() + << "\n\t Q.nr(): " << Q.nr() + << "\n\t min(alpha-lower): " << min(alpha-lower) + << "\n\t max(upper-alpha): " << max(upper-alpha) + << "\n\t eps: " << eps + << "\n\t max_iter: " << max_iter + ); + + + // Compute f'(alpha) (i.e. the gradient of f(alpha)) for the current alpha. + matrix df = Q*alpha + b; + matrix QQ = reciprocal_max(diag(Q)); + + // First we use a coordinate descent method to initialize alpha. + double max_df = 0; + for (unsigned long iter = 0; iter < alpha.size()*2; ++iter) + { + max_df = 0; + long best_r =0; + // find the best alpha to optimize. + for (long r = 0; r < Q.nr(); ++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) + { + best_r = r; + max_df = std::abs(df(r)); + } + } + + // now optimize alpha(best_r) + const long r = best_r; + const T old_alpha = alpha(r); + alpha(r) = -(df(r)-Q(r,r)*alpha(r))*QQ(r); + if (alpha(r) < lower(r)) + alpha(r) = lower(r); + else if (alpha(r) > upper(r)) + alpha(r) = upper(r); + + const T delta = old_alpha-alpha(r); + + // Now update the gradient. We will perform the equivalent of: df = Q*alpha + b; + for(long k = 0; k < df.nr(); ++k) + df(k) -= Q(r,k)*delta; + } + //cout << "max_df: " << max_df << endl; + //cout << "objective value: " << 0.5*trans(alpha)*Q*alpha + trans(b)*alpha << endl; + + + + // Now do the main iteration block of this solver. The coordinate descent method + // we used above can improve the objective rapidly in the beginning. However, + // Nesterov's method has more rapid convergence once it gets going so this is what + // we use for the main iteration. + matrix v, v_old; + v = alpha; + // We need to get an upper bound on the Lipschitz constant for this QP. Since that + // is just the max eigenvalue of Q we can do it using Gershgorin disks. + const T lipschitz_bound = max(diag(Q) + (sum_cols(abs(Q)) - abs(diag(Q)))); + double lambda = 0; + unsigned long iter; + for (iter = 0; iter < max_iter; ++iter) + { + const double next_lambda = (1 + std::sqrt(1+4*lambda*lambda))/2; + const double gamma = (1-lambda)/next_lambda; + lambda = next_lambda; + + 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 = clamp((1-gamma)*v + gamma*v_old, lower, upper); + + + // check for convergence every 10 iterations + if (iter%10 == 0) + { + max_df = 0; + for (long r = 0; r < Q.nr(); ++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)); + } + if (max_df < eps) + break; + } + } + + //cout << "max_df: " << max_df << endl; + //cout << "objective value: " << 0.5*trans(alpha)*Q*alpha + trans(b)*alpha << endl; + return iter+1; + } + // ---------------------------------------------------------------------------------------- } diff --git a/dlib/optimization/optimization_solve_qp_using_smo_abstract.h b/dlib/optimization/optimization_solve_qp_using_smo_abstract.h index ad3c95baa..212ca612e 100644 --- a/dlib/optimization/optimization_solve_qp_using_smo_abstract.h +++ b/dlib/optimization/optimization_solve_qp_using_smo_abstract.h @@ -113,6 +113,55 @@ namespace dlib converge to eps accuracy then the number returned will be max_iter+1. !*/ +// ---------------------------------------------------------------------------------------- + + template < + typename EXP1, + typename EXP2, + typename T, long NR, long NC, typename MM, typename L + > + unsigned long solve_qp_box_constrained ( + const matrix_exp& Q, + const matrix_exp& b, + matrix& alpha, + matrix& lower, + matrix& upper, + T eps, + unsigned long max_iter + ); + /*! + requires + - Q.nr() == Q.nc() + - alpha.size() == lower.size() == upper.size() + - is_col_vector(b) == true + - is_col_vector(alpha) == true + - is_col_vector(lower) == true + - is_col_vector(upper) == true + - b.size() == alpha.size() == Q.nr() + - alpha.size() > 0 + - 0 <= min(alpha-lower) + - 0 <= max(upper-alpha) + - eps > 0 + - max_iter > 0 + ensures + - This function solves the following quadratic program: + Minimize: f(alpha) == 0.5*trans(alpha)*Q*alpha + trans(b)*alpha + subject to the following box constraints on alpha: + - 0 <= min(alpha-lower) + - 0 <= max(upper-alpha) + Where f is convex. This means that Q should be positive-semidefinite. + - 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. + - 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. + !*/ + // ---------------------------------------------------------------------------------------- } diff --git a/dlib/test/mpc.cpp b/dlib/test/mpc.cpp index 36307be59..723a40bd2 100644 --- a/dlib/test/mpc.cpp +++ b/dlib/test/mpc.cpp @@ -6,6 +6,7 @@ #include #include +#include #include "tester.h" namespace @@ -314,6 +315,27 @@ namespace initial_state = A*initial_state + B*control + C; //cout << control(0) << "\t" << trans(initial_state); } + + { + // also just generally test our QP solver. + matrix Q = gaussian_randm(20,20,5); + Q = Q*trans(Q); + + matrix b = randm(20,1)-0.5; + matrix alpha, lower, upper, alpha2; + alpha = 0; + alpha2 = 0; + lower = -4; + upper = 3; + + solve_qp_box_using_smo(Q,b,alpha,lower, upper, 0.00000001, 100000); + solve_qp_box_constrained(Q,b,alpha2,lower, upper, 0.00000001, 10000); + dlog << LINFO << trans(alpha); + dlog << LINFO << trans(alpha2); + dlog << LINFO << "objective value: " << 0.5*trans(alpha)*Q*alpha + trans(b)*alpha; + dlog << LINFO << "objective value2: " << 0.5*trans(alpha2)*Q*alpha + trans(b)*alpha2; + DLIB_TEST_MSG(max(abs(alpha-alpha2)) < 1e-7, max(abs(alpha-alpha2))); + } } } a; From 0af762c5e803b263f99c54bf461d5f3c299be1da Mon Sep 17 00:00:00 2001 From: Davis King Date: Sat, 5 Mar 2016 10:46:15 -0500 Subject: [PATCH 4/5] updated docs --- docs/docs/optimization.xml | 21 +++++++++++++++++++++ docs/docs/term_index.xml | 1 + 2 files changed, 22 insertions(+) diff --git a/docs/docs/optimization.xml b/docs/docs/optimization.xml index cd71818fd..8c18cc25d 100644 --- a/docs/docs/optimization.xml +++ b/docs/docs/optimization.xml @@ -41,6 +41,7 @@
Special Purpose Optimizers + solve_qp_box_constrained solve_qp_using_smo solve_qp2_using_smo solve_qp3_using_smo @@ -429,6 +430,26 @@ subject to the following constraint: + + + + solve_qp_box_constrained + dlib/optimization.h + dlib/optimization/optimization_solve_qp_using_smo_abstract.h + + This function solves the following quadratic program: +
+   Minimize: f(alpha) == 0.5*trans(alpha)*Q*alpha + trans(b)*alpha 
+   subject to the following box constraints on alpha:
+      0 <= min(alpha-lower)
+      0 <= max(upper-alpha)
+   Where f is convex.  This means that Q should be positive-semidefinite.
+
+ + + + + diff --git a/docs/docs/term_index.xml b/docs/docs/term_index.xml index e3a70124e..5c4943a29 100644 --- a/docs/docs/term_index.xml +++ b/docs/docs/term_index.xml @@ -153,6 +153,7 @@ + From 6e6f37e512eab374478f9ffc45fef85084e0110a Mon Sep 17 00:00:00 2001 From: Davis King Date: Mon, 7 Mar 2016 07:05:25 -0500 Subject: [PATCH 5/5] Fixed some warnings and made the unit tests more robust. --- dlib/optimization/optimization_solve_qp_using_smo.h | 2 +- dlib/test/mpc.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dlib/optimization/optimization_solve_qp_using_smo.h b/dlib/optimization/optimization_solve_qp_using_smo.h index 8a37bb6f4..f20a56c6c 100644 --- a/dlib/optimization/optimization_solve_qp_using_smo.h +++ b/dlib/optimization/optimization_solve_qp_using_smo.h @@ -465,7 +465,7 @@ namespace dlib // First we use a coordinate descent method to initialize alpha. double max_df = 0; - for (unsigned long iter = 0; iter < alpha.size()*2; ++iter) + for (long iter = 0; iter < alpha.size()*2; ++iter) { max_df = 0; long best_r =0; diff --git a/dlib/test/mpc.cpp b/dlib/test/mpc.cpp index 723a40bd2..c0a98dd17 100644 --- a/dlib/test/mpc.cpp +++ b/dlib/test/mpc.cpp @@ -328,8 +328,8 @@ namespace lower = -4; upper = 3; - solve_qp_box_using_smo(Q,b,alpha,lower, upper, 0.00000001, 100000); - solve_qp_box_constrained(Q,b,alpha2,lower, upper, 0.00000001, 10000); + 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); dlog << LINFO << trans(alpha); dlog << LINFO << trans(alpha2); dlog << LINFO << "objective value: " << 0.5*trans(alpha)*Q*alpha + trans(b)*alpha;