1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Lanczos Smoothing Derivatives</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="../quadrature.html" title="Chapter 13. Quadrature and Differentiation"> 9<link rel="prev" href="autodiff.html" title="Automatic Differentiation"> 10<link rel="next" href="../filters.html" title="Chapter 14. Filters"> 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="autodiff.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="../filters.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.diff0"></a><a class="link" href="diff0.html" title="Lanczos Smoothing Derivatives">Lanczos Smoothing Derivatives</a> 28</h2></div></div></div> 29<h4> 30<a name="math_toolkit.diff0.h0"></a> 31 <span class="phrase"><a name="math_toolkit.diff0.synopsis"></a></span><a class="link" href="diff0.html#math_toolkit.diff0.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">differentiation</span><span class="special">/</span><span class="identifier">lanczos_smoothing</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span> 34 35<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">differentiation</span> <span class="special">{</span> 36 37 <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="identifier">size_t</span> <span class="identifier">order</span><span class="special">=</span><span class="number">1</span><span class="special">></span> 38 <span class="keyword">class</span> <span class="identifier">discrete_lanczos_derivative</span> <span class="special">{</span> 39 <span class="keyword">public</span><span class="special">:</span> 40 <span class="identifier">discrete_lanczos_derivative</span><span class="special">(</span><span class="identifier">Real</span> <span class="identifier">spacing</span><span class="special">,</span> 41 <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">18</span><span class="special">,</span> 42 <span class="identifier">size_t</span> <span class="identifier">approximation_order</span> <span class="special">=</span> <span class="number">3</span><span class="special">);</span> 43 44 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">RandomAccessContainer</span><span class="special">></span> 45 <span class="identifier">Real</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">RandomAccessContainer</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 46 47 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">RandomAccessContainer</span><span class="special">></span> 48 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">RandomAccessContainer</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">v</span><span class="special">,</span> <span class="identifier">RandomAccessContainer</span> <span class="special">&</span> <span class="identifier">dvdt</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 49 50 <span class="keyword">template</span><span class="special"><</span><span class="keyword">class</span> <span class="identifier">RandomAccessContainer</span><span class="special">></span> 51 <span class="identifier">RandomAccessContainer</span> <span class="keyword">operator</span><span class="special">()(</span><span class="identifier">RandomAccessContainer</span> <span class="keyword">const</span> <span class="special">&</span> <span class="identifier">v</span><span class="special">)</span> <span class="keyword">const</span><span class="special">;</span> 52 53 <span class="identifier">Real</span> <span class="identifier">get_spacing</span><span class="special">()</span> <span class="keyword">const</span><span class="special">;</span> 54 <span class="special">};</span> 55 56<span class="special">}</span> <span class="comment">// namespaces</span> 57</pre> 58<h4> 59<a name="math_toolkit.diff0.h1"></a> 60 <span class="phrase"><a name="math_toolkit.diff0.description"></a></span><a class="link" href="diff0.html#math_toolkit.diff0.description">Description</a> 61 </h4> 62<p> 63 The <code class="computeroutput"><span class="identifier">discrete_lanczos_derivative</span></code> 64 class calculates a finite-difference approximation to the derivative of a noisy 65 sequence of equally-spaced values <span class="emphasis"><em>v</em></span>. A basic usage is 66 </p> 67<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 68<span class="comment">// fill v with noisy data.</span> 69<span class="keyword">double</span> <span class="identifier">spacing</span> <span class="special">=</span> <span class="number">0.001</span><span class="special">;</span> 70<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">differentiation</span><span class="special">::</span><span class="identifier">discrete_lanczos_derivative</span><span class="special">;</span> 71<span class="keyword">auto</span> <span class="identifier">lanczos</span> <span class="special">=</span> <span class="identifier">discrete_lanczos_derivative</span><span class="special">(</span><span class="identifier">spacing</span><span class="special">);</span> 72<span class="comment">// Compute dvdt at index 30:</span> 73<span class="keyword">double</span> <span class="identifier">dvdt30</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="number">30</span><span class="special">);</span> 74<span class="comment">// Compute derivative of entire vector:</span> 75<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">dvdt</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">);</span> 76</pre> 77<p> 78 Noise-suppressing second derivatives can also be computed. The syntax is as 79 follows: 80 </p> 81<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 82<span class="comment">// fill v with noisy data.</span> 83<span class="keyword">auto</span> <span class="identifier">lanczos</span> <span class="special">=</span> <span class="identifier">lanczos_derivative</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="number">2</span><span class="special">>(</span><span class="identifier">spacing</span><span class="special">);</span> 84<span class="comment">// evaluate second derivative at a point:</span> 85<span class="keyword">double</span> <span class="identifier">d2vdt2</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="number">18</span><span class="special">);</span> 86<span class="comment">// evaluate second derivative of entire vector:</span> 87<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">d2vdt2</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">);</span> 88</pre> 89<p> 90 For memory conscious programmers, you can reuse the memory space for the derivative 91 as follows: 92 </p> 93<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 94<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">dvdt</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 95<span class="comment">// . . . define spacing, create and instance of discrete_lanczos_derivative</span> 96 97<span class="comment">// populate dvdt, perhaps in a loop:</span> 98<span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">dvdt</span><span class="special">);</span> 99</pre> 100<p> 101 If the data has variance σ<sup>2</sup>, then the variance of the computed derivative 102 is roughly σ<sup>2</sup><span class="emphasis"><em>p</em></span><sup>3</sup> <span class="emphasis"><em>n</em></span><sup>-3</sup> Δ 103 <span class="emphasis"><em>t</em></span><sup>-2</sup>, i.e., it increases cubically with the approximation 104 order <span class="emphasis"><em>p</em></span>, linearly with the data variance, and decreases 105 at the cube of the filter length <span class="emphasis"><em>n</em></span>. In addition, we must 106 not forget the discretization error which is <span class="emphasis"><em>O</em></span>(Δ 107 <span class="emphasis"><em>t</em></span><sup><span class="emphasis"><em>p</em></span></sup>). You can play around with the 108 approximation order <span class="emphasis"><em>p</em></span> and the filter length <span class="emphasis"><em>n</em></span>: 109 </p> 110<pre class="programlisting"><span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">12</span><span class="special">;</span> 111<span class="identifier">size_t</span> <span class="identifier">p</span> <span class="special">=</span> <span class="number">2</span><span class="special">;</span> 112<span class="keyword">auto</span> <span class="identifier">lanczos</span> <span class="special">=</span> <span class="identifier">lanczos_derivative</span><span class="special">(</span><span class="identifier">spacing</span><span class="special">,</span> <span class="identifier">n</span><span class="special">,</span> <span class="identifier">p</span><span class="special">);</span> 113<span class="keyword">double</span> <span class="identifier">dvdt</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v</span><span class="special">,</span> <span class="identifier">i</span><span class="special">);</span> 114</pre> 115<p> 116 If <span class="emphasis"><em>p=2n</em></span>, then the discrete Lanczos derivative is not smoothing: 117 It reduces to the standard <span class="emphasis"><em>2n+1</em></span>-point finite-difference 118 formula. For <span class="emphasis"><em>p>2n</em></span>, an assertion is hit as the filter 119 is undefined. 120 </p> 121<p> 122 In our tests with AWGN, we have found the error decreases monotonically with 123 <span class="emphasis"><em>n</em></span>, as is expected from the theory discussed above. So 124 the choice of <span class="emphasis"><em>n</em></span> is simple: As high as possible given your 125 speed requirements (larger <span class="emphasis"><em>n</em></span> implies a longer filter and 126 hence more compute), balanced against the danger of overfitting and averaging 127 over non-stationarity. 128 </p> 129<p> 130 The choice of approximation order <span class="emphasis"><em>p</em></span> for a given <span class="emphasis"><em>n</em></span> 131 is more difficult. If your signal is believed to be a polynomial, it does not 132 make sense to set <span class="emphasis"><em>p</em></span> to larger than the polynomial degree- 133 though it may be sensible to take <span class="emphasis"><em>p</em></span> less than this. 134 </p> 135<p> 136 For a sinusoidal signal contaminated with AWGN, we ran a few tests showing 137 that for SNR = 1, p = n/8 gave the best results, for SNR = 10, p = n/7 was 138 the best, and for SNR = 100, p = n/6 was the most reasonable choice. For SNR 139 = 0.1, the method appears to be useless. The user is urged to use these results 140 with caution: they have no theoretical backing and are extrapolated from a 141 single case. 142 </p> 143<p> 144 The filters are (regrettably) computed at runtime-the vast number of combinations 145 of approximation order and filter length makes the number of filters that must 146 be stored excessive for compile-time data. The constructor call computes the 147 filters. Since each filter has length <span class="emphasis"><em>2n+1</em></span> and there are 148 <span class="emphasis"><em>n</em></span> filters, whose element each consist of <span class="emphasis"><em>p</em></span> 149 summands, the complexity of the constructor call is O(<span class="emphasis"><em>n</em></span><sup>2</sup><span class="emphasis"><em>p</em></span>). 150 This is not cheap-though for most cases small <span class="emphasis"><em>p</em></span> and <span class="emphasis"><em>n</em></span> 151 not too large (< 20) is desired. However, for concreteness, on the author's 152 2.7GHz Intel laptop CPU, the <span class="emphasis"><em>n=16</em></span>, <span class="emphasis"><em>p=3</em></span> 153 filter takes 9 microseconds to compute. This is far from negligible, and as 154 such the filters may be used with multiple data, and even shared between threads: 155 </p> 156<pre class="programlisting"><span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v1</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 157<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">v2</span><span class="special">(</span><span class="number">500</span><span class="special">);</span> 158<span class="comment">// fill v1, v2 with noisy data.</span> 159<span class="keyword">auto</span> <span class="identifier">lanczos</span> <span class="special">=</span> <span class="identifier">lanczos_derivative</span><span class="special">(</span><span class="identifier">spacing</span><span class="special">);</span> 160<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">dv1dt</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v1</span><span class="special">);</span> <span class="comment">// threadsafe</span> 161<span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">dv2dt</span> <span class="special">=</span> <span class="identifier">lanczos</span><span class="special">(</span><span class="identifier">v2</span><span class="special">);</span> <span class="comment">// threadsafe</span> 162<span class="comment">// need to use a different spacing?</span> 163<span class="identifier">lanczos</span><span class="special">.</span><span class="identifier">reset_spacing</span><span class="special">(</span><span class="number">0.02</span><span class="special">);</span> <span class="comment">// not threadsafe</span> 164</pre> 165<p> 166 The implementation follows <a href="https://doi.org/10.1080/00207160.2012.666348" target="_top">McDevitt, 167 2012</a>, who vastly expanded the ideas of Lanczos to create a very general 168 framework for numerically differentiating noisy equispaced data. 169 </p> 170<h4> 171<a name="math_toolkit.diff0.h2"></a> 172 <span class="phrase"><a name="math_toolkit.diff0.example"></a></span><a class="link" href="diff0.html#math_toolkit.diff0.example">Example</a> 173 </h4> 174<p> 175 We have extracted some data from the <a href="https://www.gw-openscience.org/data/" target="_top">LIGO 176 signal</a> and differentiated it using the (<span class="emphasis"><em>n</em></span>, <span class="emphasis"><em>p</em></span>) 177 = (60, 4) Lanczos smoothing derivative, as well as using the (<span class="emphasis"><em>n</em></span>, 178 <span class="emphasis"><em>p</em></span>) = (4, 8) (nonsmoothing) derivative. 179 </p> 180<div class="blockquote"><blockquote class="blockquote"><p> 181 <span class="inlinemediaobject"><img src="../../graphs/ligo_derivative.svg" align="middle"></span> 182 183 </p></blockquote></div> 184<p> 185 The original data is in orange, the smoothing derivative in blue, and the non-smoothing 186 standard finite difference formula is in gray. (Each time series has been rescaled 187 to fit in the same graph.) We can see that the smoothing derivative tracks 188 the increase and decrease in the trend well, whereas the standard finite difference 189 formula produces nonsense and amplifies noise. 190 </p> 191<h4> 192<a name="math_toolkit.diff0.h3"></a> 193 <span class="phrase"><a name="math_toolkit.diff0.caveats"></a></span><a class="link" href="diff0.html#math_toolkit.diff0.caveats">Caveats</a> 194 </h4> 195<p> 196 The computation of the filters is ill-conditioned for large <span class="emphasis"><em>p</em></span>. 197 In order to mitigate this problem, we have computed the filters in higher precision 198 and cast the results to the desired type. However, this simply pushes the problem 199 to larger <span class="emphasis"><em>p</em></span>. In practice, this is not a problem, as large 200 <span class="emphasis"><em>p</em></span> corresponds to less powerful denoising, but keep it 201 in mind. 202 </p> 203<p> 204 In addition, the <code class="computeroutput"><span class="special">-</span><span class="identifier">ffast</span><span class="special">-</span><span class="identifier">math</span></code> flag 205 has a very large effect on the speed of these functions. In our benchmarks, 206 they were 50% faster with the flag enabled, which is much larger than the usual 207 performance increases we see by turning on this flag. Hence, if the default 208 performance is not sufficient, this flag is available, though it comes with 209 its own problems. 210 </p> 211<p> 212 This requires C++17 <code class="computeroutput"><span class="keyword">if</span> <span class="keyword">constexpr</span></code>. 213 </p> 214<h4> 215<a name="math_toolkit.diff0.h4"></a> 216 <span class="phrase"><a name="math_toolkit.diff0.references"></a></span><a class="link" href="diff0.html#math_toolkit.diff0.references">References</a> 217 </h4> 218<div class="itemizedlist"><ul class="itemizedlist" style="list-style-type: disc; "> 219<li class="listitem"> 220 Corless, Robert M., and Nicolas Fillion. <span class="emphasis"><em>A graduate introduction 221 to numerical methods.</em></span> AMC 10 (2013): 12. 222 </li> 223<li class="listitem"> 224 Lanczos, Cornelius. <span class="emphasis"><em>Applied analysis.</em></span> Courier Corporation, 225 1988. 226 </li> 227<li class="listitem"> 228 Timothy J. McDevitt (2012): <span class="emphasis"><em>Discrete Lanczos derivatives of noisy 229 data</em></span>, International Journal of Computer Mathematics, 89:7, 916-931 230 </li> 231</ul></div> 232</div> 233<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 234<td align="left"></td> 235<td align="right"><div class="copyright-footer">Copyright © 2006-2019 Nikhar 236 Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos, 237 Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan 238 Råde, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg, 239 Daryle Walker and Xiaogang Zhang<p> 240 Distributed under the Boost Software License, Version 1.0. (See accompanying 241 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>) 242 </p> 243</div></td> 244</tr></table> 245<hr> 246<div class="spirit-nav"> 247<a accesskey="p" href="autodiff.html"><img src="../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../quadrature.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="../filters.html"><img src="../../../../../doc/src/images/next.png" alt="Next"></a> 248</div> 249</body> 250</html> 251
