Past Numerical Analysis Group Internal Seminar

I will describe a novel algorithm for computing the Walsh Hadamard Transform (WHT) which consists entirely of Haar wavelet transforms. The algorithm shares precisely the same serial complexity as the popular divide-and-conquer algorithm for the WHT. There is also a natural way to parallelize the algorithm which appears to have a number of attractive features.

Medical imaging is a key diagnostic tool, and is paramount for disease detection and for patient monitoring during ongoing care. Often, to reduce the amount of radiation that a patient is subjected to, there is a strong incentive to consider image reconstruction from incomplete sets of measurements, and so the imaging process is formulated as a compressed sensing problem.

In this talk, we will focus on compressed sensing for digital tomosynthesis (DTS), in which three-dimensional images are reconstructed from a set of two-dimensional X-ray projections. We first discuss a reconstruction approach for static bodies, with a particular interest in the choice of basis for the image representation. We will then focus on the need for accurate image reconstructions when the body of interest is not stationary, but is undergoing simple motion, discussing two different approaches for tackling this dynamic problem.

The design of shapes that are in some sense optimal is a task faced by engineers in a wide range of disciplines. In shape optimisation one aims to improve a given initial shape by iteratively deforming it - if the shape is represented by a mesh, then this means that the mesh has to deformed. This is a delicate problem as overlapping or highly stretched meshes lead to poor accuracy of numerical methods.

In the presented work we consider a novel mesh deformation method motivated by the Riemannian mapping theorem and based on conformal mappings.

Time parallelisation techniques provide an additional direction for the parallelisation of the solution of time-dependent PDEs or of systems of ODEs. In particular, the Parareal algorithm has imposed itself as the canonical choice to achieve parallelisation in time, also because of its simplicity and flexibility. The algorithm works by splitting the time domain in chunks, and iteratively alternating a prediction step (parallel), in which a "fine" solver is employed to achieve a high-accuracy solution within each chunk, to a correction step (serial) where a "coarse" solver is used to quickly propagate the update between the chunks. However, the stability of the method has proven to be highly sensitive to the choice of fine and coarse solver, even more so when applied to chaotic systems or advection-dominated problems.

In this presentation, an alternative formulation of Parareal is discussed. This aims to conduct the update by estimating directly the sensitivity of the solution of the integration with respect to the initial conditions, thus eliminating altogether the necessity of choosing the most apt coarse solver, and potentially boosting its convergence properties.

 

Classical algorithms for numerical integration (quadrature/cubature) proceed by approximating the integrand with a simple function (e.g. a polynomial), and integrate the approximant exactly. In high-dimensional integration, such methods quickly become infeasible due to the curse of dimensionality.

A common alternative is the Monte Carlo method (MC), which simply takes the average of random samples, improving the estimate as more and more samples are taken. The main issue with MC is its slow "sqrt(variance/#samples)" convergence, and various techniques have been proposed to reduce the variance.

In this work we reveal a numerical analyst's interpretation of MC: it approximates the integrand with a simple(st) function, and integrates that function exactly. This observation leads naturally to MC-like methods that combines MC with function approximation theory, including polynomial approximation and sparse grids. The resulting method can be regarded as another variance reduction technique for Monte Carlo.

We develop a general purpose solver for quadratic programs based on operator splitting. We introduce a novel splitting that requires the solution of a quasi-definite linear system with the same coefficient matrix in each iteration. The resulting algorithm is very robust, and once the initial factorization is carried out, division free; it also eliminates requirements on the problem data such as positive definiteness of the objective function or linear independence of the constraint functions. Moreover, it is able to detect primal or dual infeasible problems providing infeasibility certificates. The method supports caching the factorization of the quasi-definite system and warm starting, making it efficient for solving parametrized problems arising in finance, control, and machine learning. Our open-source C implementation OSQP has a small footprint and is library-free. Numerical benchmarks on problems arising from several application domains show that OSQP is typically 10x faster than interior-point methods, especially when factorization caching or warm start is used.

