From 88778d2c0e8b1b940a5b099140479d4c2a71e4a7 Mon Sep 17 00:00:00 2001 From: Davis King Date: Wed, 15 Feb 2017 17:06:03 -0500 Subject: [PATCH] Added find_gap_between_convex_hulls() --- .../optimization_solve_qp_using_smo.h | 56 +++++++ ...optimization_solve_qp_using_smo_abstract.h | 41 +++++ dlib/test/opt_qp_solver.cpp | 154 ++++++++++++++++++ 3 files changed, 251 insertions(+) diff --git a/dlib/optimization/optimization_solve_qp_using_smo.h b/dlib/optimization/optimization_solve_qp_using_smo.h index f20a56c6c..2da3208b8 100644 --- a/dlib/optimization/optimization_solve_qp_using_smo.h +++ b/dlib/optimization/optimization_solve_qp_using_smo.h @@ -551,6 +551,62 @@ namespace dlib return iter+1; } +// ---------------------------------------------------------------------------------------- + + template < + typename EXP1, + typename EXP2, + typename T, long NRa, long NRb + > + unsigned long find_gap_between_convex_hulls ( + const matrix_exp& A, + const matrix_exp& B, + matrix& cA, + matrix& cB, + const double eps, + const unsigned long max_iter = 1000 + ) + { + DLIB_CASSERT(A.size() != 0); + DLIB_CASSERT(B.size() != 0); + DLIB_CASSERT(A.nr() == B.nr(), "The dimensionality of the points in both convex hull sets must match"); + DLIB_CASSERT(eps > 0); + DLIB_CASSERT(max_iter > 0); + + cA.set_size(A.nc()); + cB.set_size(B.nc()); + + // initialize to the centroids of A and B respectively. + cA = 1.0/cA.size(); + cB = 1.0/cB.size(); + + + matrix AA, BB, AB, ABb, ABa; + + AA = trans(A)*A; + BB = trans(B)*B; + AB = trans(A)*B; + + unsigned long iter = 0; + for (iter = 0; iter < max_iter; ++iter) + { + // find the convex combination of A that is nearest to B*cB + ABb = AB*cB; + const auto smo_iter1 = solve_qp_using_smo(AA, ABb, cA, eps, cA.size()); + + // now find the convex combination of B that is nearest to A*cA + ABa = trans(AB)*cA; + const auto smo_iter2 = solve_qp_using_smo(BB, ABa, cB, eps, cB.size()); + + // stop if the QP solvers failed to improve + if (smo_iter1 == 1 && smo_iter2 == 1) + break; + } + + + 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 212ca612e..2e8578530 100644 --- a/dlib/optimization/optimization_solve_qp_using_smo_abstract.h +++ b/dlib/optimization/optimization_solve_qp_using_smo_abstract.h @@ -162,6 +162,47 @@ namespace dlib converge to eps accuracy then the number returned will be max_iter+1. !*/ +// ---------------------------------------------------------------------------------------- + + template < + typename EXP1, + typename EXP2, + typename T, long NRa, long NRb + > + unsigned long find_gap_between_convex_hulls ( + const matrix_exp& A, + const matrix_exp& B, + matrix& cA, + matrix& cB, + const double eps, + const unsigned long max_iter = 1000 + ); + /*! + requires + - A.nr() == B.nr() + - A.size() != 0 + - B.size() != 0 + - eps > 0 + - max_iter > 0 + ensures + - If you think of A and B as sets of column vectors, then we can identify the + convex sets hullA and hullB, which are the convex hulls of A and B + respectively. This function finds the pair of points in hullA and hullB that + are nearest to each other. To be precise, this function solves the following + quadratic program: + Minimize: f(cA,cB) == length_squared(A*cA - B*cB) + subject to the following constraints on cA and cB: + - is_col_vector(cA) == true && cA.size() == A.nc() + - is_col_vector(cB) == true && cB.size() == B.nc() + - sum(cA) == 1 && min(cA) >= 0 + - sum(cB) == 1 && min(cB) >= 0 + - This function uses an iterative block coordinate descent algorithm to solve + the QP. It runs until either max_iter iterations have been performed or the + QP is solved to at least eps accuracy. + - 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/opt_qp_solver.cpp b/dlib/test/opt_qp_solver.cpp index deeabe45c..92bfeda60 100644 --- a/dlib/test/opt_qp_solver.cpp +++ b/dlib/test/opt_qp_solver.cpp @@ -356,6 +356,157 @@ namespace test_qp4_test7(); } +// ---------------------------------------------------------------------------------------- + + double max_distance_to( + const std::vector>& a, + const std::vector>& b + ) + { + double best_dist = 0; + for (auto&& aa : a) + { + for (auto&& bb : b) + { + double dist = length(aa-bb); + if (dist > best_dist) + best_dist = dist; + } + } + return best_dist; + } + + double min_distance_to( + const std::vector>& a, + const std::vector>& b + ) + { + double best_dist = std::numeric_limits::infinity(); + for (auto&& aa : a) + { + for (auto&& bb : b) + { + double dist = length(aa-bb); + if (dist < best_dist) + best_dist = dist; + } + } + return best_dist; + } + + double min_distance_to( + const std::vector>& s, + const matrix& v + ) + { + double best_dist = std::numeric_limits::infinity(); + for (auto& x : s) + { + double dist = length(v-x); + if (dist < best_dist) + { + best_dist = dist; + } + } + return best_dist; + } + + double max_distance_to( + const std::vector>& s, + const matrix& v + ) + { + double best_dist = 0; + for (auto& x : s) + { + double dist = length(v-x); + if (dist > best_dist) + { + best_dist = dist; + } + } + return best_dist; + } + + void test_find_gap_between_convex_hulls() + { + print_spinner(); + std::vector> set1, set2; + + const double dist_thresh = 5.47723; + + // generate two groups of points that are pairwise close within each set and + // pairwise far apart between each set, according to dist_thresh distance threshold. + bool which = true; + for (size_t i = 0; i < 10000; ++i) + { + matrix v = gaussian_randm(15,2,i); + const auto min_dist1 = min_distance_to(set1,v); + const auto min_dist2 = min_distance_to(set2,v); + const auto max_dist1 = max_distance_to(set1,v); + const auto max_dist2 = max_distance_to(set2,v); + if (which) + { + if ((set1.size()==0 || max_dist1 < dist_thresh) && min_dist2 > dist_thresh ) + { + set1.push_back(v); + which = !which; + } + } + else + { + if ((set2.size()==0 || max_dist2 < dist_thresh) && min_dist1 > dist_thresh) + { + set2.push_back(v); + which = !which; + } + } + } + + dlog << LINFO << "set1.size(): "<< set1.size(); + dlog << LINFO << "set2.size(): "<< set2.size(); + + + // make sure we generated the points correctly. + dlog << LINFO << "dist_thresh: "<< dist_thresh; + dlog << LINFO << "max distance between set1 and set1: "<< max_distance_to(set1,set1); + dlog << LINFO << "max distance between set2 and set2: "<< max_distance_to(set2,set2); + DLIB_TEST(max_distance_to(set1,set1) < dist_thresh); + DLIB_TEST(max_distance_to(set2,set2) < dist_thresh); + dlog << LINFO << "min distance between set2 and set1: "<< min_distance_to(set2,set1); + DLIB_TEST(min_distance_to(set2,set1) > dist_thresh); + + + // It is slightly counterintuitive but true that points picked using the above procedure + // will have elements of their convex hulls that are much closer together than + // dist_thresh, even though none of the vertices of the hulls are that close + // together. This is especially true in high dimensions. So let's use this to + // test find_gap_between_convex_hulls(). It should be able to find a pair of + // points in the convex hulls of our sets that are a lot closer together than + // dist_thresh. + + // First we need to convert the vectors to matrices. + matrix A, B; + A.set_size(set1[0].size(), set1.size()); + B.set_size(set2[0].size(), set2.size()); + for (long c = 0; c < A.nc(); ++c) + set_colm(A,c) = set1[c]; + for (long c = 0; c < B.nc(); ++c) + set_colm(B,c) = set2[c]; + + matrix c1, c2; + find_gap_between_convex_hulls(A, B, c1, c2, 0.0001); + // make sure c1 and c2 are convex combinations. + DLIB_TEST(abs(sum(c1)-1) < 1e-8); + DLIB_TEST(abs(sum(c2)-1) < 1e-8); + DLIB_TEST(min(c1) >= 0); + DLIB_TEST(min(c2) >= 0); + + // now test that the points found are close together. + dlog << LINFO << "dist: "<< length(A*c1 - B*c2); + DLIB_TEST(length(A*c1 - B*c2) < 4); + } + // ---------------------------------------------------------------------------------------- class opt_qp_solver_tester : public tester @@ -412,6 +563,9 @@ namespace dlog << LINFO << "disagreement stddev: " << rs.stddev(); DLIB_TEST_MSG(rs.mean() < 0.001, rs.mean()); DLIB_TEST_MSG(rs.stddev() < 0.001, rs.stddev()); + + + test_find_gap_between_convex_hulls(); } double do_the_test (