Source-linked AI summary
The AAA algorithm for rational approximation
Yuji Nakatsukasa, Olivier Sète, Lloyd N. Trefethen
TL;DR
Rational approximation faces numerical difficulties including spurious poles and the ill-posedness of analytic continuation. The AAA algorithm combines barycentric rational representations with greedily selected support points, providing a fast and flexible approach across varied point sets; in a zeta-function example, its pole and zero matched known values in all but the last two digits.
Problem
Spurious poles, numerical Froissart doublets, and the ill-posedness of analytic continuation make rational approximation difficult.
Method
AAA uses barycentric rational representations with interpolation at selected support points and grows degree greedily to avoid exponential instabilities.
Results
The zeta-function approximant’s pole and zero matched known values in each case in all digits but the last two.
Takeaways & Limitations
AAA provides a fast and flexible way to compute rational functions on point sets including disconnected, irregular, and possibly unbounded domains.
Takeaways & Limitations
For approximating |x|, odd iteration steps can produce a real pole near x = 0, and the authors offer no solution.
Abstract
from arXiv · showhide
We introduce a new algorithm for approximation by rational functions on a real or complex set of points, implementable in 40 lines of Matlab and requiring no user input parameters. Even on a disk or interval the algorithm may outperform existing methods, and on more complicated domains it is especially competitive. The core ideas are (1) representation of the rational approximant in barycentric form with interpolation at certain support points and (2) greedy selection of the support points to avoid exponential instabilities. The name AAA stands for "adaptive Antoulas--Anderson" in honor of the authors who introduced a scheme based on (1). We present the core algorithm with a Matlab code and nine applications and describe variants targeted at problems of different kinds. Comparisons are made with vector fitting, RKFIT, and other existing methods for rational approximation.
1. Introduction.
Rational approximation is useful for compact representation and extrapolation, but spurious poles and the ill-posedness of analytic continuation make computation difficult. AAA addresses these challenges with barycentric representations and greedy support-point selection, while remaining flexible across domains.
- Spurious poles, or Froissart doublets, are pole-zero pairs that nearly cancel or poles with very small residues.
- Analytic continuation is ill-posed, and floating-point arithmetic makes numerical Froissart doublets more common.
- AAA combines barycentric rational representations with greedy support-point selection to avoid exponential instabilities.
- AAA is not tied to a particular approximation domain and can handle disconnected, irregular, possibly unbounded point sets, including samples containing poles.
- AAA does not claim optimality in a particular norm, for which specialized methods may be preferable.
2. Rational barycentric representations.
The barycentric representation expresses a rational interpolant as a quotient of partial fractions built from support points, data values, and weights. Its apparent complexity cancels to type (m −1, m −1), while its basis can provide strong numerical conditioning.
- The representation uses distinct real or complex support points zj, data values fj, and nonzero weights wj in numerator and denominator partial fractions.
- Although the numerator and denominator partial fractions appear to have type (m −1, m), cancellation makes their quotient type (m −1, m −1).
- The barycentric form can place poles away from the support points, despite its partial-fraction appearance.
- The resulting function is a type (m −1, m −1) rational interpolant, with removable singularities at support points and values r(zj) = fj.
- The representation spans type (m −1, m −1) rational functions with no poles at the support points as the data and nonzero weights vary.
- For suitable support points, the functions 1/(z −zj) can be better conditioned than polynomial quotients and vector-fitting partial fractions.
3. Core AAA algorithm.
The core AAA iteration greedily adds the sample point with the largest current residual, then computes barycentric weights by least squares and SVD. Its residual norm is nonincreasing with degree, though strict decrease is not guaranteed.
- The new barycentric weights are computed by solving a linear least-squares problem over unselected sample points.
- The step-m approximant is generically type (m −1, m −1) and interpolates the selected data values.
- The least-squares problem uses the discrete 2-norm and requires at least m unselected points, so m ≤M/2.
- At each step, AAA selects the next support point where the current nonlinear residual has maximum absolute value.
- AAA computes the least-squares weights with the final right singular vector from a reduced SVD of the problem matrix.
- The linearized residual norm σmin(A(m)) is nonincreasing as m grows, but strict decrease is not guaranteed.
- The SVD-based iteration has complexity O(Mm3), usually modest because applications typically use small m.
4. Matlab code.
The Matlab implementation returns a callable rational approximant and supports real or complex scalar, vector, and matrix evaluation. Examples show accurate recovery of poles and residues for gamma and tangent functions, while safeguards and cleanup remain relevant in practice.
- The aaa.m implementation returns r as a function handle, with optional cleanup code separated from the core algorithm.
- The gamma-function example uses 100 samples to construct a type (9, 9) approximant on a larger interval.
- The gamma approximant matches its first four poles at 0, −1, −2, and −3 to 15, 15, 7, and 3 digits, respectively.
- The listed code lacks all safeguards of fully developed software, although the Chebfun realization includes some of them.
- The tangent example approximates 1000 spiral points and converges in 12 steps, with the reported error reaching 1.30e-13.
- For tangent, the first pole pair matches ±1 to 15 digits, the next pair matches ±3 to 7 digits, and zeros show similar agreement.
5. Removing numerical Froissart doublets.
Numerical Froissart doublets can arise in rational approximation and are addressed through a cleanup procedure. Across nine applications, cleanup was needed in only one application and reduced the detected doublets there.
- Without cleanup, the example produced 58 numerical Froissart doublets near the unit circle at m = 100.With cleanup enabled, only one Froissart doublet remained after the first 58 were eliminated.
- Cleanup reduced the example's detected doublets to one, which was identified as an artifact introduced when the first 58 were removed.
- The related literature primarily studies Froissart doublets in Padé approximation and p/q representations, while one cited work treats barycentric representations.
- Among nine applications computed with aaa, numerical Froissart doublets appeared only in Application 6.6, where six were removed by cleanup.
6. Applications.
The applications show AAA rational approximation handling analytic, branch-point, and meromorphic functions, often efficiently, while exposing limits of extrapolating boundary accuracy throughout an enclosed domain.
- 6.1. Analytic functions in the unit disk.: Rational type (7, 7) approximates tan(z) more accurately than a degree-52 polynomial on 128 unit-circle points.Polynomials require degree n ≥52 for 10-digit accuracy because tan has poles at ±π/2, whereas rational approximants can capture exterior poles.
- 6.2. Analytic functions with nearby branch points; comparison with ratdisk.: AAA approximates log(1.1 −z) with rate Oε(ρ−n), ρ ≈9.3, versus degree-2n polynomials at Oε(1.1−2n).The function has a branch point outside the disk rather than poles, yet rational approximation remains effective.
- 6.2. Analytic functions with nearby branch points; comparison with ratdisk.: Zeros and poles of the branch-point approximant interlace along a curve near the branch cut [1, ∞).The poles tend to line up along branch cuts, while the core AAA algorithm need not preserve real symmetry.
- 6.3. Meromorphic functions in the unit disk from boundary values.: AAA achieves 13-digit accuracy on the unit circle for tan(βz) using types (14, 14), (28, 28), (49, 49), and (62, 62) as β increases.The corresponding functions have 2, 10, 40, and 162 poles in the disk; the β = 256 approximant is accurate near the circle but not throughout the disk.
- 6.3. Meromorphic functions in the unit disk from boundary values.: Boundary accuracy does not ensure accuracy throughout the enclosed disk or reliable recovery of interior zeros and poles.For β = 16 the approximation is accurate throughout the disk, whereas for β = 64 it is accurate on the circle but not throughout the disk.
- 6.4. Meromorphic functions in the unit disk from boundary and interior values.: On non-circular domains, barycentric AAA avoids the domain-dependent polynomial bases and their potential exponential ill-conditioning.The method is designed for point sets including disconnected, irregular, and possibly unbounded regions.
6.5. Approximation in other connected domains.
The applications show AAA approximating functions on complex, disconnected, real, unbounded, and model-reduction domains, with strong accuracy but domain- and symmetry-related caveats.
- Complex domains: AAA approximates 1/J0(z) accurately on 2000 random points in a complex rectangle, recovering poles that nearly match the function's real poles.The algorithm exits at m = 13 with a type (12, 12) approximant.
- Disconnected domains: AAA efficiently approximates sign(Re(z)) on a disconnected square-and-circle set, though six numerical Froissart doublets require cleanup.Convergence occurs at m = 51, yielding a type (43, 43) approximant after cleanup.
- Real intervals: For |x| on [−1, 1], core AAA develops poles in the interval at odd degrees and eventually at all degrees from n = 8 onward.These poles are Froissart doublets that arise even in exact arithmetic, so sample-set errors can conceal infinite interval errors.
- Real intervals: Exploiting symmetry by transforming |x| to √x enables AAA approximations through n = 80 with better than 10 digits of accuracy and no interval poles.The paper presents this as effective for the |x| problem but not universally applicable.
- Unbounded domains and model reduction: AAA handles exp(x) on (−∞, 0] over seven orders of magnitude and approximates a clamped-beam transfer function by capturing pole pairs successively.The exp(x) computation takes about 1 second, while all clamped-beam approximant poles lie in the left half-plane.
7. Modified algorithm to treat large and small data symmetrically.
The modified algorithm seeks symmetry between large and small function values, extending AAA to data that may include infinities while introducing trade-offs in other invariance properties.
- Invariance trade-offs: The modification restores reciprocal symmetry for r and 1/r but breaks translation- and scale-invariance in f because it privileges the constant 1.Scale-invariance could instead use a threshold C that scales with f, such as the median of |f(Z)|.
- Large-small symmetry: A modified AAA formulation can treat function values equal to ∞ like other sample values, supporting approximation problems with infinite data.This is motivated by symmetry between large and small values on the Riemann sphere.
- Related symmetry variants: The core AAA algorithm can break even, odd, or Hermitian symmetries, so variants may choose support points to preserve them.The paper notes that these symmetry violations can have practical and cosmetic consequences.
8. Modified algorithms to impose even, odd, or real symmetry.
Symmetry-aware variants address cases where core AAA introduces undesirable real poles or otherwise fails to preserve even, odd, or Hermitian structure.
- Even and odd symmetry: Even and odd functions can often be handled by changing variables from z to z^1/2 rather than modifying the AAA algorithm.The paper presents this as an effective strategy for such cases, including |x| on [−1, 1].
- Hermitian symmetry: For Hermitian-symmetric real data, core AAA still increases rational type one step at a time, forcing real poles at odd steps near singular points such as x = 0.The paper states that it does not offer a solution to this problem.
- Unequal rational types: General type (µ, ν) approximation requires restricting barycentric weights to a subspace that constrains numerator and denominator degrees.This extends the representation beyond equal-degree rational functions and adds linear-algebra features.
9. Modified algorithm for approximations of type (µ, ν).
Unequal-degree rational approximation requires generalized barycentric representations and raises practical questions about controlling denominator degree.
- Generalized representations: For type (µ, ν) with µ ≠ ν, barycentric weights must be restricted to an appropriate subspace to enforce the desired numerator-versus-denominator degree constraint.The extension applies the barycentric representation while adding new linear-algebra features.
- Denominator degree: A key practical issue is finding representations with minimal denominator degree for goals such as recovering information from boundary values of meromorphic functions.The passage identifies analogous extensions of AAA as useful future development.
- Algorithm variants: The paper frames unequal-degree approximation as one context where variants of the core AAA algorithm are appropriate.This discussion is part of a broader set of algorithmic variants.
10. Other variants.
The paper describes AAA variants for weighting, continuous sample sets, confluent data, non-interpolatory approximation, iterative reweighting, and vector or matrix-valued problems.
- Weighted and continuous problems: Weighted least squares can emphasize selected regions while preventing densely sampled regions from receiving disproportionate weight.Weights scale rows of the least-squares matrix, and may also be applied in the greedy nonlinear step.
- Weighted and continuous problems: Continuous AAA replaces discrete support-point selection and least-squares computation with continuous optimization problems.The barycentric representation and support points remain discrete even when the sample set is continuous.
- Weighted and continuous problems: For continuously defined problems, enlarging a discrete sample set during computation is suggested, but success is not guaranteed.
- Confluent and non-interpolatory approximation: Confluent sample points extend AAA to matching derivative data alongside function values, requiring additional Cauchy-matrix columns and Loewner-matrix rows.Derivative conditions also require decisions about relative row weighting.
- Confluent and non-interpolatory approximation: Non-interpolatory AAA separates numerator and denominator parameters, so approximants need not interpolate sample data.Experiments suggest similar performance to interpolatory approximants in many applications, but the expanded Loewner-matrix SVD takes about four times as long.
- Further extensions: The paper also discusses Lawson-style iterative reweighting and extensions from scalar functions and variables to vector or matrix approximation.
11. Comparisons with vector fitting and other algorithms.
The comparison emphasizes AAA’s adaptive barycentric basis and contrasts it with vector fitting, RKFIT, and other representations that use different parameterizations or initialization requirements.
- AAA’s representation: AAA’s two essential features—barycentric representation and adaptive support-point selection—avoid exponential ill-conditioning.
- Vector fitting: Vector fitting uses partial fractions, fixes the pole count in advance, requires initial pole estimates, and iteratively adjusts poles rather than support points.
- RKFIT: RKFIT minimizes a least-squares norm using orthogonal rational functions and an iterative adjustment process distinct from AAA’s representation.
- Vector fitting: In the clamped beam experiment, vector fitting failed to converge satisfactorily, apparently because of initial pole guesses, while a quadrature-based variant had some success.
- Conditioning: Partial-fraction bases can exhibit exponential condition-number growth, whereas AAA’s barycentric bases avoid comparable instability in the illustrated application.The figure compares monomial, partial-fraction, and AAA barycentric bases.
- Conditioning: Across the examples, most Cauchy matrices had condition numbers near 1 and at most 40, with larger values linked to closely spaced support points.Reported exceptions include approximately 10^4 and 10^8.
- AAA’s representation: Adaptive support points produce exceptionally well-conditioned bases for representing numerators and denominators.The paper presents barycentric rational representation as flexible and robust for floating-point computation.
12. Conclusion.
The conclusion presents adaptively supported barycentric quotients as an effective computational approach, illustrates AAA on the Riemann zeta function, and leaves broader extensions open.
- Conclusion: The paper concludes that barycentric quotients with adaptively selected support points are a strikingly effective way to compute with rational functions.
- Conclusion: The authors explicitly state that AAA is not the last word in rational approximation.
- Conclusion: AAA evaluates the Riemann zeta function at 100 complex samples and constructs a type (29, 29) approximant in slightly less than one second on the authors’ desktop computer.
- Conclusion: The approximant identifies a pole and zero matching the corresponding zeta-function values in all but the last two digits.The reported pole is near 1 with residue near 1, and the zero is near 0.5 + 14.1347i.