Past Numerical Analysis Group Internal Seminar

When solving stochastic partial differential equations (SPDEs) driven by additive spatial white noise the efficient sampling of white noise realizations can be challenging. In this talk we present a novel sampling technique that can be used to efficiently compute white noise samples in a finite element and multilevel Monte Carlo (MLMC) setting.
After discretization, the action of white noise on a test function yields a Gaussian vector with the FEM mass matrix as covariance. Sampling such a vector requires an expensive Cholesky factorization and for this reason P0 representations, for which the mass matrix is diagonal, are generally preferred in the literature. This however has other disadvantages. In this talk we introduce an alternative factorization that is naturally parallelizable and has linear cost and memory complexity (in the number of mesh elements).
Moreover, in a MLMC framework the white noise samples must be coupled between subsequent levels so as to respect the telescoping sum. We show how our technique can be used to enforce this coupling even in the case in which the hierarchy is non-nested via a supermesh construction. We conclude the talk with numerical experiments that demonstrate the efficacy of our method. We observe optimal convergence rates for the finite element solution of the elliptic SPDEs of interest. In a MLMC setting, a good coupling is enforced and the telescoping sum is respected.
 

(Joint work with Coralia Cartis) The problem of finding the most extreme value of a function, also known as global optimization, is a challenging task. The difficulty is associated with the exponential increase in the computational time for a linear increase in the dimension. This is known as the ``curse of dimensionality''. In this talk, we demonstrate that such challenges can be overcome for functions with low effective dimensionality --- functions which are constant along certain linear subspaces. Such functions can often be found in applications, for example, in hyper-parameter optimization for neural networks, heuristic algorithms for combinatorial optimization problems and complex engineering simulations.
We propose the use of random subspace embeddings within a(ny) global minimisation algorithm, extending the approach in Wang et al. (2013). We introduce a new framework, called REGO (Random Embeddings for GO), which transforms the high-dimensional optimization problem into a low-dimensional one. In REGO, a new low-dimensional problem is formulated with bound constraints in the reduced space and solved with any GO solver. Using random matrix theory, we provide probabilistic bounds for the success of REGO, which indicate that this is dependent upon the dimension of the embedded subspace and the intrinsic dimension of the function, but independent of the ambient dimension. Numerical results demonstrate that high success rates can be achieved with only one embedding and that rates are for the most part invariant with respect to the ambient dimension of the problem.
 

In this talk, first we  address the convergence issues of a standard finite volume element method (FVEM) applied to simple elliptic problems. Then, we discuss discontinuous finite volume element methods (DFVEM) for elliptic problems  with emphasis on  computational and theoretical  advantages over the standard FVEM. Further, we present a natural extension of DFVEM employed for the elliptic problem to the Stokes problems. We also discuss suitability of these methods for the approximation of incompressible miscible displacement problems.
 

The vast majority of the world's electricity is generated by converting thermal energy into electric energy by use of a steam turbine. Siemens are one of the worlds leading manufacturers of such
turbines, and aim to design theirs to be as efficient as possible. Using an internally built software, Siemens can estimate the efficiency which would result from a turbine design. In this presentation, we present the approaches that have been taken to improve turbine design using mathematical optimisation software. In particular, we focus on the failings of the approach currently taken, the obstacles in place which make solving this problem difficult, and the approach we intend to take to find a locally optimal solution.

This talk will review the main Tikhonov and Bayesian smoothing formulations of inverse problems for dynamical systems with partially observed variables and parameters. The main contenders: strong-constraint, weak-constraint and penalty function formulations will be described. The relationship between these formulations and associated optimisation problems will be revealed.  To close we will indicate techniques for maintaining sparsity and for quantifying uncertainty.

In this talk we address the numerical approximation of Mean Field Games with local couplings. For finite difference discretizations of the Mean Field Game system, we follow a variational approach, proving that the schemes can be obtained as the optimality system of suitably defined optimization problems. In order to prove the existence of solutions of the scheme with a variational argument, the monotonicity of the coupling term is not used, which allow us to recover general existence results. Next, assuming next that the coupling term is monotone, the variational problem is cast as a convex optimization problem for which we study and compare several proximal type methods. These algorithms have several interesting features, such as global convergence and stability with respect to the viscosity parameter. We conclude by presenting numerical experiments assessing the performance of the proposed methods. In collaboration with L. Briceno-Arias (Valparaiso, CL) and F. J. Silva (Limoges, FR).

