dlib C++ Library - linearly_independent_subset_finder.h

    // Copyright (C) 2008  Davis E. King (davisking@users.sourceforge.net)
    // License: Boost Software License   See LICENSE.txt for the full license.
    #ifndef DLIB_LISf__
    #define DLIB_LISf__

    #include <vector>

    #include "linearly_independent_subset_finder_abstract.h"
    #include "../matrix.h"
    #include "function.h"
    #include "../std_allocator.h"
    #include "../algs.h"
    #include "../serialize.h"

    namespace dlib
    {

    // ----------------------------------------------------------------------------------------

        template <typename kernel_type>
        class linearly_independent_subset_finder
        {
            /*!
                INITIAL VALUE
                    - min_strength == 0
                    - min_vect_idx == 0
                    - K_inv.size() == 0
                    - K.size() == 0
                    - dictionary.size() == 0

                CONVENTION
                    - max_dictionary_size() == my_max_dictionary_size
                    - get_kernel() == kernel
                    - minimum_tolerance() == min_tolerance
                    - dictionary_size() == dictionary.size()
                    - get_dictionary() == vector_to_matrix(dictionary)
                    - K.nr() == dictionary.size()
                    - K.nc() == dictionary.size()
                    - for all valid r,c:
                        - K(r,c) == kernel(dictionary[r], dictionary[c])
                    - K_inv == inv(K)

                    - if (dictionary.size() == my_max_dictionary_size) then
                        - for all valid 0 < i < dictionary.size():
                            - Let STRENGTHS[i] == the delta you would get for dictionary[i] (i.e. Approximately
                              Linearly Dependent value) if you removed dictionary[i] from this object and then
                              tried to add it back in.
                            - min_strength == the minimum value from STRENGTHS
                            - min_vect_idx == the index of the element in STRENGTHS with the smallest value
            !*/

        public:
            typedef typename kernel_type::scalar_type scalar_type;
            typedef typename kernel_type::sample_type sample_type;
            typedef typename kernel_type::mem_manager_type mem_manager_type;

            linearly_independent_subset_finder (
                const kernel_type& kernel_,
                unsigned long max_dictionary_size_,
                scalar_type min_tolerance_ = 0.001
            ) :
                kernel(kernel_),
                my_max_dictionary_size(max_dictionary_size_),
                min_tolerance(min_tolerance_)
            {
                // make sure requires clause is not broken
                DLIB_ASSERT(min_tolerance_ > 0,
                    "\tlinearly_independent_subset_finder()"
                    << "\n\tinvalid argument to constructor"
                    << "\n\tmin_tolerance_: " << min_tolerance_
                    << "\n\tthis:           " << this
                    );
                clear_dictionary();
            }

            unsigned long max_dictionary_size() const
            {
                return my_max_dictionary_size;
            }

            const kernel_type& get_kernel (
            ) const
            {
                return kernel;
            }

            scalar_type minimum_tolerance(
            ) const
            {
                return min_tolerance;
            }

            void clear_dictionary ()
            {
                dictionary.clear();
                min_strength = 0;
                min_vect_idx = 0;

                K_inv.set_size(0,0);
                K.set_size(0,0);
            }

            void add (
                const sample_type& x
            )
            {
                const scalar_type kx = kernel(x,x);
                if (dictionary.size() == 0)
                {
                    // set initial state since this is the first sample we have seen
                    K_inv.set_size(1,1);
                    K_inv(0,0) = 1/kx;

                    K.set_size(1,1);
                    K(0,0) = kx;

                    dictionary.push_back(x);
                }
                else
                {
                    // fill in k
                    k.set_size(dictionary.size());
                    for (long r = 0; r < k.nr(); ++r)
                        k(r) = kernel(x,dictionary[r]);

                    // compute the error we would have if we approximated the new x sample
                    // with the dictionary.  That is, do the ALD test from the KRLS paper.
                    a = K_inv*k;
                    scalar_type delta = kx - trans(k)*a;

                    // if this new vector is approximately linearly independent of the vectors
                    // in our dictionary.  Or if our dictionary just isn't full yet.
                    if (delta > min_strength && delta > min_tolerance)
                    {
                        if (dictionary.size() == my_max_dictionary_size)
                        {

                            const long i = min_vect_idx;

                            // replace the min strength vector with x.  Put the new vector onto the end of
                            // dictionary and remove the vector at position i.
                            dictionary.erase(dictionary.begin()+i);
                            dictionary.push_back(x);

                            // compute reduced K_inv.
                            // Remove the i'th vector from the inverse kernel matrix.  This formula is basically
                            // just the reverse of the way K_inv is updated by equation 3.14 below.
                            temp = removerc(K_inv,i,i) - remove_row(colm(K_inv,i)/K_inv(i,i),i)*remove_col(rowm(K_inv,i),i);

                            // recompute these guys since they were computed with the old
                            // kernel matrix
                            k2 = remove_row(k,i);
                            a2 = temp*k2;
                            delta = kx - trans(k2)*a2;

                            // now update temp with the new dictionary vector
                            // update the middle part of the matrix
                            set_subm(K_inv, get_rect(temp)) = temp + a2*trans(a2)/delta;
                            // update the right column of the matrix
                            set_subm(K_inv, 0, temp.nr(),temp.nr(),1) = -a2/delta;
                            // update the bottom row of the matrix
                            set_subm(K_inv, temp.nr(), 0, 1, temp.nr()) = trans(-a2/delta);
                            // update the bottom right corner of the matrix
                            K_inv(temp.nr(), temp.nc()) = 1/delta;

                            // now update the kernel matrix K
                            set_subm(K,get_rect(temp)) = removerc(K, i,i);
                            set_subm(K, 0, K.nr()-1,K.nr()-1,1) = k2;
                            // update the bottom row of the matrix
                            set_subm(K, K.nr()-1, 0, 1, K.nr()-1) = trans(k2);
                            K(K.nr()-1, K.nc()-1) = kx;

                            // now we have to recompute the min_strength in this case
                            recompute_min_strength();
                        }
                        else
                        {
                            // update K_inv by computing the new one in the temp matrix (equation 3.14 from Engel)
                            temp.set_size(K_inv.nr()+1, K_inv.nc()+1);
                            // update the middle part of the matrix
                            set_subm(temp, get_rect(K_inv)) = K_inv + a*trans(a)/delta;
                            // update the right column of the matrix
                            set_subm(temp, 0, K_inv.nr(),K_inv.nr(),1) = -a/delta;
                            // update the bottom row of the matrix
                            set_subm(temp, K_inv.nr(), 0, 1, K_inv.nr()) = trans(-a/delta);
                            // update the bottom right corner of the matrix
                            temp(K_inv.nr(), K_inv.nc()) = 1/delta;
                            // put temp into K_inv
                            temp.swap(K_inv);

                            // update K (the kernel matrix)
                            temp.set_size(K.nr()+1, K.nc()+1);
                            set_subm(temp, get_rect(K)) = K;
                            // update the right column of the matrix
                            set_subm(temp, 0, K.nr(),K.nr(),1) = k;
                            // update the bottom row of the matrix
                            set_subm(temp, K.nr(), 0, 1, K.nr()) = trans(k);
                            temp(K.nr(), K.nc()) = kx;
                            // put temp into K
                            temp.swap(K);

                            // add x to the dictionary
                            dictionary.push_back(x);

                        }
                    }
                }
            }

            void swap (
                linearly_independent_subset_finder& item
            )
            {
                exchange(kernel, item.kernel);
                dictionary.swap(item.dictionary);
                exchange(min_strength, item.min_strength);
                exchange(min_vect_idx, item.min_vect_idx);
                K_inv.swap(item.K_inv);
                K.swap(item.K);
                exchange(my_max_dictionary_size, item.my_max_dictionary_size);
                exchange(min_tolerance, item.min_tolerance);

                // non-state temp members
                a.swap(item.a);
                k.swap(item.k);
                a2.swap(item.a2);
                k2.swap(item.k2);
                temp.swap(item.temp);
            }

            unsigned long dictionary_size (
            ) const { return dictionary.size(); }

            const matrix<sample_type,0,1,mem_manager_type> get_dictionary (
            ) const
            {
                return vector_to_matrix(dictionary);
            }

            friend void serialize(const linearly_independent_subset_finder& item, std::ostream& out)
            {
                serialize(item.kernel, out);
                serialize(item.dictionary, out);
                serialize(item.min_strength, out);
                serialize(item.min_vect_idx, out);
                serialize(item.K_inv, out);
                serialize(item.K, out);
                serialize(item.my_max_dictionary_size, out);
                serialize(item.min_tolerance, out);
            }

            friend void deserialize(linearly_independent_subset_finder& item, std::istream& in)
            {
                deserialize(item.kernel, in);
                deserialize(item.dictionary, in);
                deserialize(item.min_strength, in);
                deserialize(item.min_vect_idx, in);
                deserialize(item.K_inv, in);
                deserialize(item.K, in);
                deserialize(item.my_max_dictionary_size, in);
                deserialize(item.min_tolerance, in);
            }

            const sample_type& operator[] (
                unsigned long index
            ) const
            {
                return dictionary[index];
            }

        private:

            typedef std_allocator<sample_type, mem_manager_type> alloc_sample_type;
            typedef std_allocator<scalar_type, mem_manager_type> alloc_scalar_type;
            typedef std::vector<sample_type,alloc_sample_type> dictionary_vector_type;
            typedef std::vector<scalar_type,alloc_scalar_type> scalar_vector_type;

            void recompute_min_strength (
            )
            /*!
                ensures
                    - recomputes the min_strength and min_vect_idx values
                      so that they are correct with respect to the CONVENTION
            !*/
            {
                min_strength = std::numeric_limits<scalar_type>::max();

                // here we loop over each dictionary vector and compute what its delta would be if
                // we were to remove it from the dictionary and then try to add it back in.
                for (unsigned long i = 0; i < dictionary.size(); ++i)
                {
                    // compute a2 = K_inv*k but where dictionary vector i has been removed
                    a2 = (removerc(K_inv,i,i) - remove_row(colm(K_inv,i)/K_inv(i,i),i)*remove_col(rowm(K_inv,i),i)) *
                        (remove_row(colm(K,i),i));
                    scalar_type delta = K(i,i) - trans(remove_row(colm(K,i),i))*a2;

                    if (delta < min_strength)
                    {
                        min_strength = delta;
                        min_vect_idx = i;
                    }
                }
            }

            kernel_type kernel;
            dictionary_vector_type dictionary;
            scalar_type min_strength;
            unsigned long min_vect_idx;

            matrix<scalar_type,0,0,mem_manager_type> K_inv;
            matrix<scalar_type,0,0,mem_manager_type> K;

            unsigned long my_max_dictionary_size;
            scalar_type min_tolerance;

            // temp variables here just so we don't have to reconstruct them over and over.  Thus,
            // they aren't really part of the state of this object.
            matrix<scalar_type,0,1,mem_manager_type> a, a2;
            matrix<scalar_type,0,1,mem_manager_type> k, k2;
            matrix<scalar_type,0,0,mem_manager_type> temp;

        };

    // ----------------------------------------------------------------------------------------

        template <typename kernel_type>
        void swap(linearly_independent_subset_finder<kernel_type>& a, linearly_independent_subset_finder<kernel_type>& b)
        { a.swap(b); }

    // ----------------------------------------------------------------------------------------

    }

    #endif // DLIB_LISf__