gaussian_process.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <stdexcept>
#include <vector>
 
#include "pulsatrix/linear_algebra.hpp"
 
namespace pulsatrix {
 
class GaussianProcessRegressor {
public:
    GaussianProcessRegressor(float sigma_f, float length_scale, float noise_variance)
        : sigma_f_(sigma_f), length_scale_(length_scale), noise_variance_(noise_variance) {
        if (sigma_f <= 0.0f) {
            throw std::invalid_argument("GaussianProcessRegressor: sigma_f must be positive");
        }
        if (length_scale <= 0.0f) {
            throw std::invalid_argument("GaussianProcessRegressor: length_scale must be positive");
        }
        if (noise_variance < 0.0f) {
            throw std::invalid_argument("GaussianProcessRegressor: noise_variance must be non-negative");
        }
    }
 
    void Fit(std::vector<std::vector<float>> inputs, std::vector<float> targets) {
        if (inputs.empty()) {
            throw std::invalid_argument("GaussianProcessRegressor::Fit: inputs must not be empty");
        }
        if (inputs.size() != targets.size()) {
            throw std::invalid_argument("GaussianProcessRegressor::Fit: inputs and targets must be the same size");
        }
        const size_t dim = inputs[0].size();
        for (const auto& row : inputs) {
            if (row.size() != dim) {
                throw std::invalid_argument("GaussianProcessRegressor::Fit: every input row must be the same dimension");
            }
        }
 
        const size_t n = inputs.size();
        std::vector<std::vector<float>> k(n, std::vector<float>(n, 0.0f));
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < n; ++j) {
                k[i][j] = Kernel(inputs[i], inputs[j]);
            }
            k[i][i] += noise_variance_;
        }
 
        alpha_ = SolveLinearSystem(k, targets);
        training_kernel_ = std::move(k);
        inputs_ = std::move(inputs);
        targets_ = std::move(targets);
    }
 
    struct Posterior {
        double mean;
        double variance;
    };
 
    [[nodiscard]] Posterior Predict(const std::vector<float>& x) const {
        if (inputs_.empty()) {
            throw std::runtime_error("GaussianProcessRegressor::Predict: called before Fit");
        }
        if (x.size() != inputs_[0].size()) {
            throw std::invalid_argument("GaussianProcessRegressor::Predict: x dimension must match training inputs");
        }
 
        const size_t n = inputs_.size();
        std::vector<float> k_star(n);
        for (size_t i = 0; i < n; ++i) {
            k_star[i] = Kernel(x, inputs_[i]);
        }
 
        double mean = 0.0;
        for (size_t i = 0; i < n; ++i) {
            mean += static_cast<double>(k_star[i]) * static_cast<double>(alpha_[i]);
        }
 
        auto v = SolveLinearSystem(training_kernel_, k_star);
        double explained_variance = 0.0;
        for (size_t i = 0; i < n; ++i) {
            explained_variance += static_cast<double>(k_star[i]) * static_cast<double>(v[i]);
        }
        double variance = static_cast<double>(Kernel(x, x)) - explained_variance;
        if (variance < 0.0) {
            variance = 0.0;
        }
 
        return Posterior{mean, variance};
    }
 
private:
    [[nodiscard]] float Kernel(const std::vector<float>& a, const std::vector<float>& b) const {
        float squared_distance = 0.0f;
        for (size_t i = 0; i < a.size(); ++i) {
            float d = a[i] - b[i];
            squared_distance += d * d;
        }
        return sigma_f_ * sigma_f_ *
               std::exp(-squared_distance / (2.0f * length_scale_ * length_scale_));
    }
 
    float sigma_f_;
    float length_scale_;
    float noise_variance_;
    std::vector<std::vector<float>> inputs_;
    std::vector<float> targets_;
    std::vector<std::vector<float>> training_kernel_;
    std::vector<float> alpha_;
};
 
}  // namespace pulsatrix