Past Computational Mathematics and Applications Seminar

Solving the Stokes equation by an optimal domain decomposition method derived algebraically involves the use of non standard interface conditions whose discretisation is not trivial. For this reason the use of approximation methods such as hybrid discontinuous Galerkin appears as an appropriate strategy: on the one hand they provide the best compromise in terms of the number of degrees of freedom in between standard continuous and discontinuous Galerkin methods, and on the other hand the degrees of freedom used in the non standard interface conditions are naturally defined at the boundary between elements. In this work we introduce the coupling between a well chosen discretisation method (hybrid discontinuous Galerkin) and a novel and efficient domain decomposition method to solve the Stokes system. We present the detailed analysis of the hybrid discontinuous Galerkin method for the Stokes problem with non standard boundary conditions. This analysis is supported by numerical evidence. In addition, the advantage of the new preconditioners over more classical choices is also supported by numerical experiments.

This work was done in collaboration with G. Barrenechea, M. Bosy (Univ. Strathclyde) and F. Nataf, P-H Tournier (Univ of Paris VI)

In this talk we describe an approach to approximate the truncated singular value decomposition of a large matrix by first decomposing the matrix into a sum of Kronecker products. Our approach can be used to more efficiently approximate a large number of singular values and vectors than other well known schemes, such as iterative algorithms based on the Golub-Kahan bidiagonalization or randomized matrix algorithms. We provide theoretical results and numerical experiments to demonstrate accuracy of our approximation, and show how the approximation can be used to solve large scale ill-posed inverse problems, either as an approximate filtering method, or as a preconditioner to accelerate iterative algorithms.
 

We analyze a fully discrete numerical scheme for solving a parabolic PDE on a moving surface. The method is based on a diffuse interface approach that involves a level set description of the moving surface. Under suitable conditions on the spatial grid size, the time step and the interface width we obtain stability and error bounds with respect to natural norms. Test calculations are presented that confirm our analysis.

Spline curves represent a simple and efficient tool for data interpolation in Euclidean space. During the past decades, however, more and more applications have emerged that require interpolation in (often high-dimensional) nonlinear spaces such as Riemannian manifolds. An example is the generation of motion sequences in computer graphics, where the animated figure represents a curve in a Riemannian space of shapes. Two particularly useful spline interpolation methods derive from a variational principle: linear splines minimize the average squared velocity and cubic splines minimize the average squared acceleration among all interpolating curves. Those variational principles and their discrete analogues can be used to define continuous and discretized spline curves on (possibly infinite-dimensional) Riemannian manifolds. However, it turns out that well-posedness of cubic splines is much more intricate on nonlinear and high-dimensional spaces and requires quite strong conditions on the underlying manifold. We will analyse and discuss linear and cubic splines as well as their discrete counterparts on Riemannian manifolds and show a few applications.

Maintenance activities help prevent costly power generator breakdowns but because generators under maintenance are typically unavailable, the impact of maintenance schedules is significant and their cost must be accounted for when planning maintenance. In this paper we address the generator maintenance scheduling problem in hydropower systems. While this problem has been widely studied, specific operating conditions of hydroelectric systems have received less attention. We present a mixed-integer linear programming model that considers the time windows of the maintenance activities, as well as the nonlinearities and disjunctions of the hydroelectric production functions. Because the resulting model is hard to solve, we also propose an extended formulation, a set reduction approach that uses logical conditions for excluding unnecessary set elements from the model, and valid inequalities. Computational experiments using a variety of instances adapted from a real hydropower system in Canada support the conclusion that the extended formulation with set reduction achieves the best results in terms of computational time and optimality gap. This is joint work with Jesus Rodriguez, Pascal Cote and Guy Desaulniers.

Because of atmospheric turbulence, images of objects in outer space acquired via ground-based telescopes are usually blurry.  One way to estimate the blurring kernel or point spread function (PSF) is to make use of the aberration of wavefront received at the telescope, i.e., the phase. However only the low-resolution wavefront gradients can be collected by wavefront sensors. In this talk, I will discuss how to use regularization methods to reconstruct high-resolution phase gradients and then use them to recover the phase and the PSF in high accuracy. I will end by relating the problem to high-resolution image reconstruction and methods for solving it.
Joint work with Rui Zhao and research supported by HKRGC.

