BZPEERPapers, people, evidence and research paths
Menu

Discover

SearchFind papers, people, topics, institutions and opportunities from one place.GraphSee how papers, authors, topics and evidence connect.

Research tools

HomeContinue from the papers, people and topics that matter now.CollectionsTurn saved papers into reading lists, evidence maps and shareable context.ResearchReconstruct the path that led you from search to paper to question.CollaborationMove from a paper, topic or expert request into a focused thread.PublishCreate full-text research articles with metadata, formulas and versions.

Network

ExpertsFind specialists whose work and reviews explain why they can help.CommunitiesJoin research groups organized around papers, topics and shared evidence.

Opportunities

ListingsBrowse roles, calls and collaborations connected to real research context.

Account

Log inReturn to your saved papers, graph views and active threads.Create accountStart saving papers, building a research trail and joining teams.
Log inCreate account

Source author record

Alex Townsend

Alex Townsend appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record
math.NANumerical AnalysisMachine LearningArtificial IntelligenceInformation Theorymath.APmath.CVmath.DSmath.ITmath.PRmath.RAmath.STnlin.AOnlin.CDStatistics Theory

Catalog footprint

What is connected

22works
15topics
4close collaborators

Actions

Connect this record

Log in to claim
Open graphBrowse works

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

22 published item(s)

preprint2022arXiv

A generalization of the randomized singular value decomposition

The randomized singular value decomposition (SVD) is a popular and effective algorithm for computing a near-best rank $k$ approximation of a matrix $A$ using matrix-vector products with standard Gaussian vectors. Here, we generalize the randomized SVD to multivariate Gaussian vectors, allowing one to incorporate prior knowledge of $A$ into the algorithm. This enables us to explore the continuous analogue of the randomized SVD for Hilbert--Schmidt (HS) operators using operator-function products with functions drawn from a Gaussian process (GP). We then construct a new covariance kernel for GPs, based on weighted Jacobi polynomials, which allows us to rapidly sample the GP and control the smoothness of the randomly generated functions. Numerical examples on matrices and HS operators demonstrate the applicability of the algorithm.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Data-driven discovery of Green's functions with human-understandable deep learning

There is an opportunity for deep learning to revolutionize science and technology by revealing its findings in a human interpretable manner. To do this, we develop a novel data-driven approach for creating a human-machine partnership to accelerate scientific discovery. By collecting physical system responses under excitations drawn from a Gaussian process, we train rational neural networks to learn Green's functions of hidden linear partial differential equations. These functions reveal human-understandable properties and features, such as linear conservation laws and symmetries, along with shock and singularity locations, boundary effects, and dominant modes. We illustrate the technique on several examples and capture a range of physics, including advection-diffusion, viscous shocks, and Stokes flow in a lid-driven cavity.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Learning elliptic partial differential equations with randomized linear algebra

Given input-output pairs of an elliptic partial differential equation (PDE) in three dimensions, we derive the first theoretically-rigorous scheme for learning the associated Green's function $G$. By exploiting the hierarchical low-rank structure of $G$, we show that one can construct an approximant to $G$ that converges almost surely and achieves a relative error of $\mathcal{O}(Γ_ε^{-1/2}\log^3(1/ε)ε)$ using at most $\mathcal{O}(ε^{-6}\log^4(1/ε))$ input-output training pairs with high probability, for any $0<ε<1$. The quantity $0<Γ_ε\leq 1$ characterizes the quality of the training dataset. Along the way, we extend the randomized singular value decomposition algorithm for learning matrices to Hilbert--Schmidt operators and characterize the quality of covariance kernels for PDE learning.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Learning Green's functions associated with time-dependent partial differential equations

Neural operators are a popular technique in scientific machine learning to learn a mathematical model of the behavior of unknown physical systems from data. Neural operators are especially useful to learn solution operators associated with partial differential equations (PDEs) from pairs of forcing functions and solutions when numerical solvers are not available or the underlying physics is poorly understood. In this work, we attempt to provide theoretical foundations to understand the amount of training data needed to learn time-dependent PDEs. Given input-output pairs from a parabolic PDE in any spatial dimension $n\geq 1$, we derive the first theoretically rigorous scheme for learning the associated solution operator, which takes the form of a convolution with a Green's function $G$. Until now, rigorously learning Green's functions associated with time-dependent PDEs has been a major challenge in the field of scientific machine learning because $G$ may not be square-integrable when $n>1$, and time-dependent PDEs have transient dynamics. By combining the hierarchical low-rank structure of $G$ together with randomized numerical linear algebra, we construct an approximant to $G$ that achieves a relative error of $\smash{\mathcal{O}(Γ_ε^{-1/2}ε)}$ in the $L^1$-norm with high probability by using at most $\smash{\mathcal{O}(ε^{-\frac{n+2}{2}}\log(1/ε))}$ input-output training pairs, where $Γ_ε$ is a measure of the quality of the training dataset for learning $G$, and $ε>0$ is sufficiently small.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Bounding Zolotarev numbers using Faber rational functions