Semidefinite convex optimization problems often have low-rank solutions that can be represented with O(p)-storage. However, semidefinite programming methods require us to store the matrix decision variable with size O(p^2), which prevents the application of virtually all convex methods at large scale.

Indeed, storage, not arithmetic computation, is now the obstacle that prevents us from solving large- scale optimization problems. A grand challenge in contemporary optimization is therefore to design storage-optimal algorithms that provably and reliably solve large-scale optimization problems in key scientific and engineering applications. An algorithm is called storage optimal if its working storage is within a constant factor of the memory required to specify a generic problem instance and its solution.

So far, convex methods have completely failed to satisfy storage optimality. As a result, the literature has largely focused on storage optimal non-convex methods to obtain numerical solutions. Unfortunately, these algorithms have been shown to be provably correct only under unverifiable and unrealistic statistical assumptions on the problem template. They can also sacrifice the key benefits of convexity, as they do not use key convex geometric properties in their cost functions.

To this end, my talk introduces a new convex optimization algebra to obtain numerical solutions to semidefinite programs with a low-rank matrix streaming model. This streaming model provides us an opportunity to integrate sketching as a new tool for developing storage optimal convex optimization methods that go beyond semidefinite programming to more general convex templates. The resulting algorithms are expected to achieve unparalleled results for scalable matrix optimization problems in signal processing, machine learning, and computer science.

The Schrödinger equation with a time-dependent potential occurs in a wide range of applications in theoretical chemistry, quantum physics and quantum computing. In this talk I will discuss a variety of Magnus expansion based schemes that have been found to be highly effective for numerically solving these equations since the pioneering work of Tal Ezer and Kosloff in the early 90s. Recent developments in the field focus on approximation of the exponential of the Magnus expansion via exponential splittings including some asymptotic splittings and commutator-free splittings that are designed specifically for this task.

I will also present a very recently developed methodology for the case of laser-matter interaction. This methodology allows us to extend any fourth-order scheme for Schrödinger equation with time-independent potential to a fourth-order method for Schrödinger equation with laser potential with little to no additional cost. These fourth-order methods improve upon many leading schemes of order six due to their low costs and small error constants.

Many problems that involve the propagation of time-harmonic waves are naturally posed in unbounded domains. For instance, a common problem in the are a of acoustic scattering is the determination of the sound field that is generated when an incoming time-harmonic wave (which is assumed to arrive ``from infinity'') impinges onto a solid body (the scatterer). The boundary
conditions to be applied on the surface of the scatterer (most often of Dirichlet, Neumann or Robin type) tend to be easy to enforce in most numerical solution schemes. Conversely, the imposition of a suitable decay condition (typically a variant of the Sommerfeld radiation condition), which is required to ensure the well-posedness of the solution, is considerably more involved. As a result, many numerical schemes generate spurious reflections from the outer boundary of the finite computational domain.

Perfectly matched layers (PMLs) are in this context a versatile alternative to the usage of classical approaches such as employing absorbing boundary conditions or Dirichlet-to-Neumann mappings, but unfortunately most PML formulations contain adjustable parameters which have to be optimised to give the best possible performance for a particular problem. In this talk I will present a parameter-free PML formulation for the case of the two-dimensional Helmholtz equation. The performance of the proposed method is demonstrated via extensive numerical experiments, involving domains with smooth and polygonal boundaries, different solution types (smooth and singular, planar and non-planar waves), and a wide range of wavenumbers (R. Cimpeanu, A. Martinsson and M.Heil, J. Comp. Phys., 296, 329-347 (2015)). Possible extensions and generalisations will also be touched upon.

We consider numerical methods for solving  time dependent partial differential equations with convection-diffusion terms and anti-diffusive fractional operator of order $\alpha \in (1,2)$. These equations are motivated by two distinct applications: a dune morphodynamics model and a signal filtering method. 
We propose numerical schemes based on local discontinuous Galerkin methods to approximate the solutions of these equations. Numerical stability and convergence of these schemes are investigated. 
Finally numerical experiments are given to illustrate qualitative behaviors of solutions for both applications and to confirme the convergence results. 

We introduce and analyze a discontinuous Galerkin method for the Oseen equations in two dimension spaces. The boundary conditions are mixed and they are assumed to be of three different types:
the vorticity  and the normal component of the velocity are given on a first part of the boundary, the pressure and the tangential component of the velocity are given on a second part of the boundary and the Dirichlet condition is given on the remainder part . We establish a priori error estimates in the energy norm for the velocity and in the L2 norm for the pressure. An a posteriori error estimate is also carried out yielding optimal convergence rate. The analysis is based on rewriting the method in a non-consistent manner using lifting operators in the spirit of Arnold, Brezzi, Cockburn and Marini.