This is joint work with Goran Banjac, Paul Goulart, Alberto Bemporad and Stephen Boyd
 

We provide the rate of convergence of general monotone numerical schemes for parabolic Hamilton-Jacobi-Bellman equations in bounded domains with Dirichlet boundary conditions. The so-called "shaking coefficients" technique introduced by Krylov is used. This technique is based on a perturbation of the dynamics followed by a regularization step by convolution. When restricting the equation to a domain, the perturbed problem may not satisfy such a restriction, so that a special treatment near the boundary is necessary. 

The phenomenon of poor algorithmic scalability is a critical problem in large-scale machine learning and data science. This has led to a resurgence in the use of first-order (Hessian-free) algorithms from classical optimisation. One major drawback is that first-order methods tend to converge extremely slowly. However, there exist techniques for efficiently accelerating them.
    
The topic of this talk is the Dual Regularisation Nonlinear Acceleration algorithm (DRNA) (Geleta, 2017) for nonconvex optimisation. Numerical studies using the CUTEst optimisation problem set show the method to accelerate several nonconvex optimisation algorithms, including quasi-Newton BFGS and steepest descent methods. DRNA compares favourably with a number of existing accelerators in these studies.
    
DRNA extends to the nonconvex setting a recent acceleration algorithm due to Scieur et al. (Advances in Neural Information Processing Systems 29, 2016). We have proven theorems relating DRNA to the Kylov subspace method GMRES, as well as to Anderson's acceleration method and family of multi-secant quasi-Newton methods.
 

In this talk we introduce a novel dynamic programming (DP) approximation that exploits the inherent network structure present in revenue management problems. In particular, our approximation provides a new lower bound on the value function for the DP, which enables conservative revenue forecasts to be made. Existing state of the art approximations of the revenue management DP neglect the network structure, apportioning the prices of each product, whereas our proposed method does not: we partition the network of products into clusters by apportioning the capacities of resources. Our proposed approach allows, in principle, for better approximations of the DP to be made than the decomposition methods currently implemented in industry and we see it as an important stepping stone towards better approximate DP methods in practice.

One of the key challenges in revenue management is unconstraining demand data. Existing state of the art single-class unconstraining methods make restrictive assumptions about the form of the underlying demand and can perform poorly when applied to data which breaks these assumptions. In this talk, we propose a novel unconstraining method that uses Gaussian process (GP) regression. We develop a novel GP model by constructing and implementing a new non-stationary covariance function for the GP which enables it to learn and extrapolate the underlying demand trend. We show that this method can cope with important features of realistic demand data, including nonlinear demand trends, variations in total demand, lengthy periods of constraining, non-exponential inter-arrival times, and discontinuities/changepoints in demand data. In all such circumstances, our results indicate that GPs outperform existing single-class unconstraining methods.

In this talk we describe a new approach that enables the use of elliptic PDEs with white noise forcing to sample Matérn fields within the multilevel Monte Carlo (MLMC) framework.

When MLMC is used to quantify the uncertainty in the solution of PDEs with random coefficients, two key ingredients are needed: 1) a sampling technique for the coefficients that satisfies the MLMC telescopic sum and 2) a numerical solver for the forward PDE problem.

When the dimensionality of the uncertainty in the problem is infinite (i.e. coefficients are random fields), the sampling techniques commonly used in the literature are Karhunen–Loève expansions or circulant embeddings. In the specific case in which the coefficients are Gaussian fields of Mat ́ern covariance structure another sampling technique available relies on the solution of a linear elliptic PDE with white noise forcing.

When the finite element method (FEM) is used for the forward problem, the latter option can become advantageous as elliptic PDEs can be quickly and efficiently solved with the FEM, the sampling can be performed in parallel and the same FEM software can be used without the need for external packages. However, it is unclear how to enforce a good stochastic coupling of white noise between MLMC levels so as to respect the MLMC telescopic sum. In this talk we show how this coupling can be enforced in theory and in practice.

