PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // Matern1ARDKernel.cc 00004 // 00005 // Copyright (C) 2009 Nicolas Chapados 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 // Authors: Nicolas Chapados 00036 00040 #include "Matern1ARDKernel.h" 00041 #include <plearn/math/pl_math.h> 00042 00043 namespace PLearn { 00044 using namespace std; 00045 00046 PLEARN_IMPLEMENT_OBJECT( 00047 Matern1ARDKernel, 00048 "Matern kernel with nu=1/2 that can be used for Automatic Relevance Determination.", 00049 "With nu=1/2, the Matern kernel corresponds to the Ornstein-Uhlenbeck\n" 00050 "process. This function is specified as:\n" 00051 "\n" 00052 " k(x,y) = (sf / (2*a)) * exp(-a sum_i |x_i - y_i|/w_i) * k_kron(x,y)\n" 00053 "\n" 00054 "where sf = softplus(isp_signal_sigma), a = softplus(isp_persistence), w_i =\n" 00055 "softplus(isp_global_sigma + isp_input_sigma[i]), and k_kron(x,y) is the\n" 00056 "result of the KroneckerBaseKernel evaluation, or 1.0 if there are no\n" 00057 "Kronecker terms. Note that since the Kronecker terms are incorporated\n" 00058 "multiplicatively, the very presence of the term associated to this kernel\n" 00059 "can be gated by the value of some input variable(s) (that are incorporated\n" 00060 "within one or more Kronecker terms).\n" 00061 "\n" 00062 "Note that to make its operations more robust when used with unconstrained\n" 00063 "optimization of hyperparameters, all hyperparameters of this kernel are\n" 00064 "specified in the inverse softplus domain. See IIDNoiseKernel for more\n" 00065 "explanations.\n" 00066 ); 00067 00068 00069 Matern1ARDKernel::Matern1ARDKernel() 00070 : m_isp_persistence(pl_log(exp(1.0) - 1.)) // inverse-softplus(1.0) 00071 { } 00072 00073 00074 //##### declareOptions ###################################################### 00075 00076 void Matern1ARDKernel::declareOptions(OptionList& ol) 00077 { 00078 declareOption( 00079 ol, "isp_persistence", 00080 &Matern1ARDKernel::m_isp_persistence, 00081 OptionBase::buildoption, 00082 "Inverse softplus of the O-U persistence parameter. Default value =\n" 00083 "isp(1.0)."); 00084 00085 // Now call the parent class' declareOptions 00086 inherited::declareOptions(ol); 00087 } 00088 00089 00090 //##### build ############################################################### 00091 00092 void Matern1ARDKernel::build() 00093 { 00094 // ### Nothing to add here, simply calls build_ 00095 inherited::build(); 00096 build_(); 00097 } 00098 00099 00100 //##### build_ ############################################################## 00101 00102 void Matern1ARDKernel::build_() 00103 { 00104 // Ensure that we multiply in Kronecker terms 00105 inherited::m_default_value = 1.0; 00106 } 00107 00108 00109 //##### evaluate ############################################################ 00110 00111 real Matern1ARDKernel::evaluate(const Vec& x1, const Vec& x2) const 00112 { 00113 PLASSERT( x1.size() == x2.size() ); 00114 PLASSERT( !m_isp_input_sigma.size() || x1.size() == m_isp_input_sigma.size() ); 00115 00116 real gating_term = inherited::evaluate(x1,x2); 00117 if (fast_is_equal(gating_term, 0.0)) 00118 return 0.0; 00119 00120 if (x1.size() == 0) 00121 return softplus(m_isp_signal_sigma) / 00122 (2*softplus(m_isp_persistence)) * gating_term; 00123 00124 const real* px1 = x1.data(); 00125 const real* px2 = x2.data(); 00126 real sf = softplus(m_isp_signal_sigma); 00127 real persistence= softplus(m_isp_persistence); 00128 real expval = 0.0; 00129 00130 // Case where we have real ARD 00131 if (m_isp_input_sigma.size() > 0) { 00132 const real* pinpsig = m_isp_input_sigma.data(); 00133 for (int i=0, n=x1.size() ; i<n ; ++i) { 00134 real diff = *px1++ - *px2++; 00135 real absdiff = fabs(diff); 00136 expval += absdiff / softplus(m_isp_global_sigma + *pinpsig++); 00137 } 00138 } 00139 // No ARD 00140 else { 00141 real global_sigma = softplus(m_isp_global_sigma); 00142 for (int i=0, n=x1.size() ; i<n ; ++i) { 00143 real diff = *px1++ - *px2++; 00144 real absdiff = fabs(diff); 00145 expval += absdiff / global_sigma; 00146 } 00147 } 00148 00149 // Gate by Kronecker term 00150 return sf / (2. * persistence) * exp(-persistence * expval) * gating_term; 00151 } 00152 00153 00154 //##### computeGramMatrix ################################################### 00155 00156 void Matern1ARDKernel::computeGramMatrix(Mat K) const 00157 { 00158 PLASSERT( !m_isp_input_sigma.size() || dataInputsize() == m_isp_input_sigma.size() ); 00159 PLASSERT( K.size() == 0 || m_data_cache.size() > 0 ); // Ensure data cached OK 00160 00161 // Compute Kronecker gram matrix 00162 inherited::computeGramMatrix(K); 00163 00164 // Precompute some terms. Make sure that the input sigmas don't get too 00165 // small 00166 real sf = softplus(m_isp_signal_sigma); 00167 real persistence = softplus(m_isp_persistence); 00168 m_input_sigma.resize(dataInputsize()); 00169 softplusFloor(m_isp_global_sigma, 1e-6); 00170 m_input_sigma.fill(m_isp_global_sigma); // Still in ISP domain 00171 for (int i=0, n=m_input_sigma.size() ; i<n ; ++i) { 00172 if (m_isp_input_sigma.size() > 0) { 00173 softplusFloor(m_isp_input_sigma[i], 1e-6); 00174 m_input_sigma[i] += m_isp_input_sigma[i]; 00175 } 00176 m_input_sigma[i] = softplus(m_input_sigma[i]); 00177 } 00178 00179 // Compute Gram Matrix 00180 int l = data->length(); 00181 int m = K.mod(); 00182 int n = dataInputsize(); 00183 int cache_mod = m_data_cache.mod(); 00184 00185 real *data_start = &m_data_cache(0,0); 00186 real *Ki = K[0]; // Start of current row 00187 real *Kij; // Current element along row 00188 real *input_sigma_data = m_input_sigma.data(); 00189 real *xi = data_start; 00190 00191 for (int i=0 ; i<l ; ++i, xi += cache_mod, Ki+=m) 00192 { 00193 Kij = Ki; 00194 real *xj = data_start; 00195 00196 for (int j=0; j<=i; ++j, xj += cache_mod) { 00197 // Kernel evaluation per se 00198 real *x1 = xi; 00199 real *x2 = xj; 00200 real *p_inpsigma = input_sigma_data; 00201 real sum_wt = 0.0; 00202 int k = n; 00203 00204 // Use Duff's device to unroll the following loop: 00205 // while (k--) { 00206 // real diff = *x1++ - *x2++; 00207 // sum_wt += fabs(diff) / *p_inpsigma++; 00208 // } 00209 real diff; 00210 switch (k % 8) { 00211 case 0: do { diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00212 case 7: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00213 case 6: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00214 case 5: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00215 case 4: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00216 case 3: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00217 case 2: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00218 case 1: diff = *x1++ - *x2++; sum_wt += fabs(diff) / *p_inpsigma++; 00219 } while((k -= 8) > 0); 00220 } 00221 00222 // Multiplicatively update kernel matrix (already pre-filled with 00223 // Kronecker terms, or 1.0 if no Kronecker terms, as per build_). 00224 real Kij_cur = *Kij * sf / (2.*persistence) * exp(-persistence * sum_wt); 00225 *Kij++ = Kij_cur; 00226 } 00227 } 00228 if (cache_gram_matrix) { 00229 gram_matrix.resize(l,l); 00230 gram_matrix << K; 00231 gram_matrix_is_cached = true; 00232 } 00233 } 00234 00235 00236 //##### computeGramMatrixDerivative ######################################### 00237 00238 void Matern1ARDKernel::computeGramMatrixDerivative( 00239 Mat& KD, const string& kernel_param, real epsilon) const 00240 { 00241 static const string ISS("isp_signal_sigma"); 00242 static const string IGS("isp_global_sigma"); 00243 static const string IIS("isp_input_sigma["); 00244 static const string IPe("isp_persistence"); 00245 00246 if (kernel_param == ISS) { 00247 computeGramMatrixDerivIspSignalSigma(KD); 00248 00249 // computeGramMatrixDerivNV< 00250 // Matern1ARDKernel, 00251 // &Matern1ARDKernel::derivIspSignalSigma>(KD, this, -1); 00252 } 00253 /* 00254 else if (kernel_param == IGS) { 00255 computeGramMatrixDerivNV< 00256 Matern1ARDKernel, 00257 &Matern1ARDKernel::derivIspGlobalSigma>(KD, this, -1); 00258 } 00259 else if (string_begins_with(kernel_param, IIS) && 00260 kernel_param[kernel_param.size()-1] == ']') 00261 { 00262 int arg = tolong(kernel_param.substr( 00263 IIS.size(), kernel_param.size() - IIS.size() - 1)); 00264 PLASSERT( arg < m_isp_input_sigma.size() ); 00265 00266 computeGramMatrixDerivIspInputSigma(KD, arg); 00267 00268 } 00269 */ 00270 else 00271 inherited::computeGramMatrixDerivative(KD, kernel_param, epsilon); 00272 } 00273 00274 00275 //##### evaluate_all_i_x #################################################### 00276 00277 void Matern1ARDKernel::evaluate_all_i_x(const Vec& x, const Vec& k_xi_x, 00278 real squared_norm_of_x, int istart) const 00279 { 00280 evaluateAllIXNV<Matern1ARDKernel>(x, k_xi_x, istart); 00281 } 00282 00283 00284 //##### derivIspSignalSigma ################################################# 00285 00286 real Matern1ARDKernel::derivIspSignalSigma(int i, int j, int arg, real K) const 00287 { 00288 // (No longer used; see computeGramMatrixDerivIspInputSigma below) 00289 return K*sigmoid(m_isp_signal_sigma)/softplus(m_isp_signal_sigma); 00290 } 00291 00292 00293 //##### derivIspGlobalSigma ################################################# 00294 00295 real Matern1ARDKernel::derivIspGlobalSigma(int i, int j, int arg, real K) const 00296 { 00297 if (fast_is_equal(K,0.)) 00298 return 0.; 00299 00300 // The norm term inside the exponential may be accessed as Log(K/sf) 00301 real inner = pl_log(K / softplus(m_isp_signal_sigma)); 00302 return - K * inner * sigmoid(m_isp_global_sigma) / softplus(m_isp_global_sigma); 00303 00304 // Note: in the above expression for 'inner' there is the implicit 00305 // assumption that the input_sigma[i] are zero, which allows the 00306 // sigmoid/softplus term to be factored out of the norm summation. 00307 } 00308 00309 00310 //##### computeGramMatrixDerivIspSignalSigma ################################ 00311 00312 void Matern1ARDKernel::computeGramMatrixDerivIspSignalSigma(Mat& KD) const 00313 { 00314 int l = data->length(); 00315 KD.resize(l,l); 00316 PLASSERT_MSG( 00317 gram_matrix.width() == l && gram_matrix.length() == l, 00318 "To compute the derivative with respect to 'isp_signal_sigma', the\n" 00319 "Gram matrix must be precomputed and cached in Matern1ARDKernel."); 00320 00321 KD << gram_matrix; 00322 KD *= sigmoid(m_isp_signal_sigma)/softplus(m_isp_signal_sigma); 00323 } 00324 00325 00326 //##### computeGramMatrixDerivIspInputSigma ################################# 00327 00328 void Matern1ARDKernel::computeGramMatrixDerivIspInputSigma(Mat& KD, 00329 int arg) const 00330 { 00331 // Precompute some terms 00332 real input_sigma_arg = m_input_sigma[arg]; 00333 real input_sigma_sq = input_sigma_arg * input_sigma_arg; 00334 real input_sigmoid = sigmoid(m_isp_global_sigma + m_isp_input_sigma[arg]); 00335 00336 // Compute Gram Matrix derivative w.r.t. isp_input_sigma[arg] 00337 int l = data->length(); 00338 PLASSERT_MSG( 00339 gram_matrix.width() == l && gram_matrix.length() == l, 00340 "To compute the derivative with respect to 'isp_input_sigma[i]', the\n" 00341 "Gram matrix must be precomputed and cached in Matern1ARDKernel."); 00342 00343 // Variables that walk over the data matrix 00344 int cache_mod = m_data_cache.mod(); 00345 real *data_start = &m_data_cache(0,0); 00346 real *xi = data_start+arg; // Iterator on data rows 00347 00348 // Variables that walk over the gram cache 00349 int gram_cache_mod = gram_matrix.mod(); 00350 real *gram_cache_row = gram_matrix.data(); 00351 real *gram_cache_cur; 00352 00353 // Variables that walk over the kernel derivative matrix (KD) 00354 KD.resize(l,l); 00355 real* KDi = KD.data(); // Start of row i 00356 real* KDij; // Current element on row i 00357 int KD_mod = KD.mod(); 00358 00359 // Iterate on rows of derivative matrix 00360 for (int i=0 ; i<l ; ++i, xi += cache_mod, KDi += KD_mod, 00361 gram_cache_row += gram_cache_mod) 00362 { 00363 KDij = KDi; 00364 real *xj = data_start+arg; // Inner iterator on data rows 00365 gram_cache_cur = gram_cache_row; 00366 00367 // Iterate on columns of derivative matrix 00368 for (int j=0 ; j <= i 00369 ; ++j, xj += cache_mod, ++gram_cache_cur) 00370 { 00371 real diff = *xi - *xj; 00372 real sq_diff = diff * diff; 00373 real KD_cur = 0.5 * *gram_cache_cur * 00374 input_sigmoid * sq_diff / input_sigma_sq; 00375 00376 // Set into derivative matrix 00377 *KDij++ = KD_cur; 00378 } 00379 } 00380 } 00381 00382 00383 //##### makeDeepCopyFromShallowCopy ######################################### 00384 00385 void Matern1ARDKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies) 00386 { 00387 inherited::makeDeepCopyFromShallowCopy(copies); 00388 } 00389 00390 } // end of namespace PLearn 00391 00392 00393 /* 00394 Local Variables: 00395 mode:c++ 00396 c-basic-offset:4 00397 c-file-style:"stroustrup" 00398 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00399 indent-tabs-mode:nil 00400 fill-column:79 00401 End: 00402 */ 00403 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :