Past Computational Mathematics and Applications Seminar

Functions are usually approximated numerically in a basis, a non-redundant and complete set of functions that span a certain space. In this talk we highlight a number of benefits of using overcomplete sets, in particular using the more general notion of a "frame". The main 

benefit is that frames are easily constructed even for functions of several variables on domains with irregular shapes. On the other hand, allowing for possible linear depencies naturally leads to ill-conditioning of approximation algorithms. The ill-conditioning is 

potentially severe. We give some useful examples of frames and we first address the numerical stability of best approximations in a frame. Next, we briefly describe special point sets in which interpolation turns out to be stable. Finally, we review so-called Fourier extensions and an efficient algorithm to approximate functions with spectral accuracy on domains without structure.

When assigned with the task of extracting information from given image data the first challenge one faces is the derivation of a truthful model for both the information and the data. Such a model can be determined by the a-priori knowledge about the image (information), the data and their relation to each other. The source of this knowledge is either our understanding of the type of images we want to reconstruct and of the physics behind the acquisition of the data or we can thrive to learn parametric models from the data itself. The common question arises: how can we customise our model choice to a particular application? Or better how can we make our model adaptive to the given data?

Starting from the first modelling strategy this talk will lead us from nonlinear diffusion equations and subdifferential inclusions of total variation type functionals as the most successful image modeltoday to non-smooth second- and third-order variational models, with data models for Gaussian and Poisson distributed data as well as impulse noise. These models exhibit solution-dependent adaptivities in form of nonlinearities or non-smooth terms in the PDE or the variational problem, respectively. Applications for image denoising, inpainting and surface reconstruction are given. After a critical discussion of these different image and data models we will turn towards the second modelling strategy and propose to combine it with the first one using a PDE constrained optimisation method that customises a parametrised form of the model by learning from examples. In particular, we will consider optimal parameter derivation for total variation denoising with multiple noise distributions and optimising total generalised variation regularisation for its application in photography.

Equations of quantum mechanics in the semiclassical regime present an enduring challenge for numerical analysts, because their solution is highly oscillatory and evolves on two scales. Standard computational approaches to the semiclassical Schrödinger equation do not allow for long time integration as required, for example, in quantum control of atoms by short laser bursts. This has motivated our approach of asymptotic splittings. Combining techniques from Lie-algebra theory and numerical algebra, we present a new computational paradigm of symmetric Zassenhaus splittings, which lends itself to a very precise discretisation in long time intervals, at very little cost. We will illustrate our talk by examples of quantum phenomena – quantum tunnelling and quantum scattering – and their computation and, time allowing, discuss an extension of this methodology to time-dependent semiclassical systems using Magnus expansions

This talk focuses on the direct search method, arguably one of the simplest optimization algorithms. The algorithm minimizes an objective function by iteratively evaluating it along a number of (polling) directions, which are typically taken from so-called positive spanning sets. It does not use derivatives.

We first introduce the worst case complexity theory of direct search, and discuss how to choose the positive spanning set to minimize the complexity bound. The discussion leads us to a long-standing open
problem in Discrete Geometry. A recent result on this problem enables us to establish the optimal order for the worst case complexity of direct search.

We then show how to achieve even lower complexity bound by using random polling directions. It turns out that polling along two random directions at each iteration is sufficient to guarantee the convergence
of direct search for any dimension, and the resultant algorithm enjoys lower complexity both in theory and in practice.

The last part of the talk is devoted to direct search based on inaccurate function values. We address three questions:
i) what kind of solution can we obtain by direct search if the function values are inaccurate? 
ii) what is the worst case complexity to attain such a solution? iii) given
the inaccuracy in the function values, when to stop the algorithm in order
to guarantee the quality of the solution and also avoid “over-optimization”?

This talk is based on joint works with F. Delbos, M. Dodangeh, S. Gratton, B. Pauwels, C. W. Royer, and L. N. Vicente.