We propose and analyze a multilevel weighted least squares polynomial approximation method. Weighted least squares polynomial approximation uses random samples to determine projections of functions onto spaces of polynomials. It has been shown that using an optimal distribution of sample locations, the number of samples required to achieve quasi-optimal approximation in a given polynomial subspace scales, up to a logarithmic factor, linearly in the dimension of this space. However, in many applications, the computation of samples includes a numerical discretization error. Thus, obtaining polynomial approximations with a single level method can become prohibitively expensive, as it requires a sufficiently large number of samples, each computed with a sufficiently small discretization error. As a solution to this problem, we propose a multilevel method, which employs samples with different accuracies and is able to match the accuracy of single level approximations at reduced computational work. We prove complexity bounds under certain assumptions on polynomial approximability and sample work. Furthermore, we propose an adaptive
algorithm for situations where such assumptions cannot be verified a priori. Numerical experiments underline the practical applicability of our method.

In this talk, a novel discontinuous Galerkin (DG) method is introduced by utilising the principle of discrete least squares. The key idea is to build polynomial approximations by the method of  (weighted) discrete least squares instead of usual interpolation or (discrete) $L^2$ projections. The resulting method hence uses more information of the underlying function and provides a more robust alternative to common DG methods. As a result, we are able to construct high-order schemes which are conservative as well as linear stable on any set of collocation points. Several numerical tests highlight the new discontinuous Galerkin discrete least squares (DG-DLS) method to significantly outperform present-day DG methods.

High-order methods for conservation laws can be highly efficient if their stability is ensured. A suitable means mimicking estimates of the continuous level is provided by summation-by-parts (SBP) operators and the weak enforcement of boundary conditions. Recently, there has been an increasing interest in generalised SBP operators both in the finite difference and the discontinuous Galerkin spectral element framework.

However, if generalised SBP operators are used, the treatment of boundaries becomes more difficult since some properties of the continuous level are no longer mimicked discretely —interpolating the product of two functions will in general result in a value different from the product of the interpolations. Thus, desired properties such as conservation and stability are more difficult to obtain.

In this talk, the concept of generalised SBP operators and their application to entropy stable semidiscretisations will be presented. Several recent ideas extending the range of possible methods are discussed, presenting both advantages and several shortcomings.

We consider a generalization of low-rank matrix completion to the case where the data belongs to an algebraic variety, i.e., each data point is a solution to a system of polynomial equations. In this case, the original matrix is possibly high-rank, but it becomes low-rank after mapping each column to a higher dimensional space of monomial features. Many well-studied extensions of linear models, including affine subspaces and their union, can be described by a variety model. We study the sampling requirements for matrix completion under a variety model with a focus on a union of subspaces. We also propose an efficient matrix completion algorithm that minimizes a surrogate of the rank of the matrix of monomial features, which is able to recover synthetically generated data up to the predicted sampling complexity bounds. The proposed algorithm also outperforms standard low-rank matrix completion and subspace clustering techniques in experiments with real data.

Starting from an example in quantum chemistry, we explain the techniques of Numerical Tensor Calculus with particular emphasis on the convolution operation. The tensorisation technique also applies to one-dimensional grid functions and allows to perform the convolution with a cost which may be much cheaper than the fast Fourier transform.

Magnus expansion based methods are an efficient class of integrators for solving Schrödinger equations that feature time dependent potentials such as lasers. These methods have been found to be highly effective in computational quantum chemistry since the pioneering work of Tal Ezer and Kosloff in the early 90s. The convergence of the Magnus expansion, however, is usually understood only for ODEs and traditional analysis suggests a much poorer performance of these methods than observed experimentally. It was not till the work of Hochbruck and Lubich in 2003 that a rigorous analysis justifying the application to PDEs with unbounded operators, such as the Schrödinger equation, was presented. In this talk we will extend this analysis to the semiclassical regime, where the highly oscillatory solution conventionally suggests large errors and a requirement for very small time steps.
 

We present some dimensionality reduction techniques for global optimization algorithms, in order to increase their scalability. Inspired by ideas in machine learning, and extending the approach of random projections in Zhang et al (2016), we present some new algorithmic approaches for global optimisation with theoretical guarantees of good behaviour and encouraging numerical results.
 

