New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 977 - Add stable versions of normalize() and normalized()
Add stable versions of normalize() and normalized()
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Core - general
3.3 (current stable)
All All
: Low Feature Request
Assigned To: deanna.m.hood
: accuracy, Documentation, JuniorJob
Depends on:
Blocks: 3.4
  Show dependency treegraph
 
Reported: 2015-03-03 22:56 UTC by Christoph Hertzberg
Modified: 2016-01-23 21:41 UTC (History)
5 users (show)



Attachments

Description Christoph Hertzberg 2015-03-03 22:56:35 UTC
Our current implementation of normalize() and normalized() fails if squaredNorm() of the vector overflows or underflows. 
We could add stable variants of these methods, which use stableNorm() instead of norm().
Special handling for a zero-vector and for the case that the vector contains infinities (or NaNs) could also be implemented. I'd suggest returning UnitX() in the first case. The latter case is debatable. Also see the discussion here:
https://bitbucket.org/eigen/eigen/pull-request/99
Comment 1 deanna.m.hood 2015-03-20 19:58:28 UTC
Isn't it possible that a user might want the behaviour of a non-NaN return value for zero-vectors, but not want to be forced into using stableNorm? These properties (the over/underflow prevention, and the handling of special cases) seem independent.. Perhaps either the latter is applied to the existing normalize methods, or the name of the proposed methods is not simply stableNormalize.

Also I think Eigen would be the first package I've come across to return a unit vector for a zero-vector input. Common practice seems to be zero vector return and/or warning, e.g. http://www.wolframalpha.com/input/?i=normalize%28%5B0%2C0%2C0%5D%29
The direction of a vector with no direction itself should not be set to something seemingly arbitrary, in my opinion.
Comment 2 Christoph Hertzberg 2015-03-27 18:03:04 UTC
We currently generate a unit vector when converting the identity Quaternion to AngleAxis, so we would at least be consistent with ourselves.
But I agree that normalizing a zero vector is generally undefined, and the expected behavior really depends on what the user wants.
So, unless we want to provide methods for every possible use case, we should specify what our implementation of stableNormalized() returns, and leave the implementation of other special cases to the user.

I personally would expect that these conditions hold for every (finite) vector:
  v.normalized().norm()==1 and v.normalized()*v.norm()==v;
Comment 3 Gael Guennebaud 2015-03-27 22:37:06 UTC
I see normalize just like operator/ for which it is up to user to not call it with a zero denominator, and since there will never be a consensus on what should be returned for a zero-input, returning NaN/INF for the default and fast version is good to me.

Regarding the behavior of a hypothetical stableNormalized variant for zero-vector input, let me just add that returning zero would be consistent with how sign/signum/sgn functions are usually defined in 1D. That said, I do not have a strong opinion for either options.
Comment 4 Christoph Hertzberg 2015-03-28 11:49:53 UTC
Yes, I think we agreed that the default normalize[d]() methods shall keep the current fast implementation. However, we should make clear in the documentation that they can return NaNs if squaredNorm() over/underflows.

For the stable variant, the analogy of sign(0)==0 is indeed a good argument for returning a zero vector. In the end, I also don't really mind what we end up implementing.
Actually, the fact that there exist a number of different "best" behaviors, indicates that implementing this might be better left to the user?
Or do we want variants for every combination of
 - Use stableNorm instead of norm
 - vectors with stableNorm (or norm) ==0 return Zero() or UnitX()
 - vectors with stableNorm == Inf return a vector as suggested by Gael in 
   https://bitbucket.org/eigen/eigen/pull-request/99
Comment 5 Rasmus Munk Larsen 2016-01-08 00:03:47 UTC
I would like to add my vote to changing the default behavior of normalize[d] to return a zero-vector unchanged. I keep finding myself (and colleagues) hacking around this, and I have a hard time imagining a case where a NaN return value is actually helpful.

Performance-wise I doubt the real-world effect of changing

*this /= this.norm()

to 

typedef typename NumTraits<typename traits<Derived>::Scalar>::Real RealScalar;
const RealScalar norm = this->norm();
if (norm > RealScalar(0)) {
  *this /= norm;
}

will be dramatic. Even if you are normalizing very short vectors in a tight loop, I would think that branch prediction is going to make the effect rather minor. I don't have benchmark numbers to back this up, but if there is interest in moving in this direction, I can provide them. 

The issue of a hypothetical stableNormalized is orthogonal to this and I think they should have the same semantics regarding zero-vectors.

Rasmus
Comment 6 Gael Guennebaud 2016-01-09 10:15:11 UTC
I benchmarked both using Vector4f:

Here is with gcc:

0.000418993 vs 0.000460624

for which if branch has a non negligible overhead.

And with clang:

0.00037505 0.000376875

it's a tie.

Looking at the generated ASM, the slowdown in gcc (compared to clang) comes from an additional branch to handle corner cases in sqrtss by calling _sqrtf to handle inf/nan. So, replacing:  (*this) /= std::sqrt(sqnorm)
by:

(*this) /= internal::pfirst(_mm_sqrt_ss(_mm_set_ss(sqnorm)));

and we get a tie :)
Comment 7 Rasmus Munk Larsen 2016-01-11 23:21:24 UTC
(In reply to Gael Guennebaud from comment #6)
> I benchmarked both using Vector4f:
> 
> Here is with gcc:
> 
> 0.000418993 vs 0.000460624
> 
> for which if branch has a non negligible overhead.
> 
> And with clang:
> 
> 0.00037505 0.000376875
> 
> it's a tie.
> 
> Looking at the generated ASM, the slowdown in gcc (compared to clang) comes
> from an additional branch to handle corner cases in sqrtss by calling _sqrtf
> to handle inf/nan. So, replacing:  (*this) /= std::sqrt(sqnorm)
> by:
> 
> (*this) /= internal::pfirst(_mm_sqrt_ss(_mm_set_ss(sqnorm)));
> 
> and we get a tie :)

Cool. I guess you can make even make your fix portable by doing 

internal::pfirst(internal::psqrt(internal::pset1(sqnorm)))

__mm_sqrt_ps and __mm_sqrt_ss have the same throughput and latency.
Comment 8 Gael Guennebaud 2016-01-21 21:37:25 UTC
https://bitbucket.org/eigen/eigen/commits/14f468dba4d3/
Summary:     Add numext::sqrt function to enable custom optimized implementation.

https://bitbucket.org/eigen/eigen/commits/12f866a74661/
Summary:     Bug 977: avoid division by 0 in normalize() and normalized().
Comment 9 Gael Guennebaud 2016-01-23 21:41:20 UTC
https://bitbucket.org/eigen/eigen/commits/4808e7f22405/
Summary:     Bug 977: add stableNormalize[d] methods: they are analogues to normalize[d] but with carefull handling of under/over-flow

Note You need to log in before you can comment on or make changes to this bug.