dimforge
Wrong result of nalgebra::base::Matrix::min() and ::max()
The methods nalgebra::base::Matrix::min() and nalgebra::base::Matrix::max() are documented as
Returns the component with the smallest value.
and
Returns the component with the largest value.
However, A.min() actually returns N::zero().min(A[..]) (if this was a valid expression) and A.max() returns N::zero().max(A[..]). For example, this corresponds to
na::Vector3::new(-1.0, -2.0, -3.0).max()returning 0.0na::Vector3::new(1.0, 2.0, 3.0).min()returning 0.0
This is due to the fact that both rely on xcmp which always initializes the comparison value with N::zero() instead of e.g. the first value of the matrix. This is either a bug in the code or in the documentation.
Currently I have to use A[A.imin()] instead of A.min().
Thank you for reporting this! This is definitely a bug in the code. This should be fixed by replacing the initial value with the first iterator element. Something like:
fn xcmp<N2>(&self, abs: impl Fn(N) -> N2, cmp: impl Fn(N2, N2) -> bool) -> N2
where N2: Scalar + PartialOrd + Zero {
let mut iter = self.iter();
let mut max = abs(it.next().cloned().unwrap_or(N::zero()));
for e in iter {
let ae = abs(*e);
if cmp(ae, max) {
max = ae;
}
}
max
}
@sebcrozet: Note that this would give a min/max of T::zero() for matrices whose col or row count is 0 (which might frequently happen for edge cases of algorithms that need to dynamically allocate matrices). We might want to think about the desired behavior for this edge case.
I've had exactly this issue with Eigen (C++) before: they didn't properly/consistently handle edge cases for 0x0 matrices, which led to some very painful workarounds for edge cases in various algorithms that I was writing.
Another, possibly more important point, is the handling of NaNs. It seems to me that it if you have a vector [1, NaN, 3], it would be desirable to have min/max/ both return NaN, otherwise you might easily introduce silent errors in code (as you would never observe the NaN if all you ever do is call min/max).
Of course, this is part of the reason for why e.g. Iterator::min() and Iterator::max() require Ord in the first place... (I realize the reason that nalgebra only requires PartialOrd is pragmatism, but still)
Just some things to think about!
Wrt to the NaN problem I would simply propose
fn xcmp<N2>(&self, abs: impl Fn(N) -> N2, cmp: impl Fn(N2, N2) -> bool) -> N2
where N2: Scalar + PartialOrd + Zero {
let mut iter = self.iter();
let mut max = abs(it.next().cloned().unwrap_or(N::zero()));
for e in iter {
if e.is_nan() {
max = *e;
break;
}
let ae = abs(*e);
if cmp(ae, max) {
max = ae;
}
}
max
}
Wrt the zero sized matrices I would also suggest to keep returning 0 but additionally document this. Of course this should be the same behaviour for all operations that reduce matrices to a scalar, such as norm and dot. The reason that I would keep it this way is that currently all matrix -> matrix operations on zero sized matrices simply return a zero sized matrix. Returning 0 from functions returning a scalar allows code like
a = DMatrix::zeros(0,0);
b = (&a * &a).scale(a.max());
to just work. If this could break your algorithm I guess it would be the most sensible choice to check for zero sized matrices in the beginning anyway. What is your opinion on this?
btw. in Numpy there is the additional nanmax, etc. function which specifically does not propagate NaN
I feel like a panic when performing an operation that requires elements on an empty matrix would be less prone to unpleasant surprises. Note that Iterator::min returns Option.
What about try_min() returning an Option, and min panicking if the matrix has no elements? That way, you can choose a sensible value when the matrix is empty by e.g. .try_min().unwrap_or(0.0).
As for the NaN question: PartialOrd is designed such that it gives an Option<Ordering>. Thus, one could consider something like
for e in iter {
if let Some(ordering) = e.cmp(min) {
if ordering == Ordering::Less {
min = e;
}
} else if e.is_nan() {
// This way we don't check `is_nan` unless the two elements are "uncomparable",
// which for e.g. floating point numbers I believe only happens with NaNs?
return e;
}
}
As mentioned in the above comment in the code, this way you would only check explicitly for NaN when you have reason to suspect that the element might be NaN. How this actually ends up optimizing in the end is anybody's guess though I suppose.
Ok, panicking with empty matrices is sensible. But I guess this should be changed/documented for all scalar returning methods.
I just noticed that we can't actually call is_nan in the comparison function because N might be an integer. I guess it would be ok to just drop the if e.is_nan()?
Ah, yes, I also didn't consider the fact that we of course don't have a RealField bound here. Hmm, I think this could have very unintended consequences. Consider something like a Rational-like type, and consider comparing 2/4 with 1/2. I think that there are two valid choices here:
- define
2/4 == 1/2 - define the comparison to be "undefined", i.e.
partial_cmpreturnsNone. That is, they are neither <, > nor ==.
I'm not 100% sure that the second option here would be a reasonable implementation for PartialOrd, but if it is, then clearly one gets very unintended consequences when used in e.g. min or max.
On a separate note, changing the semantics of min, max would be a breaking change. It might be worthwhile to considering a separate PR that only fixes the initial bug report, which is clearly the biggest problem here, with empty matrix handling and/or NaN-handling being somewhat more minor in nature.
@Andlon Wrt to the PartialOrd problem: it's definitely something to consider but do you think there is a proper solution to this other than deciding on how to indicate this case as a failure (e.g. using an Option or by returning the first offending value)? Maybe it makes sense to split min/max in a function for PartialOrd scalars that returns an Option and one for Ord scalars?
@w1th0utnam3: I think the Ord/PartialOrd split basically is the solution here. I know nalgebra implements max/min with PartialOrd alone for pragmatic reasons, but as we see from this discussion, this is not without (arguably somewhat severe) issues. Of course, this is not really satisfying as users typically just want to compute the damn min/max value regardless.
Ideally we'd have a trait separate from both Ord and PartialOrd that would admit max and min by basically assuming a total ordering for any well-defined element (which is the case for floating point numbers and integers alike), but letting the number type admit "NaN" type values, in which case min/max of any matrix containing NaN would be defined explicitly to be NaN.
I think possibly alga's MeetSemilattice and JoinSemilattice would fit that the requirement if they would consistently return NaN for meet(a, b) and join(a, b) when either a or b is NaN. @sebcrozet: I don't know if you have any particular opinion on how alga should define meet/join for NaN values? (As far as I can see there's nothing in the documentation/implementation to suggest that this has been considered with particular care?)
I think the Ord/PartialOrd split basically is the solution here. I know nalgebra implements max/min with PartialOrd alone for pragmatic reasons, but as we see from this discussion, this is not without (arguably somewhat severe) issues. Of course, this is not really satisfying as users typically just want to compute the damn min/max value regardless.
We could for example add other methods that requires Ord instead of PartialOrd. A possibility on the line of your try_min suggestion is to have this try_min return an Result, where the error is an enum that can either be EmptyMatrix, or InvalidComparison.
@sebcrozet: I don't know if you have any particular opinion on how alga should define meet/join for NaN values? (As far as I can see there's nothing in the documentation/implementation to suggest that this has been considered with particular care?)
I don't believe it is valid for a single value to be both the supremum and the infimum of a set. So it would likely make no sense to return NaN for both meet and join if either is NaN. Though we could say that NaN is the only infimum for example (and document it).
Also I'm not sure what to do with NaN with different bit patterns.
In general I like the idea of try_* returning a Result. But that still leaves the question on how to handle the discussed problems in the non-try case
- empty matrix problem with solutions
a)
assert!(self.len() > 0)b)return N::zero() PartialOrdproblem with the solution approaches a) requiringOrd, which makes them unavailable forRealFieldwhich is a very common use case b) returning the first value that returnsNoneforpartial_cmpc) ignoringNoneresults, which would lead to a similar problem as a), as a matrix of onlyNaNvalues would then essentially be empty...
I would be fine with options 1. a) and 2. b). Personally, I don't like "removing" the methods for PartialOrd types because I'm used to this behavior from other math libraries, however I guess this might be the more rusty approach?
@w1th0utnam3 : In the non-try case, I would suggest that it panics if it encounters None, and that it's documented that it's the user's responsibility that every element can be ordered. In practice, for something like floats, this means that the user needs to guarantee that no element is a NaN. I think your option 2b) can lead to some serious bugs for unsuspecting users.
@sebcrozet: I believe the mathematical properties of MeetSemilattice are already nonsensical when comparing two NaNs, so I'm not sure if explicitly defining NaNs to be permissible as both supremum and infimum makes matters much worse.
The result option is in some sense a good solution, but the ergonomics of this solution are terrible. But I also don't really see a way around that.
An alternate option would perhaps be to introduce a brand new trait that would represent something like a quasi-total ordering of elements. Something in-between PartialOrd and Ord. You could also give a blanket implementation of Ord in this case. Something like:
/// For types that form a total order of all non-exceptional values.
///
/// If `T` implement `QuasiOrd`, then the values belonging to the set of all values in `T`
/// can be partitioned into two sets: one which admits a total ordering,
/// and one which contain exceptional values which admit no ordering whatsoever.
/// If a set of values of `T` contain only elements of the former set, then it admits a total order.
/// Otherwise, it does not. In particular, the minimum of a set of `T` is well-defined if
/// the set only contains totally orderable values, otherwise it is by definition ill-defined
/// and may take any exceptional value.
trait QuasiOrd {
fn quasi_cmp(a: &T, b: &T) -> Option<Ordering>;
}
impl<T: Ord> QuasiOrd for T {
/// ...
}
The above draft is totally rushed and unpolished. I just wanted to float the idea to see if it makes any sense.
EDIT: One could also perhaps require QuasiOrd: PartialOrd and use the same implementation, except that now QuasiOrd functions merely as a marker trait?
EDIT2: To make it clear, with the above trait, try_min() could once again simply return an Option, as the only failure case would be empty matrices.
EDIT3: I just realized that the above is not sufficient, as when comparing e.g. a NaN and 3.4, you wouldn't be able to tell which one is the actual offending value. There are multiple options to deal with this though, all of whom require some changes to the trait presented above.
