Today’s floating-point arithmetic landscape is broader than ever. While scientific computing has traditionally used single precision and double precision floating-point arithmetics, half precision is increasingly available in hardware and quadruple precision is supported in software. Lower precision arithmetic brings increased speed and reduced communication and energy costs, but it produces results of correspondingly low accuracy. Higher precisions are more expensive but can potentially provide great benefits, even if used sparingly. A variety of mixed precision algorithms have been developed that combine the superior performance of lower precisions with the better accuracy of higher precisions. Some of these algorithms aim to provide results of the same quality as algorithms running in a fixed precision but at a much lower cost; others use a little higher precision to improve the accuracy of an algorithm. This survey treats a broad range of mixed precision algorithms in numerical linear algebra, both direct and iterative, for problems including matrix multiplication, matrix factorization, linear systems, least squares, eigenvalue decomposition and singular value decomposition. We identify key algorithmic ideas, such as iterative refinement, adapting the precision to the data, and exploiting mixed precision block fused multiply–add operations. We also describe the possible performance benefits and explain what is known about the numerical stability of the algorithms. This survey should be useful to a wide community of researchers and practitioners who wish to develop or benefit from mixed precision numerical linear algebra algorithms.

Automatic, or algorithmic, differentiation addresses the need for the accurate and efficient calculation of derivative values in scientific computing. To this end procedural programs for the evaluation of problem-specific functions are transformed into programs that also compute the required derivative values at the same numerical arguments in floating point arithmetic. Disregarding many important implementation issues, we examine in this article complexity bounds and other more mathematical aspects of the program transformation task sketched above.

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.

This is the second part of an article which was started in the previous volume of Acta Numerica. References in the text to Section 1 refer to the preceding article.

We present the field of computational chemistry from the standpoint of numerical analysis. We introduce the most commonly used models and comment on their applicability. We briefly outline the results of mathematical analysis and then mostly concentrate on the main issues raised by numerical simulations. A special emphasis is laid on recent results in numerical analysis, recent developments of new methods and challenging open issues.

Numerical homogenization is a methodology for the computational solution of multiscale partial differential equations. It aims at reducing complex large-scale problems to simplified numerical models valid on some target scale of interest, thereby accounting for the impact of features on smaller scales that are otherwise not resolved. While constructive approaches in the mathematical theory of homogenization are restricted to problems with a clear scale separation, modern numerical homogenization methods can accurately handle problems with a continuum of scales. This paper reviews such approaches embedded in a historical context and provides a unified variational framework for their design and numerical analysis. Apart from prototypical elliptic model problems, the class of partial differential equations covered here includes wave scattering in heterogeneous media and serves as a template for more general multi-physics problems.

We consider a system whose state is given by the solution y to a Partial Differential Equation (PDE) of evolution, and which contains control functions, denoted by v.

The semiclassically scaled time-dependent multi-particle Schrödinger equation describes, inter alia, quantum dynamics of nuclei in a molecule. It poses the combined computational challenges of high oscillations and high dimensions. This paper reviews and studies numerical approaches that are robust to the small semiclassical parameter. We present and analyse variationally evolving Gaussian wave packets, Hagedorn’s semiclassical wave packets, continuous superpositions of both thawed and frozen Gaussians, and Wigner function approaches to the direct computation of expectation values of observables. Making good use of classical mechanics is essential for all these approaches. The arising aspects of time integration and high-dimensional quadrature are also discussed.

This paper surveys in detail the relations between numerical integration and the Hamiltonian (or hybrid) Monte Carlo method (HMC). Since the computational cost of HMC mainly lies in the numerical integrations, these should be performed as efficiently as possible. However, HMC requires methods that have the geometric properties of being volume-preserving and reversible, and this limits the number of integrators that may be used. On the other hand, these geometric properties have important quantitative implications for the integration error, which in turn have an impact on the acceptance rate of the proposal. While at present the velocity Verlet algorithm is the method of choice for good reasons, we argue that Verlet can be improved upon. We also discuss in detail the behaviour of HMC as the dimensionality of the target distribution increases.

