PLearn 0.1
ReconstructionWeightsKernel.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // ReconstructionWeightsKernel.cc
00004 //
00005 // Copyright (C) 2004 Olivier Delalleau 
00006 // 
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 // 
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 // 
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 // 
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 // 
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 // 
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 /* *******************************************************      
00036  * $Id: ReconstructionWeightsKernel.cc 6508 2006-12-15 02:35:49Z lamblin $ 
00037  ******************************************************* */
00038 
00039 // Authors: Olivier Delalleau
00040 
00044 #include "DistanceKernel.h"
00045 #include "DotProductKernel.h"
00046 #include <plearn/math/plapack.h>          
00047 #include "ReconstructionWeightsKernel.h"
00048 
00049 namespace PLearn {
00050 using namespace std;
00051 
00053 // ReconstructionWeightsKernel //
00055 ReconstructionWeightsKernel::ReconstructionWeightsKernel() 
00056     : build_in_progress(false),
00057       new_data(true),
00058       ignore_nearest(1),
00059       knn(5),
00060       regularizer(1e-6)
00061 {
00062     is_symmetric = false;
00063     sub_data = new SelectRowsVMatrix();
00064 }
00065 
00066 PLEARN_IMPLEMENT_OBJECT(ReconstructionWeightsKernel,
00067                         "Computes the reconstruction weights of a point given its neighbors.",
00068                         "K(x, x_i) = the weight of x_i in the reconstruction of x by its knn\n"
00069                         "nearest neighbors. More precisely, we compute weights W_i such that\n"
00070                         "|| x - \\sum_j W_i x_i ||^2 is minimized, and K(x,x_i) = W_i.\n"
00071                         "If the second argument is not in the training set, K(x,y) will be 0.\n"
00072                         "In order not to compute K(x_i, x_j) = delta_{ij} when applied on\n"
00073                         "training points, one can set the 'ignore_nearest' option to 1 (or more),\n"
00074                         "which will ensure we do not use x_i itself in its reconstruction by its\n"
00075                         "nearest neighbors (however, the total number of neighbors computed,\n"
00076                         "including x_i itself, will always stay equal to knn).\n"
00077                         "Note that this is NOT a symmetric kernel!\n"
00078     );
00079 
00081 // declareOptions //
00083 void ReconstructionWeightsKernel::declareOptions(OptionList& ol)
00084 {
00085     // Build options.
00086 
00087     declareOption(ol, "knn", &ReconstructionWeightsKernel::knn, OptionBase::buildoption,
00088                   "The number of nearest neighbors considered (including the point itself).");
00089 
00090     declareOption(ol, "regularizer", &ReconstructionWeightsKernel::regularizer, OptionBase::buildoption,
00091                   "A factor multiplied by the trace of the local Gram matrix and added to\n"
00092                   "the diagonal to ensure stability when solving the linear system.");
00093 
00094     declareOption(ol, "ignore_nearest", &ReconstructionWeightsKernel::ignore_nearest, OptionBase::buildoption,
00095                   "The number of nearest neighbors to ignore when computing the reconstruction weights.");
00096 
00097     declareOption(ol, "distance_kernel", &ReconstructionWeightsKernel::distance_kernel, OptionBase::buildoption,
00098                   "The kernel used to compute the distances.\n"
00099                   "If not specified, then the usual Euclidean distance will be used.");
00100 
00101     declareOption(ol, "dot_product_kernel", &ReconstructionWeightsKernel::dot_product_kernel, OptionBase::buildoption,
00102                   "The kernel used to compute dot products in the neighborhood of each data point.\n"
00103                   "If not specified, then the usual Euclidean dot product will be used.");
00104 
00105     // Now call the parent class' declareOptions
00106     inherited::declareOptions(ol);
00107 }
00108 
00110 // build //
00112 void ReconstructionWeightsKernel::build()
00113 {
00114     build_in_progress = true;
00115     inherited::build();
00116     build_();
00117 }
00118 
00120 // build_ //
00122 void ReconstructionWeightsKernel::build_()
00123 {
00124     if (distance_kernel) {
00125         dist_ker = distance_kernel;
00126     } else {
00127         dist_ker = new DistanceKernel(2);
00128         dist_ker->report_progress = this->report_progress;
00129     }
00130 
00131     if (dot_product_kernel) {
00132         dp_ker = dot_product_kernel;
00133     } else {
00134         dp_ker = new DotProductKernel();
00135         dp_ker->build();
00136     }
00137     // Safety check.
00138     if (ignore_nearest > knn)
00139         PLERROR("In ReconstructionWeightsKernel::build_ - You can't ignore more than 'knn' neighbors");
00140     build_in_progress = false;
00141     // This code, normally executed in Kernel::build_, can only be executed
00142     // now beause the kernels 'dist_ker' and 'dp_ker' have to be initialized.
00143     if (specify_dataset) {
00144         this->setDataForKernelMatrix(specify_dataset);
00145     }
00146 }
00147 
00149 // computeLLEMatrix //
00151 void ReconstructionWeightsKernel::computeLLEMatrix(const Mat& lle_mat) const {
00152     if (lle_mat.length() != n_examples || lle_mat.width() != n_examples)
00153         PLERROR("In ReconstructionWeightsKernel::computeLLEMatrix - Wrong size for 'lle_mat'");
00154     lle_mat.clear();
00155     PP<ProgressBar> pb;
00156     if (report_progress)
00157         pb = new ProgressBar("Computing LLE matrix", n_examples);
00158     int neighb_j, neighb_k;
00159     real w_ij;
00160     for (int i = 0; i < n_examples; i++) {
00161         for (int j = 0; j < knn - 1; j++) {
00162             neighb_j = neighbors(i, j + 1);
00163             w_ij = weights(i, j);
00164             lle_mat(i, neighb_j) += w_ij;
00165             lle_mat(neighb_j, i) += w_ij;
00166             for (int k = 0; k < knn - 1; k++) {
00167                 neighb_k = neighbors(i, k + 1);
00168                 lle_mat(neighb_j, neighb_k) -= w_ij * weights(i, k);
00169             }
00170         }
00171         if (report_progress)
00172             pb->update(i + 1);
00173     }
00174 }
00175 
00177 // computeWeights //
00179 void ReconstructionWeightsKernel::computeWeights() {
00180     static Vec point_i;
00181     static Vec weights_i;
00182     if (!data)
00183         PLERROR("In ReconstructionWeightsKernel::computeWeights - Can only be called if 'data' has been set");
00184     point_i.resize(data_inputsize);
00185     weights.resize(n_examples, knn - ignore_nearest); // Allocate memory for the weights.
00186     // First compute the nearest neighbors.
00187     Mat distances(n_examples, n_examples);
00188     dist_ker->computeGramMatrix(distances);
00189     neighbors =
00190         computeKNNeighbourMatrixFromDistanceMatrix(distances, knn, true,
00191                                                                                                    report_progress != 0);
00192     distances = Mat(); // Free memory.
00193     // Fill the 'is_neighbor_of' vector.
00194     is_neighbor_of.resize(n_examples);
00195     TVec<int> row(2);
00196     for (int i = 0; i < n_examples; i++)
00197         is_neighbor_of[i].resize(0, 2);
00198     for (int i = 0; i < n_examples; i++) {
00199         row[0] = i;
00200         for (int j = ignore_nearest; j < knn; j++) {
00201             row[1] = j;
00202             is_neighbor_of[neighbors(i,j)].appendRow(row);
00203         }
00204     }
00205     for (int i = 0; i < n_examples; i++)
00206         sortRows(is_neighbor_of[i]);
00207     // Then compute the weights for each point i.
00208     TVec<int> neighbors_of_i;
00209     PP<ProgressBar> pb;
00210     if (report_progress)
00211         pb = new ProgressBar("Computing reconstruction weights", n_examples);
00212     for (int i = 0; i < n_examples; i++) {
00213         // Isolate the neighbors.
00214         neighbors_of_i = neighbors(i).subVec(ignore_nearest, knn - ignore_nearest);
00215         weights_i = weights(i);
00216         data->getSubRow(i, 0, point_i);
00217         reconstruct(point_i, neighbors_of_i, weights_i);
00218         if (report_progress)
00219             pb->update(i+1);
00220     }
00221 }
00222 
00224 // evaluate //
00226 real ReconstructionWeightsKernel::evaluate(const Vec& x1, const Vec& x2) const {
00227     static int j;
00228     if (isInData(x2, &j)) {
00229         // x2 is in the training set, thus it makes sense to use it in the reconstruction.
00230         return evaluate_x_i(x1, j);
00231     } else {
00232         // x2 is not in the training set, thus its weight is 0.
00233         return 0;
00234     }
00235 }
00236 
00238 // evaluate_i_j //
00240 real ReconstructionWeightsKernel::evaluate_i_j(int i, int j) const {
00241     static TVec<int> neighbors_of_i;
00242     if (ignore_nearest == 0) {
00243         // We do not ignore the nearest neighbor, which is i itself. Thus the
00244         // weight is \delta_{ij}, since i is reconstructed exactly by itself.
00245         if (i == j)
00246             return 1.0;
00247         else
00248             return 0;
00249     } else {
00250 #ifdef BOUNDCHECK
00251         if (ignore_nearest != knn - weights.width())
00252             PLERROR("In ReconstructionWeightsKernel::evaluate_i_j - You must recompute the weights after modifying the 'ignore_nearest' option");
00253 #endif
00254         neighbors_of_i = neighbors(i);
00255         for (int k = ignore_nearest; k < knn; k++) {
00256             if (neighbors_of_i[k] == j) {
00257                 return weights(i, k - ignore_nearest);
00258             }
00259         }
00260         return 0;
00261     }
00262 }
00263 
00265 // evaluate_i_x //
00267 real ReconstructionWeightsKernel::evaluate_i_x(int i, const Vec& x, real squared_norm_of_x) const {
00268     static int j;
00269     if (isInData(x, &j))
00270         return evaluate_i_j(i,j);
00271     else
00272         return 0;
00273 }
00274 
00276 // evaluate_sum_k_i_k_j //
00278 real ReconstructionWeightsKernel::evaluate_sum_k_i_k_j(int i, int j) const {
00279     static TMat<int> i_is_neighb_of, j_is_neighb_of;
00280     i_is_neighb_of = is_neighbor_of[i];
00281     j_is_neighb_of = is_neighbor_of[j];
00282     int test_n;
00283     int k_i = 0;
00284     int k_j = 0;
00285     real sum = 0;
00286     // Safety check
00287     if (ignore_nearest != knn - weights.width())
00288         PLERROR("In ReconstructionWeightsKernel::evaluate_sum_k_i_k_j - You must recompute the weights after modifying 'ignore_nearest'");
00289     while (k_i < i_is_neighb_of.length()) {
00290         test_n = i_is_neighb_of(k_i, 0);
00291         while (k_j < j_is_neighb_of.length() && test_n > j_is_neighb_of(k_j, 0))
00292             k_j++;
00293         if (k_j < j_is_neighb_of.length()) {
00294             if (test_n == j_is_neighb_of(k_j, 0)) {
00295                 // Found a common k.
00296                 sum += weights(test_n, i_is_neighb_of(k_i, 1) - ignore_nearest) * weights(test_n, j_is_neighb_of(k_j, 1) - ignore_nearest);
00297                 k_i++;
00298                 k_j++;
00299             } else {
00300                 // Increase k_i.
00301                 test_n = j_is_neighb_of(k_j, 0);
00302                 while (k_i < i_is_neighb_of.length() && test_n > i_is_neighb_of(k_i, 0))
00303                     k_i++;
00304             }
00305         } else {
00306             // No more common point.
00307             return sum;
00308         }
00309     }
00310     return sum;
00311 }
00312 
00314 // evaluate_x_i //
00316 real ReconstructionWeightsKernel::evaluate_x_i(const Vec& x, int i, real squared_norm_of_x) const {
00317     return evaluate_x_i_again(x, i, squared_norm_of_x, true);
00318 }
00319 
00321 // evaluate_x_i_again //
00323 real ReconstructionWeightsKernel::evaluate_x_i_again(const Vec& x, int i, real squared_norm_of_x, bool first_time) const {
00324     if (first_time) {
00325         neighbors_of_x.resize(knn);
00326         // Find nearest neighbors of x.
00327         dist_ker->computeNearestNeighbors(x, k_xi_x_sorted, knn);
00328         neighbors_of_x << k_xi_x_sorted.subMat(ignore_nearest, 1, knn, 1);
00329         // Find reconstruction weights.
00330         reconstruct(x, neighbors_of_x, weights_x);
00331     }
00332     int n_j = neighbors_of_x.find(i);
00333     if (n_j == -1)
00334         // The point i is not a neighbor of x.
00335         return 0;
00336     return weights_x[n_j];
00337 }
00338 
00340 // makeDeepCopyFromShallowCopy //
00342 void ReconstructionWeightsKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00343 {
00344     inherited::makeDeepCopyFromShallowCopy(copies);
00345 
00346     // ### Call deepCopyField on all "pointer-like" fields 
00347     // ### that you wish to be deepCopied rather than 
00348     // ### shallow-copied.
00349     // ### ex:
00350     // deepCopyField(trainvec, copies);
00351 
00352     // ### Remove this line when you have fully implemented this method.
00353     PLERROR("ReconstructionWeightsKernel::makeDeepCopyFromShallowCopy not fully (correctly) implemented yet!");
00354 }
00355 
00357 // reconstruct //
00359 void ReconstructionWeightsKernel::reconstruct(const Vec& x, const TVec<int>& neighbors, Vec& w) const {
00360     static bool need_init;
00361     need_init = new_data;
00362     int k_neighb = neighbors.length();
00363     if (ones.length() != k_neighb) {
00364         // 'ones' does not have the right size.
00365         need_init = true;
00366         ones.resize(k_neighb);
00367         ones.fill(1);
00368     }
00369     w.resize(k_neighb);
00370     if (need_init) {
00371         // This is the first execution.
00372         local_gram.resize(k_neighb, k_neighb);
00373         centered_neighborhood = new ShiftAndRescaleVMatrix();
00374         centered_neighborhood->no_scale = true;
00375         centered_neighborhood->negate_shift = true;
00376         centered_neighborhood->automatic = false;
00377         centered_neighborhood->source = (SelectRowsVMatrix*) sub_data;
00378         new_data = false;
00379     }
00380     // Center data on x.
00381     sub_data->indices = neighbors;
00382     sub_data->build();
00383     centered_neighborhood->shift = x;
00384     centered_neighborhood->build();
00385     // TODO Get rid of this expensive build.
00386     // Compute the local Gram matrix.
00387     dp_ker->setDataForKernelMatrix((ShiftAndRescaleVMatrix*) centered_neighborhood);
00388     dp_ker->computeGramMatrix(local_gram);
00389     // Add regularization on the diagonal.
00390     regularizeMatrix(local_gram, regularizer);
00391     // Solve linear system.
00392     Vec weights_x = solveLinearSystem(local_gram, ones);
00393     // TODO Avoid the copy of the weights.
00394     w << weights_x;
00395     // Ensure the sum of weights is 1 to get final solution.
00396     w /= sum(w);
00397 }
00398 
00400 // setDataForKernelMatrix //
00402 void ReconstructionWeightsKernel::setDataForKernelMatrix(VMat the_data) {
00403     if (build_in_progress)
00404         return;
00405     inherited::setDataForKernelMatrix(the_data);
00406     dist_ker->setDataForKernelMatrix(the_data);
00407     sub_data->source = the_data;
00408     new_data = true;
00409     computeWeights();
00410 }
00411 
00412 } // end of namespace PLearn
00413 
00414 
00415 /*
00416   Local Variables:
00417   mode:c++
00418   c-basic-offset:4
00419   c-file-style:"stroustrup"
00420   c-file-offsets:((innamespace . 0)(inline-open . 0))
00421   indent-tabs-mode:nil
00422   fill-column:79
00423   End:
00424 */
00425 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines