Automatic finite differences

Question

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

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

Answer 2

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

Answer 3

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.

Answer 4

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.