SHOGUN  4.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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>
48 
49 using namespace Eigen;
50 
51 namespace shogun
52 {
53 
54 CNumericalVGLikelihood::CNumericalVGLikelihood()
56 {
57  init();
58 }
59 
61 {
62 }
63 
64 void CNumericalVGLikelihood::init()
65 {
66  SG_ADD(&m_log_lam, "log_lam",
67  "The result of used for computing variational expection\n",
69 
70  SG_ADD(&m_xgh, "xgh",
71  "Gaussian-Hermite quadrature base points (abscissas)\n",
73 
74  SG_ADD(&m_wgh, "wgh",
75  "Gaussian-Hermite quadrature weight factors\n",
77 
78  SG_ADD(&m_GHQ_N, "GHQ_N",
79  "The number of Gaussian-Hermite quadrature point\n",
81 
82  SG_ADD(&m_is_init_GHQ, "is_init_GHQ",
83  "Whether Gaussian-Hermite quadrature points are initialized or not\n",
85  m_GHQ_N=20;
86  m_is_init_GHQ=false;
87 
88 }
89 
91 {
92  REQUIRE(n>0, "The number (%d) of Gaussian Hermite point should be positive\n",n);
93  if (m_GHQ_N!=n)
94  {
95  m_GHQ_N=n;
96  m_is_init_GHQ=false;
97  }
98 }
99 
101  const TParameter* param) const
102 {
103  REQUIRE(param, "Param is required (param should not be NULL)\n");
104  REQUIRE(param->m_name, "Param name is required (param->m_name should not be NULL)\n");
105  if (!(strcmp(param->m_name, "mu") && strcmp(param->m_name, "sigma2")))
106  return SGVector<float64_t> ();
107 
109  res.zero();
110  Map<VectorXd> eigen_res(res.vector, res.vlen);
111 
112  //ll = ll + w(i)*lp;
113  CLabels* lab=NULL;
114 
115  if (supports_binary())
116  lab=new CBinaryLabels(m_lab);
117  else if (supports_regression())
118  lab=new CRegressionLabels(m_lab);
119 
120  for (index_t cidx = 0; cidx < m_log_lam.num_cols; cidx++)
121  {
122  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
123  SGVector<float64_t> lp=get_first_derivative(lab, tmp, param);
124  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
125  eigen_res+=eigen_lp*m_wgh[cidx];
126  }
127 
128  SG_UNREF(lab);
129 
130  return res;
131 }
132 
134 {
136  res.zero();
137  Map<VectorXd> eigen_res(res.vector, res.vlen);
138 
139  //ll = ll + w(i)*lp;
140  CLabels* lab=NULL;
141 
142  if (supports_binary())
143  lab=new CBinaryLabels(m_lab);
144  else if (supports_regression())
145  lab=new CRegressionLabels(m_lab);
146 
147  for (index_t cidx = 0; cidx < m_log_lam.num_cols; cidx++)
148  {
149  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
151  Map<VectorXd> eigen_lp(lp.vector, lp.vlen);
152  eigen_res+=eigen_lp*m_wgh[cidx];
153  }
154 
155  SG_UNREF(lab);
156 
157  return res;
158 }
159 
161  const TParameter* param) const
162 {
163  //based on the likKL(v, lik, varargin) function in infKL.m
164 
165  //compute gradient using numerical integration
166  REQUIRE(param, "Param is required (param should not be NULL)\n");
167  REQUIRE(param->m_name, "Param name is required (param->m_name should not be NULL)\n");
168  //We take the derivative wrt to param. Only mu or sigma2 can be the param
169  REQUIRE(!(strcmp(param->m_name, "mu") && strcmp(param->m_name, "sigma2")),
170  "Can't compute derivative of the variational expection ",
171  "of log LogitLikelihood using numerical integration ",
172  "wrt %s.%s parameter. The function only accepts mu and sigma2 as parameter\n",
173  get_name(), param->m_name);
174 
176  res.zero();
177  Map<VectorXd> eigen_res(res.vector, res.vlen);
178 
179  Map<VectorXd> eigen_v(m_s2.vector, m_s2.vlen);
180 
181  CLabels* lab=NULL;
182 
183  if (supports_binary())
184  lab=new CBinaryLabels(m_lab);
185  else if (supports_regression())
186  lab=new CRegressionLabels(m_lab);
187 
188  if (strcmp(param->m_name, "mu")==0)
189  {
190  //Compute the derivative wrt mu
191 
192  //df = df + w(i)*(dlp);
193  for (index_t cidx=0; cidx<m_log_lam.num_cols; cidx++)
194  {
195  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
197  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
198  eigen_res+=eigen_dlp*m_wgh[cidx];
199  }
200  }
201  else
202  {
203  //Compute the derivative wrt sigma2
204 
205  //ai = t(i)./(2*sv+eps); dvi = dlp.*ai; dv = dv + w(i)*dvi;
206  VectorXd eigen_sv=eigen_v.array().sqrt().matrix();
207  const float64_t EPS=2.2204e-16;
208 
209  for (index_t cidx=0; cidx<m_log_lam.num_cols; cidx++)
210  {
211  SGVector<float64_t> tmp(m_log_lam.get_column_vector(cidx), m_log_lam.num_rows, false);
213  Map<VectorXd> eigen_dlp(dlp.vector, dlp.vlen);
214  eigen_res+=((m_wgh[cidx]*0.5*m_xgh[cidx])*eigen_dlp.array()/(eigen_sv.array()+EPS)).matrix();
215  }
216  }
217 
218  SG_UNREF(lab);
219 
220  return res;
221 }
222 
223 
225  SGVector<float64_t> s2, const CLabels* lab)
226 {
227  bool status = true;
229 
230  if (status)
231  {
232  if (supports_binary())
233  {
235  "Labels must be type of CBinaryLabels\n");
236  }
237  else
238  {
239  if (supports_regression())
240  {
242  "Labels must be type of CRegressionLabels\n");
243  }
244  else
245  SG_ERROR("Unsupported Label type\n");
246  }
247 
248  if (supports_binary())
249  m_lab=((CBinaryLabels*)lab)->get_labels();
250  else
251  m_lab=((CRegressionLabels*)lab)->get_labels();
252 
253  if (!m_is_init_GHQ)
254  {
255  m_xgh=SGVector<float64_t>(m_GHQ_N);
256  m_wgh=SGVector<float64_t>(m_GHQ_N);
257  CIntegration::generate_gauher(m_xgh, m_wgh);
258  m_is_init_GHQ=true;
259  }
260 
261  precompute();
262 
263  }
264 
265  return status;
266 }
267 
268 void CNumericalVGLikelihood::precompute()
269 {
270  //samples-by-abscissas
271  m_log_lam=SGMatrix<float64_t>(m_s2.vlen, m_xgh.vlen);
272 
273  Map<MatrixXd> eigen_log_lam(m_log_lam.matrix, m_log_lam.num_rows, m_log_lam.num_cols);
274  Map<VectorXd> eigen_v(m_s2.vector, m_s2.vlen);
275  Map<VectorXd> eigen_f(m_mu.vector, m_mu.vlen);
276  Map<VectorXd> eigen_t(m_xgh.vector, m_xgh.vlen);
277 
278  VectorXd eigen_sv=eigen_v.array().sqrt().matrix();
279  //varargin{3} = f + sv*t(i); % coordinate transform of the quadrature points
280  eigen_log_lam = (
281  (eigen_t.replicate(1, eigen_v.rows()).array().transpose().colwise()
282  *eigen_sv.array()).colwise()+eigen_f.array()
283  ).matrix();
284 }
285 
286 } /* namespace shogun */
287 
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:62
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:20
parameter struct
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
index_t num_cols
Definition: SGMatrix.h:376
static void generate_gauher(SGVector< float64_t > xgh, SGVector< float64_t > wgh)
virtual SGVector< float64_t > get_first_derivative(const CLabels *lab, SGVector< float64_t > func, const TParameter *param) const
index_t num_rows
Definition: SGMatrix.h:374
index_t vlen
Definition: SGVector.h:492
double float64_t
Definition: common.h:50
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()
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
T * get_column_vector(index_t col) const
Definition: SGMatrix.h:113
#define SG_UNREF(x)
Definition: SGObject.h:55
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:84
virtual SGVector< float64_t > get_first_derivative_wrt_hyperparameter(const TParameter *param) const
virtual SGVector< float64_t > get_variational_first_derivative(const TParameter *param) const

SHOGUN Machine Learning Toolbox - Documentation