math/tools/condition_numbers.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
0007 #define BOOST_MATH_TOOLS_CONDITION_NUMBERS_HPP
0008 #include <cmath>
0009 #include <limits>
0010 #include <boost/math/differentiation/finite_difference.hpp>
0011 #include <boost/math/tools/config.hpp>
0012
0013 namespace boost { namespace math { namespace tools {
0014
0015 template<class Real, bool kahan=true>
0016 class summation_condition_number {
0017 public:
0018     summation_condition_number(Real const x = 0)
0019     {
0020         using std::abs;
0021         m_l1 = abs(x);
0022         m_sum = x;
0023         m_c = 0;
0024     }
0025
0026     void operator+=(Real const & x)
0027     {
0028         using std::abs;
0029         // No need to Kahan the l1 calc; it's well conditioned:
0030         m_l1 += abs(x);
0031         BOOST_IF_CONSTEXPR (kahan)
0032         {
0033             Real y = x - m_c;
0034             Real t = m_sum + y;
0035             m_c = (t-m_sum) -y;
0036             m_sum = t;
0037         }
0038         else
0039         {
0040             m_sum += x;
0041         }
0042     }
0043
0044     inline void operator-=(Real const & x)
0045     {
0046         this->operator+=(-x);
0047     }
0048
0049     // Is operator*= relevant? Presumably everything gets rescaled,
0050     // (m_sum -> k*m_sum, m_l1->k*m_l1, m_c->k*m_c),
0051     // but is this sensible? More important is it useful?
0052     // In addition, it might change the condition number.
0053
0054     Real operator()() const
0055     {
0056         using std::abs;
0057         if (m_sum == Real(0) && m_l1 != Real(0))
0058         {
0059             return std::numeric_limits<Real>::infinity();
0060         }
0061         return m_l1/abs(m_sum);
0062     }
0063
0064     Real sum() const
0065     {
0066         // Higham, 1993, "The Accuracy of Floating Point Summation":
0067         // "In [17] and [18], Kahan describes a variation of compensated summation in which the final sum is also corrected
0068         // thus s=s+e is appended to the algorithm above)."
0069         return m_sum + m_c;
0070     }
0071
0072     Real l1_norm() const
0073     {
0074         return m_l1;
0075     }
0076
0077 private:
0078     Real m_l1;
0079     Real m_sum;
0080     Real m_c;
0081 };
0082
0083 template<class F, class Real>
0084 Real evaluation_condition_number(F const & f, Real const & x)
0085 {
0086     using std::abs;
0087     using std::isnan;
0088     using std::sqrt;
0089     using boost::math::differentiation::finite_difference_derivative;
0090
0091     Real fx = f(x);
0092     if (isnan(fx))
0093     {
0094         return std::numeric_limits<Real>::quiet_NaN();
0095     }
0096     bool caught_exception = false;
0097     Real fp;
0098 #ifndef BOOST_NO_EXCEPTIONS
0099     try
0100     {
0101 #endif
0102         fp = finite_difference_derivative(f, x);
0103 #ifndef BOOST_NO_EXCEPTIONS
0104     }
0105     catch(...)
0106     {
0107         caught_exception = true;
0108     }
0109 #endif
0110     if (isnan(fp) || caught_exception)
0111     {
0112         // Check if the right derivative exists:
0113         fp = finite_difference_derivative<decltype(f), Real, 1>(f, x);
0114         if (isnan(fp))
0115         {
0116             // Check if a left derivative exists:
0117             const Real eps = (std::numeric_limits<Real>::epsilon)();
0118             Real h = - 2 * sqrt(eps);
0119             h = boost::math::differentiation::detail::make_xph_representable(x, h);
0120             Real yh = f(x + h);
0121             Real y0 = f(x);
0122             Real diff = yh - y0;
0123             fp = diff / h;
0124             if (isnan(fp))
0125             {
0126                 return std::numeric_limits<Real>::quiet_NaN();
0127             }
0128         }
0129     }
0130
0131     if (fx == 0)
0132     {
0133         if (x==0 || fp==0)
0134         {
0135             return std::numeric_limits<Real>::quiet_NaN();
0136         }
0137         return std::numeric_limits<Real>::infinity();
0138     }
0139
0140     return abs(x*fp/fx);
0141 }
0142
0143 }}} // Namespaces
0144 #endif