Infinite loop in Hypergeometric::sample with large (>2^50) parameters

[mstoeckl]

Summary

With large parameters, sampling from the Hypergeometric distribution implementation can hang in an infinite or very slow loop somewhere. I'd prefer success or an explicit error.

This is probably a consequence of limited precision floating point -- the inputs that trigger this are inching toward the 2^53 boundary . I haven't looked into the cause in detail yet, and will update when I have time to find out more. (Edit: this may be a fundamental issue with large parameters passed to algorithm "HIN", making it calculate a ratio of binomials incredibly slowly.)

Code sample

use rand::SeedableRng;
use rand::rngs::Xoshiro128PlusPlus;
use rand_distr::{Distribution, Hypergeometric};

fn main() {
    let mut rng = Xoshiro128PlusPlus::seed_from_u64(0x1);
    let dist = Hypergeometric::new(5277655813324800, 527765581332480, 52776558133248).unwrap();
    println!("{}", dist.sample(&mut rng));
}

Notes

While I encountered this while trying to write some additional property tests, it's possible to systematically look for hang and crash issues of this type with a fuzzer. As testing methods go, fuzzing is somewhat awkward and can require a lot of manual work to explore a complicated system; but most of the rand_distr sampler code doesn't have too complicated control flow or internal checks. My current fuzzer script is in https://github.com/mstoeckl/rand_distr/commits/fuzz-params/ ; if you'd like I can make a pull request out of it (or start reporting issues that come from it, or just PRs to fix them -- there are a few more distributions that can crash or hang on unusual inputs.)

[dhardy]

I haven't looked into the cause in detail yet, and will update when I have time to find out more.

Let us know if you do narrow this down.

[mstoeckl]

I haven't looked into the cause in detail yet, and will update when I have time to find out more.

Let us know if you do narrow this down.

The root cause appears to be that fraction_of_products_of_factorials(), invoked for algorithm HIN when the peak of the distribution is narrow and close to zero (well, more precisely, M - max(0, k-n2) <= 10, so there are other cases), scales slowly on large inputs; it can run for 2^50 iterations, which is slow enough to appear infinite.

Right now I think the best solution would be to write an efficient and accurate implementation for this function. Calculating the ratio naively (using ln_gamma()) runs into limited precision issues when subtracting large approximately calculated values, but it should be possible to modify existing calculation methods for ln_gamma to directly implement fraction_of_products_of_factorials() and only subtract u64 values. Or maybe just implement ln_gamma using a higher precision type (e.g. https://docs.rs/twofloat/latest/twofloat/). I'm looking into this approach in particular, because having a way to explicitly and accurately evaluate the probability of the hypergeometric distribution at any point would also be very useful when testing the distribution sampler.

Edit: implementing fraction_of_products_of_factorials() with a Stirling-like series to compute log gamma terms, plus some algebra to simplify the formula so significant subtractions are only done with u64 values, plus direct calculation of small factors, seems to work relatively well, although I haven't confirmed that it actually is accurate enough on all possible inputs (the logs of ratios may need improvement).

Example, skipping small factor handling logic

/// Given l > h > s, estimate ln of (h! / (h-s)!) / (l! (l - s)!) using Stirling's formula
/// expanded out to the -9th power terms. This may be a bit inaccurate when h-s is small,
/// so it's best to run it with (h - s) > Ω(1) so that 1/(h - s)^9 is negligible.
fn falling_factorial_log_ratio(h: u64, l: u64, s: u64) -> f64 {
    assert!(l > h);
    assert!(h > s);

    let l_h = (l - h) as f64;
    let h_s = (h - (s-1)) as f64;
    let l_s = (l - (s-1)) as f64;

    let lp1 = l.saturating_add(1) as f64;
    let hp1 = (h + 1) as f64;

    let log_lhs_hls = f64::ln(lp1 * h_s / (hp1 * l_s));
    let log_h_hs = f64::ln(hp1 / h_s);
    let log_ls_hs = f64::ln(l_s / h_s);
    let t0 = lp1 * log_lhs_hls + l_h * log_h_hs + (s as f64) * log_ls_hs;

    let t1 = -0.5 * log_lhs_hls;

    let t2 = (1.0 / 12.0) * (1.0 / lp1 - 1.0 / l_s - 1.0 / hp1 + 1.0 / h_s);

    let t3 = -(1.0 / 360.0)
        * (1.0 / lp1.powi(3) - 1.0 / l_s.powi(3) - 1.0 / hp1.powi(3) + 1.0 / h_s.powi(3));

    let t4 = (1.0 / 1260.0)
        * (1.0 / lp1.powi(5) - 1.0 / l_s.powi(5) - 1.0 / hp1.powi(5) + 1.0 / h_s.powi(5));

    let t5 = -(1.0 / 1680.0)
        * (1.0 / lp1.powi(7) - 1.0 / l_s.powi(7) - 1.0 / hp1.powi(7) + 1.0 / h_s.powi(7));

    let t6 = (1.0 / 1188.0)
        * (1.0 / lp1.powi(9) - 1.0 / l_s.powi(9) - 1.0 / hp1.powi(9) + 1.0 / h_s.powi(9));

    let ret = -(t0 + (t1 + (t2 + (t3 + (t4 + (t5 + t6))))));
    assert!(ret.is_finite() && ret <= 0.0);
    ret
}

Edit 2: actually, the InverseTransform method itself has a loop that can on certain inputs run for a huge number of iterations, so maybe more significant changes are needed.