exoplanet-dev
Implement exact Matern terms
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]
Indeed, have not gotten around publishing those explicit state-space representations. Should do so...
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]]