How to find simple pole expansion that maintains high acuracy?

[kbarros]

The script below builds a AAA approximation with very high accuracy. Great!

using LinearAlgebra
using RationalFunctionApproximation
import DoubleFloats: Double64
import ComplexRegions: Segment

f(x) = 1 / (1 + exp(1_000 * x))
bounds = (-10, 10)
segment = Segment{Double64}(bounds...)
# Large refinement and stagnation seem to be required. Are there heuristics
# to pick them or to otherwise mark the near-singular point at x=0?
r = approximate(f, segment, refinement=100, stagnation=100, tol=1e-10)

xs = range(bounds..., length=10_000)
norm(r.(xs) - f.(xs), Inf) # 5.9e-11

My ultimate goal is to find poles/residues for a simple pole expansion that maintains this high level of accuracy. Unfortunately, when I try to get that with residues(r) plus an inferred constant shift, there is significant degradation in accuracy (several orders of magnitude).

With the hope of recovering that lost accuracy, I tried building a continued fraction approximation using fixed poles. When f(x) is not too singular, the procedure can work well. However, in this example:

r0 = approximate(f, segment, poles(r), degree=0)
norm(r0.(xs) - f.(xs), Inf) # 1.2e4

r1 = approximate(f, segment, poles(r), degree=0, init=1_000)
norm(r1.(xs) - f.(xs), Inf) # 3.1e0

r2 = approximate(f, segment, poles(r), degree=0, init=10_000)
norm(r2.(xs) - f.(xs), Inf) # 1.0e-6

# ERROR: BoundsError: attempt to access 20000×4 Matrix{Double64} at index [20001, 1]
r3 = approximate(f, segment, poles(r), degree=0, init=50_000)

Here the error depends sensitively on the choice of init. The best I could achieve is order 1e-6, which is far worse than the 1e-10 achieved by AAA. Pushing the init value past a threshold generates a BoundsError.

Is there a way to get a simple pole expansion that keeps the high accuracy of the original barycentric approximation?

Another curious phenomenon is that if I rescale x as it appears in both f and bounds, I can get an overall much better approximation. Specifically, the alternative definitions,

(x) = 1 / (1 + exp(10_000 * x)) # Previously exp(1_000 * x)
bounds = (-1, 1) # Previously (-10, 10)

yield a significantly improved accuracy. Naively, I would have expected no mathematical difference under this rescaling.

Thanks in advance. This is a very nice package!

[kbarros]

I isolated the numerical stability problem to the "catch" fallback in the poles calculation:

RationalFunctionApproximation.jl/src/barycentric.jl

Lines 211 to 217 in f2878f5

	try
	pol = filter( isfinite, eigvals(E, B) )
	catch
	# generalized eigen not available in extended precision, so:
	λ = filter( z->abs(z)>1e-13, eigvals(E\B) )
	pol = 1 ./ λ
	end

Instead of eigvals(E\B), this is much more stable:

catch
    (; α, β) = schur(complex(E), complex(B))
    pol = filter( isfinite, α ./ β )
end

With this change, I'm able to extract good residues(r) from the original barycentric approximation at Double64 precision. Here is the test case:

using LinearAlgebra
using RationalFunctionApproximation
import DoubleFloats: Double64
import ComplexRegions: Segment

f(x) = 1 / (1 + exp(1_000 * x))
bounds = (-10, 10)
segment = Segment{Double64}(bounds...)
r = approximate(f, segment, refinement=100, stagnation=100, tol=1e-10)
(ζs, γs) = residues(r)

# Find constant shift empirically (is there a better way?)
x_mid = (bounds[1] + bounds[2]) / 2
c = real(r(x_mid) - sum(γ / (x_mid - ζ) for (ζ, γ) in zip(ζs, γs)))

simple_pole_approx(x) = c + sum(γ / (x - ζ) for (ζ, γ) in zip(ζs, γs))

xs = range(bounds..., length=10_000)

# Ground truth: Error for original barycentric model
norm(r.(xs) - f.(xs), Inf) # 5.9e-11

# Accuracy of simple pole expansion. Ideally, this should match the above.
norm(simple_pole_approx.(xs) - f.(xs), Inf) # With Schur: 5.9e-11, With eigvals(E\B): 2.66e-1

[tobydriscoll]

That's a great catch, thanks. There is still a deeper issue that the eigenvalue problem is not scale-invariant for nodes and values. Pending any mathematical insight to that problem, I may add one or two Newton iterations to the roots and poles methods to improve the accuracy in such cases.

[kbarros]

That’s great! Looking forward to the next release.

Regarding the broken scale invariance — are there any hidden parameters in the algorithm that are implicitly assuming that the size of the domain is order 1? If so, they could be rescaled according to some characteristic size of the domain.