diff options
author | Gael Guennebaud <g.gael@free.fr> | 2016-01-23 22:40:11 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2016-01-23 22:40:11 +0100 |
commit | 1cf85bd875ecbcfa1240b4ec08122d40d79101fd (patch) | |
tree | 64dacc9ef8e8f8fd6d4de1a707fff0c3d73c917f | |
parent | 369d6d1ae31c3e1a0f03196ccb9c792c6913ed76 (diff) |
bug #977: add stableNormalize[d] methods: they are analogues to normalize[d] but with carefull handling of under/over-flow
-rw-r--r-- | Eigen/src/Core/Dot.h | 46 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 2 | ||||
-rw-r--r-- | test/stable_norm.cpp | 14 |
3 files changed, 62 insertions, 0 deletions
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 221fc3224..82d58fc0b 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -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 { diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f3935802d..338879c73 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -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(); diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index 7561ae8be..9f12320e0 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp @@ -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() |