We construct a space of vector fields that are normal to differentiable curves in the plane. Its basis functions are defined via saddle point variational problems in reproducing kernel Hilbert spaces (RKHSs). First, we study the properties of these basis vector fields and show how to approximate them. Then, we employ this basis to discretise shape Newton methods and investigate the impact of this discretisation on convergence rates.

Structure from Motion (SfM) is a problem which asks: given photos of an object from different angles, can we reconstruct the object in 3D? This problem is important in computer vision, with applications including urban planning and autonomous navigation. A key part of SfM is bundle adjustment, where initial estimates of 3D points and camera locations are refined to match the images. This results in a high-dimensional nonlinear least-squares problem, which is typically solved using the Gauss-Newton method. In this talk, I will discuss how dimensionality reduction methods such as block coordinates and randomised sketching can be used to improve the scalability of Gauss-Newton for bundle adjustment problems.

Over the past decade, the randomized singular value decomposition (RSVD) algorithm has proven to be an efficient, reliable alternative to classical algorithms for computing low-rank approximations in a number of applications. However, in cases where no information is available on the singular value decay of the data matrix or the data matrix is known to be close to full-rank, the RSVD is ineffective. In recent years, there has been great interest in randomized algorithms for computing full factorizations that excel in this regime.  In this talk, we will give a brief overview of some key ideas in randomized numerical linear algebra and introduce a new randomized algorithm for computing a full, rank-revealing URV factorization.

In oil and gas reservoir simulation, standard preconditioners involve solving a restricted pressure system with AMG. Initially designed for isothermal models, this approach is often used in the thermal case. However, it does not incorporate heat diffusion or the effects of temperature changes on fluid flow through viscosity and density. We seek to develop preconditioners which consider this cross-coupling between pressure and temperature. In order to study the effects of both pressure and temperature on fluid and heat flow, we first consider a model of non-isothermal single phase flow through porous media. By focusing on single phase flow, we are able to isolate the properties of the pressure-temperature subsystem. We present a numerical comparison of different preconditioning approaches including block preconditioners.

The restricted isometry property is arguably the most prominent tool in the theory of compressive sensing. In its classical version, it features l_2 norms as inner and outer norms. The modified version considered in this talk features the l_1 norm as the inner norm, while the outer norm depends a priori on the distribution of the random entries populating the measurement matrix.  The modified version holds for a wider class of random matrices and still accounts for the success of sparse recovery via basis pursuit and via iterative hard thresholding. In the special case of Gaussian matrices, the outer norm actually reduces to an l_2 norm. This fact allows one to retrieve results from the theory of one-bit compressive sensing in a very simple way. Extensions to one-bit matrix recovery are then straightforward.
 

Forecasting a diagnosis of Alzheimer’s disease is a promising means of selection for clinical trials of Alzheimer’s disease therapies. A positive PET scan is commonly used as part of the inclusion criteria for clinical trials, but PET imaging is expensive, so when a positive scan is one of the trial inclusion criteria it is desirable to avoid screening failures. In this talk I will describe a scheme for pre-selecting participants using statistical learning methods, and investigate how brain regions change as the disease progresses.  As a means of generating features I apply the Chen path signature. This is a systematic way of providing feature sets for multimodal data that can probe the nonlinear interactions in the data as an extension of the usual linear features. While it can easily perform a traditional analysis, it can also probe second and higher order events for their predictive value. Combined with Lasso regularisation one can auto detect situations where the observed data has nonlinear information.

The latest CPUs by Intel and ARM support vectorised operations, where a single set of instructions (e.g. add, multiple, bit shift, XOR, etc.) are performed in parallel for small batches of data. This can provide great performance improvements if each parallel instruction performs the same operation, but carries the risk of performance loss if each needs to perform different tasks (e.g. if else conditions). I will present the work I have done so far looking into how to recover the full performance of the hardware, and some of the challenges faced when trading off between ever larger parallel tasks, risks of tasks diverging, and how certain coding styles might be modified for memory bandwidth limited applications. Examples will be taken from finance and Monte Carlo applications, inspecting some standard maths library functions and possibly random number generation.

