Address numerical stability in exact sampling of k DPPs for large matrices

[DhruvaKashyap]

Problem

The current k-DPP sampling algorithm using GS is not numerically stable and frequently runs into division warnings when dealing with large matrices (~10kx10k) with large condition numbers (~10^6) and low stable rank (~10).

Solution

I believe this can be traced to the computation of the elementary symmetric polynomials and this check

I have implemented a numerically stable recursion by working in the log space of the elementary symmetric polynomials and changing the check mentioned above. This can also be moved to GPUs using torch instead of numpy

Would be happy to contribute to this project.

[DhruvaKashyap]

@guilgautier Just a friendly bump in case it slipped through.