Derivative-free optimisation (DFO) algorithms are a category of optimisation methods for situations when one is unable to compute or estimate derivatives of the objective, such as when the objective has noise or is very expensive to evaluate. In this talk I will present a flexible DFO framework for unconstrained nonlinear least-squares problems, and in particular discuss its performance on noisy problems.

Each year Prof. Trefethen gives the Problem Solving Squad a sequence of problems with no hints, one a week, where the solution of each problem is a single real number to be computed by any method available.  We will present this year's three problems, involving (1) an S-shaped bifurcation curve, (2) shortest path around a random web, and (3) switching a time-varying system to maximize a matrix norm.

The 14 students this year are Simon Vary plus InFoMM cohort 2: Matteo Croci, Davin Lunz, Michael McPhail, Tori Pereira, Lindon Roberts, Caoimhe Rooney, Ian Roper, Thomas Roy, Tino Sulzer, Bogdan Toader, Florian Wechsung, Jess Williams, and Fabian Ying.  The presentations will be by (1) Lindon Roberts, (2) Florian Wechsung, and (3) Thomas Roy.

Delsarte-Goethals frames are a popular choice for deterministic measurement matrices in compressive sensing. I will show that it is possible to construct extremely sparse matrices which share precisely the same row space as Delsarte-Goethals frames. I will also describe the combinatorial block design underlying the construction and make a connection to Steiner equiangular tight frames.
 

What's the continuous analog of randn?  In recent months I've been exploring such questions with Abdul-Lateef Haji-Ali and other members of the Chebfun team, and Chebfun now has commands randnfun, randnfun2, randnfunsphere, and randnfundisk.  These are based on finite Fourier series with random coefficients, and interesting questions arise in the "white noise" limit as the lengths of the series approaches infinity and as random ODEs become stochastic DEs.    This work is at an early stage and we are grateful for input from stochastic experts at Oxford and elsewhere.

We analyse the numerical approximation of functions using radial basis functions in the context of frames. Frames generalize the notion of a basis by allowing redundancy, while being restricted by a so-called frame condition. The theory of numerical frame approximations allows the study of ill-conditioning, inherently due to their redundancy, and suggests discretization techniques that still offer numerical stability to machine precision. We apply the theory to radial basis functions.

The joint spectral radius (JSR) of a set of matrices characterizes the maximum growth rate that can be achieved by multiplying them in arbitrary order. This concept, which essentially generalizes the notion of the "largest eigenvalue" from one matrix to many, was introduced by Rota and Strang in the early 60s and has since emerged in many areas of application such as stability of switched linear systems, computation of the capacity of codes, convergence of consensus algorithms, tracability of graphs, and many others. The JSR is a very difficult quantity to compute even for a pair of matrices. In this talk, we present optimization-based algorithms (e.g., via semidefinite programming or dynamic programming) that can either compute the JSR exactly in special cases or approximate it with arbitrary prescribed accuracy in the general case.

Based on joint work (in different subsets) with Raphael Jungers, Pablo Parrilo, and Mardavij Roozbehani.
 

The past few years have seen a surge of interest in nonconvex reformulations of convex optimizations using nonlinear reparameterizations of the optimization variables. Compared with the convex formulations, the nonconvex ones typically involve many fewer variables, allowing them to scale to scenarios with millions of variables. However, one pays the price of solving nonconvex optimizations to global optimality, which is generally believed to be impossible. In this talk, I will characterize the nonconvex geometries of several low-rank matrix optimizations. In particular, I will argue that under reasonable assumptions, each critical point of the nonconvex problems either corresponds to the global optimum of the original convex optimizations, or is a strict saddle point where the Hessian matrix has a negative eigenvalue. Such a geometric structure ensures that many local search algorithms can converge to the global optimum with random initializations. Our analysis is based on studying how the convex geometries are transformed under nonlinear parameterizations.

