SHOGUN  3.2.1
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 宏定义  
MatrixOperations.cpp
浏览该文件的文档.
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 #ifdef HAVE_EIGEN3
38 #include <shogun/lib/SGVector.h>
39 #include <shogun/lib/SGMatrix.h>
40 
41 using namespace Eigen;
42 
43 namespace shogun
44 {
45 SGMatrix<float64_t> CMatrixOperations::get_choleksy(VectorXd eigen_W,
46  VectorXd eigen_sW, MatrixXd eigen_kernel, float64_t scale)
47 {
48  REQUIRE(eigen_W.rows()==eigen_sW.rows(),
49  "The length of W (%d) and sW (%d) should be the same\n",
50  eigen_W.rows(), eigen_sW.rows());
51  REQUIRE(eigen_kernel.rows()==eigen_kernel.cols(),
52  "Kernel should a sqaure matrix, row (%d) col (%d)\n",
53  eigen_kernel.rows(), eigen_kernel.cols());
54  REQUIRE(eigen_kernel.rows()==eigen_W.rows(),
55  "The dimension between kernel (%d) and W (%d) should be matched\n",
56  eigen_kernel.rows(), eigen_W.rows());
57 
58  SGMatrix<float64_t> L(eigen_W.rows(), eigen_W.rows());
59 
60  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
61 
62  if (eigen_W.minCoeff()<0)
63  {
64  // compute inverse of diagonal noise: iW = 1/W
65  VectorXd eigen_iW=(VectorXd::Ones(eigen_W.rows())).cwiseQuotient(eigen_W);
66 
67  FullPivLU<MatrixXd> lu(
68  eigen_kernel*CMath::sq(scale)+MatrixXd(eigen_iW.asDiagonal()));
69 
70  // compute cholesky: L = -(K + iW)^-1
71  eigen_L=-lu.inverse();
72  }
73  else
74  {
75  // compute cholesky: L = chol(sW * sW' .* K + I)
76  LLT<MatrixXd> llt(
77  (eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_kernel*CMath::sq(scale))+
78  MatrixXd::Identity(eigen_kernel.rows(), eigen_kernel.cols()));
79 
80  eigen_L=llt.matrixU();
81  }
82 
83  return L;
84 }
85 
86 SGMatrix<float64_t> CMatrixOperations::get_choleksy(SGVector<float64_t> W,
88 {
89 
90  Map<VectorXd> eigen_W(W.vector, W.vlen);
91  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
92  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
93 
94  return get_choleksy(eigen_W, eigen_sW, eigen_kernel, scale);
95 }
96 
97 SGMatrix<float64_t> CMatrixOperations::get_inverse(SGMatrix<float64_t> L, SGMatrix<float64_t> kernel,
99 {
100  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
101  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
102  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
103 
104  return get_inverse(eigen_L, eigen_kernel, eigen_sW, scale);
105 }
106 
107 SGMatrix<float64_t> CMatrixOperations::get_inverse(Eigen::MatrixXd eigen_L, Eigen::MatrixXd eigen_kernel,
108  Eigen::VectorXd eigen_sW, float64_t scale)
109 {
110  REQUIRE(eigen_kernel.rows()==eigen_kernel.cols(),
111  "Kernel should a sqaure matrix, row (%d) col (%d)\n",
112  eigen_kernel.rows(), eigen_kernel.cols());
113  REQUIRE(eigen_L.rows()==eigen_L.cols(),
114  "L should a sqaure matrix, row (%d) col (%d)\n",
115  eigen_L.rows(), eigen_L.cols());
116  REQUIRE(eigen_kernel.rows()==eigen_sW.rows(),
117  "The dimension between kernel (%d) and sW (%d) should be matched\n",
118  eigen_kernel.rows(), eigen_sW.rows());
119  REQUIRE(eigen_L.rows()==eigen_sW.rows(),
120  "The dimension between L (%d) and sW (%d) should be matched\n",
121  eigen_L.rows(), eigen_sW.rows());
122 
123  // compute V = L^(-1) * W^(1/2) * K, using upper triangular factor L^T
124  MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
125  eigen_sW.asDiagonal()*eigen_kernel*CMath::sq(scale));
126 
127  return get_inverse(eigen_L, eigen_kernel, eigen_sW, eigen_V, scale);
128 }
129 
130 SGMatrix<float64_t> CMatrixOperations::get_inverse(MatrixXd eigen_L,
131  MatrixXd eigen_kernel, VectorXd eigen_sW, MatrixXd eigen_V,
132  float64_t scale)
133 {
134  REQUIRE(eigen_kernel.rows()==eigen_kernel.cols(),
135  "Kernel should a sqaure matrix, row (%d) col (%d)\n",
136  eigen_kernel.rows(), eigen_kernel.cols());
137  REQUIRE(eigen_L.rows()==eigen_L.cols(),
138  "L should a sqaure matrix, row (%d) col (%d)\n",
139  eigen_L.rows(), eigen_L.cols());
140  REQUIRE(eigen_kernel.rows()==eigen_sW.rows(),
141  "The dimension between kernel (%d) and sW (%d) should be matched\n",
142  eigen_kernel.rows(), eigen_sW.rows());
143  REQUIRE(eigen_L.rows()==eigen_sW.rows(),
144  "The dimension between L (%d) and sW (%d) should be matched\n",
145  eigen_L.rows(), eigen_sW.rows());
146 
147 
148  SGMatrix<float64_t> tmp(eigen_L.rows(), eigen_L.cols());
149  Map<MatrixXd> eigen_tmp(tmp.matrix, tmp.num_rows, tmp.num_cols);
150 
151  // compute covariance matrix of the posterior:
152  // Sigma = K - K * W^(1/2) * (L * L^T)^(-1) * W^(1/2) * K =
153  // K - (K * W^(1/2)) * (L^T)^(-1) * L^(-1) * W^(1/2) * K =
154  // K - (W^(1/2) * K)^T * (L^(-1))^T * L^(-1) * W^(1/2) * K = K - V^T * V
155  eigen_tmp=eigen_kernel*CMath::sq(scale)-eigen_V.adjoint()*eigen_V;
156  return tmp;
157 }
158 
159 SGMatrix<float64_t> CMatrixOperations::get_inverse(SGMatrix<float64_t> L,
161  float64_t scale)
162 {
163  Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
164  Map<MatrixXd> eigen_kernel(kernel.matrix, kernel.num_rows, kernel.num_cols);
165  Map<VectorXd> eigen_sW(sW.vector, sW.vlen);
166  Map<MatrixXd> eigen_V(V.matrix, V.num_rows, V.num_cols);
167 
168  return get_inverse(eigen_L, eigen_kernel, eigen_sW, eigen_V, scale);
169 }
170 
171 } /* namespace shogun */
172 #endif /* HAVE_EIGEN3 */
#define REQUIRE(x,...)
Definition: SGIO.h:207
double float64_t
Definition: common.h:50
index_t num_rows
Definition: SGMatrix.h:298
index_t num_cols
Definition: SGMatrix.h:300
index_t vlen
Definition: SGVector.h:707

SHOGUN 机器学习工具包 - 项目文档