SHOGUN  6.1.3
NumericalVGLikelihood.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) The Shogun Machine Learning Toolbox
3  * Written (w) 2014 Wu Lin
4  * All rights reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright notice, this
10  * list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright notice,
12  * this list of conditions and the following disclaimer in the documentation
13  * and/or other materials provided with the distribution.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
19  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  *
26  * The views and conclusions contained in the software and documentation are those
27  * of the authors and should not be interpreted as representing official policies,
28  * either expressed or implied, of the Shogun Development Team.
29  *
30  * Code adapted from
31  * and Gaussian Process Machine Learning Toolbox
32  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
33  *
34  * and the reference paper is
35  * Nickisch, Hannes, and Carl Edward Rasmussen.
36  * "Approximations for Binary Gaussian Process Classification."
37  * Journal of Machine Learning Research 9.10 (2008).
38  *
39  * This code specifically adapted from function in infKL.m
40  */
41 
43 
45 #include <shogun/lib/SGVector.h>
46 #include <shogun/lib/SGMatrix.h>
47 #ifdef USE_GPL_SHOGUN
48 #include <shogun/mathematics/Integration.h>
49 #endif //USE_GPL_SHOGUN
50 
51 using namespace Eigen;
52 
53 namespace shogun
54 {
55 
56 CNumericalVGLikelihood::CNumericalVGLikelihood()
58 {
59  init();
60 }
61 
63 {
64 }
65 
66 void CNumericalVGLikelihood::init()
67 {
68  SG_ADD(&m_log_lam, "log_lam",
69  "The result of used for computing variational expection\n",
71 
72  SG_ADD(&m_xgh, "xgh",
73  "Gaussian-Hermite quadrature base points (abscissas)\n",
75 
76  SG_ADD(&m_wgh, "wgh",
77  "Gaussian-Hermite quadrature weight factors\n",
79 
80  SG_ADD(&m_GHQ_N, "GHQ_N",
81  "The number of Gaussian-Hermite quadrature point\n",
83 
84  SG_ADD(&m_is_init_GHQ, "is_init_GHQ",
85  "Whether Gaussian-Hermite quadrature points are initialized or not\n",
87  m_GHQ_N=20;
88  m_is_init_GHQ=false;
89 
90 }
91 
93 {
94  REQUIRE(n>0, "The number (%d) of Gaussian Hermite point should be positive\n",n);
95  if (m_GHQ_N!=n)
96  {
97  m_GHQ_N=n;
98  m_is_init_GHQ=false;
99  }
100 }
101 
103  const TParameter* param) const
104 {
105  REQUIRE(param, "Param is required (param should not be NULL)\n");
106  REQUIRE(param->m_name, "Param name is required (param->m_name should not be NULL)\n");
107  if (!(strcmp(param->m_name, "mu") && strcmp(param->m_name, "sigma2")))
108  return SGVector<float64_t> ();
109 
111  res.zero();
112  Map<VectorXd> eigen_res(res.vector, res.vlen);
113 
114  //ll = ll + w(i)*lp;
115  CLabels* lab=NULL;
116 
117  if (supports_binary())
118  lab=new CBinaryLabels(m_lab);
119  else if (supports_regression())
120  lab=new CRegressionLabels(m_lab);
121 
122  for (index_t cidx = 0; cidx < m_log_lam.num_cols; cidx++)
123  {
124  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
125  SGVector<float64_t> lp=get_first_derivative(lab, tmp, param);
126  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
127  eigen_res+=eigen_lp*m_wgh[cidx];
128  }
129 
130  SG_UNREF(lab);
131 
132  return res;
133 }
134 
136 {
138  res.zero();
139  Map<VectorXd> eigen_res(res.vector, res.vlen);
140 
141  //ll = ll + w(i)*lp;
142  CLabels* lab=NULL;
143 
144  if (supports_binary())
145  lab=new CBinaryLabels(m_lab);
146  else if (supports_regression())
147  lab=new CRegressionLabels(m_lab);
148 
149  for (index_t cidx = 0; cidx < m_log_lam.num_cols; cidx++)
150  {
151  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
153  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
154  eigen_res+=eigen_lp*m_wgh[cidx];
155  }
156 
157  SG_UNREF(lab);
158 
159  return res;
160 }
161 
163  const TParameter* param) const
164 {
165  //based on the likKL(v, lik, varargin) function in infKL.m
166 
167  //compute gradient using numerical integration
168  REQUIRE(param, "Param is required (param should not be NULL)\n");
169  REQUIRE(param->m_name, "Param name is required (param->m_name should not be NULL)\n");
170  //We take the derivative wrt to param. Only mu or sigma2 can be the param
171  REQUIRE(!(strcmp(param->m_name, "mu") && strcmp(param->m_name, "sigma2")),
172  "Can't compute derivative of the variational expection ",
173  "of log LogitLikelihood using numerical integration ",
174  "wrt %s.%s parameter. The function only accepts mu and sigma2 as parameter\n",
175  get_name(), param->m_name);
176 
178  res.zero();
179  Map<VectorXd> eigen_res(res.vector, res.vlen);
180 
181  Map<VectorXd> eigen_v(m_s2.vector, m_s2.vlen);
182 
183  CLabels* lab=NULL;
184 
185  if (supports_binary())
186  lab=new CBinaryLabels(m_lab);
187  else if (supports_regression())
188  lab=new CRegressionLabels(m_lab);
189 
190  if (strcmp(param->m_name, "mu")==0)
191  {
192  //Compute the derivative wrt mu
193 
194  //df = df + w(i)*(dlp);
195  for (index_t cidx=0; cidx<m_log_lam.num_cols; cidx++)
196  {
197  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
199  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
200  eigen_res+=eigen_dlp*m_wgh[cidx];
201  }
202  }
203  else
204  {
205  //Compute the derivative wrt sigma2
206 
207  //ai = t(i)./(2*sv+eps); dvi = dlp.*ai; dv = dv + w(i)*dvi;
208  VectorXd eigen_sv=eigen_v.array().sqrt().matrix();
209  const float64_t EPS=2.2204e-16;
210 
211  for (index_t cidx=0; cidx<m_log_lam.num_cols; cidx++)
212  {
213  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
215  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
216  eigen_res+=((m_wgh[cidx]*0.5*m_xgh[cidx])*eigen_dlp.array()/(eigen_sv.array()+EPS)).matrix();
217  }
218  }
219 
220  SG_UNREF(lab);
221 
222  return res;
223 }
224 
225 
227  SGVector<float64_t> s2, const CLabels* lab)
228 {
229  bool status = true;
231 
232  if (status)
233  {
234  if (supports_binary())
235  {
237  "Labels must be type of CBinaryLabels\n");
238  }
239  else
240  {
241  if (supports_regression())
242  {
244  "Labels must be type of CRegressionLabels\n");
245  }
246  else
247  SG_ERROR("Unsupported Label type\n");
248  }
249 
250  if (supports_binary())
251  m_lab=((CBinaryLabels*)lab)->get_labels();
252  else
253  m_lab=((CRegressionLabels*)lab)->get_labels();
254 
255  if (!m_is_init_GHQ)
256  {
257  m_xgh=SGVector<float64_t>(m_GHQ_N);
258  m_wgh=SGVector<float64_t>(m_GHQ_N);
259 #ifdef USE_GPL_SHOGUN
260  CIntegration::generate_gauher(m_xgh, m_wgh);
261 #else
263 #endif //USE_GPL_SHOGUN
264  m_is_init_GHQ=true;
265  }
266 
267  precompute();
268 
269  }
270 
271  return status;
272 }
273 
274 void CNumericalVGLikelihood::precompute()
275 {
276  //samples-by-abscissas
277  m_log_lam=SGMatrix<float64_t>(m_s2.vlen, m_xgh.vlen);
278 
279  Map<MatrixXd> eigen_log_lam(m_log_lam.matrix, m_log_lam.num_rows, m_log_lam.num_cols);
280  Map<VectorXd> eigen_v(m_s2.vector, m_s2.vlen);
281  Map<VectorXd> eigen_f(m_mu.vector, m_mu.vlen);
282  Map<VectorXd> eigen_t(m_xgh.vector, m_xgh.vlen);
283 
284  VectorXd eigen_sv=eigen_v.array().sqrt().matrix();
285  //varargin{3} = f + sv*t(i); % coordinate transform of the quadrature points
286  eigen_log_lam = (
287  (eigen_t.replicate(1, eigen_v.rows()).array().transpose().colwise()
288  *eigen_sv.array()).colwise()+eigen_f.array()
289  ).matrix();
290 }
291 
292 } /* namespace shogun */
293 
virtual bool set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
virtual ELabelType get_label_type() const =0
binary labels +1/-1
Definition: LabelTypes.h:18
virtual const char * get_name() const
Real Labels are real-valued labels.
int32_t index_t
Definition: common.h:72
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
real valued labels (e.g. for regression, classifier outputs)
Definition: LabelTypes.h:22
The variational Gaussian Likelihood base class. The variational distribution is Gaussian.
Definition: SGMatrix.h:25
parameter struct
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
virtual SGVector< float64_t > get_first_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
#define SG_GPL_ONLY
Definition: SGIO.h:139
double float64_t
Definition: common.h:60
virtual bool set_variational_distribution(SGVector< float64_t > mu, SGVector< float64_t > s2, const CLabels *lab)
virtual void set_GHQ_number(index_t n)
virtual SGVector< float64_t > get_variational_expection()
index_t num_rows
Definition: SGMatrix.h:495
virtual SGVector< float64_t > get_log_probability_f(const CLabels *lab, SGVector< float64_t > func) const
virtual SGVector< float64_t > get_log_probability_derivative_f(const CLabels *lab, SGVector< float64_t > func, index_t i) const
index_t num_cols
Definition: SGMatrix.h:497
#define SG_UNREF(x)
Definition: SGObject.h:53
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Binary Labels for binary classification.
Definition: BinaryLabels.h:37
#define SG_ADD(...)
Definition: SGObject.h:93
virtual SGVector< float64_t > get_first_derivative_wrt_hyperparameter(const TParameter *param) const
T * get_column_vector(index_t col) const
Definition: SGMatrix.h:144
virtual SGVector< float64_t > get_variational_first_derivative(const TParameter *param) const
index_t vlen
Definition: SGVector.h:571

SHOGUN Machine Learning Toolbox - Documentation