I start a new thread to continue the discussion from here Gaussian Process Predictive Posterior Function

I think this is specially for @Bob_Carpenter

Assume we have data `x`

as a NxD matrix, with N the number of observations and D the number of dimensions.

We would have a user defined function

```
real my_gp_covfun(vector x1, vector x2, vector params)
```

where vector x1, and vector x2 have length D, and they would usually be two rows of `x`

. vector params includes parameters of the covariance function.

We would have a function added in math and exposed in language

```
matrix gp_cov(function covfun, matrix x, vector params)
```

which we would call as

```
C = gp_cov(my_gp_covfun, x, params)
```

or

```
L = gp_cov_chol(my_gp_covfun, x, params)
```

(or something like this)

Inside `gp_cov`

and `gp_cov_chol`

we would call the user defined function, e.g. `my_gp_covfun`

to construct the covariance matrix by calling it N(N+1)/2 times (need to compute only triangular matrix) with different row combinations of matrix `x`

.

Question: If we would write other parts of gradient computation in C++ similar to what now is done inside `cov_exp_quad`

, can we use autodiff to compute just the gradients of user defined function `my_gp_covfun`

when we need the derivative of the covariance matrix with respect to one covariance function parameter in `params`

? If this works, then we can continue with hand written matrix derivative operations to get the scalar gradient for that parameter, and then repeat for the the next covariance function parameter in `params`

. This could lead to less memory use and possibility to use user defined covariance functions would give a lot of flexibility in defining models.