SHOGUN  6.1.3
MatrixOperations.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  * Gaussian Process Machine Learning Toolbox
32  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
33  */
34 
36 
37 #include <shogun/lib/SGVector.h>
38 #include <shogun/lib/SGMatrix.h>
40 
41 using namespace Eigen;
42 
43 namespace shogun
44 {
45 
46 SGMatrix<float64_t> CMatrixOperations::get_choleksy(SGVector<float64_t> W,
48 {
49  Map<VectorXd> eigen_W(W.vector, W.vlen);
50  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
51  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
52 
53  REQUIRE(eigen_W.rows()==eigen_sW.rows(),
54  "The length of W (%d) and sW (%d) should be the same\n",
55  eigen_W.rows(), eigen_sW.rows());
56  REQUIRE(eigen_kernel.rows()==eigen_kernel.cols(),
57  "Kernel should a sqaure matrix, row (%d) col (%d)\n",
58  eigen_kernel.rows(), eigen_kernel.cols());
59  REQUIRE(eigen_kernel.rows()==eigen_W.rows(),
60  "The dimension between kernel (%d) and W (%d) should be matched\n",
61  eigen_kernel.rows(), eigen_W.rows());
62 
63  SGMatrix<float64_t> L(eigen_W.rows(), eigen_W.rows());
64 
65  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
66 
67  if (eigen_W.minCoeff()<0)
68  {
69  // compute inverse of diagonal noise: iW = 1/W
70  VectorXd eigen_iW=(VectorXd::Ones(eigen_W.rows())).cwiseQuotient(eigen_W);
71 
72  FullPivLU<MatrixXd> lu(
73  eigen_kernel*CMath::sq(scale)+MatrixXd(eigen_iW.asDiagonal()));
74 
75  // compute cholesky: L = -(K + iW)^-1
76  eigen_L=-lu.inverse();
77  }
78  else
79  {
80  // compute cholesky: L = chol(sW * sW' .* K + I)
81  LLT<MatrixXd> llt(
82  (eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_kernel*CMath::sq(scale))+
83  MatrixXd::Identity(eigen_kernel.rows(), eigen_kernel.cols()));
84 
85  eigen_L=llt.matrixU();
86  }
87 
88  return L;
89 }
90 
91 SGMatrix<float64_t> CMatrixOperations::get_inverse(SGMatrix<float64_t> L, SGMatrix<float64_t> kernel,
93 {
94  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
95  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
96  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
98  Map<MatrixXd> eigen_V(V.matrix, V.num_rows, V.num_cols);
99 
100  // compute V = L^(-1) * W^(1/2) * K, using upper triangular factor L^T
101  eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
102  eigen_sW.asDiagonal()*eigen_kernel*CMath::sq(scale));
103 
104  return get_inverse(L, kernel, sW, V, scale);
105 }
106 
107 SGMatrix<float64_t> CMatrixOperations::get_inverse(SGMatrix<float64_t> L,
110 {
111  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
112  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
113  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
114  Map<MatrixXd> eigen_V(V.matrix, V.num_rows, V.num_cols);
115 
116  REQUIRE(eigen_kernel.rows()==eigen_kernel.cols(),
117  "Kernel should a sqaure matrix, row (%d) col (%d)\n",
118  eigen_kernel.rows(), eigen_kernel.cols());
119  REQUIRE(eigen_L.rows()==eigen_L.cols(),
120  "L should a sqaure matrix, row (%d) col (%d)\n",
121  eigen_L.rows(), eigen_L.cols());
122  REQUIRE(eigen_kernel.rows()==eigen_sW.rows(),
123  "The dimension between kernel (%d) and sW (%d) should be matched\n",
124  eigen_kernel.rows(), eigen_sW.rows());
125  REQUIRE(eigen_L.rows()==eigen_sW.rows(),
126  "The dimension between L (%d) and sW (%d) should be matched\n",
127  eigen_L.rows(), eigen_sW.rows());
128 
129 
130  SGMatrix<float64_t> tmp(eigen_L.rows(), eigen_L.cols());
131  Map<MatrixXd> eigen_tmp(tmp.matrix, tmp.num_rows, tmp.num_cols);
132 
133  // compute covariance matrix of the posterior:
134  // Sigma = K - K * W^(1/2) * (L * L^T)^(-1) * W^(1/2) * K =
135  // K - (K * W^(1/2)) * (L^T)^(-1) * L^(-1) * W^(1/2) * K =
136  // K - (W^(1/2) * K)^T * (L^(-1))^T * L^(-1) * W^(1/2) * K = K - V^T * V
137  eigen_tmp=eigen_kernel*CMath::sq(scale)-eigen_V.adjoint()*eigen_V;
138  return tmp;
139 }
140 
141 } /* namespace shogun */
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
Definition: SGMatrix.h:25
#define REQUIRE(x,...)
Definition: SGIO.h:181
double float64_t
Definition: common.h:60
index_t num_rows
Definition: SGMatrix.h:495
Matrix< float64_t,-1,-1, 0,-1,-1 > MatrixXd
Definition: KLInference.h:52
index_t num_cols
Definition: SGMatrix.h:497
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
index_t vlen
Definition: SGVector.h:571

SHOGUN Machine Learning Toolbox - Documentation