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.

In the last few decades the concepts of structure-preserving discretization, geometric integration and compatible discretizations have emerged as subfields in the numerical approximation of ordinary and partial differential equations. The article discusses certain selected topics within these areas; discretization techniques both in space and time are considered. Lie group integrators are discussed with particular focus on the application to partial differential equations, followed by a discussion of how time integrators can be designed to preserve first integrals in the differential equation using discrete gradients and discrete variational derivatives.

Lie group integrators depend crucially on fast and structure-preserving algorithms for computing matrix exponentials. Preservation of domain symmetries is of particular interest in the application of Lie group integrators to PDEs. The equivariance of linear operators and Fourier transforms on non-commutative groups is used to construct fast structure-preserving algorithms for computing exponentials. The theory of Weyl groups is employed in the construction of high-order spectral element discretizations, based on multivariate Chebyshev polynomials on triangles, simplexes and simplicial complexes.

The theory of mixed finite elements is developed in terms of special inverse systems of complexes of differential forms, where the inclusion of cells corresponds to pullback of forms. The theory covers, for instance, composite piecewise polynomial finite elements of variable order over polyhedral grids. Under natural algebraic and metric conditions, interpolators and smoothers are constructed, which commute with the exterior derivative and whose product is uniformly stable in Lebesgue spaces. As a consequence we obtain not only eigenpair approximation for the Hodge–Laplacian in mixed form, but also variants of Sobolev injections and translation estimates adapted to variational discretizations.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Phase retrieval, i.e. the problem of recovering a function from the squared magnitude of its Fourier transform, arises in many applications, such as X-ray crystallography, diffraction imaging, optics, quantum mechanics and astronomy. This problem has confounded engineers, physicists, and mathematicians for many decades. Recently, phase retrieval has seen a resurgence in research activity, ignited by new imaging modalities and novel mathematical concepts. As our scientific experiments produce larger and larger datasets and we aim for faster and faster throughput, it is becoming increasingly important to study the involved numerical algorithms in a systematic and principled manner. Indeed, the past decade has witnessed a surge in the systematic study of computational algorithms for phase retrieval. In this paper we will review these recent advances from a numerical viewpoint.

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

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.

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.