Among the iterative methods for solving large linear systems with a sparse (or, possibly, structured) nonsymmetric matrix, those that are based on the Lanczos process feature short recurrences for the generation of the Krylov space. This means low cost and low memory requirement. This review article introduces the reader not only to the basic forms of the Lanczos process and some of the related theory, but also describes in detail a number of solvers that are based on it, including those that are considered to be the most efficient ones. Possible breakdowns of the algorithms and ways to cure them by look-ahead are also discussed.

This article addresses the general problem of establishing upper bounds for the norms of the nth powers of square matrices. The focus is on upper bounds that grow only moderately (or stay constant) where n, or the order of the matrices, increases. The so-called resolvant condition, occuring in the famous Kreiss matrix theorem, is a classical tool for deriving such bounds.

Recently the classical upper bounds known to be valid under Kreiss's resolvant condition have been improved. Moreover, generalizations of this resolvant condition have been considered so as to widen the range of applications. The main purpose of this article is to review and extend some of these new developments.

The upper bounds for the powers of matrices discussed in this article are intimately connected with the stability analysis of numerical processes for solving initial(-boundary) value problems in ordinary and partial linear differential equations. The article highlights this connection.

The article concludes with numerical illustrations in the solution of a simple initial-boundary value problem for a partial differential equation.

Numerical simulation of parametrized differential equations is of crucial importance in the study of real-world phenomena in applied science and engineering. Computational methods for real-time and many-query simulation of such problems often require prohibitively high computational costs to achieve sufficiently accurate numerical solutions. During the last few decades, model order reduction has proved successful in providing low-complexity high-fidelity surrogate models that allow rapid and accurate simulations under parameter variation, thus enabling the numerical simulation of increasingly complex problems. However, many challenges remain to secure the robustness and efficiency needed for the numerical simulation of nonlinear time-dependent problems. The purpose of this article is to survey the state of the art of reduced basis methods for time-dependent problems and draw together recent advances in three main directions. First, we discuss structure-preserving reduced order models designed to retain key physical properties of the continuous problem. Second, we survey localized and adaptive methods based on nonlinear approximations of the solution space. Finally, we consider data-driven techniques based on non-intrusive reduced order models in which an approximation of the map between parameter space and coefficients of the reduced basis is learned. Within each class of methods, we describe different approaches and provide a comparative discussion that lends insights to advantages, disadvantages and potential open questions.

In this article, the basic constructs of algebraic kinematics (links, joints, and mechanism spaces) are introduced. This provides a common schema for many kinds of problems that are of interest in kinematic studies. Once the problems are cast in this algebraic framework, they can be attacked by tools from algebraic geometry. In particular, we review the techniques of numerical algebraic geometry, which are primarily based on homotopy methods. We include a review of the main developments of recent years and outline some of the frontiers where further research is occurring. While numerical algebraic geometry applies broadly to any system of polynomial equations, algebraic kinematics provides a body of interesting examples for testing algorithms and for inspiring new avenues of work.

A classical problem in electrostatics is the determination of the effective electrical conductivity in a composite material consisting of a collection of piecewise homogeneous inclusions embedded in a uniform background. We discuss recently developed fast algorithms for the evaluation of the potential and electrostatic fields induced in multiphase composites by an applied potential, from which the desired effective properties may be easily obtained. The schemes are based on combining a suitable boundary integral equation with the Fast Multipole Method and the GMRES iterative method; the CPU time required grows linearly with the number of points in the discretization of the interface between the inclusions and the background material.

A variety of other questions in electrostatics, magnetostatics and diffusion can be formulated in terms of interface problems. These include the evaluation of electrostatic fields in the presence of dielectric inclusions, the determination of magnetostatic fields in media with variable magnetic permeability, and the calculation of the effective thermal conductivity of a composite material. The methods presented here apply with minor modification to these other situations as well.

We survey and unify recent results on the existence of accurate algorithms for evaluating multivariate polynomials, and more generally for accurate numerical linear algebra with structured matrices. By ‘accurate’ we mean that the computed answer has relative error less than 1, i.e., has some correct leading digits. We also address efficiency, by which we mean algorithms that run in polynomial time in the size of the input. Our results will depend strongly on the model of arithmetic: most of our results will use the so-called traditional model (TM), where the computed result of op(a, b), a binary operation like a+b, is given by op(a, b) * (1+δ) where all we know is that |δ| ≤ ε ≪ 1. Here ε is a constant also known as machine epsilon.

We will see a common reason for the following disparate problems to permit accurate and efficient algorithms using only the four basic arithmetic operations: finding the eigenvalues of a suitably discretized scalar elliptic PDE, finding eigenvalues of arbitrary products, inverses, or Schur complements of totally non-negative matrices (such as Cauchy and Vandermonde), and evaluating the Motzkin polynomial. Furthermore, in all these cases the high accuracy is ‘deserved’, i.e., the answer is determined much more accurately by the data than the conventional condition number would suggest.

In contrast, we will see that evaluating even the simple polynomial x + y + z accurately is impossible in the TM, using only the basic arithmetic operations. We give a set of necessary and sufficient conditions to decide whether a high accuracy algorithm exists in the TM, and describe progress toward a decision procedure that will take any problem and provide either a high-accuracy algorithm or a proof that none exists.

When no accurate algorithm exists in the TM, it is natural to extend the set of available accurate operations by a library of additional operations, such as x + y + z, dot products, or indeed any enumerable set which could then be used to build further accurate algorithms. We show how our accurate algorithms and decision procedure for finding them extend to this case.

Finally, we address other models of arithmetic, and the relationship between (im)possibility in the TM and (in)efficient algorithms operating on numbers represented as bit strings.

The use of computational models and simulations to predict events that take place in our physical universe, or to predict the behaviour of engineered systems, has significantly advanced the pace of scientific discovery and the creation of new technologies for the benefit of humankind over recent decades, at least up to a point. That ‘point’ in recent history occurred around the time that the scientific community began to realize that true predictive science must deal with many formidable obstacles, including the determination of the reliability of the models in the presence of many uncertainties. To develop meaningful predictions one needs relevant data, itself possessing uncertainty due to experimental noise; in addition, one must determine model parameters, and concomitantly, there is the overriding need to select and validate models given the data and the goals of the simulation.

This article provides a broad overview of predictive computational science within the framework of what is often called the science of uncertainty quantification. The exposition is divided into three major parts. In Part 1, philosophical and statistical foundations of predictive science are developed within a Bayesian framework. There the case is made that the Bayesian framework provides, perhaps, a unique setting for handling all of the uncertainties encountered in scientific prediction. In Part 2, general frameworks and procedures for the calculation and validation of mathematical models of physical realities are given, all in a Bayesian setting. But beyond Bayes, an introduction to information theory, the maximum entropy principle, model sensitivity analysis and sampling methods such as MCMC are presented. In Part 3, the central problem of predictive computational science is addressed: the selection, adaptive control and validation of mathematical and computational models of complex systems. The Occam Plausibility Algorithm, OPAL, is introduced as a framework for model selection, calibration and validation. Applications to complex models of tumour growth are discussed.

We review the construction and analysis of numerical methods for strongly nonlinear PDEs, with an emphasis on convex and non-convex fully nonlinear equations and the convergence to viscosity solutions. We begin by describing a fundamental result in this area which states that stable, consistent and monotone schemes converge as the discretization parameter tends to zero. We review methodologies to construct finite difference, finite element and semi-Lagrangian schemes that satisfy these criteria, and, in addition, discuss some rather novel tools that have paved the way to derive rates of convergence within this framework.

This article is about numerical methods for the solution of nonlinear equations. We consider both the fixed-point form $\mathbf{x}=\mathbf{G}(\mathbf{x})$ and the equations form $\mathbf{F}(\mathbf{x})=0$ and explain why both versions are necessary to understand the solvers. We include the classical methods to make the presentation complete and discuss less familiar topics such as Anderson acceleration, semi-smooth Newton’s method, and pseudo-arclength and pseudo-transient continuation methods.

It used to be good enough to bound absolute of matrix eigenvalues and singular values. Not any more. Now it is fashionable to bound relative errors. We present a collection of relative perturbation results which have emerged during the past ten years.

No need to throw away all those absolute error bound, though. Deep down, the derivation of many relative bounds can be based on absolute bounds. This means that relative bounds are not always better. They may just be better sometimes – and exactly when depends on the perturbation.

Many of the boundary value problems traditionally cast as partial differential equations can be reformulated as integral equations over the boundary. After an introduction to boundary integral equations, this review describes some of the methods which have been proposed for their approximate solution. It discusses, as simply as possible, some of the techniques used in their error analysis, and points to areas in which the theory is still unsatisfactory.