1[/ 2Copyright (c) 2019 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:cond Condition Numbers] 9 10[heading Synopsis] 11 12`` 13#include <boost/math/tools/condition_numbers.hpp> 14 15 16namespace boost::math::tools { 17 18template<class Real, bool kahan=true> 19class summation_condition_number { 20public: 21 summation_condition_number(Real const x = 0); 22 23 void operator+=(Real const & x); 24 25 inline void operator-=(Real const & x); 26 27 [[nodiscard]] Real operator()() const; 28 29 [[nodiscard]] Real sum() const; 30 31 [[nodiscard]] Real l1_norm() const; 32}; 33 34template<class F, class Real> 35auto evaluation_condition_number(F const & f, Real const & x); 36 37} // namespaces 38`` 39 40[heading Summation Condition Number] 41 42Here we compute the condition number of the alternating harmonic sum: 43 44 using boost::math::tools::summation_condition_number; 45 auto cond = summation_condition_number<float, /* kahan = */ false>(); 46 float max_n = 10000000; 47 for (float n = 1; n < max_n; n += 2) 48 { 49 cond += 1/n; 50 cond -= 1/(n+1); 51 } 52 std::cout << std::setprecision(std::numeric_limits<float>::digits10); 53 std::cout << "ln(2) = " << boost::math::constants::ln_two<float>() << "\n"; 54 std::cout << "Sum = " << cond.sum() << "\n"; 55 std::cout << "Condition number = " << cond() << "\n"; 56 57Output: 58 59 ln(2) = 0.693147 60 Sum = 0.693137 61 Condition number = 22.22282 62 63We see that we have lost roughly two digits of accuracy, 64consistent with the heuristic that if the condition number is 10[super /k/], 65then we lose /k/ significant digits in the sum. 66 67Our guess it that if you're worried about whether your sum is ill-conditioned, 68the /last/ thing you want is for your condition number estimate to be inaccurate as well. 69Since the condition number estimate relies on computing the (perhaps ill-conditioned) sum, 70we have defaulted the accumulation to use Kahan summation: 71 72 auto cond = boost::math::tools::summation_condition_number<float>(); // will use Kahan summation. 73 // ... 74 75Output: 76 77 ln(2) = 0.693147 78 Kahan sum = 0.693147 79 Condition number = 22.2228 80 81If you are interested, the L[sub 1] norm is also generated by this computation, so you may query it if you like: 82 83 float l1 = cond.l1_norm(); 84 // l1 = 15.4 85 86[heading Condition Number of Function Evaluation] 87 88The [@https://en.wikipedia.org/wiki/Condition_number condition number] of function evaluation is defined as the absolute value of /xf/'(/x/)/f/(/x/)[super -1]. 89It is computed as follows: 90 91 using boost::math::tools::evaluation_condition_number; 92 auto f = [](double x)->double { return std::log(x); }; 93 double x = 1.125; 94 double cond = evaluation_condition_number(f, 1.125); 95 // cond = 1/log(x) 96 97[heading Caveats] 98 99The condition number of function evaluation gives us a measure of how sensitive our function is to roundoff error. 100Unfortunately, evaluating the condition number requires evaluating the function and its derivative, 101and this calculation is itself inaccurate whenever the condition number of function evaluation is large. 102Sadly, this is also the regime when you are most interested in evaluating a condition number! 103 104This seems to be a fundamental problem. 105However, it should not be necessary to evaluate the condition number to high accuracy, 106valuable insights can be obtained simply by looking at the change in condition number as the function 107evolves over its domain. 108 109 110[heading References] 111 112* Gautschi, Walter. ['Orthogonal polynomials: computation and approximation] Oxford University Press on Demand, 2004. 113 114* Higham, Nicholas J. ['The accuracy of floating point summation.] SIAM Journal on Scientific Computing 14.4 (1993): 783-799. 115 116* Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Vol. 80. Siam, 2002. 117 118[endsect] 119