Abstract. A spectral method is developed for the direct solution of linear ordinary differential equations with variable coefficients. The method leads to matrices which are almost banded, and a numerical solver is presented that takes Ø(m 2 n) operations, where m is the number of Chebyshev points needed to resolve the coefficients of the differential operator and n is the number of Chebyshev coefficients needed to resolve the solution to the differential equation. We prove stability of the method by relating it to a diagonally preconditioned system which has a bounded condition number, in a suitable norm. For Dirichlet boundary conditions, this implies stability in the standard 2-norm. An adaptive QR factorization is developed to efficiently solve the resulting linear system and automatically choose the optimal number of Chebyshev coefficients needed to represent the solution. The resulting algorithm can efficiently and reliably solve for solutions that require as many as a million unknowns.

Abstract. An efficient algorithm for the accurate computation of Gauss-Legendre and GaussJacobi quadrature nodes and weights is presented. The algorithm is based on Newton's root-finding method with initial guesses and function evaluations computed via asymptotic formulae. The n-point quadrature rule is computed in O(n) operations to an accuracy of essentially double precision for any n ≥ 100.

Abstract. An object-oriented Matlab system is described that extends the capabilities of Chebfun to smooth functions of two variables defined on rectangles. Functions are approximated to essentially machine precision by using iterative Gaussian elimination with complete pivoting to form "chebfun2" objects representing low rank approximations. Operations such as integration, differentiation, function evaluation, and transforms are particularly efficient. Global optimization, the singular value decomposition, and rootfinding are also extended to chebfun2 objects. Numerical applications are presented.

Matrices with displacement structure such as Pick, Vandermonde, and Hankel matrices appear in a diverse range of applications. In this paper, we use an extremal problem involving rational functions to derive explicit bounds on the singular values of such matrices. For example, we show that the kth singular value of a real n × n positive definite Hankel matrix, Hn, is bounded by Cρ −k/ log n H 2 with explicitly given constants C > 0 and ρ > 1, where Hn 2 is the spectral norm. This means that a real n × n positive definite Hankel matrix can be approximated, up to an accuracy of ǫ Hn 2 with 0 < ǫ < 1, by a rank O(log n log(1/ǫ)) matrix. Analogous results are obtained for Pick, Cauchy, real Vandermonde, Löwner, and certain Krylov matrices.

Matrices of (approximate) low rank are pervasive in data science, appearing in recommender systems, movie preferences, topic models, medical records, and genomics. While there is a vast literature on how to exploit low rank structure in these datasets, there is less attention on explaining why the low rank structure appears in the first place. Here, we explain the effectiveness of low rank models in data science by considering a simple generative model for these matrices: we suppose that each row or column is associated to a (possibly high dimensional) bounded latent variable, and entries of the matrix are generated by applying a piecewise analytic function to these latent variables. These matrices are in general full rank. However, we show that we can approximate every entry of an m × n matrix drawn from this model to within a fixed absolute error by a low rank matrix whose rank grows as O(log(m + n)). Hence any sufficiently large matrix from such a latent variable model can be approximated, up to a small entrywise error, by a low rank matrix.

Abstract. A collection of algorithms is described for numerically computing with smooth functions defined on the unit sphere. Functions are approximated to essentially machine precision by using a structure-preserving iterative variant of Gaussian elimination together with the double Fourier sphere method. We show that this procedure allows for stable differentiation, reduces the oversampling of functions near the poles, and converges for certain analytic functions. Operations such as function evaluation, differentiation, and integration are particularly efficient and can be computed by essentially one-dimensional algorithms. A highlight is an optimal complexity direct solver for Poisson's equation on the sphere using a spectral method. Without parallelization, we solve Poisson's equation with 100 million degrees of freedom in 1 minute on a standard laptop. Numerical results are presented throughout. In a companion paper (part II) we extend the ideas presented here to computing with functions on the disk.Key words. low rank approximation, Gaussian elimination, functions, approximation theory AMS subject classification. 65D05 DOI. 10.1137/15M10458551. Introduction. Spherical geometries are universal in computational science and engineering, arising in weather forecasting and climate modeling [11,13,17,24,28,31,35], geophysics [16,48], and astrophysics [3,8,39]. At various levels these applications all require the approximation of functions defined on the surface of the unit sphere. For such computational tasks, a standard approach is to use longitude-latitude coordinates (λ, θ) ∈ [−π, π] × [0, π], where λ and θ denote the azimuthal and polar angles, respectively. Thus, computations with functions on the sphere can be conveniently related to analogous tasks involving functions defined on a rectangular domain. This is a useful observation that, unfortunately, also has many severe disadvantages due to artificial pole singularities introduced by the coordinate transform.In this paper, we synthesize a classic technique known as the double Fourier sphere (DFS) method [6,18,28,31,50] together with new algorithmic techniques in low rank function approximation [5,42]. This alleviates many of the drawbacks inherent with standard coordinate transforms. Our approximants have several attractive properties: (1) no artificial pole singularities, (2) a representation that allows for fast algorithms, (3) a structure so that differentiation is stable, and (4) an underlying interpolation grid that rarely oversamples functions near the poles.

Abstract. A spectral method for solving linear partial differential equations (PDEs) with variable coefficients and general boundary conditions defined on rectangular domains is described, based on separable representations of partial differential operators and the one-dimensional ultraspherical spectral method. If a partial differential operator is of splitting rank 2, such as the operator associated with Poisson or Helmholtz, the corresponding PDE is solved via a generalized Sylvester matrix equation, and a bivariate polynomial approximation of the solution of degree (nx, ny) is computed in O((nxny) 3/2 ) operations. Partial differential operators of splitting rank ≥ 3 are solved via a linear system involving a block-banded matrix in O(min(n 3 x ny, nxn 3 y )) operations. Numerical examples demonstrate the applicability of our 2D spectral method to a broad class of PDEs, which includes elliptic and dispersive time-evolution equations. The resulting PDE solver is written in Matlab and is publicly available as part of Chebfun. It can resolve solutions requiring over a million degrees of freedom in under 60 seconds. An experimental implementation in the Julia language can currently perform the same solve in 10 seconds.

Analogues of singular value decomposition (SVD), QR, LU and Cholesky factorizations are presented for problems in which the usual discrete matrix is replaced by a 'quasimatrix', continuous in one dimension, or a 'cmatrix', continuous in both dimensions. Two challenges arise: the generalization of the notions of triangular structure and row and column pivoting to continuous variables (required in all cases except the SVD, and far from obvious), and the convergence of the infinite series that define the cmatrix factorizations. Our generalizations of triangularity and pivoting are based on a new notion of a 'triangular quasimatrix'. Concerning convergence of the series, we prove theorems asserting convergence provided the functions involved are sufficiently smooth.