mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #977: add stableNormalize[d] methods: they are analogues to normalize[d] but with carefull handling of under/over-flow

Browse Source
This commit is contained in:
Gael Guennebaud 2016-01-23 22:40:11 +01:00
parent 369d6d1ae3
commit 1cf85bd875
3 changed files with 62 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
Eigen/src/Core/Dot.h
View File

@ -142,6 +142,52 @@ inline void MatrixBase<Derived>::normalize()
derived() /= numext::sqrt(z);
}
/** \returns an expression of the quotient of \c *this by its own norm while avoiding underflow and overflow.
*
* \only_for_vectors
*
* This method is analogue to the normalized() method, but it reduces the risk of
* underflow and overflow when computing the norm.
*
* \warning If the input vector is too small (i.e., this->norm()==0),
* then this function returns a copy of the input.
*
* \sa stableNorm(), stableNormalize(), normalized()
*/
template<typename Derived>
inline const typename MatrixBase<Derived>::PlainObject
MatrixBase<Derived>::stableNormalized() const
{
typedef typename internal::nested_eval<Derived,3>::type _Nested;
_Nested n(derived());
RealScalar w = n.cwiseAbs().maxCoeff();
RealScalar z = (n/w).squaredNorm();
if(z>RealScalar(0))
return n / (numext::sqrt(z)*w);
else
return n;
}
/** Normalizes the vector while avoid underflow and overflow
*
* \only_for_vectors
*
* This method is analogue to the normalize() method, but it reduces the risk of
* underflow and overflow when computing the norm.
*
* \warning If the input vector is too small (i.e., this->norm()==0), then \c *this is left unchanged.
*
* \sa stableNorm(), stableNormalized(), normalize()
*/
template<typename Derived>
inline void MatrixBase<Derived>::stableNormalize()
{
RealScalar w = cwiseAbs().maxCoeff();
RealScalar z = (derived()/w).squaredNorm();
if(z>RealScalar(0))
derived() /= numext::sqrt(z)*w;
}
//---------- implementation of other norms ----------
namespace internal {

2
Eigen/src/Core/MatrixBase.h
View File

@ -204,7 +204,9 @@ template<typename Derived> class MatrixBase
RealScalar blueNorm() const;
RealScalar hypotNorm() const;
EIGEN_DEVICE_FUNC const PlainObject normalized() const;
EIGEN_DEVICE_FUNC const PlainObject stableNormalized() const;
EIGEN_DEVICE_FUNC void normalize();
EIGEN_DEVICE_FUNC void stableNormalize();
EIGEN_DEVICE_FUNC const AdjointReturnType adjoint() const;
EIGEN_DEVICE_FUNC void adjointInPlace();

14
test/stable_norm.cpp
View File

@ -163,6 +163,20 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
VERIFY(!(numext::isfinite)(v.blueNorm())); VERIFY((numext::isnan)(v.blueNorm()));
VERIFY(!(numext::isfinite)(v.hypotNorm())); VERIFY((numext::isnan)(v.hypotNorm()));
}
// stableNormalize[d]
{
VERIFY_IS_APPROX(vrand.stableNormalized(), vrand.normalized());
MatrixType vcopy(vrand);
vcopy.stableNormalize();
VERIFY_IS_APPROX(vcopy, vrand.normalized());
VERIFY_IS_APPROX((vrand.stableNormalized()).norm(), RealScalar(1));
VERIFY_IS_APPROX(vcopy.norm(), RealScalar(1));
VERIFY_IS_APPROX((vbig.stableNormalized()).norm(), RealScalar(1));
VERIFY_IS_APPROX((vsmall.stableNormalized()).norm(), RealScalar(1));
VERIFY_IS_APPROX(vbig, vbig.stableNorm() * vbig.stableNormalized());
VERIFY_IS_APPROX(vsmall, vsmall.stableNorm() * vsmall.stableNormalized());
}
}
void test_stable_norm()