Feature request: GP Functor specification

Description

This issue describes the skeleton of the GP functor specification that was hashed out in Aalto before StanCon Helsinki. Please comment as necessary.

The basic idea is that we can build a functor (modelled after the ODE solver) that can automatically build all of the internals of the GP density from a given covariance function f.

The basic idea is that we can autodiff through f and use it to evaluate both the GP covariance matrix C and it's derivatives without storing these to the autodiff stack. This is a huge storage saving. The derivative with respect to the hyperparameters theta are given in equation (5.9) of Rasmussen and Williams.

There are some easy computational savings here if there is more than one hyperparameter. To compute the derivatives you need to:

• Populate the covariance matrix C
• Compute the Cholesky Decomposition of C
• Populate the derivative matrix for C (computed with the autodiff on f)
• Compute some products, solves, and traces.
  The first two steps only need to be performed once.

In order to speed up non-centering, we would need a second functor that could compute from a covariance function f the matrix-vector product Lu for a vector u and where C = LL^T. We need to work out what the appropriate partials of this expression are. They should be available without too much work. Once again, this avoids an O(n^2) addition to the autodiff expression tree.

Example

The code will likely look something like this.

template <bool propto, typename F, typename T_y, typename T_X, typename T_theta>
struct gp_lpdf {

  typedef typename return_type<T_y, T_X, T_theta> T_return_type;
  typedef typename partials_return_type<T_y, T_X, T_theta> T_partials_return_type;

  const F& f_;
  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y_;
  const Eigen::Matrix<T_X, Eigen::Dynamic, Eigen::Dynamic>& X_;
  const std::vector<T_theta> theta_;
  Eigen::Matrix< partials_return_type, Eigen::Dynamic, Eigen::Dynamic> C_dbl_;
  Eigen::Matrix< partials_return_type, Eigen::Dynamic, 1> alpha_; // See R&W equation (5.9)
 
   // Constructor 
  gp_lpdf(const F& f, const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,  
                const Eigen::Matrix<T_X, Eigen::Dynamic, Eigen::Dynamic>& X,
               const std::vector<T_theta> theta_) {
    // Do something
  }

  T_return_type operator()() {
    // Evaluate lpdf + gradient if necessary (using propto)
  }
};

A GP would then be specified as something like

GP<true, F, double, double, var> gp_instance(f, y, X, theta);

and the lpdf and its derivatives would be evaluated as, for example, var lp = gp_instance();.

A similar class would need to be built for the non-centred parameterisation.

Current Math Version

v2.18.0