Follow up to #2
we should probably store
C_S = cholesky(W.S)
C_B
for
B = C_S.U' \ (W.A * cholesky(W.D).U')
C_B = cholesky(Symmetric(I + B'B))
as well as the original matrix.
As these tend to be needed for many operations, and don't need to be recomputed as long as the matrix isn't changed.
One question is making sure this works OK with AD.