I will review some recent work on the problem of reliable automatic detection of blow-up behaviour for nonlinear parabolic PDEs. The adaptive algorithms developed are based on rigorous conditional a posteriori error bounds. The use of space-time adaptivity is crucial in making the problem computationally tractable. The results presented are applicable to quite general spatial operators, rendering the approach potentially useful in informing respective PDE theory. The new adaptive algorithm is shown to accurately estimate the blow-up time of a number of problems, including ones exhibiting regional blow-up. 

The talk will introduce the concepts of multilevel optimization and motivate them in the context of problems arising from the discretization of infinite dimensional applications. It will be shown how optimization methods can accomodate a number of useful (and classical) ideas from the multigrid
community, and thereby produce substantial efficiency improvements compared to existing large-scale minimization techniques.  Two different classes of multilevel methods will be discussed: trust-region and linesearch algorithms.
The first class will be described in the context of a multilevel generalization of the well-known trust-region-Newton method.  The second will focus on limited-memory quasi-Newton algorithms.  Preliminary numerical results will be presented which indicate that both types of multilevel algorithms may be practically very advantageous.

In numerical atmosphere models, values of relevant physical parameters are often uncertain by more than 100% and weather forecast skill is significantly reduced after a couple of days. Still, numerical operations are typically calculated in double precision with 15 significant decimal digits. If we reduce numerical precision, we can reduce power consumption and increase computational performance significantly. If savings are reinvested to build larger supercomputers, this would allow an increase in resolution in weather and climate models and might lead to better predictions of future weather and climate. 
I will discuss approaches to reduce numerical precision beyond single precision in high performance computing and in particular in weather and climate modelling. I will present results that show that precision can be reduced significantly in atmosphere models and that potential savings can be huge. I will also discuss how rounding errors will impact model dynamics and interact with model uncertainty and predictability.

We propose the use of a constraint preconditioner for the iterative solution of the linear system arising from the finite element discretization of the coupled Stokes-Darcy system. The Stokes-Darcy system is a set of coupled PDEs that can be used to model a freely flowing fluid over porous media flow. The fully coupled system matrix is large, sparse, non-symmetric, and of saddle point form. We provide for exact versions of the constraint preconditioner spectral and field-of-values bounds that are independent of the underlying mesh width. We present several numerical experiments, using the deal.II finite element library, that illustrate our results in both two and three dimensions. We compare exact and inexact versions of the constraint preconditioner against standard block diagonal and block lower triangular preconditioners to illustrate its favorable properties.

When formulated appropriately, the broad families of sequential quadratic programming, augmented Lagrangian and interior-point methods all require the solution of symmetric saddle-point linear systems. When regularization is employed, the systems become symmetric and quasi definite. The latter are
indefinite but their rich structure and strong relationships with definite systems enable specialized linear algebra, and make them prime candidates for matrix-free implementations of optimization methods. In this talk, I explore various formulations of the step equations in optimization and corresponding
iterative methods that exploit their structure.

Security Constrained Optimal Power Flow is an increasingly important problem for power systems operation both in its own right and as a subproblem for more complex problems such as transmission switching or
unit commitment.

The structure of the problem resembles stochastic programming problems in that one aims to find a cost optimal operation schedule that is feasible for all possible equipment outage scenarios
(contingencies). Due to the presence of power flow constraints (in their "DC" or "AC" version), the resulting problem is a large scale linear or nonlinear programming problem.

However it is known that only a small subset of the contingencies is active at the solution. We show how Interior Point methods can exploit this structure both by simplifying the linear algebra operations as
well as generating necessary contingencies on the fly and integrating them into the algorithm using IPM warmstarting techniques. The final problem solved by this scheme is significantly smaller than the full
contingency constrained problem, resulting in substantial speed gains.

Numerical and theoretical results of our algorithm will be presented.

Can we extend the FEM to general polytopic, i.e. polygonal and polyhedral, meshes while retaining 
the ease of implementation and computational cost comparable to that of standard FEM? Within this talk, I present two approaches that achieve just  that (and much more): the Virtual Element Method (VEM) and an hp-version discontinuous Galerkin (dG) method.

