41 using namespace Eigen;
45 SGMatrix<float64_t> CMatrixOperations::get_choleksy(VectorXd eigen_W,
46 VectorXd eigen_sW, MatrixXd eigen_kernel,
float64_t scale)
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());
58 SGMatrix<float64_t> L(eigen_W.rows(), eigen_W.rows());
60 Map<MatrixXd> eigen_L(L.matrix, L.num_rows, L.num_cols);
62 if (eigen_W.minCoeff()<0)
65 VectorXd eigen_iW=(VectorXd::Ones(eigen_W.rows())).cwiseQuotient(eigen_W);
67 FullPivLU<MatrixXd> lu(
68 eigen_kernel*CMath::sq(scale)+MatrixXd(eigen_iW.asDiagonal()));
71 eigen_L=-lu.inverse();
77 (eigen_sW*eigen_sW.transpose()).cwiseProduct(eigen_kernel*CMath::sq(scale))+
78 MatrixXd::Identity(eigen_kernel.rows(), eigen_kernel.cols()));
80 eigen_L=llt.matrixU();
94 return get_choleksy(eigen_W, eigen_sW, eigen_kernel, scale);
104 return get_inverse(eigen_L, eigen_kernel, eigen_sW, scale);
108 Eigen::VectorXd eigen_sW,
float64_t scale)
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());
124 MatrixXd eigen_V=eigen_L.triangularView<Upper>().adjoint().solve(
125 eigen_sW.asDiagonal()*eigen_kernel*CMath::sq(scale));
127 return get_inverse(eigen_L, eigen_kernel, eigen_sW, eigen_V, scale);
131 MatrixXd eigen_kernel, VectorXd eigen_sW, MatrixXd eigen_V,
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());
149 Map<MatrixXd> eigen_tmp(tmp.matrix, tmp.num_rows, tmp.num_cols);
155 eigen_tmp=eigen_kernel*CMath::sq(scale)-eigen_V.adjoint()*eigen_V;
168 return get_inverse(eigen_L, eigen_kernel, eigen_sW, eigen_V, scale);