FIX: compare to epsilon during rref #90
Conversation
src/structure/matrix.rs
Outdated
| /// | ||
| /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form) | ||
| fn rref(&self) -> Matrix { | ||
| let epsilon = 1e-10; |
There was a problem hiding this comment.
This particular value of epsilon feels a bit arbitrary when applied universally here.
Would you be interested in considering floating-point constants like machine epsilon or smallest positive normal number? I also assume the choice of a particular threshold constant should be justified formally if possible
There was a problem hiding this comment.
I agree with @djmaxus, but if we use f64::EPSILON, then since it's too strict, so it may not be practical. How's your opinion? @developing-human @djmaxus.
Furthermore, I think it's better to use relative tolerance as follows:
let max_abs = self.data.iter().fold(0f64, |acc, &x| acc.max(x.abs()));
let epsilon = (max_abs * 1e-12).max(1e-15);There was a problem hiding this comment.
This is good feedback, thank you.
I had been thinking of having the function accept a value for epsilon (sympy does similar) but this would be a breaking change.
I think computing epsilon based on the values in the matrix is a good solution. I confirmed it works for my use case and pushed an update.
3 participants
I found a matrix that was numerically unstable while calculating reduced row echelon form. I've updated the comparisons to zero so they instead compare to epsilon, which has been more stable for my use case.
I added a test case for the unstable scenario, as well as a more basic test since one wasn't present yet for rref.
Updating numerical crates is outside my usual wheelhouse, so I'm very open to feedback here!