PLearn 0.1
|
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 :