In various research fields such as machine learning, compressed sensing and operations research, optimization problems which seek sparsity of solutions by the cardinality constraint or rank constraint are studied. We formulate such problems as DC (Difference of two Convex functions) optimization problems and apply DC Algorithm (DCA) to them. While a subproblem needs to be solved in each DCA iteration, its closed-form solution can be easily obtained by soft-thresholding operation. Numerical experiments demonstrate the efficiency of the proposed DCA in comparison with existing methods.
This is a joint work with J. Gotoh (Chuo Univ.) and K. Tono (U. Tokyo). 

Joint work with Wolfgang Hackbusch and Daniel Kressner

We consider the classic problem of establishing a statistical ranking of a set of n items given a set of inconsistent and incomplete pairwise comparisons between such items. Instantiations of this problem occur in numerous applications in data analysis (e.g., ranking teams in sports data), computer vision, and machine learning. We formulate the above problem of ranking with incomplete noisy information as an instance of the group synchronization problem over the group SO(2) of planar rotations, whose usefulness has been demonstrated in numerous applications in recent years. Its least squares solution can be approximated by either a spectral or a semidefinite programming (SDP) relaxation, followed by a rounding procedure. We perform extensive numerical simulations on both synthetic and real-world data sets (Premier League soccer games, a Halo 2 game tournament and NCAA College Basketball games) showing that our proposed method compares favorably to other algorithms from the recent literature.

We propose a similar synchronization-based algorithm for the rank-aggregation problem, which integrates in a globally consistent ranking pairwise comparisons given by different rating systems on the same set of items. We also discuss the problem of semi-supervised ranking when there is available information on the ground truth rank of a subset of players, and propose an algorithm based on SDP which recovers the ranks of the remaining players. Finally, synchronization-based ranking, combined with a spectral technique for the densest subgraph problem, allows one to extract locally-consistent partial rankings, in other words, to identify the rank of a small subset of players whose pairwise comparisons are less noisy than the rest of the data, which other methods are not able to identify. 
 

The trigonometric interpolants to a periodic function f in equispaced points converge if f is Dini-continuous, and the associated quadrature formula, the trapezoidal rule, converges if f is continuous.  What if the points are perturbed?  Amazingly little has been done on this problem, or on its algebraic (i.e. nonperiodic) analogue.  I will present new results joint with Anthony Austin which show some surprises.

Using a variational approach applied to generalized Rayleigh functionals, we extend the concepts of singular values and singular functions to trivariate functions defined on a rectangular parallelepiped. We also consider eigenvalues and eigenfunctions for trivariate functions whose domain is a cube. For a general finite-rank trivariate function, we describe an algorithm for computing the canonical polyadic (CP) decomposition, provided that the CP factors are linearly independent in two variables. All these notions are computed using Chebfun. Application in finding the best rank-1 approximation of trivariate functions is investigated. We also prove that if the function is analytic and two-way orthogonally decomposable (odeco), then the CP values decay geometrically, and optimal finite-rank approximants converge at the same rate.
 

Joint work with Jen Pestana.

There are several classes of random function that appear naturally in mathematical physics, probability, number theory, and other areas of mathematics. I will give a brief overview of some of these random functions and explain what they are and why they are important. Finally, I will explain how I use chebfun to study these functions.
 

Stochastic Hamiltonian systems with multiplicative noise are a mathematical model for many physical systems with uncertainty. For example, they can be used to describe synchrotron oscillations of a particle in a storage ring. Just like their deterministic counterparts, stochastic Hamiltonian systems possess several important geometric features; for instance, their phase flows preserve the canonical symplectic form. When simulating these systems numerically, it is therefore advisable that the numerical scheme also preserves such geometric structures. In this talk we propose a variational principle for stochastic Hamiltonian systems and use it to construct stochastic Galerkin variational integrators. We show that such integrators are indeed symplectic, preserve integrals of motion related to Lie group symmetries, demonstrate superior long-time energy behavior compared to nonsymplectic methods, and they include stochastic symplectic Runge-Kutta methods as a special case. We also analyze their convergence properties and present the results of several numerical experiments.