119 lines
3.8 KiB
Plaintext
119 lines
3.8 KiB
Plaintext
[/
|
|
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]
|