Open pnevyk opened 10 months ago
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
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?
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)
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
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:
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)
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.
yeah, looks good to me!
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?
Few notes:
Display
andZero
traits for complex numbers was necessary due toDisplay
bound on ElementwiseComparator::Error type andZero
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 formax_relative
as implemented by approx crate). Ideally, we might move these features to thematrixcompare
crate.