123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118 |
- [/
- Copyright (c) 2019 Nick Thompson
- Use, modification and distribution are subject to the
- Boost Software License, Version 1.0. (See accompanying file
- LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- ]
- [section:cond Condition Numbers]
- [heading Synopsis]
- ``
- #include <boost/math/tools/condition_numbers.hpp>
- namespace boost::math::tools {
- template<class Real, bool kahan=true>
- class summation_condition_number {
- public:
- summation_condition_number(Real const x = 0);
- void operator+=(Real const & x);
- inline void operator-=(Real const & x);
- [[nodiscard]] Real operator()() const;
- [[nodiscard]] Real sum() const;
- [[nodiscard]] Real l1_norm() const;
- };
- template<class F, class Real>
- auto evaluation_condition_number(F const & f, Real const & x);
- } // namespaces
- ``
- [heading Summation Condition Number]
- Here we compute the condition number of the alternating harmonic sum:
- using boost::math::tools::summation_condition_number;
- auto cond = summation_condition_number<float, /* kahan = */ false>();
- float max_n = 10000000;
- for (float n = 1; n < max_n; n += 2)
- {
- cond += 1/n;
- cond -= 1/(n+1);
- }
- std::cout << std::setprecision(std::numeric_limits<float>::digits10);
- std::cout << "ln(2) = " << boost::math::constants::ln_two<float>() << "\n";
- std::cout << "Sum = " << cond.sum() << "\n";
- std::cout << "Condition number = " << cond() << "\n";
- Output:
- ln(2) = 0.693147
- Sum = 0.693137
- Condition number = 22.22282
- We see that we have lost roughly two digits of accuracy,
- consistent with the heuristic that if the condition number is 10[super /k/],
- then we lose /k/ significant digits in the sum.
- Our guess it that if you're worried about whether your sum is ill-conditioned,
- the /last/ thing you want is for your condition number estimate to be inaccurate as well.
- Since the condition number estimate relies on computing the (perhaps ill-conditioned) sum,
- we have defaulted the accumulation to use Kahan summation:
- auto cond = boost::math::tools::summation_condition_number<float>(); // will use Kahan summation.
- // ...
- Output:
- ln(2) = 0.693147
- Kahan sum = 0.693147
- Condition number = 22.2228
- If you are interested, the L[sub 1] norm is also generated by this computation, so you may query it if you like:
- float l1 = cond.l1_norm();
- // l1 = 15.4
- [heading Condition Number of Function Evaluation]
- The [@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].
- It is computed as follows:
- using boost::math::tools::evaluation_condition_number;
- auto f = [](double x)->double { return std::log(x); };
- double x = 1.125;
- double cond = evaluation_condition_number(f, 1.125);
- // cond = 1/log(x)
- [heading Caveats]
- The condition number of function evaluation gives us a measure of how sensitive our function is to roundoff error.
- Unfortunately, evaluating the condition number requires evaluating the function and its derivative,
- and this calculation is itself inaccurate whenever the condition number of function evaluation is large.
- Sadly, this is also the regime when you are most interested in evaluating a condition number!
- This seems to be a fundamental problem.
- However, it should not be necessary to evaluate the condition number to high accuracy,
- valuable insights can be obtained simply by looking at the change in condition number as the function
- evolves over its domain.
- [heading References]
- * Gautschi, Walter. ['Orthogonal polynomials: computation and approximation] Oxford University Press on Demand, 2004.
- * Higham, Nicholas J. ['The accuracy of floating point summation.] SIAM Journal on Scientific Computing 14.4 (1993): 783-799.
- * Higham, Nicholas J. ['Accuracy and stability of numerical algorithms.] Vol. 80. Siam, 2002.
- [endsect]
|