1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Condition Numbers</title>
<link rel="stylesheet" href="../math.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
<link rel="home" href="../index.html" title="Math Toolkit 4.2.1">
<link rel="up" href="../utils.html" title="Chapter 2. Floating Point Utilities">
<link rel="prev" href="float_comparison.html" title="Floating-point Comparison">
<link rel="next" href="ulps_plots.html" title="ULPs Plots">
<meta name="viewport" content="width=device-width, initial-scale=1">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td>
<td align="center"><a href="../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="float_comparison.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../utils.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="ulps_plots.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="math_toolkit.cond"></a><a class="link" href="cond.html" title="Condition Numbers">Condition Numbers</a>
</h2></div></div></div>
<h4>
<a name="math_toolkit.cond.h0"></a>
<span class="phrase"><a name="math_toolkit.cond.synopsis"></a></span><a class="link" href="cond.html#math_toolkit.cond.synopsis">Synopsis</a>
</h4>
<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">tools</span><span class="special">/</span><span class="identifier">condition_numbers</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
<span class="keyword">namespace</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span> <span class="special">{</span>
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">Real</span><span class="special">,</span> <span class="keyword">bool</span> <span class="identifier">kahan</span><span class="special">=</span><span class="keyword">true</span><span class="special">></span>
<span class="keyword">class</span> <span class="identifier">summation_condition_number</span> <span class="special">{</span>
<span class="keyword">public</span><span class="special">:</span>
<span class="identifier">summation_condition_number</span><span class="special">(</span><span class="identifier">Real</span> <span class="keyword">const</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">0</span><span class="special">);</span>
<span class="keyword">void</span> <span class="keyword">operator</span><span class="special">+=(</span><span class="identifier">Real</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">x</span><span class="special">);</span>
<span class="keyword">inline</span> <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">-=(</span><span class="identifier">Real</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">x</span><span class="special">);</span>
<span class="special">[[</span><span class="identifier">nodiscard</span><span class="special">]]</span> <span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="special">[[</span><span class="identifier">nodiscard</span><span class="special">]]</span> <span class="identifier">Real</span> <span class="identifier">sum</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="special">[[</span><span class="identifier">nodiscard</span><span class="special">]]</span> <span class="identifier">Real</span> <span class="identifier">l1_norm</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span>
<span class="special">};</span>
<span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">F</span><span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Real</span><span class="special">></span>
<span class="keyword">auto</span> <span class="identifier">evaluation_condition_number</span><span class="special">(</span><span class="identifier">F</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">f</span><span class="special">,</span> <span class="identifier">Real</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">x</span><span class="special">);</span>
<span class="special">}</span> <span class="comment">// namespaces</span>
</pre>
<h4>
<a name="math_toolkit.cond.h1"></a>
<span class="phrase"><a name="math_toolkit.cond.summation_condition_number"></a></span><a class="link" href="cond.html#math_toolkit.cond.summation_condition_number">Summation
Condition Number</a>
</h4>
<p>
Here we compute the condition number of the alternating harmonic sum:
</p>
<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">summation_condition_number</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">cond</span> <span class="special">=</span> <span class="identifier">summation_condition_number</span><span class="special"><</span><span class="keyword">float</span><span class="special">,</span> <span class="comment">/* kahan = */</span> <span class="keyword">false</span><span class="special">>();</span>
<span class="keyword">float</span> <span class="identifier">max_n</span> <span class="special">=</span> <span class="number">10000000</span><span class="special">;</span>
<span class="keyword">for</span> <span class="special">(</span><span class="keyword">float</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">1</span><span class="special">;</span> <span class="identifier">n</span> <span class="special"><</span> <span class="identifier">max_n</span><span class="special">;</span> <span class="identifier">n</span> <span class="special">+=</span> <span class="number">2</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">cond</span> <span class="special">+=</span> <span class="number">1</span><span class="special">/</span><span class="identifier">n</span><span class="special">;</span>
<span class="identifier">cond</span> <span class="special">-=</span> <span class="number">1</span><span class="special">/(</span><span class="identifier">n</span><span class="special">+</span><span class="number">1</span><span class="special">);</span>
<span class="special">}</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">(</span><span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">float</span><span class="special">>::</span><span class="identifier">digits10</span><span class="special">);</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"ln(2) = "</span> <span class="special"><<</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">constants</span><span class="special">::</span><span class="identifier">ln_two</span><span class="special"><</span><span class="keyword">float</span><span class="special">>()</span> <span class="special"><<</span> <span class="string">"\n"</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Sum = "</span> <span class="special"><<</span> <span class="identifier">cond</span><span class="special">.</span><span class="identifier">sum</span><span class="special">()</span> <span class="special"><<</span> <span class="string">"\n"</span><span class="special">;</span>
<span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Condition number = "</span> <span class="special"><<</span> <span class="identifier">cond</span><span class="special">()</span> <span class="special"><<</span> <span class="string">"\n"</span><span class="special">;</span>
</pre>
<p>
Output:
</p>
<pre class="programlisting"><span class="identifier">ln</span><span class="special">(</span><span class="number">2</span><span class="special">)</span> <span class="special">=</span> <span class="number">0.693147</span>
<span class="identifier">Sum</span> <span class="special">=</span> <span class="number">0.693137</span>
<span class="identifier">Condition</span> <span class="identifier">number</span> <span class="special">=</span> <span class="number">22.22282</span>
</pre>
<p>
We see that we have lost roughly two digits of accuracy, consistent with the
heuristic that if the condition number is 10<sup><span class="emphasis"><em>k</em></span></sup>, then we
lose <span class="emphasis"><em>k</em></span> significant digits in the sum.
</p>
<p>
Our guess it that if you're worried about whether your sum is ill-conditioned,
the <span class="emphasis"><em>last</em></span> 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:
</p>
<pre class="programlisting"><span class="keyword">auto</span> <span class="identifier">cond</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">summation_condition_number</span><span class="special"><</span><span class="keyword">float</span><span class="special">>();</span> <span class="comment">// will use Kahan summation.</span>
<span class="comment">// ...</span>
</pre>
<p>
Output:
</p>
<pre class="programlisting"><span class="identifier">ln</span><span class="special">(</span><span class="number">2</span><span class="special">)</span> <span class="special">=</span> <span class="number">0.693147</span>
<span class="identifier">Kahan</span> <span class="identifier">sum</span> <span class="special">=</span> <span class="number">0.693147</span>
<span class="identifier">Condition</span> <span class="identifier">number</span> <span class="special">=</span> <span class="number">22.2228</span>
</pre>
<p>
If you are interested, the L<sub>1</sub> norm is also generated by this computation, so
you may query it if you like:
</p>
<pre class="programlisting"><span class="keyword">float</span> <span class="identifier">l1</span> <span class="special">=</span> <span class="identifier">cond</span><span class="special">.</span><span class="identifier">l1_norm</span><span class="special">();</span>
<span class="comment">// l1 = 15.4</span>
</pre>
<h4>
<a name="math_toolkit.cond.h2"></a>
<span class="phrase"><a name="math_toolkit.cond.condition_number_of_function_eva"></a></span><a class="link" href="cond.html#math_toolkit.cond.condition_number_of_function_eva">Condition
Number of Function Evaluation</a>
</h4>
<p>
The <a href="https://en.wikipedia.org/wiki/Condition_number" target="_top">condition number</a>
of function evaluation is defined as the absolute value of <span class="emphasis"><em>xf</em></span>'(<span class="emphasis"><em>x</em></span>)<span class="emphasis"><em>f</em></span>(<span class="emphasis"><em>x</em></span>)<sup>-1</sup>.
It is computed as follows:
</p>
<pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">tools</span><span class="special">::</span><span class="identifier">evaluation_condition_number</span><span class="special">;</span>
<span class="keyword">auto</span> <span class="identifier">f</span> <span class="special">=</span> <span class="special">[](</span><span class="keyword">double</span> <span class="identifier">x</span><span class="special">)-></span><span class="keyword">double</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">log</span><span class="special">(</span><span class="identifier">x</span><span class="special">);</span> <span class="special">};</span>
<span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">1.125</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">cond</span> <span class="special">=</span> <span class="identifier">evaluation_condition_number</span><span class="special">(</span><span class="identifier">f</span><span class="special">,</span> <span class="number">1.125</span><span class="special">);</span>
<span class="comment">// cond = 1/log(x)</span>
</pre>
<h4>
<a name="math_toolkit.cond.h3"></a>
<span class="phrase"><a name="math_toolkit.cond.caveats"></a></span><a class="link" href="cond.html#math_toolkit.cond.caveats">Caveats</a>
</h4>
<p>
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!
</p>
<p>
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.
</p>
<h4>
<a name="math_toolkit.cond.h4"></a>
<span class="phrase"><a name="math_toolkit.cond.references"></a></span><a class="link" href="cond.html#math_toolkit.cond.references">References</a>
</h4>
<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; ">
<li class="listitem">
Gautschi, Walter. <span class="emphasis"><em>Orthogonal polynomials: computation and approximation</em></span>
Oxford University Press on Demand, 2004.
</li>
<li class="listitem">
Higham, Nicholas J. <span class="emphasis"><em>The accuracy of floating point summation.</em></span>
SIAM Journal on Scientific Computing 14.4 (1993): 783-799.
</li>
<li class="listitem">
Higham, Nicholas J. <span class="emphasis"><em>Accuracy and stability of numerical algorithms.</em></span>
Vol. 80. Siam, 2002.
</li>
</ul></div>
</div>
<div class="copyright-footer">Copyright © 2006-2021 Nikhar Agrawal, Anton Bikineev, Matthew Borland,
Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, Hubert Holin, Bruno
Lalande, John Maddock, Evan Miller, Jeremy Murphy, Matthew Pulver, Johan Råde,
Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, Daryle
Walker and Xiaogang Zhang<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
</p>
</div>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="float_comparison.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../utils.html"><img src="../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="ulps_plots.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>
|
