PLearn 0.1
RationalQuadraticARDKernel.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // RationalQuadraticARDKernel.cc
00004 //
00005 // Copyright (C) 2006-2007 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 "RationalQuadraticARDKernel.h"
00041 
00042 namespace PLearn {
00043 using namespace std;
00044 
00045 PLEARN_IMPLEMENT_OBJECT(
00046     RationalQuadraticARDKernel,
00047     "Rational-Quadratic kernel that can be used for Automatic Relevance Determination",
00048     "This kernel can be interpreted as an infinite mixture of\n"
00049     "SquaredExponentialARDKernel (with different characteristic length-scales),\n"
00050     "allowing a greater variety of \"interesting\" functions to be generated.\n"
00051     "Similar to C.E. Rasmussen's GPML code (see http://www.gaussianprocess.org),\n"
00052     "this kernel is specified as:\n"
00053     "\n"
00054     "  k(x,y) = sf * [1 + (sum_i (x_i - y_i)^2 / w_i)/(2*alpha)]^(-alpha) * k_kron(x,y)\n"
00055     "\n"
00056     "where sf is softplus(isp_signal_sigma), w_i is softplus(isp_global_sigma +\n"
00057     "isp_input_sigma[i]), and k_kron(x,y) is the result of the\n"
00058     "KroneckerBaseKernel evaluation, or 1.0 if there are no Kronecker terms.\n"
00059     "Note that since the Kronecker terms are incorporated multiplicatively, the\n"
00060     "very presence of the term associated to this kernel can be gated by the\n"
00061     "value of some input variable(s) (that are incorporated within one or more\n"
00062     "Kronecker terms).\n"
00063     "\n"
00064     "Note that contrarily to previous versions that incorporated IID noise and\n"
00065     "Kronecker terms ADDITIVELY, this version does not add any noise at all (and\n"
00066     "as explained above incorporates the Kronecker terms multiplicatively).  For\n"
00067     "best results, especially with moderately noisy data, IT IS IMPERATIVE to\n"
00068     "use whis kernel within a SummationKernel in conjunction with an\n"
00069     "IIDNoiseKernel, as follows (e.g. within a GaussianProcessRegressor):\n"
00070     "\n"
00071     "    kernel = SummationKernel(terms = [ RationalQuadraticARDKernel(),\n"
00072     "                                       IIDNoiseKernel() ] )\n"
00073     "\n"
00074     "In order to make its operations more robust when used with unconstrained\n"
00075     "optimization of hyperparameters, all hyperparameters of this kernel are\n"
00076     "specified in the inverse softplus domain.  See IIDNoiseKernel for more\n"
00077     "explanations.\n"
00078     );
00079 
00080 
00081 RationalQuadraticARDKernel::RationalQuadraticARDKernel()
00082     : m_isp_alpha(0.0)
00083 { }
00084 
00085 
00086 //#####  declareOptions  ######################################################
00087 
00088 void RationalQuadraticARDKernel::declareOptions(OptionList& ol)
00089 {
00090     declareOption(
00091         ol, "isp_alpha",
00092         &RationalQuadraticARDKernel::m_isp_alpha,
00093         OptionBase::buildoption,
00094         "Inverse softplus of the alpha parameter in the rational-quadratic kernel.\n"
00095         "Default value=0.0");
00096 
00097     // Now call the parent class' declareOptions
00098     inherited::declareOptions(ol);
00099 }
00100 
00101 
00102 //#####  build  ###############################################################
00103 
00104 void RationalQuadraticARDKernel::build()
00105 {
00106     // ### Nothing to add here, simply calls build_
00107     inherited::build();
00108     build_();
00109 }
00110 
00111 
00112 //#####  build_  ##############################################################
00113 
00114 void RationalQuadraticARDKernel::build_()
00115 {
00116     // Ensure that we multiply in Kronecker terms
00117     inherited::m_default_value = 1.0;
00118 }
00119 
00120 
00121 //#####  makeDeepCopyFromShallowCopy  #########################################
00122 
00123 void RationalQuadraticARDKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00124 {
00125     inherited::makeDeepCopyFromShallowCopy(copies);
00126     deepCopyField(m_pow_minus_alpha_minus_1, copies);
00127 }
00128 
00129 
00130 //#####  evaluate  ############################################################
00131 
00132 real RationalQuadraticARDKernel::evaluate(const Vec& x1, const Vec& x2) const
00133 {
00134     PLASSERT( x1.size() == x2.size() );
00135     PLASSERT( !m_isp_input_sigma.size() || x1.size() == m_isp_input_sigma.size() );
00136 
00137     real gating_term = inherited::evaluate(x1,x2);
00138     if (fast_is_equal(gating_term, 0.0))
00139         return 0.0;
00140     
00141     if (x1.size() == 0)
00142         return softplus(m_isp_signal_sigma) * gating_term;
00143     
00144     const real* px1 = x1.data();
00145     const real* px2 = x2.data();
00146     real sf         = softplus(m_isp_signal_sigma);
00147     real alpha      = softplus(m_isp_alpha);
00148     real sum_wt     = 0.0;
00149     
00150     if (m_isp_input_sigma.size() > 0) {
00151         const real* pinpsig = m_isp_input_sigma.data();
00152         for (int i=0, n=x1.size() ; i<n ; ++i) {
00153             real diff   = *px1++ - *px2++;
00154             real sqdiff = diff * diff;
00155             sum_wt     += sqdiff / softplus(m_isp_global_sigma + *pinpsig++);
00156         }
00157     }
00158     else {
00159         real global_sigma = softplus(m_isp_global_sigma);
00160         for (int i=0, n=x1.size() ; i<n ; ++i) {
00161             real diff   = *px1++ - *px2++;
00162             real sqdiff = diff * diff;
00163             sum_wt     += sqdiff / global_sigma;
00164         }
00165     }
00166 
00167     // Gate by Kronecker term
00168     return sf * pow(1 + sum_wt / (real(2.)*alpha), -alpha) * gating_term;
00169 }
00170 
00171 
00172 //#####  evaluate_all_i_x  ####################################################
00173 
00174 void RationalQuadraticARDKernel::evaluate_all_i_x(const Vec& x1, const Vec& k_xi_x,
00175                                                   real , int istart) const
00176 {
00177     if (x1.size() == 0) {
00178         k_xi_x.fill(0.0);
00179         return;
00180     }
00181  
00182     // Precompute some terms
00183     real sf    = softplus(m_isp_signal_sigma);
00184     real alpha = softplus(m_isp_alpha);
00185     m_input_sigma.resize(dataInputsize());
00186     for (int i=0, n=m_input_sigma.size() ; i<n ; ++i)
00187         m_input_sigma[i] = softplus(m_isp_global_sigma + m_isp_input_sigma[i]);
00188     
00189     const real* px1_start = x1.data();
00190     const real* pinpsig_start = m_input_sigma.data();
00191     int i_max = min(istart + k_xi_x.size(), data->length());
00192     int j = 0;
00193     for (int i=istart ; i<i_max ; ++i, ++j) {
00194         Vec* train_row = dataRow(i);
00195         const real* px2 = train_row->data();
00196         const real* px1 = px1_start;
00197     
00198         real gating_term = inherited::evaluate(x1,*train_row);
00199         if (fast_is_equal(gating_term, 0.0)) {
00200             k_xi_x[j] = 0.0;
00201             continue;
00202         }
00203     
00204         real sum_wt     = 0.0;
00205         const real* pinpsig = pinpsig_start;
00206         for (int i=0, n=x1.size() ; i<n ; ++i) {
00207             real diff   = *px1++ - *px2++;
00208             real sqdiff = diff * diff;
00209             sum_wt     += sqdiff / *pinpsig++;
00210         }
00211 
00212         // Gate by Kronecker term
00213         k_xi_x[j] = sf * pow(1 + sum_wt / (real(2.)*alpha), -alpha) * gating_term;
00214     }
00215 }
00216 
00217 
00218 //#####  computeGramMatrix  ###################################################
00219 
00220 void RationalQuadraticARDKernel::computeGramMatrix(Mat K) const
00221 {
00222     PLASSERT( !m_isp_input_sigma.size() || dataInputsize() == m_isp_input_sigma.size() );
00223     PLASSERT( K.size() == 0 || m_data_cache.size() > 0 );  // Ensure data cached OK
00224 
00225     // Compute Kronecker gram matrix.  No need to cache it.
00226     inherited::computeGramMatrix(K);
00227 
00228     // Precompute some terms. Make sure that the alpha and input sigmas don't
00229     // get too small
00230     real sf    = softplus(m_isp_signal_sigma);
00231     real alpha = softplusFloor(m_isp_alpha, 1e-6);
00232     m_input_sigma.resize(dataInputsize());
00233     softplusFloor(m_isp_global_sigma, 1e-6);
00234     m_input_sigma.fill(m_isp_global_sigma);  // Still in ISP domain
00235     for (int i=0, n=m_input_sigma.size() ; i<n ; ++i) {
00236         if (m_isp_input_sigma.size() > 0) {
00237             softplusFloor(m_isp_input_sigma[i], 1e-6);
00238             m_input_sigma[i] += m_isp_input_sigma[i];
00239         }
00240         m_input_sigma[i] = softplus(m_input_sigma[i]);
00241     }
00242     
00243     // Prepare the cache for the pow terms
00244     m_pow_minus_alpha_minus_1.resize(K.length(), K.width());
00245     int   pow_cache_mod = m_pow_minus_alpha_minus_1.mod();
00246     real* pow_cache_row = m_pow_minus_alpha_minus_1.data();
00247     
00248     // Compute Gram Matrix
00249     int  l = data->length();
00250     int  m = K.mod();
00251     int  n = dataInputsize();
00252     int  cache_mod = m_data_cache.mod();
00253 
00254     real *data_start = &m_data_cache(0,0);
00255     real *Ki = K[0];                         // Start of current row
00256     real *Kij;                               // Current element along row
00257     real *input_sigma_data = m_input_sigma.data();
00258     real *xi = data_start;
00259     
00260     for (int i=0 ; i<l
00261              ; ++i, xi += cache_mod, pow_cache_row+=pow_cache_mod, Ki+=m)
00262     {
00263         Kij = Ki;
00264         real *xj = data_start;
00265         real *pow_cache_cur = pow_cache_row;
00266 
00267         // This whole loop can be optimized further when a Kronecker term is 0
00268         for (int j=0; j<=i; ++j, xj += cache_mod) {
00269             // Kernel evaluation per se
00270             real *x1 = xi;
00271             real *x2 = xj;
00272             real *p_inpsigma = input_sigma_data;
00273             real sum_wt = 0.0;
00274             int  k = n;
00275 
00276             // Use Duff's device to unroll the following loop:
00277             //     while (k--) {
00278             //         real diff = *x1++ - *x2++;
00279             //         sum_wt += (diff * diff) / *p_inpsigma++;
00280             //     }
00281             real diff;
00282             switch (k % 8) {
00283             case 0: do { diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00284             case 7:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00285             case 6:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00286             case 5:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00287             case 4:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00288             case 3:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00289             case 2:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00290             case 1:      diff = *x1++ - *x2++; sum_wt += (diff*diff) / *p_inpsigma++;
00291                        } while((k -= 8) > 0);
00292             }
00293 
00294             // Multiplicatively update kernel matrix (already pre-filled with
00295             // Kronecker terms, or 1.0 if no Kronecker terms, as per build_).
00296             real inner_pow   = 1 + sum_wt / (2.*alpha);
00297             real pow_alpha   = pow(inner_pow, -alpha);
00298             real Kij_cur     = *Kij * sf * pow_alpha;       // Mind *Kij here
00299             *pow_cache_cur++ = Kij_cur / inner_pow;
00300             *Kij++           = Kij_cur;
00301         }
00302     }
00303     if (cache_gram_matrix) {
00304         gram_matrix.resize(l,l);
00305         gram_matrix << K;
00306         gram_matrix_is_cached = true;
00307     }
00308 }
00309 
00310 
00311 //#####  computeGramMatrixDerivative  #########################################
00312 
00313 void RationalQuadraticARDKernel::computeGramMatrixDerivative(
00314     Mat& KD, const string& kernel_param, real epsilon) const
00315 {
00316     static const string ISS("isp_signal_sigma");
00317     static const string IGS("isp_global_sigma");
00318     static const string IIS("isp_input_sigma[");
00319     static const string IAL("isp_alpha");
00320 
00321     if (kernel_param == ISS) {
00322         computeGramMatrixDerivIspSignalSigma(KD);
00323 
00324         // computeGramMatrixDerivNV<
00325         //     RationalQuadraticARDKernel,
00326         //     &RationalQuadraticARDKernel::derivIspSignalSigma>(KD, this, -1);
00327     }
00328     else if (kernel_param == IGS) {
00329         computeGramMatrixDerivNV<
00330             RationalQuadraticARDKernel,
00331             &RationalQuadraticARDKernel::derivIspGlobalSigma>(KD, this, -1);
00332     }
00333     else if (string_begins_with(kernel_param, IIS) &&
00334              kernel_param[kernel_param.size()-1] == ']')
00335     {
00336         int arg = tolong(kernel_param.substr(
00337                              IIS.size(), kernel_param.size() - IIS.size() - 1));
00338         PLASSERT( arg < m_isp_input_sigma.size() );
00339 
00340         computeGramMatrixDerivIspInputSigma(KD, arg);
00341 
00342         // computeGramMatrixDerivNV<
00343         //     RationalQuadraticARDKernel,
00344         //     &RationalQuadraticARDKernel::derivIspInputSigma>(KD, this, arg);
00345     }
00346     else if (kernel_param == IAL) {
00347         computeGramMatrixDerivIspAlpha(KD);
00348 
00349         // computeGramMatrixDerivNV<
00350         //     RationalQuadraticARDKernel,
00351         //     &RationalQuadraticARDKernel::derivIspAlpha>(KD, this, -1);
00352     }
00353     else
00354         inherited::computeGramMatrixDerivative(KD, kernel_param, epsilon);
00355 
00356     // Compare against finite differences
00357     // Mat KD1;
00358     // Kernel::computeGramMatrixDerivative(KD1, kernel_param, epsilon);
00359     // cerr << "Kernel hyperparameter: " << kernel_param << endl;
00360     // cerr << "Analytic derivative (200th row):" << endl
00361     //      << KD(200) << endl
00362     //      << "Finite differences:" << endl
00363     //      << KD1(200) << endl;
00364 }
00365 
00366 
00367 //#####  derivIspSignalSigma  #################################################
00368 
00369 real RationalQuadraticARDKernel::derivIspSignalSigma(int i, int j, int arg, real K) const
00370 {
00371     return K*sigmoid(m_isp_signal_sigma)/softplus(m_isp_signal_sigma);
00372 }
00373 
00374 
00375 //#####  derivIspGlobalSigma  #################################################
00376 
00377 real RationalQuadraticARDKernel::derivIspGlobalSigma(int i, int j, int arg, real K) const
00378 {
00379     // The rational quadratic gives us:
00380     //     K = s * k^(-alpha) * kron
00381     // where kron is 0 or 1.  Rederive the value of k == (K/s)^(-1/alpha)
00382     if (fast_is_equal(K, 0.))
00383         return 0.;
00384     real alpha = softplus(m_isp_alpha);
00385     real k     = pow(K / softplus(m_isp_signal_sigma), real(-1.) / alpha);
00386     real inner = (k - 1) * alpha * sigmoid(m_isp_global_sigma) / softplus(m_isp_global_sigma);
00387     return (K / k) * inner;
00388 
00389     // Note: in the above expression for 'inner' there is the implicit
00390     // assumption that the input_sigma[i] are zero, which allows the
00391     // sigmoid/softplus term to be factored out of the norm summation.
00392 }
00393 
00394 
00395 //#####  derivIspInputSigma  ##################################################
00396 
00397 // This function computes the derivative element-wise.  The function actually
00398 // used now is computeGramMatrixDerivIspInputSigma, which computes the whole
00399 // matrix much faster.
00400 real RationalQuadraticARDKernel::derivIspInputSigma(int i, int j, int arg, real K) const
00401 {
00402     // The rational quadratic gives us:
00403     //     K = s * k^(-alpha) * kron
00404     // where kron is 0 or 1.  Rederive the value of k == (K/s)^(-1/alpha)
00405     if (fast_is_equal(K, 0.))
00406         return 0.;
00407     real alpha   = softplus(m_isp_alpha);
00408     Vec& row_i   = *dataRow(i);
00409     Vec& row_j   = *dataRow(j);
00410     real k       = pow(K / softplus(m_isp_signal_sigma), real(-1.) / alpha);
00411     real diff    = row_i[arg] - row_j[arg];
00412     real sq_diff = diff * diff;
00413     real inner   = m_isp_global_sigma + m_isp_input_sigma[arg];
00414     real sig_inn = sigmoid(inner);
00415     real spl_inn = softplus(inner);
00416     return 0.5 * (K / k) * sig_inn * sq_diff / (spl_inn * spl_inn);
00417 }
00418 
00419 
00420 //#####  derivIspAlpha  #######################################################
00421 
00422 real RationalQuadraticARDKernel::derivIspAlpha(int i, int j, int arg, real K) const
00423 {
00424     // The rational quadratic gives us:
00425     //     K = s * k^(-alpha) * kron
00426     // where kron is 0 or 1.  Rederive the value of k == (K/s)^(-1/alpha)
00427     if (fast_is_equal(K, 0.))
00428         return 0.;
00429     real alpha = softplus(m_isp_alpha);
00430     real k     = pow(K / softplus(m_isp_signal_sigma), real(-1.) / alpha);
00431     return sigmoid(m_isp_alpha) * K * (1 - pl_log(k) - 1 / k);
00432 }
00433 
00434 
00435 //#####  computeGramMatrixDerivIspSignalSigma  ################################
00436 
00437 void RationalQuadraticARDKernel::computeGramMatrixDerivIspSignalSigma(Mat& KD) const
00438 {
00439     int l = data->length();
00440     KD.resize(l,l);
00441     PLASSERT_MSG(
00442         gram_matrix.width() == l && gram_matrix.length() == l,
00443         "To compute the derivative with respect to 'isp_signal_sigma', the\n"
00444         "Gram matrix must be precomputed and cached in SquaredExponentialARDKernel.");
00445     
00446     KD << gram_matrix;
00447     KD *= sigmoid(m_isp_signal_sigma)/softplus(m_isp_signal_sigma);
00448 }
00449 
00450 
00451 //#####  computeGramMatrixDerivIspInputSigma  #################################
00452 
00453 void RationalQuadraticARDKernel::computeGramMatrixDerivIspInputSigma(Mat& KD,
00454                                                                      int arg) const
00455 {
00456     // Precompute some terms
00457     real input_sigma_arg = m_input_sigma[arg];
00458     real input_sigma_sq  = input_sigma_arg * input_sigma_arg;
00459     real input_sigmoid   = sigmoid(m_isp_global_sigma + m_isp_input_sigma[arg]);
00460     
00461     // Compute Gram Matrix derivative w.r.t. isp_input_sigma[arg]
00462     int  l = data->length();
00463 
00464     // Variables that walk over the data matrix
00465     int  cache_mod = m_data_cache.mod();
00466     real *data_start = &m_data_cache(0,0);
00467     real *xi = data_start+arg;               // Iterator on data rows
00468 
00469     // Variables that walk over the pow cache
00470     int   pow_cache_mod = m_pow_minus_alpha_minus_1.mod();
00471     real *pow_cache_row = m_pow_minus_alpha_minus_1.data();
00472     real *pow_cache_cur;
00473     
00474     // Variables that walk over the kernel derivative matrix (KD)
00475     KD.resize(l,l);
00476     real* KDi = KD.data();                   // Start of row i
00477     real* KDij;                              // Current element on row i
00478     int   KD_mod = KD.mod();
00479 
00480     // Iterate on rows of derivative matrix
00481     for (int i=0 ; i<l ; ++i, xi += cache_mod, KDi += KD_mod,
00482              pow_cache_row += pow_cache_mod)
00483     {
00484         KDij = KDi;
00485         real *xj  = data_start+arg;           // Inner iterator on data rows
00486         pow_cache_cur = pow_cache_row;
00487 
00488         // Iterate on columns of derivative matrix
00489         for (int j=0 ; j <= i
00490                  ; ++j, xj += cache_mod, ++pow_cache_cur)
00491         {
00492             real diff    = *xi - *xj;
00493             real sq_diff = diff * diff;
00494             real KD_cur  = 0.5 * *pow_cache_cur *
00495                            input_sigmoid * sq_diff / input_sigma_sq;
00496 
00497             // Set into derivative matrix
00498             *KDij++ = KD_cur;
00499         }
00500     }
00501 }
00502 
00503 
00504 //#####  computeGramMatrixDerivIspAlpha  ######################################
00505 
00506 void RationalQuadraticARDKernel::computeGramMatrixDerivIspAlpha(Mat& KD) const
00507 {
00508     // Precompute some terms
00509     real alpha_sigmoid = sigmoid(m_isp_alpha);
00510     
00511     // Compute Gram Matrix derivative w.r.t. isp_alpha
00512     int  l     = data->length();
00513     PLASSERT_MSG(
00514         gram_matrix.width() == l && gram_matrix.length() == l,
00515         "To compute the derivative with respect to 'isp_alpha', the\n"
00516         "Gram matrix must be precomputed and cached in RationalQuadraticARDKernel.");
00517 
00518     // Variables that walk over the pre-computed kernel matrix (K) 
00519     int  k_mod = gram_matrix.mod();
00520     real *Ki = &gram_matrix(0,0);            // Current row of kernel matrix
00521     real *Kij;                               // Current element of kernel matrix
00522 
00523     // Variables that walk over the pow cache
00524     int   pow_cache_mod = m_pow_minus_alpha_minus_1.mod();
00525     real *pow_cache_row = m_pow_minus_alpha_minus_1.data();
00526     real *pow_cache_cur;
00527 
00528     // Variables that walk over the kernel derivative matrix (KD)
00529     KD.resize(l,l);
00530     real* KDi = KD.data();                   // Start of row i
00531     real* KDij;                              // Current element on row i
00532     int   KD_mod = KD.mod();
00533 
00534     // Iterate on rows of derivative matrix
00535     for (int i=0 ; i<l ; ++i, Ki += k_mod,
00536              KDi += KD_mod, pow_cache_row += pow_cache_mod)
00537     {
00538         Kij  = Ki;
00539         KDij = KDi;
00540         pow_cache_cur   = pow_cache_row;
00541 
00542         // Iterate on columns of derivative matrix
00543         for (int j=0 ; j <= i ; ++j, ++Kij, ++pow_cache_cur)
00544         {
00545             real pow_cur = *pow_cache_cur;
00546             if (fast_is_equal(pow_cur, 0)) 
00547                 *KDij++ = 0.;
00548             else {
00549                 real K      = *Kij;
00550                 real k      = K / pow_cur;
00551                 real KD_cur = alpha_sigmoid * K * (1 - pl_log(k) - 1/k);
00552             
00553                 // Set into derivative matrix
00554                 *KDij++ = KD_cur;
00555             }
00556         }
00557     }
00558 }
00559 
00560 
00561 } // end of namespace PLearn
00562 
00563 
00564 /*
00565   Local Variables:
00566   mode:c++
00567   c-basic-offset:4
00568   c-file-style:"stroustrup"
00569   c-file-offsets:((innamespace . 0)(inline-open . 0))
00570   indent-tabs-mode:nil
00571   fill-column:79
00572   End:
00573 */
00574 // 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