1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Condition Numbers</title> 5<link rel="stylesheet" href="../math.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../index.html" title="Math Toolkit 2.12.0"> 8<link rel="up" href="../utils.html" title="Chapter 2. Floating Point Utilities"> 9<link rel="prev" href="float_comparison.html" title="Floating-point Comparison"> 10<link rel="next" href="../cstdfloat.html" title="Chapter 3. Specified-width floating-point typedefs"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../boost.png"></td> 15<td align="center"><a href="../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<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="../cstdfloat.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h2 class="title" style="clear: both"> 27<a name="math_toolkit.cond"></a><a class="link" href="cond.html" title="Condition Numbers">Condition Numbers</a> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.cond.h0"></a> 31 <span class="phrase"><a name="math_toolkit.cond.synopsis"></a></span><a class="link" href="cond.html#math_toolkit.cond.synopsis">Synopsis</a> 32 </h4> 33<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> 34 35 36<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> 37 38<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> 39<span class="keyword">class</span> <span class="identifier">summation_condition_number</span> <span class="special">{</span> 40<span class="keyword">public</span><span class="special">:</span> 41 <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> 42 43 <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> 44 45 <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> 46 47 <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> 48 49 <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> 50 51 <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> 52<span class="special">};</span> 53 54<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> 55<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> 56 57<span class="special">}</span> <span class="comment">// namespaces</span> 58</pre> 59<h4> 60<a name="math_toolkit.cond.h1"></a> 61 <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 62 Condition Number</a> 63 </h4> 64<p> 65 Here we compute the condition number of the alternating harmonic sum: 66 </p> 67<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> 68<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> 69<span class="keyword">float</span> <span class="identifier">max_n</span> <span class="special">=</span> <span class="number">10000000</span><span class="special">;</span> 70<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> 71<span class="special">{</span> 72 <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> 73 <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> 74<span class="special">}</span> 75<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> 76<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> 77<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> 78<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> 79</pre> 80<p> 81 Output: 82 </p> 83<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> 84<span class="identifier">Sum</span> <span class="special">=</span> <span class="number">0.693137</span> 85<span class="identifier">Condition</span> <span class="identifier">number</span> <span class="special">=</span> <span class="number">22.22282</span> 86</pre> 87<p> 88 We see that we have lost roughly two digits of accuracy, consistent with the 89 heuristic that if the condition number is 10<sup><span class="emphasis"><em>k</em></span></sup>, then we 90 lose <span class="emphasis"><em>k</em></span> significant digits in the sum. 91 </p> 92<p> 93 Our guess it that if you're worried about whether your sum is ill-conditioned, 94 the <span class="emphasis"><em>last</em></span> thing you want is for your condition number estimate 95 to be inaccurate as well. Since the condition number estimate relies on computing 96 the (perhaps ill-conditioned) sum, we have defaulted the accumulation to use 97 Kahan summation: 98 </p> 99<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> 100<span class="comment">// ...</span> 101</pre> 102<p> 103 Output: 104 </p> 105<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> 106<span class="identifier">Kahan</span> <span class="identifier">sum</span> <span class="special">=</span> <span class="number">0.693147</span> 107<span class="identifier">Condition</span> <span class="identifier">number</span> <span class="special">=</span> <span class="number">22.2228</span> 108</pre> 109<p> 110 If you are interested, the L<sub>1</sub> norm is also generated by this computation, so 111 you may query it if you like: 112 </p> 113<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> 114<span class="comment">// l1 = 15.4</span> 115</pre> 116<h4> 117<a name="math_toolkit.cond.h2"></a> 118 <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 119 Number of Function Evaluation</a> 120 </h4> 121<p> 122 The <a href="https://en.wikipedia.org/wiki/Condition_number" target="_top">condition number</a> 123 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>. 124 It is computed as follows: 125 </p> 126<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> 127<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> 128<span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">1.125</span><span class="special">;</span> 129<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> 130<span class="comment">// cond = 1/log(x)</span> 131</pre> 132<h4> 133<a name="math_toolkit.cond.h3"></a> 134 <span class="phrase"><a name="math_toolkit.cond.caveats"></a></span><a class="link" href="cond.html#math_toolkit.cond.caveats">Caveats</a> 135 </h4> 136<p> 137 The condition number of function evaluation gives us a measure of how sensitive 138 our function is to roundoff error. Unfortunately, evaluating the condition 139 number requires evaluating the function and its derivative, and this calculation 140 is itself inaccurate whenever the condition number of function evaluation is 141 large. Sadly, this is also the regime when you are most interested in evaluating 142 a condition number! 143 </p> 144<p> 145 This seems to be a fundamental problem. However, it should not be necessary 146 to evaluate the condition number to high accuracy, valuable insights can be 147 obtained simply by looking at the change in condition number as the function 148 evolves over its domain. 149 </p> 150<h4> 151<a name="math_toolkit.cond.h4"></a> 152 <span class="phrase"><a name="math_toolkit.cond.references"></a></span><a class="link" href="cond.html#math_toolkit.cond.references">References</a> 153 </h4> 154<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 155<li class="listitem"> 156 Gautschi, Walter. <span class="emphasis"><em>Orthogonal polynomials: computation and approximation</em></span> 157 Oxford University Press on Demand, 2004. 158 </li> 159<li class="listitem"> 160 Higham, Nicholas J. <span class="emphasis"><em>The accuracy of floating point summation.</em></span> 161 SIAM Journal on Scientific Computing 14.4 (1993): 783-799. 162 </li> 163<li class="listitem"> 164 Higham, Nicholas J. <span class="emphasis"><em>Accuracy and stability of numerical algorithms.</em></span> 165 Vol. 80. Siam, 2002. 166 </li> 167</ul></div> 168</div> 169<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 170<td align="left"></td> 171<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 172 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 173 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 174 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 175 Daryle Walker and Xiaogang Zhang<p> 176 Distributed under the Boost Software License, Version 1.0. (See accompanying 177 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>) 178 </p> 179</div></td> 180</tr></table> 181<hr> 182<div class="spirit-nav"> 183<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="../cstdfloat.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 184</div> 185</body> 186</html> 187
