gp.hpp 12.9 KB
Newer Older
1
2
#ifndef LIMBO_MODEL_GP_HPP
#define LIMBO_MODEL_GP_HPP
3

4
#include <iostream>
5
#include <cassert>
6
#include <limits>
7
#include <vector>
8

9
10
11
12
#include <Eigen/Core>
#include <Eigen/LU>
#include <Eigen/Cholesky>

13
#include <limbo/model/gp/no_lf_opt.hpp>
14

15
namespace limbo {
16
    namespace model {
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
17
18
19
20
21
        /// @ingroup model
        /// A classic Gaussian process.
        /// It is parametrized by:
        /// - a mean function
        /// - [optionnal] an optimizer for the hyper-parameters
22
        template <typename Params, typename KernelFunction, typename MeanFunction, class HyperParamsOptimizer = gp::NoLFOpt<Params>>
23
        class GP {
24
        public:
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
25
            /// useful because the model might be created before knowing anything about the process
26
            GP() : _dim_in(-1), _dim_out(-1) {}
27

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
28
            /// useful because the model might be created  before having samples
29
30
31
            GP(int dim_in, int dim_out)
                : _dim_in(dim_in), _dim_out(dim_out), _kernel_function(dim_in), _mean_function(dim_out) {}

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
32
            /// Compute the GP from samples, observation, noise. [optionnal: blacklisted samples]. This call needs to be explicit!
33
            void compute(const std::vector<Eigen::VectorXd>& samples,
34
                const std::vector<Eigen::VectorXd>& observations, double noise,
35
36
                const std::vector<Eigen::VectorXd>& bl_samples = std::vector<Eigen::VectorXd>())
            {
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
                //should be checked each time! not only the first time
                assert(samples.size() != 0);
                assert(observations.size() != 0);
                assert(samples.size() == observations.size());

                if (_dim_in != samples[0].size()) {
                    _dim_in = samples[0].size();
                    _kernel_function = KernelFunction(_dim_in); // the cost of building a functor should be relatively low
                }

                if (_dim_out != observations[0].size()) {
                    _dim_out = observations[0].size();
                    _mean_function = MeanFunction(_dim_out); // the cost of building a functor should be relatively low
                }

52
                _samples = samples;
53
                _observations.resize(observations.size(), _dim_out);
54
55
56
57
58
59
60
61
62
63
64
65
66
                for (int i = 0; i < _observations.rows(); ++i)
                    _observations.row(i) = observations[i];

                _mean_observation.resize(_dim_out);
                for (int i = 0; i < _observations.cols(); i++)
                    _mean_observation(i) = _observations.col(i).sum() / _observations.rows();

                _noise = noise;

                _bl_samples = bl_samples;

                _compute_obs_mean();
                _compute_kernel();
67

68
                HyperParamsOptimizer()(*this);
69
70
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
71
72
            /**
             \rst
73
             return :math:`\mu`, :math:`\sigma^2` (unormalized). If there is no sample, return the value according to the mean function. Using this method instead of separate calls to mu() and sigma() is more efficient because some computations are shared between mu() and sigma().
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
74
75
             \endrst
            */
76
            std::tuple<Eigen::VectorXd, double> query(const Eigen::VectorXd& v) const
77
78
79
80
81
82
83
84
85
86
87
88
89
            {
                if (_samples.size() == 0 && _bl_samples.size() == 0)
                    return std::make_tuple(_mean_function(v, *this),
                        sqrt(_kernel_function(v, v)));

                if (_samples.size() == 0)
                    return std::make_tuple(_mean_function(v, *this),
                        _sigma(v, _compute_k_bl(v, _compute_k(v))));

                Eigen::VectorXd k = _compute_k(v);
                return std::make_tuple(_mu(v, k), _sigma(v, _compute_k_bl(v, k)));
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
90
91
            /**
             \rst
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
92
             return :math:`\mu` (unormalized). If there is no sample, return the value according to the mean function.
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
93
94
             \endrst
            */
95
            Eigen::VectorXd mu(const Eigen::VectorXd& v) const
96
97
98
99
100
101
            {
                if (_samples.size() == 0)
                    return _mean_function(v, *this);
                return _mu(v, _compute_k(v));
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
102
103
            /**
             \rst
104
             return :math:`\sigma^2` (unormalized). If there is no sample, return the value according to the mean function.
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
105
106
             \endrst
            */
107
108
109
110
111
112
113
            double sigma(const Eigen::VectorXd& v) const
            {
                if (_samples.size() == 0 && _bl_samples.size() == 0)
                    return sqrt(_kernel_function(v, v));
                return _sigma(v, _compute_k_bl(v, _compute_k(v)));
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
114
            /// return the number of dimensions of the input
115
116
117
118
119
120
            int dim_in() const
            {
                assert(_dim_in != -1); // need to compute first !
                return _dim_in;
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
121
            /// return the number of dimensions of the output
122
123
124
125
126
127
128
129
            int dim_out() const
            {
                assert(_dim_out != -1); // need to compute first !
                return _dim_out;
            }

            const KernelFunction& kernel_function() const { return _kernel_function; }

130
131
            KernelFunction& kernel_function() { return _kernel_function; }

132
133
            const MeanFunction& mean_function() const { return _mean_function; }

134
135
            MeanFunction& mean_function() { return _mean_function; }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
136
            /// return the maximum observation (only call this if the output of the GP is of dimension 1)
137
            Eigen::VectorXd max_observation() const
138
139
            {
                if (_observations.cols() > 1)
140
                    std::cout << "WARNING max_observation with multi dimensional "
141
142
143
144
                                 "observations doesn't make sense" << std::endl;
                return _observations.maxCoeff();
            }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
145
            /// return the mean observation (only call this if the output of the GP is of dimension 1)
146
            Eigen::VectorXd mean_observation() const
147
            {
148
                // TO-DO: Check if _dim_out is correct?!
149
                return _samples.size() > 0 ? _mean_observation
150
                                           : Eigen::VectorXd::Zero(_dim_out);
151
152
153
154
            }

            const Eigen::MatrixXd& mean_vector() const { return _mean_vector; }

155
            const Eigen::MatrixXd& obs_mean() const { return _obs_mean; }
156

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
157
            /// return the number of samples used to compute the GP
158
159
            int nb_samples() const { return _samples.size(); }

160
161
            /** return the number of blacklisted samples used to compute the GP
            \rst
162
            For the blacklist concept, see the Limbo-specific concept guide.
163
164
            \endrst
            */
165
166
            int nb_bl_samples() const { return _bl_samples.size(); }

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
167
            /// update the GP
168
169
170
171
172
173
            void update()
            {
                this->_compute_obs_mean(); // ORDER MATTERS
                this->_compute_kernel();
            }

174
            /// return the likelihood (do not compute it!)
175
            double get_lik() const { return _lik; }
176

177
            /// set the likelihood (you need to compute it from outside!)
178
            void set_lik(const double& lik) { _lik = lik; }
179

180
            /// LLT matrix (from Cholesky decomposition)
181
            const Eigen::LLT<Eigen::MatrixXd>& llt() const { return _llt; }
182

183
            const Eigen::MatrixXd& alpha() const { return _alpha; }
184

185
            /// return the list of samples that have been tested so far
186
            const std::vector<Eigen::VectorXd>& samples() const { return _samples; }
187

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
        protected:
            int _dim_in;
            int _dim_out;

            KernelFunction _kernel_function;
            MeanFunction _mean_function;

            std::vector<Eigen::VectorXd> _samples;
            Eigen::MatrixXd _observations;
            std::vector<Eigen::VectorXd> _bl_samples; // black listed samples
            Eigen::MatrixXd _mean_vector;
            Eigen::MatrixXd _obs_mean;

            double _noise;
            Eigen::MatrixXd _alpha;
203
            Eigen::VectorXd _mean_observation;
204
205
206
207
208
209
210

            Eigen::MatrixXd _kernel;
            // Eigen::MatrixXd _inverted_kernel;
            Eigen::MatrixXd _l_matrix;
            Eigen::LLT<Eigen::MatrixXd> _llt;
            Eigen::MatrixXd _inv_bl_kernel;

211
            double _lik;
212

213
214
215
216
217
218
219
220
221
222
223
224
            void _compute_obs_mean()
            {
                _mean_vector.resize(_samples.size(), _dim_out);
                for (int i = 0; i < _mean_vector.rows(); i++)
                    _mean_vector.row(i) = _mean_function(_samples[i], *this);
                _obs_mean = _observations - _mean_vector;
            }

            void _compute_kernel()
            {
                // O(n^2) [should be negligible]
                _kernel.resize(_samples.size(), _samples.size());
225
226
                for (size_t i = 0; i < _samples.size(); i++)
                    for (size_t j = 0; j < _samples.size(); ++j)
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
                        _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();

                _llt = Eigen::LLT<Eigen::MatrixXd>(_kernel);

                // alpha = K^{-1} * this->_obs_mean;
                _alpha = _llt.matrixL().solve(_obs_mean);
                _llt.matrixL().adjoint().solveInPlace(_alpha);
                if (_bl_samples.size() == 0)
                    return;

                Eigen::MatrixXd A1 = Eigen::MatrixXd::Identity(this->_samples.size(), this->_samples.size());
                _llt.matrixL().solveInPlace(A1);
                _llt.matrixL().transpose().solveInPlace(A1);
                _inv_bl_kernel.resize(_samples.size() + _bl_samples.size(),
                    _samples.size() + _bl_samples.size());

                Eigen::MatrixXd B(_samples.size(), _bl_samples.size());
                for (size_t i = 0; i < _samples.size(); i++)
                    for (size_t j = 0; j < _bl_samples.size(); ++j)
                        B(i, j) = _kernel_function(_samples[i], _bl_samples[j]);

                Eigen::MatrixXd D(_bl_samples.size(), _bl_samples.size());
                for (size_t i = 0; i < _bl_samples.size(); i++)
                    for (size_t j = 0; j < _bl_samples.size(); ++j)
                        D(i, j) = _kernel_function(_bl_samples[i], _bl_samples[j]) + ((i == j) ? _noise : 0);

                Eigen::MatrixXd comA = (D - B.transpose() * A1 * B);
                Eigen::LLT<Eigen::MatrixXd> llt_bl(comA);
                Eigen::MatrixXd comA1 = Eigen::MatrixXd::Identity(_bl_samples.size(), _bl_samples.size());
                llt_bl.matrixL().solveInPlace(comA1);
                llt_bl.matrixL().transpose().solveInPlace(comA1);

                // fill the matrix block wise
                _inv_bl_kernel.block(0, 0, _samples.size(), _samples.size()) = A1 + A1 * B * comA1 * B.transpose() * A1;
                _inv_bl_kernel.block(0, _samples.size(), _samples.size(),
                    _bl_samples.size()) = -A1 * B * comA1;
                _inv_bl_kernel.block(_samples.size(), 0, _bl_samples.size(),
                    _samples.size()) = _inv_bl_kernel.block(0, _samples.size(), _samples.size(),
                                                          _bl_samples.size()).transpose();
                _inv_bl_kernel.block(_samples.size(), _samples.size(), _bl_samples.size(),
                    _bl_samples.size()) = comA1;
            }

273
            Eigen::VectorXd _mu(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const
274
275
276
277
278
279
280
            {
                return (k.transpose() * _alpha) + _mean_function(v, *this).transpose();
            }

            double _sigma(const Eigen::VectorXd& v, const Eigen::VectorXd& k) const
            {
                double res;
281
                if (_bl_samples.size() == 0) {
282
283
284
                    Eigen::VectorXd z = _llt.matrixL().solve(k);
                    res = _kernel_function(v, v) - z.dot(z);
                }
285
                else {
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
                    res = _kernel_function(v, v) - k.transpose() * _inv_bl_kernel * k;
                }

                return (res <= std::numeric_limits<double>::epsilon()) ? 0 : res;
            }

            Eigen::VectorXd _compute_k(const Eigen::VectorXd& v) const
            {
                Eigen::VectorXd k(_samples.size());
                for (int i = 0; i < k.size(); i++)
                    k[i] = _kernel_function(_samples[i], v);
                return k;
            }

            Eigen::VectorXd _compute_k_bl(const Eigen::VectorXd& v,
                const Eigen::VectorXd& k) const
            {
303
                if (_bl_samples.size() == 0) {
304
305
306
307
308
309
310
311
312
313
314
315
                    return k;
                }

                Eigen::VectorXd k_bl(_samples.size() + _bl_samples.size());

                k_bl.head(_samples.size()) = k;
                for (size_t i = 0; i < _bl_samples.size(); i++)
                    k_bl[i + this->_samples.size()] = this->_kernel_function(_bl_samples[i], v);
                return k_bl;
            }
        };
    }
316
}
317

318
#endif