Proptest tests

[pnevyk]

Few notes:

• Implementing Display and Zero traits for complex numbers was necessary due to Display bound on ElementwiseComparator::Error type and Zero bound on compare_matrices function. If you feel strongly against, we could find workarounds.
• matrixcompare crate has support for proptest, but I was lacking support for custom comparators (or at least support for max_relative as implemented by approx crate). Ideally, we might move these features to the matrixcompare crate.
• The tolerances used for absolute and relative float comparison are up for discussion. I reused sqrt(epsilon) for absolute difference and epsilon for relative comparison (as this article suggests epsilon or a small multiply of epsilon).

[sarah-ek]

i think using epsilon for relative comparisons is too strict.

in fact, the elementwise absolute/relative comparisons might both be unsuitable in this case. what we'd ideally want is to compare the norm of the difference of vectors/matrices.

this is because matrix operations aren't going to act on elements individually. they're going to combine all the elements in different ways

[pnevyk]

I now compare the numbers using absolute difference with a threshold that is relative to the norms and dimensions of the matrices, hope my attempt is how you imagined it. I am sure you will find improvements to the code.

We still need to use custom macro and comparator if we want to support matrices of complex numbers.

The QR and eigendecomposition tests are now passing. The lower triangular still fails, is the failure valid or do I have a bug in the proptest code?

[sarah-ek]

the failure is caused by the bad conditioning of the input matrix. for solving linear systems, the error also depends on the condition number of the matrix we're inverting (in this case the lower triangular one)

given the singular values singvals of an invertible matrix A, the condition number is defined as max(singvals) / min(singvals)

[sarah-ek]

i'm trying to go through the formulas and i think this should be good

norm_l2(Ax - b) <= eps × (r + r × cond(A)) × norm_l2(b)

where r is some relaxation parameter >= 1, again we could try something like 10 or 100

[pnevyk]

Here is my implementation of eps × (r + r × cond(A)) × norm_l2(b)

pub fn relative_epsilon_cond<E>(
    a: impl AsMatRef<E>,
    b: impl AsMatRef<E>,
    threshold: E::Real,
) -> E::Real
where
    E: ComplexField,
{
    let a = a.as_mat_ref();
    let b = b.as_mat_ref();

    let nrows = a.nrows();
    if nrows == 0 {
        return E::Real::faer_zero();
    }

    let svd = Svd::new(a);

    let mut sing_min = E::Real::faer_from_f64(f64::INFINITY);
    let mut sing_max = E::Real::faer_from_f64(0.0);
    zipped!(svd.s_diagonal()).for_each(|unzipped!(sing)| {
        let sing = sing.read().faer_real();

        if sing < sing_min {
            sing_min = sing;
        }

        if sing > sing_max {
            sing_max = sing;
        }
    });

    let a_cond = sing_max.faer_div(sing_min);
    let b_norm = b.norm_l2();

    let eps = E::Real::faer_epsilon().unwrap();

    eps.faer_mul(threshold.faer_add(threshold.faer_mul(a_cond)))
        .faer_mul(b_norm)
}

Now a few questions:

1. Is my implementation for condition number correct?
2. If so, could it be somehow improved?
3. Proptest is able to generate matrices for which there is a singular value equal to 0. What to do in this case?

[sarah-ek]

the implementation is correct, though you can also use Faer::singular_values (https://docs.rs/faer/latest/faer/trait.Faer.html#tymethod.singular_values) instead of computing the full svd, since that's more efficient as for the issue with singular matrices, maybe we can skip the proptests for solving linear systems, and just merge the existing stuff that works for now (qr and eigendecompositions)

[pnevyk]

I polished the proptest strategy implementation a bit. And I commented out the test for solving linear systems. Can you have a final look? I will then resolve the conflict and we can merge it.

[sarah-ek]

yeah, looks good to me!

[pnevyk]

Rebased on the current main. Apologies for the delay.

What would be the next steps in your opinion? I could add more proptest tests (in a separate PR), maybe you would have some good candidates?