For Jean-Paul Berrut, the pioneer of numerical algorithms based on rational barycentric representations, on his 65th birthday.Abstract. 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.

The state-of-the-art algorithms for solving the trust-region subproblem (TRS) are based on an iterative process, involving solutions of many linear systems, eigenvalue problems, subspace optimization, or line search steps. A relatively underappreciated fact, due to Gander, Golub, and von Matt [Linear Algebra Appl., 114 (1989), pp. 815-839], is that TRSs can be solved by one generalized eigenvalue problem, with no outer iterations. In this paper we rediscover this fact and discover its great practicality, which exhibits good performance both in accuracy and efficiency. Moreover, we generalize the approach in various directions, namely by allowing for an ellipsoidal constraint, dealing with the so-called hard case, and obtaining approximate solutions efficiently when high accuracy is unnecessary. We demonstrate that the resulting algorithm is a general-purpose TRS solver, effective both for dense and large-sparse problems, including the so-called hard case. Our algorithm is easy to implement: its essence is a few lines of MATLAB code.

Abstract. Spectral divide and conquer algorithms solve the eigenvalue problem for all the eigenvalues and eigenvectors by recursively computing an invariant subspace for a subset of the spectrum and using it to decouple the problem into two smaller subproblems. A number of such algorithms have been developed over the last 40 years, often motivated by parallel computing and, most recently, with the aim of achieving minimal communication costs. However, none of the existing algorithms has been proved to be backward stable, and they all have a significantly higher arithmetic cost than the standard algorithms currently used. We present new spectral divide and conquer algorithms for the symmetric eigenvalue problem and the singular value decomposition that are backward stable, achieve lower bounds on communication costs recently derived by Ballard, Demmel, Holtz, and Schwartz, and have operation counts within a small constant factor of those for the standard algorithms. The new algorithms are built on the polar decomposition and exploit the recently developed QR-based dynamically weighted Halley algorithm of Nakatsukasa, Bai, and Gygi, which computes the polar decomposition using a cubically convergent iteration based on the building blocks of QR factorization and matrix multiplication. [24]).An important step in the practical development of spectral divide and conquer algorithms was the toolbox of Bai and Demmel [5], [6], which provides the building blocks for constructing algorithms via the Newton iteration for the matrix sign function. A parallel implementation based on these ideas is described in [8].

Abstract. The symmetric eigenvalue decomposition and the singular value decomposition (SVD) are fundamental matrix decompositions with many applications. Conventional algorithms for computing these decompositions are suboptimal in view of recent trends in computer architectures, which require minimizing communication together with arithmetic costs. Spectral divideand-conquer algorithms, which recursively decouple the problem into two smaller subproblems, can achieve both requirements. Such algorithms can be constructed with the polar decomposition playing two key roles: it forms a bridge between the symmetric eigendecomposition and the SVD, and its connection to the matrix sign function naturally leads to spectral-decoupling. For computing the polar decomposition, the scaled Newton and QDWH iterations are two of the most popular algorithms, as they are backward stable and converge in at most nine and six iterations, respectively. Following this framework, we develop a higher-order variant of the QDWH iteration for the polar decomposition. The key idea of this algorithm comes from approximation theory: we use the best rational approximant for the scalar sign function due to Zolotarev in 1877. The algorithm exploits the extraordinary property enjoyed by the sign function that a high-degree Zolotarev function (best rational approximant) can be obtained by appropriately composing low-degree Zolotarev functions. This lets the algorithm converge in just two iterations in double-precision arithmetic, with the whopping rate of convergence seventeen. The resulting algorithms for the symmetric eigendecompositions and the SVD have higher arithmetic costs than the QDWH-based algorithms, but are better-suited for parallel computing and exhibit excellent numerical backward stability.

We introduce a dynamically weighted Halley (DWH) iteration for computing the polar decomposition of a matrix, and we prove that the new method is globally and asymptotically cubically convergent. For matrices with condition number no greater than 10 16 , the DWH method needs at most six iterations for convergence with the tolerance 10 −16. The Halley iteration can be implemented via QR decompositions without explicit matrix inversions. Therefore, it is an inverse free communication friendly algorithm for the emerging multicore and hybrid high performance computing systems.