Based on "Reducing the Variance of Gaussian Process Hyperparameter Optimization with Preconditioning" by Wenger et al., 2021.

When using iterative methods (i.e. CG/SLQ) to compute the log determinant, the forward pass currently computes:
logdet K \approx logdet P + SLQ( P^{-1/2} K P^{-1/2} ),
where P is a preconditioner, and SLQ is a stochastic estimate of the log determinant. If the preconditioner is a good approximation of K, then this forward pass can be seen as a form of variance reduction.

In this PR, we apply this same variance reduction strategy to the backward pass. We compute the backward pass as:
d logdet(K)/dtheta \approx d logdet(P)/dtheta + d SLQ/dtheta

TODOs:

• Implement pivoted cholesky as a torch.autograd.Function, so that we can compute backward passes through it.
• Redo inv_quad_logdet function to apply variance reduction in the forward and backward passes.