EQ2EMu/EQ2/source/depends/boost_1_72_0/libs/math/doc/fp_utilities/condition_numbers.qbk

[/
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]