Several optimization problems combine nonlinear constraints with the integrality of a subset of variables. For an important class of problems  called Mixed Integer Second-Order Cone Optimization (MISOCO), with applications in facility location, robust optimization, and finance, among others, these nonlinear constraints are second-order (or Lorentz) cones.

For such problems, as for many discrete optimization problems, it is crucial to understand the properties of the union of two disjoint sets of feasible solutions. To this end, we apply the disjunctive programming paradigm to MISOCO and present conditions under which the convex hull of two disjoint sets can be obtained by intersecting the feasible set with a specially constructed second-order cone. Computational results show that such cone has a positive impact on the solution of MISOCO problems.

The mathematical analysis and numerical simulation of acoustic and electromagnetic wave scattering by planar screens is a classical topic. The standard technique involves reformulating the problem as a boundary integral equation on the screen, which can be solved numerically using a boundary element method. Theory and computation are both well-developed for the case where the screen is an open subset of the plane with smooth (e.g. Lipschitz or smoother) boundary. In this talk I will explore the case where the screen is an arbitrary subset of the plane; in particular, the screen could have fractal boundary, or itself be a fractal. Such problems are of interest in the study of fractal antennas in electrical engineering, light scattering by snowflakes/ice crystals in atmospheric physics, and in certain diffraction problems in laser optics. The roughness of the screen presents challenging questions concerning how boundary conditions should be enforced, and the appropriate function space setting. But progress is possible and there is interesting behaviour to be discovered: for example, a sound-soft screen with zero area (planar measure zero) can scatter waves provided the fractal dimension of the set is large enough. Accurate computations are also challenging because of the need to adapt the mesh to the fine structure of the fractal. As well as presenting numerical results, I will outline some of the outstanding open questions from the point of view of numerical analysis. This is joint work with Simon Chandler-Wilde (Reading) and Andrea Moiola (Pavia).
 

In this talk, I will present a novel solution strategy to efficiently and accurately compute approximate solutions to semilinear optimal control problems, focusing on the optimal control of phase field formulations of geometric evolution laws.
The optimal control of geometric evolution laws arises in a number of applications in fields including material science, image processing, tumour growth and cell motility.
Despite this, many open problems remain in the analysis and approximation of such problems.
In the current work we focus on a phase field formulation of the optimal control problem, hence exploiting the well developed mathematical theory for the optimal control of semilinear parabolic partial differential equations.
Approximation of the resulting optimal control problem is computationally challenging, requiring massive amounts of computational time and memory storage.
The main focus of this work is to propose, derive, implement and test an efficient solution method for such problems. The solver for the discretised partial differential equations is based upon a geometric multigrid method incorporating advanced techniques to deal with the nonlinearities in the problem and utilising adaptive mesh refinement.
An in-house two-grid solution strategy for the forward and adjoint problems, that significantly reduces memory requirements and CPU time, is proposed and investigated computationally.
Furthermore, parallelisation as well as an adaptive-step gradient update for the control are employed to further improve efficiency.
Along with a detailed description of our proposed solution method together with its implementation we present a number of computational results that demonstrate and evaluate our algorithms with respect to accuracy and efficiency.
A highlight of the present work is simulation results on the optimal control of phase field formulations of geometric evolution laws in 3-D which would be computationally infeasible without the solution strategies proposed in the present work.

Adjoint derivatives reveal the sensitivity of a computer program's output to changes in its inputs. These derivatives are useful as a building block for optimisation, uncertainty quantification, noise estimation, inverse design, etc., in many industrial and scientific applications that use PDE solvers or other codes.
Algorithmic differentiation (AD) is an established method to transform a given computation into its corresponding adjoint computation. One of the key challenges in this process is the efficiency of the resulting adjoint computation. This becomes especially pressing with the increasing use of shared-memory parallelism on multi- and many-core architectures, for which AD support is currently insufficient.
In this talk, I will present an overview of challenges and solutions for the differentiation of shared-memory-parallel code, using two examples: an unstructured-mesh CFD solver, and a structured-mesh stencil kernel, both parallelised with OpenMP. I will show how AD can be used to generate adjoint solvers that scale as well as their underlying original solvers on CPUs and a KNC XeonPhi. The talk will conclude with some recent efforts in using AD and formal verification tools to check the correctness of manually optimised adjoint solvers.
 

