# Automatic finite differences

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".

• Probably also have a look at the complex step derivative (which basically is automatic differentiation in disguise). Commented Jun 10, 2021 at 6:24
• @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? Commented Jun 10, 2021 at 6:32
• @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. Commented Jun 10, 2021 at 21:44

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.

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

There is an analogue of automatic differentiation for divided differences. Here is a minimal sample implementation in C++ that computes the divided difference of $$f(x)=x^3+x^2$$ but which can be extended in a number of directions. Implementations will look similar in Haskell, Python and other languages:

#include <iostream>

template<typename T>
struct FD
{
T a, b, c;

FD operator+(const FD& other)
{
return FD{a + other.a, b + other.b, c + other.c};
}

FD operator*(const FD& other)
{
return FD{a * other.a,
b * other.b,
c * other.b + a * other.c};
}
};

template<typename T>
float FiniteDifference(T (*)(T), T a, T b)
{
return f(FD{a, b, 1.f}).c;
}

template<typename T>
T f(T x)
{
return x * x * x + x * x;
}

int main()
{
std::cout << FiniteDifference(f, 3.f, 8.f) << std::endl;
}


There are a number of ways of looking at this:

(1) if we define $$M=\pmatrix{a&1\\0&b}$$ we can extend analytic functions $$f$$ to act on matrices and so compute $$f(M)$$. The top right entry is the divided difference. The type FD represents the matrix $$\pmatrix{a&c\\0&b}$$ and gives an efficient way to compute products and sums. Note that many methods for optimizing computations over reals carry over to matrices. Eg. we can compute finite differences of $$f(x)=x^n$$ by repeated squaring and the same method can be applied here. For related methods see Kahan & Fateman.

(2) Another point of view is that the variable c is using a modified chain rule to track the differences as functions like + and * are applied. The rule for * is a modified form of the Leibniz rule.

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.

• 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. Commented Jun 10, 2021 at 9:54
• Well, it seems that I am missing something. Wouldn't your ratio be equal to the derivative when $y-x$ is small?
– NNN
Commented Jun 10, 2021 at 9:56
• No; its limit is equal to the derivative, but the ratio itself is not. Commented Jun 10, 2021 at 9:57
• 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). Commented Jun 10, 2021 at 14:52
• Thanks everyone. For some reason, this reminds me of 'modified equations'. ocw.mit.edu/courses/mathematics/…
– NNN
Commented Jun 11, 2021 at 6:45