The Virtual Element spaces are like the usual (polynomial) finite element spaces with the addition of suitable non-polynomial functions. This is far from being a novel idea. The novelty of the VEM approach is that it avoids expensive evaluations of the non-polynomial "virtual" functions by basing all 
computations solely on the method's carefully chosen degrees of freedom. This way we can easily deal 
with complicated element geometries and/or higher continuity requirements (like C1, C2, etc.), while 
maintaining the computational complexity comparable to that of standard finite element computations.

As you might expect, the choice and number of the degrees of freedom depends on such continuity 
requirements. If mesh flexibility is the goal, while one is ready to  give up on global regularity, other approaches can be considered. For instance, dG approaches are naturally suited to deal with polytopic meshes. Here I present an extension of the classical Interior Penalty dG method which achieves optimal rates of convergence on polytopic meshes even under elemental edge/face degeneration. 

The next step is to exploit mesh flexibility in the efficient resolution of problems characterised by 
complicated geometries and solution features, for instance within the framework of automatic FEM 
adaptivity. I shall finally introduce ongoing work in this direction.

In this seminar I will present a semi-langrangian discretisation of the Monge-Ampère operator, which is of interest in optimal transport 
and differential geometry as well as in related fields of application.

I will discuss the proof of convergence to viscosity solutions. To address the challenge of uniqueness and convexity we draw upon the classical relationship with Hamilton-Jacobi-Bellman equations, which we extend to the viscosity setting. I will explain that the monotonicity of semi-langrangian schemes implies that they possess large stencils, which in turn requires careful treatment of the boundary conditions.

The contents of the seminar is based on current work with X Feng from the University of Tennessee.

The Singular Value Decomposition (SVD) of matrices and the related Principal Components Analysis (PCA) express a matrix in terms of singular vectors, which are linear combinations of all the input data and lack an intuitive physical interpretation. Motivated by the application of PCA and SVD in the analysis of populations genetics data, we will discuss the notion of leverage scores: a simple statistic that reveals columns/rows of a matrix that lie in the subspace spanned by the top principal components (left/right singular vectors). We will then use the leverage scores to present matrix decompositions that express the structure in a matrix in terms of actual columns (and/or rows) of the matrix. Such decompositions are easier to interpret in applications, since the selected columns and rows are subsets of the data. We will also discuss extensions of the leverage scores to reveal influential entries of a matrix.

We present a trust region algorithm for solving nonconvex optimization problems that, in the worst-case, is able to drive the norm of the gradient of the objective below a prescribed threshold $\epsilon > 0$ after at most ${\cal O}(\epsilon^{-3/2})$ function evaluations, gradient evaluations, or iterations.  Our work has been inspired by the recently proposed Adaptive Regularisation framework using Cubics (i.e., the ARC algorithm), which attains the same worst-case complexity bound.  Our algorithm is modeled after a traditional trust region algorithm, but employs modified step acceptance criteria and a novel trust region updating mechanism that allows it to achieve this desirable property.  Importantly, our method also maintains standard global and fast local convergence guarantees.

Although Toeplitz matrices are often dense, matrix-vector products with Toeplitz matrices can be quickly performed via circulant embedding and the fast Fourier transform. This makes their solution by preconditioned Krylov subspace methods attractive. 

For a wide class of symmetric Toeplitz matrices, symmetric positive definite circulant preconditioners that cluster eigenvalues have been proposed. MINRES or the conjugate gradient method can be applied to these problems and descriptive convergence theory based on eigenvalues guarantees fast convergence. 

In contrast, although circulant preconditioners have been proposed for nonsymmetric Toeplitz systems, guarantees of fast convergence are generally only available for CG for the normal equations (CGNE). This is somewhat unsatisfactory because CGNE has certain drawbacks, including slow convergence and a larger condition number. In this talk we discuss a simple alternative symmetrization of nonsymmetric Toeplitz matrices, that allows us to use MINRES to solve the resulting linear system. We show how existing circulant preconditioners for nonsymmetric Toeplitz matrices can be straightforwardly adapted to this situation and give convergence estimates similar to those in the symmetric case.

