Commit 0545f944 authored by Antoine Cully's avatar Antoine Cully
Browse files

minor fix for the nothing case

parent 8bf3cc85
......@@ -59,7 +59,7 @@ namespace limbo {
_noise = noise;
if (samples.size() != _samples.size() || !change_samples) //no changes in samples since last computation.
if (samples.size() != _samples.size() || change_samples) //no changes in samples since last computation.
{
_samples = samples;
......@@ -242,31 +242,41 @@ namespace limbo {
void _compute_kernel(bool incremental = false)
{
// Check if new kernel is the same as before
if (!incremental) //_kernel.cols()==0 || !kernel.block(0,0,_kernel.rows(),_kernel.cols()).isApprox(_kernel, 1e-10)) // incremental LLT is impossible
{
_kernel.resize(_samples.size(), _samples.size());
_kernel.resize(_samples.size(), _samples.size());
// O(n^2) [should be negligible]
for (size_t i = 0; i < _samples.size(); i++)
for (size_t j = 0; j < _samples.size(); ++j)
_kernel(i, j) = _kernel_function(_samples[i], _samples[j]) + ((i == j) ? _noise : 0); // noise only on the diagonal
// O(n^3)
// _inverted_kernel = _kernel.inverse();
_matrixL = Eigen::LLT<Eigen::MatrixXd>(_kernel).matrixL();
}
else //incremental LLT
//This part of the code is inpired from the Bayesopt Library (cholesky_add_row function). However, the mathematical fundations can be easily retrieved by detailling the equations of the extended L matrix that produces the desired kernel.
{
size_t n = _samples.size();
_kernel.conservativeResize(n, n);
size_t n = _samples.size();
_kernel.conservativeResize(n, n);
for (size_t i = 0; i < n; ++i) {
_kernel(i, n - 1) = _kernel_function(_samples[i], _samples[n - 1]) + ((i == n - 1) ? _noise : 0); // noise only on the diagonal
_kernel(n - 1, i) = _kernel(i, n - 1);
}
_matrixL.conservativeResizeLike(Eigen::MatrixXd::Zero(n, n));
_matrixL.conservativeResizeLike(Eigen::MatrixXd::Zero(n, n));
double L_j;
for (size_t j = 0; j < n - 1; ++j) {
L_j = _kernel(n - 1, j) - (_matrixL.block(j, 0, 1, j) * _matrixL.block(n - 1, 0, 1, j).transpose())(0, 0);
......@@ -275,12 +285,16 @@ namespace limbo {
L_j = _kernel(n - 1, n - 1) - (_matrixL.block(n - 1, 0, 1, n - 1) * _matrixL.block(n - 1, 0, 1, n - 1).transpose())(0, 0);
_matrixL(n - 1, n - 1) = sqrt(L_j);
}
// alpha = K^{-1} * this->_obs_mean;
_alpha = _matrixL.template triangularView<Eigen::Lower>().solve(_obs_mean);
_matrixL.template triangularView<Eigen::Lower>().adjoint().solveInPlace(_alpha); //can probably be improved by avoiding to generate the view twice
_change_kf = false;
}
void _compute_bl_kernel()
......
......@@ -260,3 +260,4 @@ BOOST_AUTO_TEST_CASE(test_gp_auto)
BOOST_CHECK(std::abs((mu(0) - 5)) < 1);
BOOST_CHECK(sigma < 1e-5);
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment