Given numbers $x, y \in \mathbb{R}$ where $$\frac{|y-x|}{|x|}$$ is small, and code that implements the function $f$ with a sequence of arithmetic operations, I would like to compute to high accuracy the finite difference $$ \frac{f(y) - f(x)}{y-x}, $$ which is plagued by cancellation. Is there a variant of automatic differentiation that lets me do it (without working in higher precision, possibly)? At a first sight it looks like the basic ideas of AD transfer to this similar problem: start from the pair $(x, y-x)$, and apply the same sequence of operations to it in "forward mode". But I can't get anything useful with a search for "automatic finite differences".
$\begingroup$
$\endgroup$
3
asked Jun 9, 2021 at 6:13
-
3$\begingroup$ Probably also have a look at the complex step derivative (which basically is automatic differentiation in disguise). $\endgroup$– davidhighJun 10, 2021 at 6:24
-
$\begingroup$ @davidhigh thanks, but I don't really see how that would help in my case. How is $f(x+iz)$ related to an incremental ratio? Could you expand, if you have something in mind? $\endgroup$– Federico PoloniJun 10, 2021 at 6:32
-
$\begingroup$ @FredericoPoloni: nope, nothing in mind. It was just a reflex based on the words AD and Finite Difference. But as you stated in the comment section of the other answer, your question is a different setup. $\endgroup$– davidhighJun 10, 2021 at 21:44
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
1
This is what Griewank et al. call "Piecewise linearization in secant mode", see for instance https://opus4.kobv.de/opus4-zib/files/6164/newton_secant_approx_paper.pdf. The aim of that research was to capture the kinks of absolute value operations with the same precision a tangent or a secant captures the local behavior of a smooth function, with an implementation as an extension of the forward mode of automatic/algorithmic differentiation (Adol-C and home-brew python classes). But as part of that it of course also computes secants of smooth functions.
The secant slopes are computed using the available tricks for the elementary operations, for instance if $v(x)=\sin(u(x))$, then $$ S_v=\frac{v(y)-v(x)}{y-x}=2\frac{\sin\frac{u(y)-u(x)}{2}}{y-x}·\cos\frac{u(y)+u(x)}2 \\ =2\cos\frac{u(y)+u(x)}2·\frac{\sin\left(\frac{S_u}2·(y-x)\right)}{y-x} $$ where $S_u$, $S_v$ are the slopes of the secants of $u$ and $v$ for the given pair of points.
The numerical evaluation of the sine close to zero is usually implemented in a way that the quotient $\frac{\sin\epsilon}{\epsilon}$ is correctly evaluated. If one is unsure about that, one could of course approximate that quotient for values below some threshold using the Taylor series $1-\frac16\epsilon^2+O(\epsilon^4)$.
This basic procedure for elementary operations gets extended to the full function by propagating it along the computational graph, the same way as the forward mode of AD. And as there, variants exist in the implementation philosophy, such as reading out the function into a tape and using a "tape machine" for the secant computation, or transforming the computational graph/AST directly into a tree data structure where the secant computation proceeds in node operations, or by propagating a secant "number" type through the original (overloaded) function without reading out the function structure.
answered Jun 9, 2021 at 7:00
$\begingroup$
$\endgroup$
Potentially related and useful (I gave these resources to my students when teaching Intro to Computational Mathematics, kinda useful pedagogically too):
"Automatic Source-to-Source Error Compensation of Floating-Point Programs" by Laurent Thévenoux, Philippe Langlois and Matthieu Martel : https://hal.archives-ouvertes.fr/hal-01158399/document
Herbie: Automatically Improving Floating Point Accuracy (https://herbie.uwplse.org/) by Programming Languages and Software Engineering group of the University of Washington (https://uwplse.org/): the main publication is "Automatically Improving Accuracy for Floating Point Expressions" by Pavel Panchekha, Alex Sachez-Stern, James R. Wilcox and Zachary Tatlock which can be found at https://herbie.uwplse.org/pldi15-paper.pdf
answered Jun 10, 2021 at 7:08
$\begingroup$
$\endgroup$
5
Perhaps what @davidhigh has in mind is the following
$$ \frac{f(y)-f(x)}{y-x} = \frac{f(x+h) - f(x)}{h} = f'(x) \,\, \text{ as } h\rightarrow{0}\,\, \text{substituting} \,\, y-x=h $$
And you can get $f'(x)$ by complex step differentiation.
-
$\begingroup$ Thanks, but I am still missing something; in this question I am not interested in $f'(x)$: I want to compute the incremental ratio exactly given $x$ and $y$ (close to each other but distinct), not its limit. This is a different setup than the classical one of algorithmic differentiation. $\endgroup$ Jun 10, 2021 at 9:54
-
$\begingroup$ Well, it seems that I am missing something. Wouldn't your ratio be equal to the derivative when $y-x$ is small? $\endgroup$– NNNJun 10, 2021 at 9:56
-
$\begingroup$ No; its limit is equal to the derivative, but the ratio itself is not. $\endgroup$ Jun 10, 2021 at 9:57
-
5$\begingroup$ Take $f(x) = \sin(x)$ for example. Its derivative is $\cos(x)$, but Herbie (mentioned in my answer) suggests a much different function to implement $(\sin(y)-\sin(x))/(y-x)$ if $|h| = |y-x|<10^{-3}$: $\left(\cos x - h \cdot \left(\sin x \cdot 0.5 + h \cdot \left(\cos x \cdot 0.16666666666666666\right)\right)\right) + 0.041666666666666664 \cdot \left(\sin x \cdot {h}^{3}\right)$. This function reduces the average error from 46.2 bits to 0 bits in the vicinity of zero. (Also if $h=0$ you get $\cos(x)$ back which is nice too). $\endgroup$ Jun 10, 2021 at 14:52
-
$\begingroup$ Thanks everyone. For some reason, this reminds me of 'modified equations'. ocw.mit.edu/courses/mathematics/… $\endgroup$– NNNJun 11, 2021 at 6:45