Convergence of AAA rational approximations

using RationalFunctionApproximation, CairoMakie

For a function that is analytic on its domain, the AAA algorithm typically converges at least root-exponentially. In order to observe the convergence, we can construct an approximation that preserves the history of its construction. For example, we approximate an oscillatory function over $[-1,1]$ via:

f = x -> cos(11x)
r = approximate(f, unit_interval, stats=true)
convergenceplot(r)

In the plot above, the markers show the estimated max-norm error of the $n$-point AAA rational interpolant over the domain as a function of the numerator/denominator degree $n-1$. The red circles indicate that a pole of the rational interpolant lies on the interval itself. We can verify this by using rewind to recover the degree-7 approximation that was found along the way to r:

julia> r7 = rewind(r, 7)Barycentric rational function of type (7,7) on the domain: Path with 1 curve
julia> poles(r7)7-element Vector{ComplexF64}:
 -0.7865385850382732 - 0.26769746910591097im
 -0.7865385850382732 + 0.2676974691059109im
 -0.3174146424804872 + 0.28894380984739326im
 -0.3174146424804871 - 0.2889438098473932im
 0.16129605360549806 + 0.2496796849023167im
  0.1612960536054981 - 0.24967968490231668im
  0.7365067113018727 + 0.0im

When the AAA iteration encounters an approximation with such undesired poles, or having less accuracy than a predecessor, the AAA iteration simply disregards that approximation and continues–-unless there have been more than a designated number of consecutive failures, at which the best interpolant ever encountered is returned. That interpolant is indicated by the gold halo in the convergence plot above.

When a singularity is very close to the approximation domain, it can cause stagnation and a large number of bad-pole failures:

f = x -> tanh(3000*(x - 1/5))
r = approximate(f, unit_interval, stats=true)
convergenceplot(r)

This effect is thought to be mainly due to roundoff and conditioning of the problem. If we use more accurate floating-point arithmetic, we can see that the AAA convergence continues steadily past the previous plateau:

using DoubleFloats
r = approximate(f, unit_interval, float_type=Double64, stats=true)
convergenceplot(r)

In the extreme case of a function with a singularity on the domain, the convergence can be substantially affected:

f = x -> abs(x - 1/8)
r = approximate(f, unit_interval, stats=true)
convergenceplot(r)

In such a case, we might get improvement by increasing the number of allowed consecutive failures via the lookahead keyword argument:

r = approximate(f, unit_interval, stats=true, lookahead=20)
convergenceplot(r)