SHOGUN  6.1.3
KLLowerTriangularInference.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  * http://hannes.nickisch.org/code/approxXX.tar.gz
32  * and Gaussian Process Machine Learning Toolbox
33  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
34  * and the reference paper is
35  * Challis, Edward, and David Barber.
36  * "Concave Gaussian variational approximations for inference in large-scale Bayesian linear models."
37  * International conference on Artificial Intelligence and Statistics. 2011.
38  *
39  * This code specifically adapted from function in approxKL.m and infKL.m
40  */
41 
43 
47 
48 using namespace Eigen;
49 
50 namespace shogun
51 {
52 
53 CKLLowerTriangularInference::CKLLowerTriangularInference() : CKLInference()
54 {
55  init();
56 }
57 
59  CFeatures* feat, CMeanFunction* m, CLabels* lab, CLikelihoodModel* mod)
60  : CKLInference(kern, feat, m, lab, mod)
61 {
62  init();
63 }
64 
65 void CKLLowerTriangularInference::init()
66 {
67  SG_ADD(&m_InvK_Sigma, "invk_Sigma",
68  "K^{-1}Sigma'",
70  SG_ADD(&m_mean_vec, "mean_vec",
71  "The mean vector generated from mean function",
73  SG_ADD(&m_log_det_Kernel, "log_det_kernel",
74  "The Log-determinant of Kernel",
76 
77  SG_ADD(&m_Kernel_LsD, "L_sqrt_D",
78  "The L*sqrt(D) matrix, where L and D are defined in LDLT factorization on Kernel*sq(m_scale)",
80  SG_ADD(&m_Kernel_P, "Permutation_P",
81  "The permutation sequence of P, where P are defined in LDLT factorization on Kernel*sq(m_scale)",
84 }
85 
87 {
88 }
89 
91 {
95  return SGVector<float64_t>();
96 }
97 
99 {
103 }
104 
106 {
108  MatrixXd Kernel_D=ldlt.vectorD();
109  MatrixXd Kernel_L=ldlt.matrixL();
111  m_Kernel_LsD.zero();
113  eigen_Kernel_LsD.triangularView<Lower>()=Kernel_L*Kernel_D.array().sqrt().matrix().asDiagonal();
114  m_log_det_Kernel=2.0*eigen_Kernel_LsD.diagonal().array().abs().log().sum();
115 
117  for (index_t i=0; i<m_Kernel_P.vlen; i++)
118  m_Kernel_P[i]=i;
119  Map<VectorXi> eigen_Kernel_P(m_Kernel_P.vector, m_Kernel_P.vlen);
120  eigen_Kernel_P=ldlt.transpositionsP()*eigen_Kernel_P;
121 
123 }
124 
126 {
127  Map<VectorXi> eigen_Kernel_P(m_Kernel_P.vector, m_Kernel_P.vlen);
129 
130  //re-construct the Permutation Matrix
131  PermutationMatrix<Dynamic> P(m_Kernel_P.vlen);
132  P.setIdentity();
134  for (index_t i=0; i<tmp.vlen; i++)
135  {
136  while(tmp[i]>i)
137  {
138  P.applyTranspositionOnTheLeft(i,tmp[i]);
139  index_t idx=tmp[i];
140  tmp[i]=tmp[idx];
141  tmp[idx]=idx;
142  }
143  }
144  P=P.transpose();
145  //(P'LDL'P)\eigen_A
146  MatrixXd tmp1=P*eigen_A;
147  MatrixXd tmp2=eigen_Kernel_LsD.triangularView<Lower>().solve(tmp1);
148  MatrixXd tmp3=eigen_Kernel_LsD.triangularView<Lower>().transpose().solve(tmp2);
149  return P.transpose()*tmp3;
150 }
151 
153 {
154  Map<MatrixXd> eigen_dK(dK.matrix, dK.num_rows, dK.num_cols);
155  Map<VectorXd> eigen_alpha(m_alpha.vector, m_mu.vlen);
156  Map<VectorXd> eigen_mu(m_mu.vector, m_mu.vlen);
158 
159  //dnlZ(j)=0.5*sum(sum(dK.*(K\((eye(n)- (invK_V+alpha*m'))')),2),1);
160  MatrixXd tmp1=eigen_InvK_Sigma+eigen_alpha*(eigen_mu.transpose());
161  MatrixXd tmp2=(MatrixXd::Identity(m_ktrtr.num_rows,m_ktrtr.num_cols)-tmp1).transpose();
162  MatrixXd tmp3=solve_inverse(tmp2);
163  return 0.5*(tmp3.array()*eigen_dK.array()).sum();
164 }
165 
167 {
171 }
172 
174 {
175  update_Sigma();
177 
181  MatrixXd tmp2=(eigen_InvK_Sigma-MatrixXd::Identity(m_ktrtr.num_rows,m_ktrtr.num_cols)).transpose();
182 
183  eigen_L=solve_inverse(tmp2);
184 }
185 
186 } /* namespace shogun */
187 
int32_t index_t
Definition: common.h:72
The class Labels models labels, i.e. class assignments of objects.
Definition: Labels.h:43
Definition: SGMatrix.h:25
virtual SGVector< float64_t > get_mean_vector(const CFeatures *features) const =0
An abstract class of the mean function.
Definition: MeanFunction.h:49
virtual Eigen::LDLT< Eigen::MatrixXd, 0x1 > update_init_helper()
CFeatures * m_features
Definition: Inference.h:473
SGMatrix< float64_t > m_ktrtr
Definition: Inference.h:488
virtual float64_t get_derivative_related_cov(SGMatrix< float64_t > dK)
CMeanFunction * m_mean
Definition: Inference.h:467
The KL approximation inference method class.
Definition: KLInference.h:75
double float64_t
Definition: common.h:60
index_t num_rows
Definition: SGMatrix.h:495
SGVector< float64_t > m_mu
Definition: KLInference.h:367
SGMatrix< float64_t > m_L
Definition: Inference.h:482
SGVector< T > clone() const
Definition: SGVector.cpp:262
index_t num_cols
Definition: SGMatrix.h:497
Eigen::MatrixXd solve_inverse(Eigen::MatrixXd A)
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
T sum(const Container< T > &a, bool no_diag=false)
The class Features is the base class of all feature objects.
Definition: Features.h:69
The Kernel base class.
virtual SGVector< float64_t > get_diagonal_vector()
#define SG_ADD(...)
Definition: SGObject.h:93
The Likelihood model base class.
index_t vlen
Definition: SGVector.h:571
SGVector< float64_t > m_alpha
Definition: Inference.h:479

SHOGUN Machine Learning Toolbox - Documentation