Gauss invented Gaussian quadrature following an approach entirely different from the one we now find in textbooks. I will describe leisurely the contents of Gauss's original memoir on quadrature, an impressive piece of mathematics, based on continued fractions, Padé approximation, generating functions, the hypergeometric series and more.

The talk will describe accelerated algorithms for computing full or partial matrix factorizations such as the eigenvalue decomposition, the QR factorization, etc. The key technical novelty is the use of  randomized projections to reduce the effective dimensionality of  intermediate steps in the computation. The resulting algorithms execute faster on modern hardware than traditional algorithms, and are particularly well suited for processing very large data sets.

The algorithms described are supported by a rigorous mathematical analysis that exploits recent work in random matrix theory. The talk will briefly review some representative theoretical results.

For many applications in science and engineering, the ability to efficiently and accurately approximate solutions to elliptic PDEs dictates what physical phenomena can be simulated numerically.  In this seminar, we present a high-order accurate discretization technique for variable coefficient PDEs with smooth coefficients.  The technique comes with a nested dissection inspired direct solver that scales linearly or nearly linearly with respect to the number of unknowns.  Unlike the application of nested dissection methods to classic discretization techniques, the constant prefactors do not grow with the order of the discretization.  The discretization is robust even for problems with highly oscillatory solutions.  For example, a problem 100 wavelengths in size can be solved to 9 digits of accuracy with 3.7 million unknowns on a desktop computer.  The precomputation of the direct solver takes 6 minutes on a desktop computer.  Then applying the computed solver takes 3 seconds.  The recent application of the algorithm to inverse media scattering also will be presented.

Structural optimization can be interpreted as the attempt to find the best mechanical structure to support specific load cases respecting some possible constraints. Within this context, topology optimization aims to obtain the connectivity, shape and location of voids inside a prescribed structural design domain. The methods for the design of stiff lightweight structures are well established and can already be used in a specific range of industries where such structures are important, e.g., in aerospace and automobile industries.

In this seminar, we will go through the basic engineering concepts used to quantify and analyze the computational models of mechanical structures. After presenting the motivation, the methods and mathematical tools used in structural topology optimization will be discussed. In our method, an implicit level set function is used to describe the structural boundaries. The optimization problem is approximated by linearization of the objective and constraint equations via Taylor’s expansion. Shape sensitivities are used to evaluate the change in the structural performance due to a shape movement and to feed the mathematical optimiser in an iterative procedure. Recent developments comprising multiscale and Multiphysics problems will be presented and a specific application proposal including acoustic-structure interaction will be discussed.

We describe a convergence acceleration technique for generic optimization problems. Our scheme computes estimates of the optimum from a nonlinear average of the iterates produced by any optimization method. The weights in this average are computed via a simple linear system, whose solution can be updated online. This acceleration scheme runs in parallel to the base algorithm, providing improved estimates of the solution on the fly, while the original optimization method is running. Numerical experiments are detailed on classical classification problems.
 

Abstract: We study nonuniform sampling in shift-invariant spaces whose generator is a totally positive function. For a subclass of such generators the sampling theorems can be formulated in analogy to the theorems of Beurling and Landau for bandlimited functions. These results are  optimal and validate  the  heuristic reasonings in the engineering literature. In contrast to the cardinal series, the reconstruction procedures for sampling in a shift-invariant space with a totally positive generator  are local and thus accessible to numerical linear algebra.

A subtle  connection between sampling in shift-invariant spaces and the theory of Gabor frames leads to new and optimal  results for Gabor frames.  We show that the set of phase-space shifts of  $g$ (totally positive with a Gaussian part) with respect to a rectangular lattice forms a frame, if and only if the density of the lattice  is strictly larger than 1. This solves an open problem going backto Daubechies in 1990 for the class of totally positive functions of Gaussian type.
 

Almost all real-world applications involve a degree of uncertainty. This may be the result of noisy measurements, restrictions on observability, or simply unforeseen events. Since many models in both engineering and the natural sciences make use of partial differential equations (PDEs), it is natural to consider PDEs with random inputs. In this context, passing from modelling and simulation to optimization or control results in stochastic PDE-constrained optimization problems. This leads to a number of theoretical, algorithmic, and numerical challenges.

 From a mathematical standpoint, the solution of the underlying PDE is a random field, which in turn makes the quantity of interest or the objective function an implicitly defined random variable. In order to minimize this distributed objective, one can use, e.g., stochastic order constraints, a distributionally robust approach, or risk measures. In this talk, we will make use of risk measures.