Low-rank plus sparse matrices arise in many data-oriented applications, most notably in a foreground-background separation from a moving camera. It is known that low-rank matrix recovery from a few entries (low-rank matrix completion) requires low coherence (Candes et al 2009) as in the extreme cases when the low-rank matrix is also sparse, where matrix completion can miss information and be unrecoverable. However, the requirement of low coherence does not suffice in the low-rank plus sparse model, as the set of low-rank plus sparse matrices is not closed. We will discuss the relation of non-closedness of the low-rank plus sparse model to the notion of matrix rigidity function in complexity theory.

The object of this talk is a class of generalised Newtonian fluids with implicit constitutive law.
Both in the steady and the unsteady case, existence of weak solutions was proven by Bul\'\i{}\v{c}ek et al. (2009, 2012) and the main challenge is the small growth exponent qq and the implicit law.
I will discuss the application of a splitting and regularising strategy to show convergence of FEM approximations to weak solutions of the flow. 
In the steady case this allows to cover the full range of growth exponents and thus generalises existing work of Diening et al. (2013). If time permits, I will also address the unsteady case.
This is joint work with Endre Suli.

We consider the problem of localising non-negative point sources, namely finding their locations and amplitudes from noisy samples which consist of the convolution of the input signal with a known kernel (e.g. Gaussian). In contrast to the existing literature, which focuses on TV-norm minimisation, we analyse the feasibility problem. In the presence of noise, we show that the localised error is proportional to the level of noise and depends on the distance between each source and the closest samples. This is achieved using duality and considering the spectrum of the associated sampling matrix.

A simple experiment in the field of electrochemistry involves  controlling the applied potential in an electrochemical cell. This  causes electron transfer to take place at the electrode surface and in turn this causes a current to flow. The current depends on parameters in  the system and the inverse problem requires us to estimate these  parameters given an experimental trace of the current. We briefly  describe recent work in this area from simple least squares approximation of the parameters, through bootstrapping to estimate the distributions of the parameters, to MCMC methods which allow us to see correlations between parameters.

Within the wide field of sparse approximation, convolutional sparse coding (CSC) has gained considerable attention in the computer vision and machine learning communities. While several works have been devoted to the practical aspects of this model, a systematic theoretical understanding of CSC seems to have been left aside. In this talk, I will present a novel analysis of the CSC problem based on the observation that, while being global, this model can be characterized and analyzed locally. By imposing only local sparsity conditions, we show that uniqueness of solutions, stability to noise contamination and success of pursuit algorithms are globally guaranteed. I will then present a Multi-Layer extension of this model and show its close relation to Convolutional Neural Networks (CNNs). This connection brings a fresh view to CNNs, as one can attribute to this architecture theoretical claims under local sparse assumptions, which shed light on ways of improving the design and implementation of these networks. Last, but not least, we will derive a learning algorithm for this model and demonstrate its applicability in unsupervised settings.

In this talk, we revisit the cubic regularization (CR) method for solving smooth non-convex optimization problems and study its local convergence behaviour. In their seminal paper, Nesterov and Polyak showed that the sequence of iterates of the CR method converges quadratically a local minimum under a non-degeneracy assumption, which implies that the local minimum is isolated. However, many optimization problems from applications such as phase retrieval and low-rank matrix recovery have non-isolated local minima. In the absence of the non-degeneracy assumption, the result was downgraded to the superlinear convergence of function values. In particular, they showed that the sequence of function values enjoys a superlinear convergence of order 4/3 (resp. 3/2) if the function is gradient dominated (resp. star-convex and globally non-degenerate). To remedy the situation, we propose a unified local error bound (EB) condition and show that the sequence of iterates of the CR method converges quadratically a local minimum under the EB condition. Furthermore, we prove that the EB condition holds if the function is gradient dominated or if it is star-convex and globally non-degenerate, thus improving the results of Nesterov and Polyak in three aspects: weaker assumption, faster rate and iterate instead of function value convergence. Finally, we apply our results to two concrete non-convex optimization problems that arise from phase retrieval and low-rank matrix recovery. For both problems, we prove that with overwhelming probability, the local EB condition is satisfied and the CR method converges quadratically to a global optimizer. We also present some numerical results on these two problems.

