incremental llt, it compiles and divides by 3 the comput time. BL test still fails. WIP

src/limbo/model/gp/kernel_lf_opt.hpp

@@ -45,7 +45,7 @@ namespace limbo {
                         // --- cholesky ---
                         // see:
                         // http://xcorr.net/2008/06/11/log-determinant-of-positive-definite-matrices-in-matlab/
-                        Eigen::MatrixXd l = gp.llt().matrixL();
+                        Eigen::MatrixXd l = gp.matrixL();
                         long double det = 2 * l.diagonal().array().log().sum();
 
                         double a = (gp.obs_mean().transpose() * gp.alpha())
@@ -58,8 +58,10 @@ namespace limbo {
 
                         // K^{-1} using Cholesky decomposition
                         Eigen::MatrixXd w = Eigen::MatrixXd::Identity(n, n);
-                        gp.llt().matrixL().solveInPlace(w);
-                        gp.llt().matrixL().transpose().solveInPlace(w);
+
+			gp.matrixL().template triangularView<Eigen::Lower>().solveInPlace(w);
+                        gp.matrixL().template triangularView<Eigen::Lower>().transpose().solveInPlace(w);
+
 
                         // alpha * alpha.transpose() - K^{-1}
                         w = gp.alpha() * gp.alpha().transpose() - w;

src/limbo/model/gp/kernel_mean_lf_opt.hpp

@@ -47,7 +47,7 @@ namespace limbo {
                         // --- cholesky ---
                         // see:
                         // http://xcorr.net/2008/06/11/log-determinant-of-positive-definite-matrices-in-matlab/
-                        Eigen::MatrixXd l = gp.llt().matrixL();
+                        Eigen::MatrixXd l = gp.matrixL();
                         long double det = 2 * l.diagonal().array().log().sum();
 
                         double a = (gp.obs_mean().transpose() * gp.alpha())
@@ -60,8 +60,10 @@ namespace limbo {
 
                         // K^{-1} using Cholesky decomposition
                         Eigen::MatrixXd K = Eigen::MatrixXd::Identity(n, n);
-                        gp.llt().matrixL().solveInPlace(K);
-                        gp.llt().matrixL().transpose().solveInPlace(K);
+
+			
+			gp.matrixL().template triangularView<Eigen::Lower>().solveInPlace(K);
+                        gp.matrixL().template triangularView<Eigen::Lower>().transpose().solveInPlace(K);
 
                         // alpha * alpha.transpose() - K^{-1}
                         Eigen::MatrixXd w = gp.alpha() * gp.alpha().transpose() - K;

src/limbo/model/gp/mean_lf_opt.hpp

@@ -45,7 +45,7 @@ namespace limbo {
                         // --- cholesky ---
                         // see:
                         // http://xcorr.net/2008/06/11/log-determinant-of-positive-definite-matrices-in-matlab/
-                        Eigen::MatrixXd l = gp.llt().matrixL();
+                        Eigen::MatrixXd l = gp.matrixL();
                         long double det = 2 * l.diagonal().array().log().sum();
 
                         double a = (gp.obs_mean().transpose() * gp.alpha())
@@ -58,8 +58,8 @@ namespace limbo {
 
                         // K^{-1} using Cholesky decomposition
                         Eigen::MatrixXd K = Eigen::MatrixXd::Identity(n, n);
-                        gp.llt().matrixL().solveInPlace(K);
-                        gp.llt().matrixL().transpose().solveInPlace(K);
+                        gp.matrixL().template triangularView<Eigen::Lower>().solveInPlace(K);
+                        gp.matrixL().template triangularView<Eigen::Lower>().transpose().solveInPlace(K);
 
                         Eigen::VectorXd grad = Eigen::VectorXd::Zero(params.size());
                         for (int i_obs = 0; i_obs < gp.dim_out(); ++i_obs)

src/tests/test_gp.cpp

@@ -115,6 +115,66 @@ BOOST_AUTO_TEST_CASE(test_gp)
     }
 }
 
+
+
+BOOST_AUTO_TEST_CASE(test_gp_bw_inversion){
+  using namespace limbo;
+
+  typedef kernel::MaternFiveHalfs<Params> KF_t;
+  typedef mean::Constant<Params> Mean_t;
+  typedef model::GP<Params, KF_t, Mean_t> GP_t;
+  
+  std::vector<Eigen::VectorXd> observations;
+  std::vector<Eigen::VectorXd> samples;
+  tools::rgen_double_t rgen(0.0, 10);                                                                                                                                                    
+  for(size_t i=0;i<100; i++)
+    {
+      observations.push_back(make_v1(rgen.rand()));
+      samples.push_back(make_v1(rgen.rand()));
+    }
+  
+  GP_t gp;
+  auto t1 = std::chrono::steady_clock::now();
+  gp.compute(samples, observations, 0.0);
+  auto time_init = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
+  std::cout.precision(17);
+  std::cout << "Time running first batch: " << time_init << "us" << std::endl<<std::endl;
+
+
+  observations.push_back(make_v1(rgen.rand()));
+  samples.push_back(make_v1(rgen.rand()));
+
+
+  t1 = std::chrono::steady_clock::now();
+  gp.compute(samples, observations, 0.0);
+  auto time_increment = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
+  std::cout << "Time running increment: " << time_increment << "us" << std::endl<<std::endl; 
+
+
+  
+  t1 = std::chrono::steady_clock::now();
+  gp.compute(samples, observations, 0.0);
+  auto time_nothing = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
+  std::cout << "Time running no change: " << time_nothing << "us" << std::endl<<std::endl;  
+
+  GP_t gp2;
+  t1 = std::chrono::steady_clock::now();
+  gp2.compute(samples, observations, 0.0);
+  auto time_full = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
+  std::cout << "Time running whole batch: " << time_full << "us" << std::endl<<std::endl;
+  
+
+
+  Eigen::VectorXd s=make_v1(rgen.rand());  
+  BOOST_CHECK((gp.mu(s)-gp2.mu(s)).norm()<1e-5);
+  BOOST_CHECK(gp.matrixL().isApprox(gp2.matrixL(),1e-5));
+  BOOST_CHECK(time_full>time_increment);
+  BOOST_CHECK(time_increment > time_nothing);
+
+}
+
+
+
 BOOST_AUTO_TEST_CASE(test_gp_no_samples_acqui_opt)
 {
     using namespace limbo;