After motivating the approach via a model for the mitigation of an airborne pollutant, we build up an analytical framework and introduce some useful risk measures. This allows us to prove the existence of solutions and derive optimality conditions. We then present several approximation schemes for handling non-smooth risk measures in order to leverage existing numerical methods from PDE-constrained optimization. Finally, we discuss solutions techniques and illustrate our results with numerical examples.

During the last decade, the progress in the computational performance of commercial mixed-integer programming solvers have been significant. Part of this success is due to faster computers and better software engineering but a more significant part of it is due to the power of the cutting planes used in these solvers.
In the first part of this talk, we will discuss main components of a MIP solver and describe some classical families of valid inequalities (Gomory mixed integer cuts, mixed integer rounding cuts, split cuts, etc.) that are routinely used in these solvers. In the second part, we will discuss recent progress in cutting plane theory that has not yet made its way to commercial solvers. In particular, we will discuss cuts from lattice-free convex sets and answer a long standing question in the affirmative by deriving a finite cutting plane algorithm for mixed-integer programming.

In variational imaging and other inverse problem modeling, regularisation plays a major role. In recent years, high order regularizers such as the total generalised variation, the mean curvature and the Gaussian curvature are increasingly studied and applied, and many improved results over the widely-used total variation model are reported.
Here we first introduce the fractional order derivatives and the total fractional-order variation which provides an alternative  regularizer and is not yet formally analysed. We demonstrate that existence and uniqueness properties of the new model can be analysed in a fractional BV space, and, equally, the new model performs as well as the high order regularizers (which do not yet have much theory). 
In the usual framework, the algorithms of a fractional order model are not fast due to dense matrices involved. Moreover, written in a Bregman framework, the resulting Sylvester equation with Toeplitz coefficients can be solved efficiently by a preconditioned solver. Further ideas based on adaptive integration can also improve the computational efficiency in a dramatic way.
 Numerical experiments will be given to illustrate the advantages of the new regulariser for both restoration and registration problems.
 

We will present a very general framework for unconstrained stochastic optimization which is based on standard trust region framework using  random models. In particular this framework retains the desirable features such step acceptance criterion, trust region adjustment and ability to utilize of second order models. We make assumptions on the stochasticity that are different from the typical assumptions of stochastic and simulation-based optimization. In particular we assume that our models and function values satisfy some good quality conditions with some probability fixed, but can be arbitrarily bad otherwise. We will analyze the convergence and convergence rates of this general framework and discuss the requirement on the models and function values. We will will contrast our results with existing results from stochastic approximation literature. We will finish with examples of applications arising the area of machine learning. 
 

Work with Jemima Tabeart, Sarah Dance, Nancy Nichols, Joanne Waller (University of Reading) and Stefano Migliorini, Fiona Smith (Met Office). 
In environmental prediction variational data assimilation (DA) is a method for using observational data to estimate the current state of the system. The DA problem is usually solved as a very large nonlinear least squares problem, in which the fit to the measurements is balanced against the fit to a previous model forecast. These two terms are weighted by matrices describing the correlations of the errors in the forecast and in the observations. Until recently most operational weather and ocean forecasting systems assumed that the errors in the observations are uncorrelated. However, as we move to higher resolution observations then it is becoming more important to specify observation error correlations. In this work we look at the effect this has on the conditioning of the optimization problem. In the context of a linear system we develop bounds on the condition number of the problem in the presence of correlated observation errors. We show that the condition number is very dependent on the minimum eigenvalue of the observation error correlation matrix. We then present results using the Met Office data assimilation system, in which different methods for reconditioning the correlation matrix are tested. We investigate the effect of these different methods on the conditioning and the final solution of the problem.
 

Inverse problems are ubiquitous in many areas of Science and Engineering and, once discretised, they lead to ill-conditioned linear systems, often of huge dimensions: regularisation consists in replacing the original system by a nearby problem with better numerical properties, in order to find meaningful approximations of its solution. In this talk we will explore the regularisation properties of many iterative methods based on Krylov subspaces. After surveying some basic methods such as CGLS and GMRES, innovative approaches based on flexible variants of CGLS and GMRES will be presented, in order to efficiently enforce nonnegativity and sparsity into the solution.