By closely following a construction by Ganelius, we construct Faber rational functions that allow us to derive tight and explicit bounds on Zolotarev numbers. We use our results to bound the singular values of matrices, including complex-valued Cauchy matrices and Vandermonde matrices with nodes inside the unit disk. We construct Faber rational functions using doubly-connected conformal maps and use their zeros and poles to supply shift parameters in the alternating direction implicit method.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Dense networks that do not synchronize and sparse ones that do

For any network of identical Kuramoto oscillators with identical positive coupling, there is a critical connectivity above which the system is guaranteed to converge to the in-phase synchronous state, for almost all initial conditions. But the precise value of this critical connectivity remains unknown. In 2018, Ling, Xu, and Bandeira proved that if each oscillator is coupled to at least 79.29 percent of all the others, global synchrony is ensured. In 2019, Lu and Steinerberger improved this bound to 78.89 percent. Here, by focusing on circulant networks, we find clues that the critical connectivity may be exactly 75 percent. Our methods yield a slight improvement on the best known lower bound on the critical connectivity, from $68.18\%$ to $68.28\%$. We also consider the opposite end of the connectivity spectrum, where the networks are sparse rather than dense. In this regime, we ask how few edges one needs to add to a ring of $n$ oscillators to turn it into a globally synchronizing network. We prove a partial result: all the twisted states in a ring of size $n=2^m$ can be destabilized by adding just $\mathcal{O}(n \log_2 n)$ edges. To finish the proof, one also needs to rule out all other candidate attractors. We have done this for $n=8$ with computational algebraic geometry, but the problem remains open for larger $n$. Thus, even for systems as simple as Kuramoto oscillators, much remains to be learned about dense networks that do not globally synchronize and sparse ones that do.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Error localization of best L1 polynomial approximants

An important observation in compressed sensing is that the $\ell_0$ minimizer of an underdetermined linear system is equal to the $\ell_1$ minimizer when there exists a sparse solution vector and a certain restricted isometry property holds. Here, we develop a continuous analogue of this observation and show that the best $L_0$ and $L_1$ polynomial approximants of a polynomial that is corrupted on a set of small measure are nearly equal. We go on to demonstrate an error localization property of best $L_1$ polynomial approximants and use our observations to develop an improved algorithm for computing best $L_1$ polynomial approximants to continuous functions.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

On the compressibility of tensors

Tensors are often compressed by expressing them in low rank tensor formats. In this paper, we develop three methodologies that bound the compressibility of a tensor: (1) Algebraic structure, (2) Smoothness, and (3) Displacement structure. For each methodology, we derive bounds on storage costs that partially explain the abundance of compressible tensors in applied mathematics. For example, we show that the solution tensor $\mathcal{X} \in \mathbb{C}^{n \times n \times n}$ of a discretized Poisson equation $-\nabla^2 u =1$ on $[-1,1]^3$ with zero Dirichlet conditions can be approximated to a relative accuracy of $0<ε<1$ in the Frobenius norm by a tensor in tensor-train format with $\mathcal{O}(n (\log n)^2 (\log(1/ε))^2)$ degrees of freedom. As this bound is constructive, we are also able to solve this equation spectrally with $\mathcal{O}(n (\log n)^3 (\log(1/ε))^3)$ complexity.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Computing with functions in spherical and polar geometries I. The sphere

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 one 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.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Fast polynomial transforms based on Toeplitz and Hankel matrices

Many standard conversion matrices between coefficients in classical orthogonal polynomial expansions can be decomposed using diagonally-scaled Hadamard products involving Toeplitz and Hankel matrices. This allows us to derive $\smash{\mathcal{O}(N(\log N)^2)}$ algorithms, based on the fast Fourier transform, for converting coefficients of a degree $N$ polynomial in one polynomial basis to coefficients in another. Numerical results show that this approach is competitive with state-of-the-art techniques, requires no precomputational cost, can be implemented in a handful of lines of code, and is easily adapted to extended precision arithmetic.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Gaussian elimination corrects pivoting mistakes

Gaussian elimination (GE) is the archetypal direct algorithm for solving linear systems of equations and this has been its primary application for thousands of years. In the last decade, GE has found another major use as an iterative algorithm for low rank approximation. In this setting, GE is often employed with complete pivoting and designed to allow for non-optimal pivoting, i.e., pivoting mistakes, that could render GE numerically unstable when implemented in floating point arithmetic. While it may appear that pivoting mistakes could accumulate and lead to a large growth factor, we show that later GE steps correct earlier pivoting mistakes, even while more are being made. In short, GE is very robust to non-optimal pivots, allowing for its iterative variant to flourish.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Numerical instability of resultant methods for multidimensional rootfinding

Hidden-variable resultant methods are a class of algorithms for solving multidimensional polynomial rootfinding problems. In two dimensions, when significant care is taken, they are competitive practical rootfinders. However, in higher dimensions they are known to miss zeros, calculate roots to low precision, and introduce spurious solutions. We show that the hidden variable resultant method based on the Cayley (Dixon or Bézout) matrix is inherently and spectacularly numerically unstable by a factor that grows exponentially with the dimension. We also show that the Sylvester matrix for solving bivariate polynomial systems can square the condition number of the problem. In other words, two popular hidden variable resultant methods are numerically unstable, and this mathematically explains the difficulties that are frequently reported by practitioners. Regardless of how the constructed polynomial eigenvalue problem is solved, severe numerical difficulties will be present. Along the way, we prove that the Cayley resultant is a generalization of Cramer's rule for solving linear systems and generalize Clenshaw's algorithm to an evaluation scheme for polynomials expressed in a degree-graded polynomial basis.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

On the singular values of matrices with displacement structure

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 $k$th singular value of a real $n\times n$ positive definite Hankel matrix, $H_n$, is bounded by $Cρ^{-k/\log n}\|H\|_2$ with explicitly given constants $C>0$ and $ρ>1$, where $\|H_n\|_2$ is the spectral norm. This means that a real $n\times n$ positive definite Hankel matrix can be approximated, up to an accuracy of $ε\|H_n\|_2$ with $0<ε<1$, by a rank $\mathcal{O}(\log n\log(1/ε) )$ matrix. Analogous results are obtained for Pick, Cauchy, real Vandermonde, Löwner, and certain Krylov matrices.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Stable extrapolation of analytic functions

This paper examines the problem of extrapolation of an analytic function for $x > 1$ given perturbed samples from an equally spaced grid on $[-1,1]$. Mathematical folklore states that extrapolation is in general hopelessly ill-conditioned, but we show that a more precise statement carries an interesting nuance. For a function $f$ on $[-1,1]$ that is analytic in a Bernstein ellipse with parameter $ρ> 1$, and for a uniform perturbation level $ε$ on the function samples, we construct an asymptotically best extrapolant $e(x)$ as a least squares polynomial approximant of degree $M^*$ given explicitly. We show that the extrapolant $e(x)$ converges to $f(x)$ pointwise in the interval $I_ρ\in[1,(ρ+ρ^{-1})/2)$ as $ε\to 0$, at a rate given by a $x$-dependent fractional power of $ε$. More precisely, for each $x \in I_ρ$ we have \[ |f(x) - e(x)| = \mathcal{O}\left( ε^{-\log r(x) / \logρ} \right), \qquad\qquad r(x) = \frac{x+\sqrt{x^2-1}}ρ, \] up to log factors, provided that the oversampling conditioning is satisfied. That is, \[ M^* \leq \frac{1}{2} \sqrt{N}, \] which is known to be needed from approximation theory. In short, extrapolation enjoys a weak form of stability, up to a fraction of the characteristic smoothness length. The number of function samples, $N+1$, does not bear on the size of the extrapolation error provided that it obeys the oversampling condition. We also show that one cannot construct an asymptotically more accurate extrapolant from $N+1$ equally spaced samples than $e(x)$, using any other linear or nonlinear procedure. The proofs involve original statements on the stability of polynomial approximation in the Chebyshev basis from equally spaced samples and these are expected to be of independent interest.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Vector spaces of linearizations for matrix polynomials: a bivariate polynomial approach

We revisit the landmark paper [D. S. Mackey, N. Mackey, C. Mehl, and V. Mehrmann, SIAM J. Matrix Anal. Appl., 28 (2006), pp.~971--1004] and, by viewing matrices as coefficients for bivariate polynomials, we provide concise proofs for key properties of linearizations for matrix polynomials. We also show that every pencil in the double ansatz space is intrinsically connected to a Bézout matrix, which we use to prove the eigenvalue exclusion theorem. In addition our exposition allows for any polynomial basis and for any field. The new viewpoint also leads to new results. We generalize the double ansatz space by exploiting its algebraic interpretation as a space of Bézout pencils to derive new linearizations with potential applications in the theory of structured matrix polynomials. Moreover, we analyze the conditioning of double ansatz space linearizations in the important practical case of a Chebyshev basis.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

A fast analysis-based discrete Hankel transform using asymptotic expansions

A fast and numerically stable algorithm is described for computing the discrete Hankel transform of order $0$ as well as evaluating Schlömilch and Fourier--Bessel expansions in $\mathcal{O}(N(\log N)^2/\log\!\log N)$ operations. The algorithm is based on an asymptotic expansion for Bessel functions of large arguments, the fast Fourier transform, and the Neumann addition formula. All the algorithmic parameters are selected from error bounds to achieve a near-optimal computational cost for any accuracy goal. Numerical results demonstrate the efficiency of the resulting algorithm.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

A fast FFT-based discrete Legendre transform

An $\mathcal{O}(N(\log N)^2/\log\!\log N)$ algorithm for computing the discrete Legendre transform and its inverse is described. The algorithm combines a recently developed fast transform for converting between Legendre and Chebyshev coefficients with a Taylor series expansion for Chebyshev polynomials about equally-spaced points in the frequency domain. Both components are based on the FFT, and as an intermediate step we obtain an $\mathcal{O}(N\log N)$ algorithm for evaluating a degree $N-1$ Chebyshev expansion at an $N$-point Legendre grid. Numerical results are given to demonstrate performance and accuracy.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

The automatic solution of partial differential equations using a global spectral method

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 $(n_x,n_y)$ is computed in $\mathcal{O}((n_x n_y)^{3/2})$ operations. Partial differential operators of splitting rank $\geq 3$ are solved via a linear system involving a block-banded matrix in $\mathcal{O}(\min(n_x^{3} n_y,n_x n_y^{3}))$ 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.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

A practical framework for infinite-dimensional linear algebra

We describe a framework for solving a broad class of infinite-dimensional linear equations, consisting of almost banded operators, which can be used to resepresent linear ordinary differential equations with general boundary conditions. The framework contains a data structure on which row operations can be performed, allowing for the solution of linear equations by the adaptive QR approach. The algorithm achieves $O(n^{\rm opt})$ complexity, where $n^{\rm opt}$ is the number of degrees of freedom required to achieve a desired accuracy, which is determined adaptively. In addition, special tensor product equations, such as partial differential equations on rectangles, can be solved by truncating the operator in the $y$-direction with $n_y$ degrees of freedom and using a generalized Schur decomposition to upper triangularize, before applying the adaptive QR approach to the $x$-direction, requiring $O(n_y^2 n_x^{\rm opt})$ operations. The framework is implemented in the ApproxFun package written in the Julia programming language, which achieves highly competitive computational costs by exploiting unique features of Julia.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

Fast computation of Gauss quadrature nodes and weights on the whole real line

A fast and accurate algorithm for the computation of Gauss-Hermite and generalized Gauss-Hermite quadrature nodes and weights is presented. The algorithm is based on Newton's method with carefully selected initial guesses for the nodes and a fast evaluation scheme for the associated orthogonal polynomial. In the Gauss-Hermite case the initial guesses and evaluation scheme rely on explicit asymptotic formulas. For generalized Gauss-Hermite, the initial guesses are furnished by sampling a certain equilibrium measure and the associated polynomial evaluated via a Riemann-Hilbert reformulation. In both cases the $n$-point quadrature rule is computed in $\mathcal{O}(n)$ operations to an accuracy that is close to machine precision. For sufficiently large $n$, some of the quadrature weights have a value less than the smallest positive normalized floating-point number in double precision and we exploit this fact to achieve a complexity as low as $\mathcal{O}(\sqrt{n})$.

0 citations0 reviewsNo code linkedOpen access
preprint2013arXiv

Fast inverse transform sampling in one and two dimensions

We develop a computationally efficient and robust algorithm for generating pseudo-random samples from a broad class of smooth probability distributions in one and two dimensions. The algorithm is based on inverse transform sampling with a polynomial approximation scheme using Chebyshev polynomials, Chebyshev grids, and low rank function approximation. Numerical experiments demonstrate that our algorithm outperforms existing approaches.

0 citations0 reviewsNo code linkedOpen access
preprint2012arXiv

A fast and well-conditioned spectral method

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 O(m^2n) 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.

0 citations0 reviewsNo code linkedOpen access