The question of what happens to the eigenvalues of a matrix after an additive perturbation has a long history, with notable contributions from Wilkinson, Sorensen, Golub, H\"ormander, Ipsen and Mehrmann, among many others. If the perturbed matrix $C \in \mathbb{C}^{n \times n}$ is given by $C = A + B$, these theorems typically consider the case where $A$ and/or $B$ are symmetric and $B$ has rank one. In this talk we will prove a theorem that bounds the number of distinct eigenvalues of $C$ in terms of the number of distinct eigenvalues of $A$, the diagonalisability of $A$, and the rank of $B$. This new theorem is more general in that it applies to arbitrary matrices $A$ perturbed by matrices of arbitrary rank $B$. We will also discuss various refinements of my bound recently developed by other authors.
 

We consider a system of nonlinear partial differential equations modelling the steady motion of an incompressible non-Newtonian fluid, which is chemically reacting. The governing system consists of a steady convection-diffusion equation for the concentration and the generalized steady Navier–Stokes equations, where the viscosity coefficient is a power-law type function of the shear-rate, and the coupling between the equations results from the concentration-dependence of the power-law index. This system of nonlinear partial differential equations arises in mathematical models of the synovial fluid found in the cavities of moving joints. We construct a finite element approximation of the model and perform the mathematical analysis of the numerical method. Key technical tools include discrete counterparts of the Bogovski operator, De Giorgi’s regularity theorem and the Acerbi–Fusco Lipschitz truncation of Sobolev functions, in function spaces with variable integrability exponents.

In this talk we discuss the convergence rate of the Newton method for finding the singularity point on vetor fields. It is well-known that the Newton Method has local quadratic convergence rate with nonsingularity and Lipschitz condition. Here we release Lipschitz condition. With only nonsingularity, the Newton Method has superlinear convergence. If we have enough time, we can quickly give the damped Newton method on finding singularity on vector fields with superlinear convergence under nonsingularity condition only.

In this talk we consider the issue of mass loss in fragmentation models due to 'shattering'. As a solution we propose a hybrid discrete/continuous model whereby the smaller particles are considered as having discrete mass, whilst above a certain cut-off, mass is taken to be a continuous variable. The talk covers the development of such a model, its initial analysis via the theory of operator semigroups and its numerical approximation using a finite volume discretisation.

We explore the use of applying multiple preconditioners for solving linear systems arising in simulations of incompressible two-phase flow. In particular, we use a selective MPGMRES algorithm, for which the search space grows linearly throughout the iterative solver, and block preconditioners based on Schur complement approximations

We introduce a discontinuous Galerkin finite element method (DGFEM) for Hamilton–Jacobi–Bellman equations on piecewise curved domains, and prove that the method is consistent, stable, and produces optimal convergence rates. Upon utilising a long standing result due to N. Krylov, we may characterise the Monge–Ampère equation as a HJB equation; in two dimensions, this HJB equation can be characterised further as uniformly elliptic HJB equation, allowing for the application of the DGFEM

We propose a new parameter estimation technique for SDEs, based on the inverse problem of finding a forward operator describing the evolution of temporal data. Nonlinear dynamical systems on a state-space can be lifted to linear dynamical systems on spaces of higher, often infinite, dimension. Recently, much work has gone into approximating these higher-dimensional systems with linear operators calculated from data, using what is called Dynamic Mode Decomposition (DMD). For SDEs, this linear system is given by a second-order differential operator, which we can quickly calculate and compare to the DMD operator.

In 2016-17, Fjordholm, Kappeli, Mishra and Tadmor developed a numerical method by which one could compute measure-valued solutions to systems of hyperbolic conservation laws with either measure-valued or deterministic initial data. In this talk I will discuss the ideas behind this method, and discuss how it can be adapted to systems of quasi-linear parabolic PDEs whose nonlinearity fails to satisfy a monotonicity condition.

The purpose of this work is to introduce a new constraint functional for shape optimization problems, which enforces the constructibility by means of additive manufacturing processes, and helps in preventing the appearance of overhang features - large regions hanging over void which are notoriously difficult to assemble using such technologies. The proposed constraint relies on a simplified model for the construction process: it involves a continuum of shapes, namely the intermediate shapes corresponding to the stages of the construction process where the total, final shape is erected only up to a certain level. The shape differentiability of this constraint functional is analyzed - which is not a standard issue because of its peculiar structure. Several numerical strategies and examples are then presented. This is a joint work with G. Allaire, R. Estevez, A. Faure and G. Michailidis.

A new generation of low cost 3D tomography systems is based on multiple emitters and sensors that partially convolve measurements. A successful approach to deconvolve the measurements is to use nonlinear compressed sensing models. We discuss such models, as well as algorithms for their solution. 

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.