How does the Pivoted Cholesky Decomposition implementation differ from Harbrecht et.al. 2012 ?

[kunalghosh]

Hi GPyTorch devs, I am trying to implement the mBCG algorithm from the original GPyTorch paper for PyMC. I am currently working on implementing Pivoted Cholesky decomposition and have implemented the one from the Original Harbrecht et.al 2012 paper

To me it seems like the implementation in GPyTorch differs from the original in three places,

1. The inner for-loop runs until max_iter or rank instead of n
2. Initial error is max(diag) instead of l1_norm(diag)
3. Error computed in the loop is divided by the Initial error

Are there any other differences ?

I wrote a simple test (with my Pivoted Cholesky implementation in NumPy) and the results differ from GPyTorch. link to the notebook

Any help / pointers would be very helpful. (tagging @gpleiss since he implemented the algorithm on GPyTorch)

[gpleiss]

The inner for-loop runs until max_iter or rank instead of n

You could fix this by setting gpytorch.settings.max_preconditioner_size to n. We have a maximum number of iterations because we are primarily using this function for preconditioning.

Initial error is max(diag) instead of l1_norm(diag)

I don't think that's true?

gpytorch/gpytorch/functions/_pivoted_cholesky.py

Line 32 in 78796f7

errors = torch.norm(matrix_diag, 1, dim=-1) / orig_error

Error computed in the loop is divided by the Initial error

Most of our numerical methods in GPyTorch use relative errors, not absolute errors.

It sounds like the only difference between your implementation and ours is the stopping conditions. If you need strict stopping conditions based on absolute error, then our implementation probably isn't suitable for your needs.