Added find_gap_between_convex_hulls()

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<EXP1>& A,
+        const matrix_exp<EXP2>& B,
+        matrix<T,NRa,1>& cA,
+        matrix<T,NRb,1>& 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<T> 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;
+    }
+
 // ----------------------------------------------------------------------------------------
 
 }

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<EXP1>& A,
+        const matrix_exp<EXP2>& B,
+        matrix<T,NRa,1>& cA,
+        matrix<T,NRb,1>& 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.
+    !*/
+
 // ----------------------------------------------------------------------------------------
 
 }

dlib/test/opt_qp_solver.cpp

@@ -356,6 +356,157 @@ namespace
         test_qp4_test7();
     }
 
+// ----------------------------------------------------------------------------------------
+
+    double max_distance_to(
+        const std::vector<matrix<double,0,1>>& a,
+        const std::vector<matrix<double,0,1>>& 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<matrix<double,0,1>>& a,
+        const std::vector<matrix<double,0,1>>& b
+    )
+    {
+        double best_dist = std::numeric_limits<double>::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<matrix<double,0,1>>& s,
+        const matrix<double,0,1>& v
+    )
+    {
+        double best_dist = std::numeric_limits<double>::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<matrix<double,0,1>>& s,
+        const matrix<double,0,1>& 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<matrix<double,0,1>> 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<double,0,1> 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<double> 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<double,0,1> 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 (