The remarkable practical success of deep learning has revealed some major surprises from a theoretical perspective. In particular, simple gradient methods easily find near-optimal solutions to non-convex optimization problems, and despite giving a near-perfect fit to training data without any explicit effort to control model complexity, these methods exhibit excellent predictive accuracy. We conjecture that specific principles underlie these phenomena: that overparametrization allows gradient methods to find interpolating solutions, that these methods implicitly impose regularization, and that overparametrization leads to benign overfitting, that is, accurate predictions despite overfitting training data. In this article, we survey recent progress in statistical learning theory that provides examples illustrating these principles in simpler settings. We first review classical uniform convergence results and why they fall short of explaining aspects of the behaviour of deep learning methods. We give examples of implicit regularization in simple settings, where gradient methods lead to minimal norm functions that perfectly fit the training data. Then we review prediction methods that exhibit benign overfitting, focusing on regression problems with quadratic loss. For these methods, we can decompose the prediction rule into a simple component that is useful for prediction and a spiky component that is useful for overfitting but, in a favourable setting, does not harm prediction accuracy. We focus specifically on the linear regime for neural networks, where the network can be approximated by a linear model. In this regime, we demonstrate the success of gradient flow, and we consider benign overfitting with two-layer networks, giving an exact asymptotic analysis that precisely demonstrates the impact of overparametrization. We conclude by highlighting the key challenges that arise in extending these insights to realistic deep learning settings.

Over the past decade considerable progress has been made towards the numerical solution of large-scale eigenvalue problems, particularly for nonsymmetric matrices. Krylov methods and variants of subspace iteration have been improved to the point that problems of the order of several million variables can be solved. The methods and software that have led to these advances are surveyed.

Atomistic-to-continuum (a/c) coupling methods are a class of computational multiscale schemes that combine the accuracy of atomistic models with the efficiency of continuum elasticity. They are increasingly being utilized in materials science to study the fundamental mechanisms of material failure such as crack propagation and plasticity, which are governed by the interaction between crystal defects and long-range elastic fields.

In the construction of a/c coupling methods, various approximation errors are committed. A rigorous numerical analysis approach that classifies and quantifies these errors can give confidence in the simulation results, as well as enable optimization of the numerical methods for accuracy and computational cost. In this article, we present such a numerical analysis framework, which is inspired by recent research activity.

Any logical procedure that is used to reason or to infer either deductively or inductively, so as to draw conclusions or make decisions, can be called, in a broad sense, a realization process. A realization process usually assumes the recursive form that one state develops into another state by following a certain specific rule. Such an action is generally formalized as a dynamical system. In mathematics, especially for existence questions, a realization process often appears in the form of an iterative procedure or a differential equation. For years researchers have taken great effort to describe, analyse, and modify realization processes for various applications.

The thrust in this exposition is to exploit the notion of dynamical systems as a special realization process for problems arising from the field of linear algebra. Several differential equations whose solutions evolve in submanifolds of matrices are cast in fairly general frameworks, of which special cases have been found to afford unified and fundamental insights into the structure and behaviour of existing discrete methods and, now and then, suggest new and improved numerical methods. In some cases, there are remarkable connections between smooth flows and discrete numerical algorithms. In other cases, the flow approach seems advantageous in tackling very difficult open problems. Various aspects of the recent development and application in this direction are discussed in this paper.

Differential algebraic equations (DAE) are special implicit ordinary differential equations (ODE)

where the partial Jacobian f′y(y, x, t) is singular for all values of its arguments.

This article addresses the inference of physics models from data, from the perspectives of inverse problems and model reduction. These fields develop formulations that integrate data into physics-based models while exploiting the fact that many mathematical models of natural and engineered systems exhibit an intrinsically low-dimensional solution manifold. In inverse problems, we seek to infer uncertain components of the inputs from observations of the outputs, while in model reduction we seek low-dimensional models that explicitly capture the salient features of the input–output map through approximation in a low-dimensional subspace. In both cases, the result is a predictive model that reflects data-driven learning yet deeply embeds the underlying physics, and thus can be used for design, control and decision-making, often with quantified uncertainties. We highlight recent developments in scalable and efficient algorithms for inverse problems and model reduction governed by large-scale models in the form of partial differential equations. Several illustrative applications to large-scale complex problems across different domains of science and engineering are provided.

The need to evaluate a function f(A) ∈ ℂn×n of a matrix A ∈ ℂn×n arises in a wide and growing number of applications, ranging from the numerical solution of differential equations to measures of the complexity of networks. We give a survey of numerical methods for evaluating matrix functions, along with a brief treatment of the underlying theory and a description of two recent applications. The survey is organized by classes of methods, which are broadly those based on similarity transformations, those employing approximation by polynomial or rational functions, and matrix iterations. Computation of the Fréchet derivative, which is important for condition number estimation, is also treated, along with the problem of computing f(A)b without computing f(A). A summary of available software completes the survey.

We are concerned here with processing discontinuous functions from their spectral information. We focus on two main aspects of processing such piecewise smooth data: detecting the edges of a piecewise smooth f, namely, the location and amplitudes of its discontinuities; and recovering with high accuracy the underlying function in between those edges. If f is a smooth function, say analytic, then classical Fourier projections recover f with exponential accuracy. However, if f contains one or more discontinuities, its global Fourier projections produce spurious Gibbs oscillations which spread throughout the smooth regions, enforcing local loss of resolution and global loss of accuracy. Our aim in the computation of the Gibbs phenomenon is to detect edges and to reconstruct piecewise smooth functions, while regaining the high accuracy encoded in the spectral data.

To detect edges, we utilize a general family of edge detectors based on concentration kernels. Each kernel forms an approximate derivative of the delta function, which detects edges by separation of scales. We show how such kernels can be adapted to detect edges with one- and two-dimensional discrete data, with noisy data, and with incomplete spectral information. The main feature is concentration kernels which enable us to convert global spectral moments into local information in physical space. To reconstruct f with high accuracy we discuss novel families of mollifiers and filters. The main feature here is making these mollifiers and filters adapted to the local region of smoothness while increasing their accuracy together with the dimension of the data. These mollifiers and filters form approximate delta functions which are properly parametrized to recover f with (root-) exponential accuracy.

The traditional metric for the efficiency of a numerical algorithm has been the number of arithmetic operations it performs. Technological trends have long been reducing the time to perform an arithmetic operation, so it is no longer the bottleneck in many algorithms; rather, communication, or moving data, is the bottleneck. This motivates us to seek algorithms that move as little data as possible, either between levels of a memory hierarchy or between parallel processors over a network. In this paper we summarize recent progress in three aspects of this problem. First we describe lower bounds on communication. Some of these generalize known lower bounds for dense classical (O(n3)) matrix multiplication to all direct methods of linear algebra, to sequential and parallel algorithms, and to dense and sparse matrices. We also present lower bounds for Strassen-like algorithms, and for iterative methods, in particular Krylov subspace methods applied to sparse matrices. Second, we compare these lower bounds to widely used versions of these algorithms, and note that these widely used algorithms usually communicate asymptotically more than is necessary. Third, we identify or invent new algorithms for most linear algebra problems that do attain these lower bounds, and demonstrate large speed-ups in theory and practice.

Homogenization is an important mathematical framework for developing effective models of differential equations with oscillations. We include in the presentation techniques for deriving effective equations, a brief discussion on analysis of related limit processes and numerical methods that are based on homogenization principles. We concentrate on first- and second-order partial differential equations and present results concerning both periodic and random media for linear as well as nonlinear problems. In the numerical sections, we comment on computations of multi-scale problems in general and then focus on projection-based numerical homogenization and the heterogeneous multi-scale method.

A number of fundamental results in modern statistical theory involve thresholding estimators. This survey paper aims at reconstructing the history of how thresholding rules came to be popular in statistics and describing, in a not overly technical way, the domain of their application. Two notions play a fundamental role in our narrative: sparsity and oracle inequalities. Sparsity is a property of the object to estimate, which seems to be characteristic of many modern problems, in statistics as well as applied mathematics and theoretical computer science, to name a few. ‘Oracle inequalities’ are a powerful decision-theoretic tool which has served to understand the optimality of thresholding rules, but which has many other potential applications, some of which we will discuss.

Our story is also the story of the dialogue between statistics and applied harmonic analysis. Starting with the work of Wiener, we will see that certain representations emerge as being optimal for estimation. A leitmotif throughout our exposition is that efficient representations lead to efficient estimation.

We give examples of problem areas in interpolation, approximation, and quadrature, that call for orthogonal polynomials not of the classical kind. We then discuss numerical methods of computing the respective Gauss-type quadrature rules and orthogonal polynomials. The basic task is to compute the coefficients in the three-term recurrence relation for the orthogonal polynomials. This can be done by methods relying either on moment information or on discretization procedures. The effect on the recurrence coefficients of multiplying the weight function by a rational function is also discussed. Similar methods are applicable to computing Sobolev orthogonal polynomials, although their recurrence relations are more complicated. The paper concludes with a brief account of available software.

Shallow-water equations are widely used to model water flow in rivers, lakes, reservoirs, coastal areas, and other situations in which the water depth is much smaller than the horizontal length scale of motion. The classical shallow-water equations, the Saint-Venant system, were originally proposed about 150 years ago and still are used in a variety of applications. For many practical purposes, it is extremely important to have an accurate, efficient and robust numerical solver for the Saint-Venant system and related models. As their solutions are typically non-smooth and even discontinuous, finite-volume schemes are among the most popular tools. In this paper, we review such schemes and focus on one of the simplest (yet highly accurate and robust) methods: central-upwind schemes. These schemes belong to the family of Godunov-type Riemann-problem-solver-free central schemes, but incorporate some upwinding information about the local speeds of propagation, which helps to reduce an excessive amount of numerical diffusion typically present in classical (staggered) non-oscillatory central schemes. Besides the classical one- and two-dimensional Saint-Venant systems, we will consider the shallow-water equations with friction terms, models with moving bottom topography, the two-layer shallow-water system as well as general non-conservative hyperbolic systems.

Complexity theory of numerical analysis is the study of the number of arithmetic operations required to pass from the input to the output of a numerical problem.

In this article we shall present a unified and axiomatized view of several theories and algorithms of image multiscale analysis (and low level vision) which have been developed in the past twenty years. We shall show that under reasonable invariance and assumptions, all image (and shape) analyses can be reduced to a single partial differential equation. In the same way, movie analysis leads to a single parabolic differential equation. We discuss some applications to image segmentation and movie restoration. The experiments show how accurate and invariant the numerical schemes must be and we compare several (old and new) algorithms by discussing how well they match the axiomatic invariance requirements.

The Delaunay triangulation of a finite point set is a central theme in computational geometry. It finds its major application in the generation of meshes used in the simulation of physical processes. This paper connects the predominantly combinatorial work in classical computational geometry with the numerical interest in mesh generation. It focuses on the two- and three-dimensional case and covers results obtained during the twentieth century.

A standard paradigm for the existence of solutions in fluid dynamics is based on the construction of sequences of approximate solutions or approximate minimizers. This approach faces serious obstacles, most notably in multi-dimensional problems, where the persistence of oscillations at ever finer scales prevents compactness. Indeed, these oscillations are an indication, consistent with recent theoretical results, of the possible lack of existence/uniqueness of solutions within the standard framework of integrable functions. It is in this context that Young measures – parametrized probability measures which can describe the limits of such oscillatory sequences – offer the more general paradigm of measure-valued solutions for these problems.

We present viable numerical algorithms to compute approximate measure-valued solutions, based on the realization of approximate measures as laws of Monte Carlo sampled random fields. We prove convergence of these algorithms to measure-valued solutions for the equations of compressible and incompressible inviscid fluid dynamics, and present a large number of numerical experiments which provide convincing evidence for the viability of the new paradigm. We also discuss the use of these algorithms, and their extensions, in uncertainty quantification and contexts other than fluid dynamics, such as non-convex variational problems in materials science.

In the past decade the mathematical theory of machine learning has lagged far behind the triumphs of deep neural networks on practical challenges. However, the gap between theory and practice is gradually starting to close. In this paper I will attempt to assemble some pieces of the remarkable and still incomplete mathematical mosaic emerging from the efforts to understand the foundations of deep learning. The two key themes will be interpolation and its sibling over-parametrization. Interpolation corresponds to fitting data, even noisy data, exactly. Over-parametrization enables interpolation and provides flexibility to select a suitable interpolating model.

As we will see, just as a physical prism separates colours mixed within a ray of light, the figurative prism of interpolation helps to disentangle generalization and optimization properties within the complex picture of modern machine learning. This article is written in the belief and hope that clearer understanding of these issues will bring us a step closer towards a general theory of deep learning and machine learning.

The development of Krylov subspace methods for the solution of operator equations has shown that two basic construction principles underlie the most commonly used algorithms: the orthogonal residual (OR) and minimal residual (MR) approaches. It is shown that these can both be formulated as techniques for solving an approximation problem on a sequence of nested subspaces of a Hilbert space, an abstract problem not necessarily related to an operator equation. Essentially all Krylov subspace algorithms result when these subspaces form a Krylov sequence. The well-known relations among the iterates and residuals of MR/OR pairs are shown to hold also in this rather general setting. We further show that a common error analysis for these methods involving the canonical angles between subspaces allows many of the known residual and error bounds to be derived in a simple and consistent manner. An application of this analysis to compact perturbations of the identity shows that MR/OR pairs of Krylov subspace methods converge q-superlinearly when applied to such operator equations.

The subject of ‘wavelets’ is expanding at such a tremendous rate that it is impossible to give, within these few pages, a complete introduction to all aspects of its theory. We hope, however, to allow the reader to become sufficiently acquainted with the subject to understand, in part, the enthusiasm of its proponents toward its potential application to various numerical problems. Furthermore, we hope that our exposition can guide the reader who wishes to make more serious excursions into the subject. Our viewpoint is biased by our experience in approximation theory and data compression; we warn the reader that there are other viewpoints that are either not represented here or discussed only briefly. For example, orthogonal wavelets were developed primarily in the context of signal processing, an application upon which we touch only indirectly. However, there are several good expositions (e.g. Daubechies (1990) and Rioul and Vetterli (1991)) of this application. A discussion of wavelet decompositions in the context of Littlewood-Paley theory can be found in the monograph of Frazier et al. (1991). We shall also not attempt to give a complete discussion of the history of wavelets. Historical accounts can be found in the book of Meyer (1990) and the introduction of the article of Daubechies (1990). We shall try to give sufficient historical commentary in the course of our presentation to provide some feeling for the subject's development.