I experimented locally with a change to make the two factorize methods consistent and could make a pull request soon.
However, I am wondering if we should modify all the cholesky calls (also in the dense factorize) to use the pivoted Cholesky decomposition, which also works on rank deficient matrices and is significantly faster than bunchkaufman in those cases. Opinions?