PLearn 0.1
Matern1ARDKernel.cc
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines