Ben Adcock

Ben Adcock, Anders C. Hansen, Clarice Poon et al.

This paper provides an extension of compressed sensing which bridges a substantial gap between existing theory and its current use in real-world applications. It introduces a mathematical framework that generalizes the three standard pillars of compressed sensing - namely, sparsity, incoherence and uniform random subsampling - to three new concepts: asymptotic sparsity, asymptotic incoherence and multilevel random sampling. The new theorems show that compressed sensing is also possible, and reveals several advantages, under these substantially relaxed conditions. The importance of this is threefold. First, inverse problems to which compressed sensing is currently applied are typically coherent. The new theory provides the first comprehensive mathematical explanation for a range of empirical usages of compressed sensing in real-world applications, such as medical imaging, microscopy, spectroscopy and others. Second, in showing that compressed sensing does not require incoherence, but instead that asymptotic incoherence is sufficient, the new theory offers markedly greater flexibility in the design of sensing mechanisms. Third, by using asymptotic incoherence and multi-level sampling to exploit not just sparsity, but also structure, i.e. asymptotic sparsity, the new theory shows that substantially improved reconstructions can be obtained from fewer measurements.

Ben Adcock, Anders C. Hansen

We introduce a generalized framework for sampling and reconstruction in separable Hilbert spaces. Specifically, we establish that it is always possible to stably reconstruct a vector in an arbitrary Riesz basis from sufficiently many of its samples in any other Riesz basis. This framework can be viewed as an extension of that of Eldar et al. However, whilst the latter imposes stringent assumptions on the reconstruction basis, and may in practice be unstable, our framework allows for recovery in any (Riesz) basis in a manner that is completely stable. Whilst the classical Shannon Sampling Theorem is a special case of our theorem, this framework allows us to exploit additional information about the approximated vector (or, in this case, function), for example sparsity or regularity, to design a reconstruction basis that is better suited. Examples are presented illustrating this procedure.

Ben Adcock, Anders C. Hansen, Clarice Poon

Generalized sampling is a recently developed linear framework for sampling and reconstruction in separable Hilbert spaces. It allows one to recover any element in any finite-dimensional subspace given finitely many of its samples with respect to an arbitrary frame. Unlike more common approaches for this problem, such as the consistent reconstruction technique of Eldar et al, it leads to completely stable numerical methods possessing both guaranteed stability and accuracy. The purpose of this paper is twofold. First, we give a complete and formal analysis of generalized sampling, the main result of which being the derivation of new, sharp bounds for the accuracy and stability of this approach. Such bounds improve those given previously, and result in a necessary and sufficient condition, the stable sampling rate, which guarantees a priori a good reconstruction. Second, we address the topic of optimality. Under some assumptions, we show that generalized sampling is an optimal, stable reconstruction. Correspondingly, whenever these assumptions hold, the stable sampling rate is a universal quantity. In the final part of the paper we illustrate our results by applying generalized sampling to the so-called uniform resampling problem.

Ben Adcock, Anders C. Hansen

We introduce a method to reconstruct an element of a Hilbert space in terms of an arbitrary finite collection of linearly independent reconstruction vectors, given a finite number of its samples with respect to any Riesz basis. As we establish, provided the dimension of the reconstruction space is chosen suitably in relation to the number of samples, this procedure can be numerically implemented in a stable manner. Moreover, the accuracy of the resulting approximation is completely determined by the choice of reconstruction basis, meaning that the reconstruction vectors can be tailored to the particular problem at hand. An important example of this approach is the accurate recovery of a piecewise analytic function from its first few Fourier coefficients. Whilst the standard Fourier projection suffers from the Gibbs phenomenon, by reconstructing in a piecewise polynomial basis, we obtain an approximation with root exponential accuracy in terms of the number of Fourier samples and exponential accuracy in terms of the degree of the reconstruction function. Numerical examples illustrate the advantage of this approach over other existing methods.

Ben Adcock, Daan Huybrechs, Jesus Martin-Vaquero

An effective means to approximate an analytic, nonperiodic function on a bounded interval is by using a Fourier series on a larger domain. When constructed appropriately, this so-called Fourier extension is known to converge geometrically fast in the truncation parameter. Unfortunately, computing a Fourier extension requires solving an ill-conditioned linear system, and hence one might expect such rapid convergence to be destroyed when carrying out computations in finite precision. The purpose of this paper is to show that this is not the case. Specifically, we show that Fourier extensions are actually numerically stable when implemented in finite arithmetic, and achieve a convergence rate that is at least superalgebraic. Thus, in this instance, ill-conditioning of the linear system does not prohibit a good approximation. In the second part of this paper we consider the issue of computing Fourier extensions from equispaced data. A result of Platte, Trefethen & Kuijlaars states that no method for this problem can be both numerically stable and exponentially convergent. We explain how Fourier extensions relate to this theoretical barrier, and demonstrate that they are particularly well suited for this problem: namely, they obtain at least superalgebraic convergence in a numerically stable manner.

Ben Adcock

We introduce and analyze a framework for function interpolation using compressed sensing. This framework - which is based on weighted $\ell^1$ minimization - does not require a priori bounds on the expansion tail in either its implementation or its theoretical guarantees, and in the absence of noise leads to genuinely interpolatory approximations. We also establish a new recovery guarantee for compressed sensing with weighted $\ell^1$ minimization based on this framework. This guarantee has several key benefits. First, unlike existing results, it is sharp (up to constants and log factors) for large classes of functions regardless of the choice of weights. Second, by examining the measurement condition in the recovery guarantee, we are able to suggest a good overall strategy for selecting the weights. In particular, when applied to the important case of multivariate approximation with orthogonal polynomials, this weighting strategy leads to provably optimal estimates on the number of measurements required, whenever the support set of the significant coefficients is a so-called lower set. Finally, this guarantee can also be used to theoretically confirm the benefits of alternative weighting strategies where the weights are chosen based on prior support information. This provides a theoretical basis for a number of recent numerical studies showing the effectiveness of such approaches.

Ben Adcock, Anders C. Hansen, Clarice Poon

In this paper we study the problem of computing wavelet coefficients of compactly supported functions from their Fourier samples. For this, we use the recently introduced framework of generalized sampling. Our first result demonstrates that using generalized sampling one obtains a stable and accurate reconstruction, provided the number of Fourier samples grows linearly in the number of wavelet coefficients recovered. For the class of Daubechies wavelets we derive the exact constant of proportionality. Our second result concerns the optimality of generalized sampling for this problem. Under some mild assumptions we show that generalized sampling cannot be outperformed in terms of approximation quality by more than a constant factor. Moreover, for the class of so-called perfect methods, any attempt to lower the sampling ratio below a certain critical threshold necessarily results in exponential ill-conditioning. Thus generalized sampling provides a nearly-optimal solution to this problem.

Ben Adcock, Anders C. Hansen, Alexei Shadrin

We prove that any stable method for resolving the Gibbs phenomenon - that is, recovering high-order accuracy from the first $m$ Fourier coefficients of an analytic and nonperiodic function - can converge at best root-exponentially fast in $m$. Any method with faster convergence must also be unstable, and in particular, exponential convergence implies exponential ill-conditioning. This result is analogous to a recent theorem of Platte, Trefethen & Kuijlaars concerning recovery from pointwise function values on an equispaced $m$-grid. The main step in our proof is an estimate for the maximal behaviour of a polynomial of degree $n$ with bounded $m$-term Fourier series, which is related to a conjecture of Hrycak & Groechenig. In the second part of the paper we discuss the implications of our main theorem to polynomial-based interpolation and least-squares approaches for overcoming the Gibbs phenomenon. Finally, we consider the use of so-called Fourier extensions as an attractive alternative for this problem. We present numerical results demonstrating rapid convergence in a stable manner.

Ben Adcock, Daan Huybrechs

Functions that are smooth but non-periodic on a certain interval possess Fourier series that lack uniform convergence and suffer from the Gibbs phenomenon. However, they can be represented accurately by a Fourier series that is periodic on a larger interval. This is commonly called a Fourier extension. When constructed in a particular manner, Fourier extensions share many of the same features of a standard Fourier series. In particular, one can compute Fourier extensions which converge spectrally fast whenever the function is smooth, and exponentially fast if the function is analytic, much the same as the Fourier series of a smooth/analytic and periodic function. With this in mind, the purpose of this paper is to describe, analyze and explain the observation that Fourier extensions, much like classical Fourier series, also have excellent resolution properties for representing oscillatory functions. The resolution power, or required number of degrees of freedom per wavelength, depends on a user-controlled parameter and, as we show, it varies between 2 and π. The former value is optimal and is achieved by classical Fourier series for periodic functions, for example. The latter value is the resolution power of algebraic polynomial approximations. Thus, Fourier extensions with an appropriate choice of parameter are eminently suitable for problems with moderate to high degrees of oscillation.

Simone Brugiapaglia, Ben Adcock

Quadratically-constrained basis pursuit has become a popular device in sparse regularization; in particular, in the context of compressed sensing. However, the majority of theoretical error estimates for this regularizer assume an a priori bound on the noise level, which is usually lacking in practice. In this paper, we develop stability and robustness estimates which remove this assumption. First, we introduce an abstract framework and show that robust instance optimality of any decoder in the noise-aware setting implies stability and robustness in the noise-blind setting. This is based on certain sup-inf constants referred to as quotients, strictly related to the quotient property of compressed sensing. We then apply this theory to prove the robustness of quadratically-constrained basis pursuit under unknown error in the cases of random Gaussian matrices and of random matrices with heavy-tailed rows, such as random sampling matrices from bounded orthonormal systems. We illustrate our results in several cases of practical importance, including subsampled Fourier measurements and recovery of sparse polynomial expansions.

Alexander Bastounis, Paolo Campodonico, Mihaela van der Schaar et al.

We introduce the Consistent Reasoning Paradox (CRP). Consistent reasoning, which lies at the core of human intelligence, is the ability to handle tasks that are equivalent, yet described by different sentences ('Tell me the time!' and 'What is the time?'). The CRP asserts that consistent reasoning implies fallibility -- in particular, human-like intelligence in AI necessarily comes with human-like fallibility. Specifically, it states that there are problems, e.g. in basic arithmetic, where any AI that always answers and strives to mimic human intelligence by reasoning consistently will hallucinate (produce wrong, yet plausible answers) infinitely often. The paradox is that there exists a non-consistently reasoning AI (which therefore cannot be on the level of human intelligence) that will be correct on the same set of problems. The CRP also shows that detecting these hallucinations, even in a probabilistic sense, is strictly harder than solving the original problems, and that there are problems that an AI may answer correctly, but it cannot provide a correct logical explanation for how it arrived at the answer. Therefore, the CRP implies that any trustworthy AI (i.e., an AI that never answers incorrectly) that also reasons consistently must be able to say 'I don't know'. Moreover, this can only be done by implicitly computing a new concept that we introduce, termed the 'I don't know' function -- something currently lacking in modern AI. In view of these insights, the CRP also provides a glimpse into the behaviour of Artificial General Intelligence (AGI). An AGI cannot be 'almost sure', nor can it always explain itself, and therefore to be trustworthy it must be able to say 'I don't know'.

Ben Adcock, Anyi Bao, John D. Jakeman et al.

The recovery of approximately sparse or compressible coefficients in a Polynomial Chaos Expansion is a common goal in modern parametric uncertainty quantification (UQ). However, relatively little effort in UQ has been directed toward theoretical and computational strategies for addressing the sparse corruptions problem, where a small number of measurements are highly corrupted. Such a situation has become pertinent today since modern computational frameworks are sufficiently complex with many interdependent components that may introduce hardware and software failures, some of which can be difficult to detect and result in a highly polluted simulation result. In this paper we present a novel compressive sampling-based theoretical analysis for a regularized $\ell^1$ minimization algorithm that aims to recover sparse expansion coefficients in the presence of measurement corruptions. Our recovery results are uniform, and prescribe algorithmic regularization parameters in terms of a user-defined a priori estimate on the ratio of measurements that are believed to be corrupted. We also propose an iteratively reweighted optimization algorithm that automatically refines the value of the regularization parameter, and empirically produces superior results. Our numerical results test our framework on several medium-to-high dimensional examples of solutions to parameterized differential equations, and demonstrate the effectiveness of our approach.

Ben Adcock

We consider the problem of approximating a smooth function from finitely-many pointwise samples using $\ell^1$ minimization techniques. In the first part of this paper, we introduce an infinite-dimensional approach to this problem. Three advantages of this approach are as follows. First, it provides interpolatory approximations in the absence of noise. Second, it does not require \textit{a priori} bounds on the expansion tail in order to be implemented. In particular, the truncation strategy we introduce as part of this framework is independent of the function being approximated, provided the function has sufficient regularity. Third, it allows one to explain the key role weights play in the minimization; namely, that of regularizing the problem and removing aliasing phenomena. In the second part of this paper we present a worst-case error analysis for this approach. We provide a general recipe for analyzing this technique for arbitrary deterministic sets of points. Finally, we use this tool to show that weighted $\ell^1$ minimization with Jacobi polynomials leads to an optimal method for approximating smooth, one-dimensional functions from scattered data.

Ben Adcock, Yi Sui

We consider the sparse polynomial approximation of a multivariate function on a tensor product domain from samples of both the function and its gradient. When only function samples are prescribed, weighted $\ell^1$ minimization has recently been shown to be an effective procedure for computing such approximations. We extend this work to the gradient-augmented case. Our main results show that for the same asymptotic sample complexity, gradient-augmented measurements achieve an approximation error bound in a stronger Sobolev norm, as opposed to the $L^2$-norm in the unaugmented case. For Chebyshev and Legendre polynomial approximations, this sample complexity estimate is algebraic in the sparsity $s$ and at most logarithmic in the dimension $d$, thus mitigating the curse of dimensionality to a substantial extent. We also present several experiments numerically illustrating the benefits of gradient information over an equivalent number of function samples only.

Ben Adcock, Rodrigo Platte, Alexei Shadrin

We consider the problem of approximating an analytic function on a compact interval from its values at $M+1$ distinct points. When the points are equispaced, a recent result (the so-called impossibility theorem) has shown that the best possible convergence rate of a stable method is root-exponential in $M$, and that any method with faster exponential convergence must also be exponentially ill-conditioned at a certain rate. This result hinges on a classical theorem of Coppersmith & Rivlin concerning the maximal behaviour of polynomials bounded on an equispaced grid. In this paper, we first generalize this theorem to arbitrary point distributions. We then present an extension of the impossibility theorem valid for general nonequispaced points, and apply it to the case of points that are equidistributed with respect to (modified) Jacobi weight functions. This leads to a necessary sampling rate for stable approximation from such points. We prove that this rate is also sufficient, and therefore exactly quantify (up to constants) the precise sampling rate for approximating analytic functions from such node distributions with stable methods. Numerical results -- based on computing the maximal polynomial via a variant of the classical Remez algorithm -- confirm our main theorems. Finally, we discuss the implications of our results for polynomial least-squares approximations. In particular, we theoretically confirm the well-known heuristic that stable least-squares approximation using polynomials of degree $N < M$ is possible only once $M$ is sufficiently large for there to be a subset of $N$ of the nodes that mimic the behaviour of the $N$th set of Chebyshev nodes.

Ben Adcock, Simone Brugiapaglia, Nick Dexter et al.

Sparse polynomial approximation has become indispensable for approximating smooth, high- or infinite-dimensional functions from limited samples. This is a key task in computational science and engineering, e.g., surrogate modelling in uncertainty quantification where the function is the solution map of a parametric or stochastic differential equation (DE). Yet, sparse polynomial approximation lacks a complete theory. On the one hand, there is a well-developed theory of best $s$-term polynomial approximation, which asserts exponential or algebraic rates of convergence for holomorphic functions. On the other, there are increasingly mature methods such as (weighted) $\ell^1$-minimization for computing such approximations. While the sample complexity of these methods has been analyzed with compressed sensing, whether they achieve best $s$-term approximation rates is not fully understood. Furthermore, these methods are not algorithms per se, as they involve exact minimizers of nonlinear optimization problems. This paper closes these gaps. Specifically, we consider the following question: are there robust, efficient algorithms for computing approximations to finite- or infinite-dimensional, holomorphic and Hilbert-valued functions from limited samples that achieve best $s$-term rates? We answer this affirmatively by introducing algorithms and theoretical guarantees that assert exponential or algebraic rates of convergence, along with robustness to sampling, algorithmic, and physical discretization errors. We tackle both scalar- and Hilbert-valued functions, this being key to parametric or stochastic DEs. Our results involve significant developments of existing techniques, including a novel restarted primal-dual iteration for solving weighted $\ell^1$-minimization problems in Hilbert spaces. Our theory is supplemented by numerical experiments demonstrating the efficacy of these algorithms.

Ben Adcock, Juan M. Cardenas, Nick Dexter

The problem of approximating smooth, multivariate functions from sample points arises in many applications in scientific computing, e.g., in computational Uncertainty Quantification (UQ) for science and engineering. In these applications, the target function may represent a desired quantity of interest of a parameterized Partial Differential Equation (PDE). Due to the large cost of solving such problems, where each sample is computed by solving a PDE, sample efficiency is a key concerning these applications. Recently, there has been increasing focus on the use of Deep Neural Networks (DNN) and Deep Learning (DL) for learning such functions from data. In this work, we propose an adaptive sampling strategy, CAS4DL (Christoffel Adaptive Sampling for Deep Learning) to increase the sample efficiency of DL for multivariate function approximation. Our novel approach is based on interpreting the second to last layer of a DNN as a dictionary of functions defined by the nodes on that layer. With this viewpoint, we then define an adaptive sampling strategy motivated by adaptive sampling schemes recently proposed for linear approximation schemes, wherein samples are drawn randomly with respect to the Christoffel function of the subspace spanned by this dictionary. We present numerical experiments comparing CAS4DL with standard Monte Carlo (MC) sampling. Our results demonstrate that CAS4DL often yields substantial savings in the number of samples required to achieve a given accuracy, particularly in the case of smooth activation functions, and it shows a better stability in comparison to MC. These results therefore are a promising step towards fully adapting DL towards scientific computing applications.

Maksym Neyra-Nesterenko, Ben Adcock

Solving inverse problems is a fundamental component of science, engineering and mathematics. With the advent of deep learning, deep neural networks have significant potential to outperform existing state-of-the-art, model-based methods for solving inverse problems. However, it is known that current data-driven approaches face several key issues, notably hallucinations, instabilities and unpredictable generalization, with potential impact in critical tasks such as medical imaging. This raises the key question of whether or not one can construct deep neural networks for inverse problems with explicit stability and accuracy guarantees. In this work, we present a novel construction of accurate, stable and efficient neural networks for inverse problems with general analysis-sparse models, termed NESTANets. To construct the network, we first unroll NESTA, an accelerated first-order method for convex optimization. The slow convergence of this method leads to deep networks with low efficiency. Therefore, to obtain shallow, and consequently more efficient, networks we combine NESTA with a novel restart scheme. We then use compressed sensing techniques to demonstrate accuracy and stability. We showcase this approach in the case of Fourier imaging, and verify its stability and performance via a series of numerical experiments. The key impact of this work is demonstrating the construction of efficient neural networks based on unrolling with guaranteed stability and accuracy.

Ben Adcock, Milana Gataric, José Luis Romero

We study the problem of recovering an unknown compactly-supported multivariate function from samples of its Fourier transform that are acquired nonuniformly, i.e. not necessarily on a uniform Cartesian grid. Reconstruction problems of this kind arise in various imaging applications, where Fourier samples are taken along radial lines or spirals for example. Specifically, we consider finite-dimensional reconstructions, where a limited number of samples is available, and investigate the rate of convergence of such approximate solutions and their numerical stability. We show that the proportion of Fourier samples that allow for stable approximations of a given numerical accuracy is independent of the specific sampling geometry and is therefore universal for different sampling scenarios. This allows us to relate both sufficient and necessary conditions for different sampling setups and to exploit several results that were previously available only for very specific sampling geometries. The results are obtained by developing: (i) a transference argument for different measures of the concentration of the Fourier transform and Fourier samples; (ii) frame bounds valid up to the critical sampling density, which depend explicitly on the sampling set and the spectrum. As an application, we identify sufficient and necessary conditions for stable and accurate reconstruction of algebraic polynomials or wavelet coefficients from nonuniform Fourier data.

Ben Adcock, Juan M. Cardenas, Nick Dexter

We introduce a general framework for active learning in regression problems. Our framework extends the standard setup by allowing for general types of data, rather than merely pointwise samples of the target function. This generalization covers many cases of practical interest, such as data acquired in transform domains (e.g., Fourier data), vector-valued data (e.g., gradient-augmented data), data acquired along continuous curves, and, multimodal data (i.e., combinations of different types of measurements). Our framework considers random sampling according to a finite number of sampling measures and arbitrary nonlinear approximation spaces (model classes). We introduce the concept of generalized Christoffel functions and show how these can be used to optimize the sampling measures. We prove that this leads to near-optimal sample complexity in various important cases. This paper focuses on applications in scientific computing, where active learning is often desirable, since it is usually expensive to generate data. We demonstrate the efficacy of our framework for gradient-augmented learning with polynomials, Magnetic Resonance Imaging (MRI) using generative models and adaptive sampling for solving PDEs using Physics-Informed Neural Networks (PINNs).

Ben Adcock, Milana Gataric, Anders C. Hansen

In this paper, we consider the problem of reconstructing piecewise smooth functions to high accuracy from nonuniform samples of their Fourier transform. We use the framework of nonuniform generalized sampling (NUGS) to do this, and to ensure high accuracy we employ reconstruction spaces consisting of splines or (piecewise) polynomials. We analyze the relation between the dimension of the reconstruction space and the bandwidth of the nonuniform samples, and show that it is linear for splines and piecewise polynomials of fixed degree, and quadratic for piecewise polynomials of varying degree.

Ben Adcock, Matthew J. Colbrook, Maksym Neyra-Nesterenko

Sharpness is an almost generic assumption in continuous optimization that bounds the distance from minima by objective function suboptimality. It facilitates the acceleration of first-order methods through restarts. However, sharpness involves problem-specific constants that are typically unknown, and restart schemes typically reduce convergence rates. Moreover, these schemes are challenging to apply in the presence of noise or with approximate model classes (e.g., in compressive imaging or learning problems), and they generally assume that the first-order method used produces feasible iterates. We consider the assumption of approximate sharpness, a generalization of sharpness that incorporates an unknown constant perturbation to the objective function error. This constant offers greater robustness (e.g., with respect to noise or relaxation of model classes) for finding approximate minimizers. By employing a new type of search over the unknown constants, we design a restart scheme that applies to general first-order methods and does not require the first-order method to produce feasible iterates. Our scheme maintains the same convergence rate as when the constants are known. The convergence rates we achieve for various first-order methods match the optimal rates or improve on previously established rates for a wide range of problems. We showcase our restart scheme in several examples and highlight potential future applications and developments of our framework and theory.

Simone Brugiapaglia, Ben Adcock, Richard K. Archibald

From a numerical analysis perspective, assessing the robustness of l1-minimization is a fundamental issue in compressed sensing and sparse regularization. Yet, the recovery guarantees available in the literature usually depend on a priori estimates of the noise, which can be very hard to obtain in practice, especially when the noise term also includes unknown discrepancies between the finite model and data. In this work, we study the performance of l1-minimization when these estimates are not available, providing robust recovery guarantees for quadratically constrained basis pursuit and random sampling in bounded orthonormal systems. Several applications of this work are approximation of high-dimensional functions, infinite-dimensional sparse regularization for inverse problems, and fast algorithms for non-Cartesian Magnetic Resonance Imaging.

Ben Adcock, Juan M. Cardenas, Nick Dexter

This work considers the fundamental problem of learning an unknown object from training data using a given model class. We introduce a unified framework that allows for objects in arbitrary Hilbert spaces, general types of (random) linear measurements as training data and general types of nonlinear model classes. We establish a series of learning guarantees for this framework. These guarantees provide explicit relations between the amount of training data and properties of the model class to ensure near-best generalization bounds. In doing so, we also introduce and develop the key notion of the variation of a model class with respect to a distribution of sampling operators. To exhibit the versatility of this framework, we show that it can accommodate many different types of well-known problems of interest. We present examples such as matrix sketching by random sampling, compressed sensing with isotropic vectors, active learning in regression and compressed sensing with generative models. In all cases, we show how known results become straightforward corollaries of our general learning guarantees. For compressed sensing with generative models, we also present a number of generalizations and improvements of recent results. In summary, our work not only introduces a unified way to study learning unknown objects from general types of data, but also establishes a series of general theoretical guarantees which consolidate and improve various known results.

Ben Adcock

Least-squares approximation is one of the most important methods for recovering an unknown function from data. While in many applications the data is fixed, in many others there is substantial freedom to choose where to sample. In this paper, we review recent progress on near-optimal random sampling strategies for (weighted) least-squares approximation in arbitrary linear spaces. We introduce the Christoffel function as a key quantity in the analysis of (weighted) least-squares approximation from random samples, then show how it can be used to construct a random sampling strategy, termed Christoffel sampling, that possesses near-optimal sample complexity: namely, the number of samples scales log-linearly in the dimension of the approximation space $n$. We discuss a series of variations, extensions and further topics, and throughout highlight connections to approximation theory, machine learning, information-based complexity and numerical linear algebra. Finally, motivated by various contemporary applications, we consider a generalization of the classical setting where the samples need not be pointwise samples of a scalar-valued function, and the approximation space need not be linear. We show that, even in this significantly more general setting, suitable generalizations of Christoffel function still determine the sample complexity. Consequently, these can be used to design enhanced, Christoffel sampling strategies in a unified way for general recovery problems. This article is largely self-contained, and intended to be accessible to nonspecialists.

Ben Adcock, Simone Brugiapaglia

This paper concerns the approximation of smooth, high-dimensional functions from limited samples using polynomials. This task lies at the heart of many applications in computational science and engineering - notably, some of those arising from parametric modelling and computational uncertainty quantification. It is common to use Monte Carlo sampling in such applications, so as not to succumb to the curse of dimensionality. However, it is well known that such a strategy is theoretically suboptimal. Specifically, there are many polynomial spaces of dimension $n$ for which the sample complexity scales log-quadratically, i.e., like $c \cdot n^2 \cdot \log(n)$ as $n \rightarrow \infty$. This well-documented phenomenon has led to a concerted effort over the last decade to design improved, and moreover, near-optimal strategies, whose sample complexities scale log-linearly, or even linearly in $n$. In this work we demonstrate that Monte Carlo is actually a perfectly good strategy in high dimensions, despite its apparent suboptimality. We first document this phenomenon empirically via a systematic set of numerical experiments. Next, we present a theoretical analysis that rigorously justifies this fact in the case of holomorphic functions of infinitely-many variables. We show that there is a least-squares approximation based on $m$ Monte Carlo samples whose error decays algebraically fast in $m/\log(m)$, with a rate that is the same as that of the best $n$-term polynomial approximation. This result is non-constructive, since it assumes knowledge of a suitable polynomial subspace in which to perform the approximation. We next present a compressed sensing-based scheme that achieves the same rate, except for a larger polylogarithmic factor. This scheme is practical, and numerically it performs as well as or better than well-known adaptive least-squares schemes.

James Rowbottom, Elizabeth L. Baker, Nick Huang et al.

Score-based diffusion models in infinite-dimensional function spaces provide a mathematically principled framework for modelling function-valued data, offering key advantages such as resolution invariance and the ability to handle irregular discretisations. However, practical implementations have struggled to fully realise these benefits. Existing backbones like Fourier neural operators are often biased towards regular grids and fail to generalise to complex domain topologies. We propose a novel architecture for function-space diffusion models that represents generalised graph convolutional kernels as finite element functions, enabling the model to naturally handle unstructured meshes and complex geometries. We demonstrate the efficacy of our network architecture through a series of unconditional and conditional sampling experiments across diverse geometries, including non-convex and multiply-connected domains. Our results show that the proposed method maintains resolution invariance and achieves high fidelity in capturing functional distributions on non-trivial geometries.

James Rowbottom, Nick Huang, Carola-Bibiane Schönlieb et al.

State estimation is a critical task in scientific, engineering and control applications. Since the reliability of reconstructions depends on the number and position of sensors, optimal sensor placement (OSP) is essential in scenarios where measurements are sparse and expensive. Classical OSP approaches rely on Gaussian assumptions and are consequently unable to account for the complex distributions encountered in many real-world systems. Generative-model-based reconstruction using sensor guided diffusion posterior sampling (DPS) has emerged as a promising technique for reconstructing states from highly complex distributions. However, existing sensor-selection methods either require unrealistically many sensors or emulate classical OSP, creating a mismatch between modern recovery models with classical OSP tools motivating the need for fundamentally new ideas towards OSP that match the recent advances made in powerful recovery models. We introduce a distribution-free sensor placement framework based on the Christoffel function: a mathematical formulation of optimal sampling and recovery guarantees for posterior sampling with arbitrary sensors and signal distributions, from which we derive a new OSP strategy with non-asymptotic bounds on the number of sensors needed for recovery. We develop Christoffel-DPS, with offline and online variants, instantiating Christoffel sampling for generative models. Christoffel-DPS outperforms Gaussian OSP baselines and existing generative-model placement methods, validating that distribution-free sensing is both theoretically principled and practically superior. The framework is model-agnostic; we demonstrate its application to a range of unconditional DPS and flow-matching models on structurally non-Gaussian benchmarks, showing the efficacy of Christoffel-DPS in low sensor budget regimes.

Ben Adcock, Khiem Can, Xuemeng Wang

Random Feature Models (RFMs) have become a powerful tool for approximating multivariate functions and solving partial differential equations efficiently. Sparse Random Feature Expansions (SRFE) improve traditional RFMs by incorporating sparsity, making it particularly effective in data-scarce settings. In this work, we integrate active learning with sparse random feature approximations to improve sampling efficiency. Specifically, we incorporate the Christoffel function to guide an adaptive sampling process, dynamically selecting informative sample points based on their contribution to the function space. This approach optimizes the distribution of sample points by leveraging the Christoffel function associated with an iteratively-chosen basis obtained by the sparse recovery solver. We conduct numerical experiments comparing adaptive and nonadaptive sampling strategies with the SRFE framework and examine their accuracy for various function approximation tasks. Overall, our results demonstrate the advantages of adaptive sampling in maintaining high accuracy while reducing sample complexity for SRFE, highlighting its potential for scientific computing tasks where data is expensive to acquire.

Ben Adcock, Simone Brugiapaglia, Nick Dexter et al.

Learning approximations to smooth target functions of many variables from finite sets of pointwise samples is an important task in scientific computing and its many applications in computational science and engineering. Despite well over half a century of research on high-dimensional approximation, this remains a challenging problem. Yet, significant advances have been made in the last decade towards efficient methods for doing this, commencing with so-called sparse polynomial approximation methods and continuing most recently with methods based on Deep Neural Networks (DNNs). In tandem, there have been substantial advances in the relevant approximation theory and analysis of these techniques. In this work, we survey this recent progress. We describe the contemporary motivations for this problem, which stem from parametric models and computational uncertainty quantification; the relevant function classes, namely, classes of infinite-dimensional, Banach-valued, holomorphic functions; fundamental limits of learnability from finite data for these classes; and finally, sparse polynomial and DNN methods for efficiently learning such functions from finite data. For the latter, there is currently a significant gap between the approximation theory of DNNs and the practical performance of deep learning. Aiming to narrow this gap, we develop the topic of practical existence theory, which asserts the existence of dimension-independent DNN architectures and training strategies that achieve provably near-optimal generalization errors in terms of the amount of training data.

Ben Adcock, Avi Gupta

A key problem in approximation theory is the recovery of high-dimensional functions from samples. In many cases, the functions of interest exhibit anisotropic smoothness, and, in many practical settings, the nature of this anisotropy may be unknown a priori. Therefore, an important question involves the development of universal algorithms, namely, algorithms that simultaneously achieve optimal or near-optimal rates of convergence across a range of different anisotropic smoothness classes. In this work, we consider universal approximation of periodic functions that belong to anisotropic Sobolev spaces and anisotropic dominating mixed smoothness Sobolev spaces. Our first result is the construction of a universal algorithm. This recasts function recovery as a sparse recovery problem for Fourier coefficients and then exploits compressed sensing to yield the desired approximation rates. Note that this algorithm is nonadaptive, as it does not seek to learn the anisotropic smoothness of the target function. We then demonstrate optimality of this algorithm up to a dimension-independent polylogarithmic factor. We do this by presenting a lower bound for the adaptive $m$-width for the unit balls of such function classes. Finally, we demonstrate the necessity of nonlinear algorithms. We show that universal linear algorithms can achieve rates that are at best suboptimal by a dimension-dependent polylogarithmic factor. In other words, they suffer from a curse of dimensionality in the rate -- a phenomenon which justifies the necessity of nonlinear algorithms for universal recovery.

Ben Adcock, Bernhard Hientzsch, Akil Narayan et al.

Motivated by the need for efficient estimation of conditional expectations, we consider a least-squares function approximation problem with heavily polluted data. Existing methods that are powerful in the small noise regime are suboptimal when large noise is present. We propose a hybrid approach that combines Christoffel sampling with certain types of optimal experimental design to address this issue. We show that the proposed algorithm enjoys appropriate optimality properties for both sample point generation and noise mollification, leading to improved computational efficiency and sample complexity compared to existing methods. We also extend the algorithm to convex-constrained settings with similar theoretical guarantees. When the target function is defined as the expectation of a random field, we extend our approach to leverage adaptive random subspaces and establish results on the approximation capacity of the adaptive procedure. Our theoretical findings are supported by numerical studies on both synthetic data and on a more challenging stochastic simulation problem in computational finance.

Ben Adcock, Michael Griebel, Gregor Maier

Operator learning, the approximation of mappings between infinite-dimensional function spaces using machine learning, has gained increasing research attention in recent years. Approximate operators, learned from data, can serve as efficient surrogate models for problems in computational science and engineering, complementing traditional methods. However, despite their empirical success, our understanding of the underlying mathematical theory is in large part still incomplete. In this paper, we study the approximation of Lipschitz operators with respect to Gaussian measures. We prove higher Gaussian Sobolev regularity of Lipschitz operators and establish lower and upper bounds on the Hermite polynomial approximation error. We then study general reconstruction strategies of Lipschitz operators from $m$ arbitrary (potentially adaptive) linear samples. As a key finding, we tightly characterize the corresponding sample complexity, that is, the smallest achievable worst-case error among all possible choices of (adaptive) sampling and reconstruction strategies in terms of $m$. As a consequence, we identify an inherent curse of sample complexity: No method to approximate Lipschitz operators based on $m$ linear samples can achieve algebraic convergence rates in $m$. On the positive side, we prove that a sufficiently fast spectral decay of the covariance operator of the underlying Gaussian measure guarantees convergence rates which are arbitrarily close to any algebraic rate. Overall, by tightly characterizing the sample complexity, our work confirms the intrinsic difficulty of learning Lipschitz operators, regardless of the data or learning technique.

Ben Adcock, Nick Huang

We study the sample complexity of Bayesian recovery for solving inverse problems with general prior, forward operator and noise distributions. We consider posterior sampling according to an approximate prior $\mathcal{P}$, and establish sufficient conditions for stable and accurate recovery with high probability. Our main result is a non-asymptotic bound that shows that the sample complexity depends on (i) the intrinsic complexity of $\mathcal{P}$, quantified by its so-called approximate covering number, and (ii) concentration bounds for the forward operator and noise distributions. As a key application, we specialize to generative priors, where $\mathcal{P}$ is the pushforward of a latent distribution via a Deep Neural Network (DNN). We show that the sample complexity scales log-linearly with the latent dimension $k$, thus establishing the efficacy of DNN-based priors. Generalizing existing results on deterministic (i.e., non-Bayesian) recovery for the important problem of random sampling with an orthogonal matrix $U$, we show how the sample complexity is determined by the coherence of $U$ with respect to the support of $\mathcal{P}$. Hence, we establish that coherence plays a fundamental role in Bayesian recovery as well. Overall, our framework unifies and extends prior work, providing rigorous guarantees for the sample complexity of solving Bayesian inverse problems with arbitrary distributions.

Ben Adcock, Nick Dexter, Sebastian Moraga

Operator learning problems arise in many key areas of scientific computing where Partial Differential Equations (PDEs) are used to model physical systems. In such scenarios, the operators map between Banach or Hilbert spaces. In this work, we tackle the problem of learning operators between Banach spaces, in contrast to the vast majority of past works considering only Hilbert spaces. We focus on learning holomorphic operators - an important class of problems with many applications. We combine arbitrary approximate encoders and decoders with standard feedforward Deep Neural Network (DNN) architectures - specifically, those with constant width exceeding the depth - under standard $\ell^2$-loss minimization. We first identify a family of DNNs such that the resulting Deep Learning (DL) procedure achieves optimal generalization bounds for such operators. For standard fully-connected architectures, we then show that there are uncountably many minimizers of the training problem that yield equivalent optimal performance. The DNN architectures we consider are `problem agnostic', with width and depth only depending on the amount of training data $m$ and not on regularity assumptions of the target operator. Next, we show that DL is optimal for this problem: no recovery procedure can surpass these generalization bounds up to log terms. Finally, we present numerical results demonstrating the practical performance on challenging problems including the parametric diffusion, Navier-Stokes-Brinkman and Boussinesq PDEs.

Ben Adcock, Simone Brugiapaglia, Nick Dexter et al.

Accurate approximation of scalar-valued functions from sample points is a key task in computational science. Recently, machine learning with Deep Neural Networks (DNNs) has emerged as a promising tool for scientific computing, with impressive results achieved on problems where the dimension of the data or problem domain is large. This work broadens this perspective, focusing on approximating functions that are Hilbert-valued, i.e. take values in a separable, but typically infinite-dimensional, Hilbert space. This arises in science and engineering problems, in particular those involving solution of parametric Partial Differential Equations (PDEs). Such problems are challenging: 1) pointwise samples are expensive to acquire, 2) the function domain is high dimensional, and 3) the range lies in a Hilbert space. Our contributions are twofold. First, we present a novel result on DNN training for holomorphic functions with so-called hidden anisotropy. This result introduces a DNN training procedure and full theoretical analysis with explicit guarantees on error and sample complexity. The error bound is explicit in three key errors occurring in the approximation procedure: the best approximation, measurement, and physical discretization errors. Our result shows that there exists a procedure (albeit non-standard) for learning Hilbert-valued functions via DNNs that performs as well as, but no better than current best-in-class schemes. It gives a benchmark lower bound for how well DNNs can perform on such problems. Second, we examine whether better performance can be achieved in practice through different types of architectures and training. We provide preliminary numerical results illustrating practical performance of DNNs on parametric PDEs. We consider different parameters, modifying the DNN architecture to achieve better and competitive results, comparing these to current best-in-class schemes.

Ben Adcock, Nick Dexter

Deep learning (DL) is transforming industry as decision-making processes are being automated by deep neural networks (DNNs) trained on real-world data. Driven partly by rapidly-expanding literature on DNN approximation theory showing they can approximate a rich variety of functions, such tools are increasingly being considered for problems in scientific computing. Yet, unlike traditional algorithms in this field, little is known about DNNs from the principles of numerical analysis, e.g., stability, accuracy, computational efficiency and sample complexity. In this paper we introduce a computational framework for examining DNNs in practice, and use it to study empirical performance with regard to these issues. We study performance of DNNs of different widths & depths on test functions in various dimensions, including smooth and piecewise smooth functions. We also compare DL against best-in-class methods for smooth function approx. based on compressed sensing (CS). Our main conclusion from these experiments is that there is a crucial gap between the approximation theory of DNNs and their practical performance, with trained DNNs performing relatively poorly on functions for which there are strong approximation results (e.g. smooth functions), yet performing well in comparison to best-in-class methods for other functions. To analyze this gap further, we provide some theoretical insights. We establish a practical existence theorem, asserting existence of a DNN architecture and training procedure that offers the same performance as CS. This establishes a key theoretical benchmark, showing the gap can be closed, albeit via a strategy guaranteed to perform as well as, but no better than, current best-in-class schemes. Nevertheless, it demonstrates the promise of practical DNN approx., by highlighting potential for better schemes through careful design of DNN architectures and training strategies.

Nina M. Gottschling, Vegard Antun, Anders C. Hansen et al.

Methods inspired by Artificial Intelligence (AI) are starting to fundamentally change computational science and engineering through breakthrough performances on challenging problems. However, reliability and trustworthiness of such techniques is a major concern. In inverse problems in imaging, the focus of this paper, there is increasing empirical evidence that methods may suffer from hallucinations, i.e., false, but realistic-looking artifacts; instability, i.e., sensitivity to perturbations in the data; and unpredictable generalization, i.e., excellent performance on some images, but significant deterioration on others. This paper provides a theoretical foundation for these phenomena. We give mathematical explanations for how and when such effects arise in arbitrary reconstruction methods, with several of our results taking the form of `no free lunch' theorems. Specifically, we show that (i) methods that overperform on a single image can wrongly transfer details from one image to another, creating a hallucination, (ii) methods that overperform on two or more images can hallucinate or be unstable, (iii) optimizing the accuracy-stability trade-off is generally difficult, (iv) hallucinations and instabilities, if they occur, are not rare events, and may be encouraged by standard training, (v) it may be impossible to construct optimal reconstruction maps for certain problems. Our results trace these effects to the kernel of the forward operator whenever it is nontrivial, but also apply to the case when the forward operator is ill-conditioned. Based on these insights, our work aims to spur research into new ways to develop robust and reliable AI-based methods for inverse problems in imaging.

Il Yong Chun, David Hong, Ben Adcock et al.

Convolutional analysis operator learning (CAOL) enables the unsupervised training of (hierarchical) convolutional sparsifying operators or autoencoders from large datasets. One can use many training images for CAOL, but a precise understanding of the impact of doing so has remained an open question. This paper presents a series of results that lend insight into the impact of dataset size on the filter update in CAOL. The first result is a general deterministic bound on errors in the estimated filters, and is followed by a bound on the expected errors as the number of training samples increases. The second result provides a high probability analogue. The bounds depend on properties of the training data, and we investigate their empirical values with real data. Taken together, these results provide evidence for the potential benefit of using more training data in CAOL.

Vegard Antun, Francesco Renna, Clarice Poon et al.

Deep learning, due to its unprecedented success in tasks such as image classification, has emerged as a new tool in image reconstruction with potential to change the field. In this paper we demonstrate a crucial phenomenon: deep learning typically yields unstablemethods for image reconstruction. The instabilities usually occur in several forms: (1) tiny, almost undetectable perturbations, both in the image and sampling domain, may result in severe artefacts in the reconstruction, (2) a small structural change, for example a tumour, may not be captured in the reconstructed image and (3) (a counterintuitive type of instability) more samples may yield poorer performance. Our new stability test with algorithms and easy to use software detects the instability phenomena. The test is aimed at researchers to test their networks for instabilities and for government agencies, such as the Food and Drug Administration (FDA), to secure safe use of deep learning methods.

Ben Adcock, Anyi Bao, Simone Brugiapaglia

We consider sparsity-based techniques for the approximation of high-dimensional functions from random pointwise evaluations. To date, almost all the works published in this field contain some a priori assumptions about the error corrupting the samples that are hard to verify in practice. In this paper, we instead focus on the scenario where the error is unknown. We study the performance of four sparsity-promoting optimization problems: weighted quadratically-constrained basis pursuit, weighted LASSO, weighted square-root LASSO, and weighted LAD-LASSO. From the theoretical perspective, we prove uniform recovery guarantees for these decoders, deriving recipes for the optimal choice of the respective tuning parameters. On the numerical side, we compare them in the pure function approximation case and in applications to uncertainty quantification of ODEs and PDEs with random inputs. Our main conclusion is that the lesser-known square-root LASSO is better suited for high-dimensional approximation than the other procedures in the case of bounded noise, since it avoids (both theoretically and numerically) the need for parameter tuning.

Ben Adcock, Simone Brugiapaglia

We show the potential of greedy recovery strategies for the sparse approximation of multivariate functions from a small dataset of pointwise evaluations by considering an extension of the orthogonal matching pursuit to the setting of weighted sparsity. The proposed recovery strategy is based on a formal derivation of the greedy index selection rule. Numerical experiments show that the proposed weighted orthogonal matching pursuit algorithm is able to reach accuracy levels similar to those of weighted $\ell^1$ minimization programs while considerably improving the computational efficiency for small values of the sparsity level.

Ben Adcock, Simone Brugiapaglia, Clayton G. Webster

In recent years, the use of sparse recovery techniques in the approximation of high-dimensional functions has garnered increasing interest. In this work we present a survey of recent progress in this emerging topic. Our main focus is on the computation of polynomial approximations of high-dimensional functions on $d$-dimensional hypercubes. We show that smooth, multivariate functions possess expansions in orthogonal polynomial bases that are not only approximately sparse, but possess a particular type of structured sparsity defined by so-called lower sets. This structure can be exploited via the use of weighted $\ell^1$ minimization techniques, and, as we demonstrate, doing so leads to sample complexity estimates that are at most logarithmically dependent on the dimension $d$. Hence the curse of dimensionality - the bane of high-dimensional approximation - is mitigated to a significant extent. We also discuss several practical issues, including unknown noise (due to truncation or numerical error), and highlight a number of open problems and challenges.

Chen Li, Ben Adcock

In compressed sensing, it is often desirable to consider signals possessing additional structure beyond sparsity. One such structured signal model - which forms the focus of this paper - is the local sparsity in levels class. This class has recently found applications in problems such as compressive imaging, multi-sensor acquisition systems and sparse regularization in inverse problems. In this paper we present uniform recovery guarantees for this class when the measurement matrix corresponds to a subsampled isometry. We do this by establishing a variant of the standard restricted isometry property for sparse in levels vectors, known as the restricted isometry property in levels. Interestingly, besides the usual log factors, our uniform recovery guarantees are simpler and less stringent than existing nonuniform recovery guarantees. For the particular case of discrete Fourier sampling with Haar wavelet sparsity, a corollary of our main theorem yields a new recovery guarantee which improves over the current state-of-the-art.

Ben Adcock, Milana Gataric, Anders C. Hansen

We provide sufficient density condition for a set of nonuniform samples to give rise to a set of sampling for multivariate bandlimited functions when the measurements consist of pointwise evaluations of a function and its first $k$ derivatives. Along with explicit estimates of corresponding frame bounds, we derive the explicit density bound and show that, as $k$ increases, it grows linearly in $k+1$ with the constant of proportionality $1/\mathrm{e}$. Seeking larger gap conditions, we also prove a multivariate perturbation result for nonuniform samples that are sufficiently close to sets of sampling, e.g. to uniform samples taken at $k+1$ times the Nyquist rate. Additionally, in the univariate setting, we consider a related problem of so-called nonuniform bunched sampling, where in each sampling interval $s+1$ bunched measurements of a function are taken and the sampling intervals are permitted to be of different length. We derive an explicit density condition which grows linearly in $s+1$ for large $s$, with the constant of proportionality depending on the width of the bunches. The width of the bunches is allowed to be arbitrarily small, and moreover, for sufficiently narrow bunches and sufficiently large $s$, we obtain the same result as in the case of univariate sampling with $s$ derivatives.

Ben Adcock, Jésus Martín-Vaquero, Mark Richardson

We introduce a numerical method for the approximation of functions which are analytic on compact intervals, except at the endpoints. This method is based on variable transforms using particular parametrized exponential and double-exponential mappings, in combination with Fourier-like approximation in a truncated domain. We show theoretically that this method is superior to variable transform techniques based on the standard exponential and double-exponential mappings. In particular, it can resolve oscillatory behaviour using near-optimal degrees of freedom, whereas the standard mappings require degrees of freedom that grow superlinearly with the frequency of oscillation. We highlight these results with several numerical experiments. Therein it is observed that near-machine epsilon accuracy is achieved using a number of degrees of freedom that is between four and ten times smaller than those of existing techniques.

Ben Adcock, Milana Gataric, Anders C. Hansen

In this paper, we consider the problem of recovering a compactly supported multivariate function from a collection of pointwise samples of its Fourier transform taken nonuniformly. We do this by using the concept of weighted Fourier frames. A seminal result of Beurling shows that sample points give rise to a classical Fourier frame provided they are relatively separated and of sufficient density. However, this result does not allow for arbitrary clustering of sample points, as is often the case in practice. Whilst keeping the density condition sharp and dimension independent, our first result removes the separation condition and shows that density alone suffices. However, this result does not lead to estimates for the frame bounds. A known result of Groechenig provides explicit estimates, but only subject to a density condition that deteriorates linearly with dimension. In our second result we improve these bounds by reducing the dimension dependence. In particular, we provide explicit frame bounds which are dimensionless for functions having compact support contained in a sphere. Next, we demonstrate how our two main results give new insight into a reconstruction algorithm---based on the existing generalized sampling framework---that allows for stable and quasi-optimal reconstruction in any particular basis from a finite collection of samples. Finally, we construct sufficiently dense sampling schemes that are often used in practice---jittered, radial and spiral sampling schemes---and provide several examples illustrating the effectiveness of our approach when tested on these schemes.

Alex D. Jones, Ben Adcock, Anders C. Hansen

Recently, it has been shown that incoherence is an unrealistic assumption for compressed sensing when applied to many inverse problems. Instead, the key property that permits efficient recovery in such problems is so-called local incoherence. Similarly, the standard notion of sparsity is also inadequate for many real world problems. In particular, in many applications, the optimal sampling strategy depends on asymptotic incoherence and the signal sparsity structure. The purpose of this paper is to study asymptotic incoherence and its implications towards the design of optimal sampling strategies and efficient sparsity bases. It is determined how fast asymptotic incoherence can decay in general for isometries. Furthermore it is shown that Fourier sampling and wavelet sparsity, whilst globally coherent, yield optimal asymptotic incoherence as a power law up to a constant factor. Sharp bounds on the asymptotic incoherence for Fourier sampling with polynomial bases are also provided. A numerical experiment is also presented to demonstrate the role of asymptotic incoherence in finding good subsampling strategies.