Numerical simulation tools for fluid and solid mechanics are often based on the discretisation of coupled systems of partial differential equations, which can easily be identified in terms of physical
conservation laws.  In contrast, much physical insight is often gained from the equivalent formulation of the relevant energy or free-energy functional, possibly subject to constraints.  Motivated by the
nonlinear static and dynamic behaviour of nematic liquid crystals and of magnetosensitive elastomers, we propose a finite-element framework for minimising these free-energy functionals, using Lagrange multipliers to enforce additional constraints.  This talk will highlight challenges, limitations, and successes, both in the formulation of these models and their use in numerical simulation.
This is joint work with PhD students Thomas Benson, David Emerson, and Dong Han, and with James Adler, Timothy Atherton, and Luis Dorfmann.

Preconditioning is of significant importance in the solution of large dimensional systems of linear equations such as those that arise from the numerical solution of partial differential equation problems. In this talk we will attempt a broad ranging review of preconditioning.

In this talk, we discuss the development of fast iterative solvers for matrix systems arising from various constrained optimization problems. In particular, we seek to exploit the saddle point structure of these problems to construct powerful preconditioners for the resulting systems, using appropriate approximations of the (1,1)-block and Schur complement.

The problems we consider arise from two well-studied subject areas within computational optimization. Specifically, we investigate the
numerical solution of PDE-constrained optimization problems, and the interior point method (IPM) solution of linear/quadratic programming
problems. Indeed a particular focus in this talk is the interior point method solution of PDE-constrained optimization problems with
additional inequality constraints on the state and control variables.

We present a range of optimization problems which we seek to solve using our methodology, and examine the theoretical and practical
convergence properties of our iterative methods for these problems.
 

Stability of the hp-Raviart-Thomas projection operator as a mapping H^1(K) -> H^1(K) on the unit cube K in R^3 has been addressed e.g. in [2], see also [1]. These results are suboptimal with respect to the polynomial degree. In this talk we present quasi-optimal stability estimates for the hp-Raviart-Thomas projection operator on the cube. The analysis involves elements of the polynomial approximation theory on an interval and the real method of Banach space interpolation.

(Joint work with Herbert Egger, TU Darmstadt)

[1] Mark Ainsworth and Katia Pinchedez. hp-approximation theory for BDFM and RT finite elements on quadrilaterals. SIAM J. Numer. Anal., 40(6):2047–2068 (electronic) (2003), 2002.

[2] Dominik Schötzau, Christoph Schwab, and Andrea Toselli. Mixed hp-DGFEM for incompressible flows. SIAM J. Numer. Anal., 40(6):2171–2194 (electronic) (2003), 2002.

Nonlinear systems of partial differential equations (PDEs) may permit several distinct solutions. The typical current approach to finding distinct solutions is to start Newton's method with many different initial guesses, hoping to find starting points that lie in different basins of attraction. In this talk, we present an infinite-dimensional deflation algorithm for systematically modifying the residual of a nonlinear PDE problem to eliminate known solutions from consideration. This enables the Newton--Kantorovitch iteration to converge to several different solutions, even starting from the same initial guess. The deflated Jacobian is dense, but an efficient preconditioning strategy is devised, and the number of Krylov iterations is observed not to grow as solutions are deflated. The technique is then applied to computing distinct solutions of nonlinear PDEs, tracing bifurcation diagrams, and to computing multiple local minima of PDE-constrained optimisation problems.

In this talk we will describe the steps towards self-consistent computer simulations of the evolution of the universe beginning soon after the Big Bang and ending with the formation of realistic stellar systems like the Milky Way. This is a multi-scale problem of vast proportions. The first step has been the Millennium Simulation, one of the largest and most successful numerical simulations of the Universe ever carried out. Now we are in the midst of the next step, where this is carried to a much higher level of physical fidelity on the latest supercomputing platforms. This talk will be illustrate how the role of mathematics is essential in this endeavor. Also computer simulations will be shown. This is joint work among others with Volker Springel.