Implement exact Matern terms

[dfm]

Demo here: https://gist.github.com/dfm/01d4b172d5cc3b18f34143aeb2f8cb8e

[dfm]

I spent some more time on this and it is possible to implement this in a numerically stable by changing the update step from diag(P) to a block diagonal matrix. For Matern-3/2 the P block is:

[[1 + f * dt,   dt],
 [-f ** 2 * dt, 1 - f * dt]]

where f = sqrt(3) / rho and U = V = [1, 0]. This result (and the Matern-5/2 result) was previously derived by @andres-jordan et al. and implemented here https://github.com/andres-jordan/gpstate (although I think that the paper is not available anywhere yet?). Would be nice to talk about the relationship to state-space models and perhaps provide an interface for general polynomials. But, this will require some updates to the backend so I'm going to leave this as a stretch goal.

Matmul lower in Python:

def get_p(dt):
    return np.array([[1 + f * dt, dt], [-f ** 2 * dt, 1 - f * dt]])

z = np.zeros_like(y)
F = np.zeros((2, y.shape[1]))
for n in range(1, N):
    dt = t[n] - t[n - 1]
    p = np.exp(-f * dt) * get_p(dt)
    F[0, :] += y[n - 1]
    F = np.dot(p, F)
    z[n] = sigma * F[0]

[andres-jordan]

Indeed, have not gotten around publishing those explicit state-space representations. Should do so...

[dfm]

Since I don't have anywhere else to put it...

For the derivatives discussed in #17, it's useful to know how to factor that matrix into a product where the left matrix only depends on tn and the right matrix only on tm. In this case, that factorization is:

L = [[1 / f + tn,  1],
     [-f * tn,    -f]]

R = [[f,       1          ],
     [-f * tm, -1 / f - tm]]