Added initial cut of manifold regularization stuff. Code needs to be cleaned up.

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

dlib/manifold_regularization.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_MANIFOLD_REGULARIzATION_HEADER
+#define DLIB_MANIFOLD_REGULARIzATION_HEADER
+
+#include "manifold_regularization/sample_pair.h"
+#include "manifold_regularization/graph_creation.h"
+#include "manifold_regularization/linear_manifold_regularizer.h"
+#include "manifold_regularization/function_objects.h"
+
+#endif // DLIB_MANIFOLD_REGULARIzATION_HEADER
+
+
+

dlib/manifold_regularization/function_objects.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_MR_FUNCTION_ObJECTS_H__
+#define DLIB_MR_FUNCTION_ObJECTS_H__
+
+#include "function_objects_abstract.h"
+#include "../matrix.h"
+#include <cmath>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    struct squared_euclidean_distance
+    {
+        template <typename sample_type>
+        double operator() (
+            const sample_type& a,
+            const sample_type& b
+        ) const
+        { 
+            return length_squared(a-b);
+        }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    struct use_weights_of_one 
+    {
+        template <typename edge_type>
+        double operator() (
+            const edge_type&
+        ) const
+        { 
+            return 1;
+        }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    struct use_gaussian_weights 
+    {
+        use_gaussian_weights (
+            double g
+        )
+        {
+            gamma = g;
+        }
+
+        double gamma;
+
+        template <typename edge_type>
+        double operator() (
+            const edge_type& e
+        ) const
+        { 
+            return std::exp(-gamma*e.distance());
+        }
+    };
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_MR_FUNCTION_ObJECTS_H__
+
+

dlib/manifold_regularization/graph_creation.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_GRAPH_CrEATION_H__
+#define DLIB_GRAPH_CrEATION_H__
+
+#include "graph_creation_abstract.h"
+#include <limits>
+#include <vector>
+#include "../string.h"
+#include "../rand.h"
+#include <algorithm>
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    namespace impl
+    {
+        template <typename iterator>
+        iterator iterator_of_worst (
+            iterator begin,
+            const iterator& end
+        ) 
+        /*!
+            ensures
+                - returns an iterator that points to the element in the given range that has the biggest
+                  distance 
+        !*/
+        {
+            float dist = begin->distance();
+            iterator worst = begin;
+            for (; begin != end; ++begin)
+            {
+                if (begin->distance() > dist)
+                {
+                    dist = begin->distance();
+                    worst = begin;
+                }
+            }
+
+            return worst;
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename vector_type,
+        typename distance_function_type,
+        typename alloc,
+        typename T
+        >
+    void find_percent_shortest_edges_randomly (
+        const vector_type& samples,
+        const distance_function_type& dist_funct,
+        const double percent,
+        const unsigned long num,
+        const T& random_seed,
+        std::vector<sample_pair, alloc>& out
+    )
+    /*!
+        requires
+            - samples.size() > 1
+        - 0 < percent <= 1
+        - num > 0
+        - random_seed must be convertible to a string by dlib::cast_to_string()
+    !*/
+    {
+        std::vector<sample_pair, alloc> edges;
+        edges.reserve(num);
+
+        dlib::rand::kernel_1a rnd;
+        rnd.set_seed(cast_to_string(random_seed));
+
+        // randomly sample a bunch of edges
+        while (edges.size() < num)
+        {
+            const unsigned long idx1 = rnd.get_random_32bit_number()%samples.size();
+            const unsigned long idx2 = rnd.get_random_32bit_number()%samples.size();
+            if (idx1 != idx2)
+            {
+                edges.push_back(sample_pair(idx1, idx2, dist_funct(samples[idx1], samples[idx2])));
+            }
+        }
+
+
+        // sort the edges so that duplicate edges will be adjacent
+        std::sort(edges.begin(), edges.end(), &order_by_index);
+
+        // now put edges into out while avoiding duplicates
+        out.clear();
+        out.reserve(edges.size());
+        out.push_back(edges[0]);
+        for (unsigned long i = 1; i < edges.size(); ++i)
+        {
+            if (edges[i] != edges[i-1])
+            {
+                out.push_back(edges[i]);
+            }
+        }
+
+
+        // now sort all the edges by distance and take the percent with the smallest distance
+        std::sort(out.begin(), out.end(), &order_by_distance);
+        out.swap(edges);
+        out.assign(edges.begin(), edges.begin() + edges.size()*percent);
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename vector_type,
+        typename distance_function_type,
+        typename alloc
+        >
+    void find_k_nearest_neighbors (
+        const vector_type& samples,
+        const distance_function_type& dist_funct,
+        const unsigned long k,
+        std::vector<sample_pair, alloc>& out
+    )
+    /*!
+        requires
+            - samples.size() > k
+        - k > 0
+    !*/
+    {
+        using namespace impl;
+        std::vector<sample_pair> edges;
+
+        edges.resize(samples.size()*k);
+
+        std::vector<float> worst_dists(samples.size(), std::numeric_limits<float>::max());
+
+        std::vector<sample_pair>::iterator begin_i, end_i, begin_j, end_j, itr;
+        begin_i = edges.begin();
+        end_i = begin_i + k;
+
+        // Loop over all combinations of samples.   We will maintain the iterator ranges so that
+        // within the inner for loop we have:
+        //   [begin_i, end_i) == the range in edges that contains neighbors of samples[i]
+        //   [begin_j, end_j) == the range in edges that contains neighbors of samples[j]
+        for (unsigned long i = 0; i < samples.size(); ++i)
+        {
+            begin_j = begin_i + 1;
+            end_j = begin_j + k;
+
+            for (unsigned long j = i+1; j < samples.size(); ++j)
+            {
+                const float dist = dist_funct(samples[i], samples[j]);
+
+                if (dist < worst_dists[i])
+                {
+                    *iterator_of_worst(begin_i, end_i) = sample_pair(i, j, dist);
+                    worst_dists[i] = iterator_of_worst(begin_i, end_i)->distance();
+                }
+
+                if (dist < worst_dists[j])
+                {
+                    *iterator_of_worst(begin_j, end_j) = sample_pair(i, j, dist);
+                    worst_dists[j] = iterator_of_worst(begin_j, end_j)->distance();
+                }
+
+                begin_j += k;
+                end_j += k;
+            }
+
+            begin_i += k;
+            end_i += k;
+        }
+
+        // sort the edges so that duplicate edges will be adjacent
+        std::sort(edges.begin(), edges.end(), &order_by_index);
+
+        // now put edges into out while avoiding duplicates
+        out.clear();
+        out.reserve(edges.size());
+        out.push_back(edges[0]);
+        for (unsigned long i = 1; i < edges.size(); ++i)
+        {
+            if (edges[i] != edges[i-1])
+            {
+                out.push_back(edges[i]);
+            }
+        }
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_GRAPH_CrEATION_H__
+
+

dlib/manifold_regularization/linear_manifold_regularizer.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_LINEAR_MANIFOLD_ReGULARIZER_H__
+#define DLIB_LINEAR_MANIFOLD_ReGULARIZER_H__
+
+#include "linear_manifold_regularizer.h"
+#include <limits>
+#include <vector>
+#include "../serialize.h"
+#include "../matrix.h"
+
+namespace dlib
+{
+    namespace impl
+    {
+        class undirected_adjacency_list
+        {
+            /*!
+                WHAT THIS OBJECT REPRESENTS
+                    This object is simply a tool for turning a vector of sample_pair objects
+                    into an adjacency list with floating point weights on each edge.  
+            !*/
+        public:
+
+            undirected_adjacency_list (
+            )
+            {
+                _size = 0;
+            }
+
+            struct neighbor 
+            {
+                neighbor(unsigned long idx, float w):index(idx), weight(w) {}
+                neighbor():index(0), weight(0) {}
+
+                unsigned long index;
+                float weight;
+            };
+
+            typedef std::vector<neighbor>::const_iterator const_iterator;
+
+            unsigned long size (
+            ) const
+            {
+                return _size;
+            }
+
+            const_iterator begin(
+                unsigned long idx
+            ) const
+            /*!
+                requires
+                    - idx < size()
+            !*/
+            {
+                return blocks[idx];
+            }
+
+            const_iterator end(
+                unsigned long idx
+            ) const
+            /*!
+                requires
+                    - idx < size()
+            !*/
+            {
+                return blocks[idx+1];
+            }
+
+
+            template <typename vector_type, typename weight_function_type>
+            void build (
+                const vector_type& edges,
+                const weight_function_type& weight_funct
+            ) 
+            /*!
+                requires
+                    - vector_type == a type with an interface compatible with std::vector and 
+                      it must in turn contain objects with an interface compatible with dlib::sample_pair
+                    - edges.size() > 0
+                    - all the elements of edges are unique.  That is:
+                        - for all valid i and j where i != j:
+                          it must be true that edges[i] != edges[j]
+                    - weight_funct(edges[i]) must be a valid expression that evaluates to a
+                      floating point number
+            !*/
+            {
+
+
+                // Figure out how many neighbors each sample ultimately has.  We do this so 
+                // we will know how much space to allocate in the data vector.
+                std::vector<unsigned long> num_neighbors;
+                num_neighbors.reserve(edges.size());
+
+                for (unsigned long i = 0; i < edges.size(); ++i)
+                {
+                    // make sure num_neighbors is always big enough 
+                    const unsigned long min_size = std::max(edges[i].index1(), edges[i].index2())+1;
+                    if (num_neighbors.size() < min_size)
+                        num_neighbors.resize(min_size,  0);
+
+                    num_neighbors[edges[i].index1()] += 1;
+                    num_neighbors[edges[i].index2()] += 1;
+                }
+
+                _size = num_neighbors.size();
+
+                // Now setup the iterators in blocks.  Also setup a version of blocks that holds
+                // non-const iterators so we can use it below when we populate data.
+                std::vector<std::vector<neighbor>::iterator> mutable_blocks;
+                data.resize(edges.size()*2); // each edge will show up twice 
+                blocks.resize(_size + 1);
+                blocks[0] = data.begin();
+                mutable_blocks.resize(_size + 1);
+                mutable_blocks[0] = data.begin();
+                for (unsigned long i = 0; i < num_neighbors.size(); ++i)
+                {
+                    blocks[i+1]         = blocks[i]         + num_neighbors[i];
+                    mutable_blocks[i+1] = mutable_blocks[i] + num_neighbors[i];
+                }
+
+                // finally, put the edges into data
+                for (unsigned long i = 0; i < edges.size(); ++i)
+                {
+                    const float weight = weight_funct(edges[i]);
+                    *mutable_blocks[edges[i].index1()]++ = neighbor(edges[i].index2(), weight);
+                    *mutable_blocks[edges[i].index2()]++ = neighbor(edges[i].index1(), weight);
+                }
+
+            }
+
+        private:
+
+            /*!
+                INITIAL VALUE
+                    - _size == 0
+                    - data.size() == 0
+                    - blocks.size() == 0
+
+                CONVENTION
+                    - size() == _size
+                    - blocks.size() == _size + 1
+                    - blocks == a vector of iterators that point into data.  
+                      For all valid i:
+                        - The iterator range [blocks[i], blocks[i+1]) contains all the edges
+                          for the i'th node in the graph
+            !*/
+
+            std::vector<neighbor> data;
+            std::vector<const_iterator> blocks; 
+            unsigned long _size;
+        };
+
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    template <
+        typename matrix_type
+        >
+    class linear_manifold_regularizer
+    {
+        /*!
+            REQUIREMENTS ON matrix_type
+                Must be some type of dlib::matrix.
+
+            WHAT THIS OBJECT REPRESENTS
+                This object computes the inv(T) matrix described in the following paper:
+                    Linear Manifold Regularization for Large Scale Semi-supervised Learning
+                    by Vikas Sindhwani, Partha Niyogi, and Mikhail Belkin
+        !*/
+
+    public:
+        typedef typename matrix_type::mem_manager_type mem_manager_type;
+        typedef typename matrix_type::type scalar_type;
+        typedef typename matrix_type::layout_type layout_type;
+        typedef matrix<scalar_type,0,0,mem_manager_type,layout_type> general_matrix;
+
+
+        template <
+            typename vector_type1, 
+            typename vector_type2, 
+            typename weight_function_type
+            >
+        void build (
+            const vector_type1& samples,
+            const vector_type2& edges,
+            const weight_function_type& weight_funct
+        )
+        {
+            impl::undirected_adjacency_list graph;
+            graph.build(edges, weight_funct);
+
+            make_mr_matrix(samples, graph);
+        }
+
+        general_matrix get_transformation_matrix (
+            scalar_type intrinsic_regularization_strength
+        ) const
+        /*!
+            requires
+                - intrinsic_regularization_strength >= 0
+        !*/
+        {
+            return inv_lower_triangular(chol(identity_matrix<scalar_type>(reg_mat.nr()) + intrinsic_regularization_strength*reg_mat));
+        }
+
+    private:
+
+        template <typename vector_type>
+        void make_mr_matrix (
+            const vector_type& samples,
+            const impl::undirected_adjacency_list& graph
+        )
+        /*!
+            requires
+                - samples.size() == graph.size()
+        !*/
+        {
+            const unsigned long dims = samples[0].size();
+            reg_mat.set_size(dims,dims);
+            reg_mat = 0;
+
+            // Compute trans(X)*lap(graph)*X where X is the data matrix 
+            // (i.e. the matrix that contains all the samples in its rows)
+            // and lap(graph) is the laplacian matrix of the graph.
+
+            typename impl::undirected_adjacency_list::const_iterator beg, end;
+
+            // loop over the columns of the X matrix
+            for (unsigned long d = 0; d < dims; ++d)
+            {
+                // loop down the row of X
+                for (unsigned long i = 0; i < graph.size(); ++i)
+                {
+                    beg = graph.begin(i);
+                    end = graph.end(i);
+
+                    // if this node in the graph has any neighbors
+                    if (beg != end)
+                    {
+                        double weight_sum = 0;
+                        double val = 0;
+                        for (; beg != end; ++beg)
+                        {
+                            val -= beg->weight * samples[beg->index](d);
+                            weight_sum += beg->weight;
+                        }
+                        val += weight_sum * samples[i](d);
+
+                        for (unsigned long j = 0; j < dims; ++j)
+                        {
+                            reg_mat(d,j) += val*samples[i](j);
+                        }
+                    }
+                }
+            }
+
+        }
+
+        general_matrix reg_mat;
+    };
+
+}
+
+// ----------------------------------------------------------------------------------------
+
+#endif // DLIB_LINEAR_MANIFOLD_ReGULARIZER_H__
+

dlib/manifold_regularization/sample_pair.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_SAMPLE_PaIR_H__
+#define DLIB_SAMPLE_PaIR_H__
+
+#include "sample_pair_abstract.h"
+#include <limits>
+#include "../serialize.h"
+
+namespace dlib
+{
+
+// ----------------------------------------------------------------------------------------
+
+    struct sample_pair 
+    {
+        sample_pair(
+        ) : 
+            _index1(0),
+            _index2(0)
+        {
+            _distance = std::numeric_limits<float>::max();
+        }
+
+        sample_pair (
+            const unsigned long idx1,
+            const unsigned long idx2,
+            const float dist
+        )
+        {
+            _distance = dist;
+            if (idx1 < idx2)
+            {
+                _index1 = idx1;
+                _index2 = idx2;
+            }
+            else
+            {
+                _index1 = idx2;
+                _index2 = idx1;
+            }
+        }
+
+        const unsigned long& index1 (
+        ) const { return _index1; }
+
+        const unsigned long& index2 (
+        ) const { return _index2; }
+
+        const float& distance (
+        ) const { return _distance; }
+
+    private:
+        unsigned long _index1;
+        unsigned long _index2;
+        float _distance;
+    };
+
+// ----------------------------------------------------------------------------------------
+
+    inline bool order_by_index (
+        const sample_pair& a,
+        const sample_pair& b
+    )
+    {
+        return a.index1() < b.index1() || (a.index1() == b.index1() && a.index2() < b.index2());
+    }
+
+    inline bool order_by_distance (
+        const sample_pair& a,
+        const sample_pair& b
+    )
+    {
+        return a.distance() < b.distance();
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    inline bool operator == (
+        const sample_pair& a,
+        const sample_pair& b
+    ) 
+    {
+        return a.index1() == b.index1() && a.index2() == b.index2();
+    }
+
+    inline bool operator != (
+        const sample_pair& a,
+        const sample_pair& b
+    ) 
+    {
+        return !(a == b); 
+    }
+
+// ----------------------------------------------------------------------------------------
+
+    inline void serialize (
+        const sample_pair& item,
+        std::ostream& out
+    )
+    {
+        try
+        {
+            serialize(item.index1(),out);
+            serialize(item.index2(),out);
+            serialize(item.distance(),out);
+        }
+        catch (serialization_error& e)
+        { 
+            throw serialization_error(e.info + "\n   while serializing object of type sample_pair"); 
+        }
+    }
+
+    inline void deserialize (
+        sample_pair& item,
+        std::istream& in 
+    )
+    {
+        try
+        {
+            unsigned long idx1, idx2;
+            float dist;
+
+            deserialize(idx1,in);
+            deserialize(idx2,in);
+            deserialize(dist,in);
+            item = sample_pair(idx1, idx2, dist);
+        }
+        catch (serialization_error& e)
+        { 
+            throw serialization_error(e.info + "\n   while deserializing object of type sample_pair"); 
+        }
+    }
+
+// ----------------------------------------------------------------------------------------
+
+}
+
+#endif // DLIB_SAMPLE_PaIR_H__
+