Past Computational Mathematics and Applications Seminar

The finite element method (FEM) is one of the great triumphs of applied mathematics, numerical analysis and software development. Recent developments in sensor and signalling technologies enable the phenomenological study of systems. The connection between sensor data and FEM is restricted to solving inverse problems placing unwarranted faith in the fidelity of the mathematical description of the system. If one concedes mis-specification between generative reality and the FEM then a framework to systematically characterise this uncertainty is required. This talk will present a statistical construction of the FEM which systematically blends mathematical description with observations.

>

The task of solving large scale linear algebraic problems such as factorising matrices or solving linear systems is of central importance in many areas of scientific computing, as well as in data analysis and computational statistics. The talk will describe how randomisation can be used to design algorithms that in many environments have both better asymptotic complexities and better practical speed than standard deterministic methods.

The talk will in particular focus on randomised algorithms for solving large systems of linear equations. Both direct solution techniques based on fast factorisations of the coefficient matrix, and techniques based on randomised preconditioners, will be covered.

Note: There is a related talk in the Random Matrix Seminar on Tuesday Feb 25, at 15:30 in L4. That talk describes randomised methods for computing low rank approximations to matrices. The two talks are independent, but the Tuesday one introduces some of the analytical framework that supports the methods described here.

In this talk I will present a Perron-Frobenius type result for nonlinear eigenvector problems which allows us to compute the global maximum of a class of constrained nonconvex optimization problems involving multihomogeneous functions.

I will structure the talk into three main parts:

First, I will motivate the optimization of homogeneous functions from a graph partitioning point of view, showing an intriguing generalization of the famous Cheeger inequality.

Second, I will define the concept of multihomogeneous function and I will state our main Perron-Frobenious theorem. This theorem exploits the connection between optimization of multihomogeneous functions and nonlinear eigenvectors to provide an optimization scheme that has global convergence guarantees.

Third, I will discuss a few example applications in network science and machine learning that require the optimization of multihomogeneous functions and that can be solved using nonlinear Perron eigenvectors.

Systems of nonlinear polynomial equations arise in a variety of fields in mathematics, science, and engineering.  Many numerical techniques for solving and analyzing solution sets of polynomial equations over the complex numbers, collectively called numerical algebraic geometry, have been developed over the past several decades.  However, since real solutions are the only solutions of interest in many applications, there is a current emphasis on developing new methods for computing and analyzing real solution sets.  This talk will summarize some numerical real algebraic geometric approaches as well as recent successes of these methods for solving a variety of problems in science and engineering.

(Joint work with: Jüri Lember, Heinrich Matzinger, Raul Kangro)

Principal component analysis is an important pattern recognition and dimensionality reduction tool in many applications and are computed as eigenvectors

of a maximum likelihood covariance that approximates a population covariance. The eigenvectors are often used to extract structural information about the variables (or attributes) of the studied population. Since PCA is based on the eigen-decomposition of the proxy covariance rather than the ground-truth, it is important to understand the approximation error in each individual eigenvector as a function of the number of available samples. The combination of recent results of Koltchinskii & Lounici [8] and Yu, Wang & Samworth [11] yields such bounds. In the presented work we sharpen these bounds and show that eigenvectors can often be reconstructed to a required accuracy from a sample of strictly smaller size order.

We discuss the design of algorithms and codes for the solution of large sparse systems of linear equations on extreme scale computers that are characterized by having many nodes with multi-core CPUs or GPUs. We first use two approaches to get good single node performance. For symmetric systems we use task-based algorithms based on an assembly tree representation of the factorization. We then use runtime systems for scheduling the computation on both multicore CPU nodes and GPU nodes [6]. In this work, we are also concerned with the efficient parallel implementation of the solve phase using the computed sparse factors, and we show impressive results relative to other state-of-the-art codes [3]. Our second approach was to design a new parallel threshold Markowitz algorithm [4] based on Luby’s method [7] for obtaining a maximal independent set in an undirected graph. This is a significant extension since our graph model is a directed graph. We then extend the scope of both these approaches to exploit distributed memory parallelism. In the first case, we base our work on the block Cimmino algorithm [1] using the ABCD software package coded by Zenadi in Toulouse [5, 8]. The kernel for this algorithm is the direct factorization of a symmetric indefinite submatrix for which we use the above symmetric code. To extend the unsymmetric code to distributed memory, we use the Zoltan code from Sandia [2] to partition the matrix to singly bordered block diagonal form and then use the above unsymmetric code on the blocks on the diagonal. In both cases, we illustrate the added parallelism obtained from combining the distributed memory parallelism with the high single-node performance and show that our codes out-perform other state-of-the-art codes. This work is joint with a number of people. We developed the algorithms and codes in an EU Horizon 2020 Project, called NLAFET, that finished on 30 April 2019. Coworkers in this were: Sebastien Cayrols, Jonathan Hogg, Florent Lopez, and Stojce ´ ∗@email 1 Nakov. Collaborators in the block Cimmino part of the project were: Philippe Leleux, Daniel Ruiz, and Sukru Torun. Our codes available on the github repository https://github.com/NLAFET.

References [1] M. ARIOLI, I. S. DUFF, J. NOAILLES, AND D. RUIZ, A block projection method for sparse matrices, SIAM J. Scientific and Statistical Computing, 13 (1992), pp. 47–70. [2] E. BOMAN, K. DEVINE, L. A. FISK, R. HEAPHY, B. HENDRICKSON, C. VAUGHAN, U. CATALYUREK, D. BOZDAG, W. MITCHELL, AND J. TERESCO, Zoltan 3.0: Parallel Partitioning, Load-balancing, and Data Management Services; User’s Guide, Sandia National Laboratories, Albuquerque, NM, 2007. Tech. Report SAND2007-4748W http://www.cs.sandia. gov/Zoltan/ug_html/ug.html. [3] S. CAYROLS, I. S. DUFF, AND F. LOPEZ, Parallelization of the solve phase in a task-based Cholesky solver using a sequential task flow model, Int. J. of High Performance Computing Applications, To appear (2019). NLAFET Working Note 20. RAL-TR-2018-008. [4] T. A. DAVIS, I. S. DUFF, AND S. NAKOV, Design and implementation of a parallel Markowitz threshold algorithm, Technical Report RAL-TR-2019-003, Rutherford Appleton Laboratory, Oxfordshire, England, 2019. NLAFET Working Note 22. Submitted to SIMAX. [5] I. S. DUFF, R. GUIVARCH, D. RUIZ, AND M. ZENADI, The augmented block Cimmino distributed method, SIAM J. Scientific Computing, 37 (2015), pp. A1248–A1269. [6] I. S. DUFF, J. HOGG, AND F. LOPEZ, A new sparse symmetric indefinite solver using a posteriori threshold pivoting, SIAM J. Scientific Computing, To appear (2019). NLAFET Working Note 21. RAL-TR-2018-012. [7] M. LUBY, A simple parallel algorithm for the maximal independent set problem, SIAM J. Computing, 15 (1986), pp. 1036–1053. [8] M. ZENADI, The solution of large sparse linear systems on parallel computers using a hybrid implementation of the block Cimmino method., These de Doctorat, ´ Institut National Polytechnique de Toulouse, Toulouse, France, decembre 2013.

Boundary integral equations are an elegant tool to model and simulate a range of physical phenomena in bounded and unbounded domains.

While mathematically well understood, the numerical implementation (e.g. via boundary element methods) still poses a number of computational challenges, from the efficient assembly of the underlying linear systems up to the fast preconditioned solution in complex applications. In this talk we provide an overview of some of these challenges and demonstrate the efficient implementation of boundary element methods on modern CPU and GPU architectures. As part of the talk we will present a number of practical examples using the Bempp-cl boundary element software, our next generation boundary element package, that has been developed in Python and supports modern vectorized CPU instruction sets and a number of GPU types.

Iterative algorithms for the solution of convex quadratic optimization problems are investigated, which exploit inaccurate matrix-vector products. Theoretical bounds on the performance of a Conjugate Gradients method are derived, the necessary quantities occurring in the theoretical bounds estimated and a new practical algorithm derived. Numerical experiments suggest that the new method has significant potential, including in the steadily more important context of multi-precision computations.

Everybody is familiar with the concept of eigenvalues of a matrix. In this talk, we consider the nonlinear eigenvalue problem. These are problems for which the eigenvalue parameter appears in a nonlinear way in the equation. In physics, the Schroedinger equation for determining the bound states in a semiconductor device, introduces terms with square roots of different shifts of the eigenvalue. In mechanical and civil engineering, new materials often have nonlinear damping properties. For the vibration analysis of such materials, this leads to nonlinear functions of the eigenvalue in the system matrix.

One particular example is the sandwhich beam problem, where a layer of damping material is sandwhiched between two layers of steel. Another example is the stability analysis of the Helmholtz equation with a noise excitation produced by burners in a combustion chamber. The burners lead to a boundary condition with delay terms (exponentials of the eigenvalue).

We often receive the question: “How can we solve a nonlinear eigenvalue problem?” This talk explains the different steps to be taken for using Krylov methods. The general approach works as follows: 1) approximate the nonlinearity by a rational function; 2) rewrite this rational eigenvalue problem as a linear eigenvalue problem and then 3) solve this by a Krylov method. We explain each of the three steps.

The fully coupled numerical simulation of different physical processes, which can typically occur
at variable time and space scales, is often a very challenging task. A common feature of such models is that
their discretization gives rise to systems of linearized equations with an inherent block structure, which
reflects the properties of the set of governing PDEs. The efficient solution of a sequence of systems with
matrices in a block form is usually one of the most time- and memory-demanding issue in a coupled simulation.
This effort can be carried out by using either iteratively coupled schemes or monolithic approaches, which
tackle the problem of the system solution as a whole.

This talk aims to discuss recent advances in the monolithic solution of coupled multi-physics problems, with
application to poromechanical simulations in fractured porous media. The problem is addressed either by proper
sparse approximations of the Schur complements or special splittings that can partially uncouple the variables
related to different physical processes. The selected approaches can be included in a more general preconditioning
framework that can help accelerate the convergence of Krylov subspace solvers. The generalized preconditioner
relies on approximately decoupling the different processes, so as to address each single-physics problem
independently of the others. The objective is to provide an algebraic framework that can be employed as a
general ``black-box'' tool and can be regarded as a common starting point to be later specialized for the
particular multi-physics problem at hand.

Numerical experiments, taken from real-world examples of poromechanical problems and fractured media, are used to
investigate the behaviour and the performance of the proposed strategies.

Domain decomposition methods are widely employed for the numerical solution of partial differential equations on parallel computers. We develop an adjoint-based a posteriori error analysis for overlapping multiplicative Schwarz domain decomposition and for overlapping additive Schwarz. In both cases the numerical error in a user-specified functional of the solution (quantity of interest), is decomposed into a component that arises due to the spatial discretization and a component that results from of the finite iteration between the subdomains. The spatial discretization error can be further decomposed in to the errors arising on each subdomain. This decomposition of the total error can then be used as part of a two-stage approach to construct a solution strategy that efficiently reduces the error in the quantity of interest.

The development of effective solvers for high frequency wave propagation problems, such as those described by the Helmholtz equation, presents significant challenges. One promising class of solvers for such problems are parallel domain decomposition methods, however, an appropriate coarse space is typically required in order to obtain robust behaviour (scalable with respect to the number of domains, weakly dependant on the wave number but also on the heterogeneity of the physical parameters). In this talk we introduce a coarse space based on generalised eigenproblems in the overlap (GenEO) for the Helmholtz equation. Numerical results within FreeFEM demonstrate convergence that is effectively independent of the wave number and contrast in the heterogeneous coefficient as well as good performance for minimal overlap.

Can we get rigorous answers when computing with real and complex numbers? There are now many applications where this is possible thanks to a combination of tools from computer algebra and traditional numerical computing. I will give an overview of such methods in the context of two projects I'm developing. The first project, Arb, is a library for arbitrary-precision ball arithmetic, a form of interval arithmetic enabling numerical computations with rigorous error bounds. The second project, Fungrim, is a database of knowledge about mathematical functions represented in symbolic form. It is intended to function both as a traditional reference work and as a software library to support symbolic-numeric methods for problems involving transcendental functions. I will explain a few central algorithmic ideas and explain the research goals of these projects.

Partial differential equations (PDEs) are among the most universal tools used in modelling problems in nature and man-made complex systems. For example, stochastic PDEs are a fundamental ingredient in models for nonlinear filtering problems in chemical engineering and weather forecasting, deterministic Schroedinger PDEs describe the wave function in a quantum physical system, deterministic Hamiltonian-Jacobi-Bellman PDEs are employed in operations research to describe optimal control problems where companys aim to minimise their costs, and deterministic Black-Scholes-type PDEs are highly employed in portfolio optimization models as well as in state-of-the-art pricing and hedging models for financial derivatives. The PDEs appearing in such models are often high-dimensional as the number of dimensions, roughly speaking, corresponds to the number of all involved interacting substances, particles, resources, agents, or assets in the model. For instance, in the case of the above mentioned financial engineering models the dimensionality of the PDE often corresponds to the number of financial assets in the involved hedging portfolio. Such PDEs can typically not be solved explicitly and it is one of the most challenging tasks in applied mathematics to develop approximation algorithms which are able to approximatively compute solutions of high-dimensional PDEs. Nearly all approximation algorithms for PDEs in the literature suffer from the so-called "curse of dimensionality" in the sense that the number of required computational operations of the approximation algorithm to achieve a given approximation accuracy grows exponentially in the dimension of the considered PDE. With such algorithms it is impossible to approximatively compute solutions of high-dimensional PDEs even when the fastest currently available computers are used. In the case of linear parabolic PDEs and approximations at a fixed space-time point, the curse of dimensionality can be overcome by means of Monte Carlo approximation algorithms and the Feynman-Kac formula. In this talk we introduce new nonlinear Monte Carlo algorithms for high-dimensional nonlinear PDEs. We prove that such algorithms do indeed overcome the curse of dimensionality in the case of a general class of semilinear parabolic PDEs and we thereby prove, for the first time, that a general semilinear parabolic PDE with a nonlinearity depending on the PDE solution can be solved approximatively without the curse of dimensionality.

The Landau-DeGennes Q-model of uniaxial nematic liquid crystals seeks a rank-one

traceless tensor Q that minimizes a Frank-type energy plus a double well potential

that confines the eigenvalues of Q to lie between -1/2 and 1. We propose a finite

element method (FEM) which preserves this basic structure and satisfies a discrete

form of the fundamental energy estimates. We prove that the discrete problem Gamma

converges to the continuous one as the meshsize tends to zero, and propose a discrete

gradient flow to compute discrete minimizers. Numerical experiments confirm the ability

of the scheme to approximate configurations with half-integer defects, and to deal with

colloidal and electric field effects. This work, joint with J.P. Borthagaray and S.

Walker, builds on our previous work for the Ericksen's model which we review briefly.

As parallel computers approach Exascale (10^18 floating point operations per second), processor failure and data corruption are of increasing concern. Numerical linear algebra solvers are at the heart of many scientific and engineering applications, and with the increasing failure rates, they may fail to compute a solution or produce an incorrect solution. It is therefore crucial to develop novel parallel linear algebra solvers capable of providing correct solutions on unreliable computing systems. The common way to mitigate failures in high performance computing systems consists of periodically saving data onto a reliable storage device such as a remote disk. But considering the increasing failure rate and the ever-growing volume of data involved in numerical simulations, the state-of-the-art fault-tolerant strategies are becoming time consuming, therefore unsuitable for large-scale simulations. In this talk, we will present a  novel class of fault-tolerant algorithms that do not require any additional resources. The key idea is to leverage the knowledge of numerical properties of solvers involved in a simulation to regenerate lost data due to system failures. We will also share the lessons learned and report on the numerical properties and the performance of the new resilience algorithms.

One of the major steps in the adaptive finite element methods (AFEM) is the adaptive selection of the next partition. The process is usually governed by a strategy based on carefully chosen local error indicators and aims at convergence results with optimal rates. One can formally relate the refinement of the partitions with growing an oriented graph or a tree. Then each node of the tree/graph corresponds to a cell of a partition and the approximation of a function on adaptive partitions can be expressed trough the local errors related to the cell, i.e., the node. The total approximation error is then calculated as the sum of the errors on the leaves (the terminal nodes) of the tree/graph and the problem of finding an optimal error for a given budget of nodes is known as tree approximation. Establishing a near-best tree approximation result is a key ingredient in proving optimal convergence rates for AFEM.

The classical tree approximation problems are usually related to the so-called h-adaptive approximation in which the improvements a due to reducing the size of the cells in the partition. This talk will consider also an extension of this framework to hp-adaptive approximation allowing different polynomial spaces to be used for the local approximations at different cells while maintaining the near-optimality in terms of the combined number of degrees of freedom used in the approximation.

The problem of conformity of the resulting partition will be discussed as well. Typically in AFEM, certain elements of the current partition are marked and subdivided together with some additional ones to maintain desired properties of the partition like conformity. This strategy is often described as “mark → subdivide → complete”. The process is very well understood for triangulations received via newest vertex bisection procedure. In particular, it is proven that the number of elements in the final partition is limited by constant times the number of marked cells. This hints at the possibility to design a marking procedure that is limited only to cells of the partition whose subdivision will result in a conforming partition and therefore no completion step would be necessary. This talk will present such a strategy together with theoretical results about its near-optimal performance.

In this talk, I am going to give an introduction to operator preconditioning as a general and robust strategy to precondition linear systems arising from Galerkin discretization of PDEs or Boundary Integral Equations. Then, in order to illustrate the applicability of this preconditioning technique, I will discuss the simple case of weakly singular and hypersingular integral equations, arising from exterior Dirichlet and Neumann BVPs for the Laplacian in 3D. Finally, I will show how we can also tackle operators with a more difficult structure, like the electric field integral equation (EFIE) on screens, which models the scattering of time-harmonic electromagnetic waves at perfectly conducting bounded infinitely thin objects, like patch antennas in 3D.

We present a novel approach to the solution of time-dependent PDEs via the so-called monolithic or all-at-once formulation.

This approach will be explained for simple parabolic problems and its utility in the context of PDE constrained optimization problems will be elucidated.

The underlying linear algebra includes circulant matrix approximations of Toeplitz-structured matrices and allows for effective parallel implementation. Simple computational results will be shown for the heat equation and the wave equation which indicate the potential as a parallel-in-time method.

This is joint work with Elle McDonald (CSIRO, Australia), Jennifer Pestana (Strathclyde University, UK) and Anthony Goddard (Durham University, UK)

ABSTRACT

We approximate the solution of the stationary Stokes equations with various conforming and nonconforming inf-sup stable pairs of finite element spaces on simplicial meshes. Based on each pair, we design a discretization that is quasi-optimal and pressure robust, in the sense that the velocity H^1-error is proportional to the best H^1-error to the analytical velocity. This shows that such a property can be achieved without using conforming and divergence-free pairs. We bound also the pressure L^2-error, only in terms of the best approximation errors to the analytical velocity and the analytical pressure. Our construction can be summarized as follows. First, a linear operator acts on discrete velocity test functions, before the application of the load functional, and maps the discrete kernel into the analytical one.

Second, in order to enforce consistency, we  possibly employ a new augmented Lagrangian formulation, inspired by Discontinuous Galerkin methods.

Small block overlapping, and non-overlapping, Schwarz methods are theoretically highly attractive as multilevel smoothers for a wide variety of problems that are not amenable to point relaxation methods.  Examples include monolithic Vanka smoothers for Stokes, overlapping vertex-patch decompositions for $H(\text{div})$ and  $H(\text{curl})$ problems, along with nearly incompressible elasticity, and augmented Lagrangian schemes.

 While it is possible to manually program these different schemes,  their use in general purpose libraries has been held back by a lack   of generic, composable interfaces. We present a new approach to the   specification and development such additive Schwarz methods in PETSc  that cleanly separates the topological space decomposition from the  discretisation and assembly of the equations. Our preconditioner is  flexible enough to support overlapping and non-overlapping additive  Schwarz methods, and can be used to formulate line, and plane smoothers, Vanka iterations, amongst others. I will illustrate these new features with some examples utilising the Firedrake finite element library, in particular how the design of an approriate computational interface enables these schemes to be used as building blocks inside block preconditioners.

This is joint work with Patrick Farrell and Florian Wechsung (Oxford), and Matt Knepley (Buffalo).

Classical methods for X-ray computed tomography (CT) are based on the assumption that the X-ray source intensity is known. In practice, however, the intensity is measured and hence uncertain. Under normal circumstances, when the exposure time is sufficiently high, this kind of uncertainty typically has a negligible effect on the reconstruction quality. However, in time- or dose-limited applications such as dynamic CT, this uncertainty may cause severe and systematic artifacts known as ring artifacts.
By modeling the measurement process and by taking uncertainties into account, it is possible to derive a convex reconstruction model that leads to improved reconstructions when the signal-to-noise ratio is low. We discuss some computational challenges associated with the model and illustrate its merits with some numerical examples based on simulated and real data.

 Multiscale methods based on coupled partial differential equations defined on bulk and embedded manifolds are still poorly explored from the theoretical standpoint, although they are successfully used in applications, such as microcirculation and flow in perforated subsurface reservoirs. This work aims at shedding light on some theoretical aspects of a multiscale method consisting of coupled partial differential equations defined on one-dimensional domains embedded into three-dimensional ones. Mathematical issues arise because the dimensionality gap between the bulk and the inclusions is larger than one, named as the high dimensionality gap case. First, we show that such model derives from a system of full three-dimensional equations, by the application of a topological model reduction approach. Secondly, we rigorously analyze the problem, showing that the averaging operators applied for the model reduction introduce a regularization effect that resolves the issues due to the singularity of solutions and to the ill-posedness of restriction operators. Then, we discretize the problem by means of the finite element method and we analyze the approximation error. Finally, we exploit the structure of the model reduction technique to analyze the modeling error. This study confirms that for infinitesimally small inclusions, the modeling error vanishes.

This is a joint work with Federica Laurino, Department of Mathematics, Politecnico di Milano.

Driven by its numerous applications in computational science, the approximation of smooth, high-dimensional functions via sparse polynomial expansions has received significant attention in the last five to ten years.  In the first part of this talk, I will give a brief survey of recent progress in this area.  In particular, I will demonstrate how the proper use of compressed sensing tools leads to new techniques for high-dimensional approximation which can mitigate the curse of dimensionality to a substantial extent.  The rest of the talk is devoted to approximating functions defined on irregular domains.  The vast majority of works on high-dimensional approximation assume the function in question is defined over a tensor-product domain.  Yet this assumption is often unrealistic.  I will introduce a method, known as polynomial frame approximation, suitable for broad classes of irregular domains and present theoretical guarantees for its approximation error, stability, and sample complexity.  These results show the suitability of this approach for high-dimensional approximation through the independence (or weak dependence) of the various guarantees on the ambient dimension d.  Time permitting, I will also discuss several extensions.

When the linear system in Newton’s method is approximately solved using an iterative method we have an inexact or truncated Newton method. The outer method is Newton’s method and the inner iterations will be the iterative method. The Inexact Newton framework is now close to 30 years old and is widely used and given names like Newton-Arnoldi, Newton-CG depending on the inner iterative method. In this talk we will explore convergence properties when the outer iterative method is Gauss-Newton, the Halley method or an interior point method for linear programming problems.

We discuss the key role that bespoke linear algebra plays in modelling PDEs with random coefficients using stochastic Galerkin approximation methods. As a specific example, we consider nearly incompressible linear elasticity problems with an uncertain spatially varying Young's modulus. The uncertainty is modelled with a finite set of parameters with prescribed probability distribution.  We introduce a novel three-field mixed variational formulation of the PDE model and and  assess the stability with respect to a weighted norm. The main focus will be  the efficient solution of the associated high-dimensional indefinite linear system of equations. Eigenvalue bounds for the preconditioned system can be  established and shown to be independent of the discretisation parameters and the Poisson ratio.  We also  discuss an associated a posteriori error estimation strategy and assess proxies for the error reduction associated with selected enrichments of the approximation spaces.  We will show by example that these proxies enable the design of efficient  adaptive solution algorithms that terminate when the estimated error falls below a user-prescribed tolerance.

This is joint work with Arbaz Khan and Catherine Powell

Just like optimization needs derivatives, shape optimization needs shape derivatives. Their definition and computation is a classical subject, at least concerning first order shape derivatives. Second derivatives have been studied as well, but some aspects of their theory still remains a bit mysterious for practitioners. As a result, most algorithms for shape optimization are first order methods.

To understand this situation better and in a general way, we consider first and second order shape sensitivities of integrals on smooth submanifolds using a variant of shape differentiation. Instead of computing the second derivative as the derivative of the first derivative, we choose a one-parameter family of perturbations  and compute first and second derivatives with respect to that parameter. The result is a  quadratic form in terms of a perturbation vector field that yields a second order quadratic model of the perturbed functional, which can be used as the basis of a second order shape optimization algorithm. We discuss the structure of this derivative, derive domain expressions and Hadamard forms in a general geometric framework, and give a detailed geometric interpretation of the arising terms.

Finally, we use our results to construct a second order SQP-algorithm for shape optimization that exhibits indeed local fast convergence.

Standard quadratic programs have numerous applications and play an important role in copositivity detection. We consider reformulating a standard quadratic program as a mixed integer linear programming (MILP) problem. We propose alternative MILP reformulations that exploit the specific structure of standard quadratic programs. We report extensive computational results on various classes of instances. Our experiments reveal that our MILP reformulations significantly outperform other global solution approaches. 
This is joint work with Jacek Gondzio.

In this talk I will present two recent findings concerning the preconditioning of elliptic problems. The first result concerns preconditioning of elliptic problems with variable coefficient K by an inverse Laplacian. Here we show that there is a close relationship between the eigenvalues of the preconditioned system and K. 
The second results concern the problem on mixed form where K approaches zero. Here, we show a uniform inf-sup condition and corresponding robust preconditioning. 

In many applications requiring the solution of a linear system Ax=b, the matrix A has been shown to have a low-rank property: its off-diagonal blocks have low numerical rank, i.e., they can be well approximated by matrices of small rank. Several matrix formats have been proposed to exploit this property depending on how the block partitioning of the matrix is computed.
In this talk, I will discuss the block low-rank (BLR) format, which partitions the matrix with a simple, flat 2D blocking. I will present the main characteristics of BLR matrices, in particular in terms of asymptotic complexity and parallel performance. I will then discuss some recent advances and ongoing research on BLR matrices: their multilevel extension, their use as preconditioners for iterative solvers, the error analysis of their factorization, and finally the use of fast matrix arithmetic to accelerate BLR matrix operations.

A posteriori error estimators are a key tool for the quality assessment of given finite element approximations to an unknown PDE solution as well as for the application of adaptive techniques. Typically, the estimators are equivalent to the error up to an additive term, the so called oscillation. It is a common believe that this is the price for the `computability' of the estimator and that the oscillation is of higher order than the error. Cohen, DeVore, and Nochetto [CoDeNo:2012], however, presented an example, where the error vanishes with the generic optimal rate, but the oscillation does not. Interestingly, in this example, the local $H^{-1}$-norms are assumed to be computed exactly and thus the computability of the estimator cannot be the reason for the asymptotic overestimation. In particular, this proves both believes wrong in general. In this talk, we present a new approach to posteriori error analysis, where the oscillation is dominated by the error. The crucial step is a new splitting of the data into oscillation and oscillation free data. Moreover, the estimator is computable if the discrete linear system can essentially be assembled exactly.
 

The FEniCS system [1] allows the description of finite element discretisations of partial differential equations using a high-level syntax, and the automated conversion of these representations to working code via automated code generation. In previous work described in [2] the high-level representation is processed automatically to derive discrete tangent-linear and adjoint models. The processing of the model code at a high level eases the technical difficulty associated with management of data in adjoint calculations, allowing the use of optimal data management strategies [3].

This previous methodology is extended to enable the calculation of higher order partial differential equation constrained derivative information. The key additional step is to treat tangent-linear
equations on an equal footing with originating forward equations, and in particular to treat these in a manner which can themselves be further processed to enable the derivation of associated adjoint information, and the derivation of higher order tangent-linear equations, to arbitrary order. This enables the calculation of higher order derivative information -- specifically the contraction of a Kth order derivative against (K - 1) directions -- while still making use of optimal data management strategies. Specific applications making use of Hessian information associated with models written using the FEniCS system are presented.

[1] "Automated solution of differential equations by the finite element method: The FEniCS book", A. Logg, K.-A. Mardal, and  G. N. Wells (editors), Springer, 2012
[2] P. E. Farrell, D. A. Ham, S. W. Funke, and M. E. Rognes, "Automated derivation of the adjoint of high-level transient finite element programs", SIAM Journal on Scientific Computing 35(4), C369--C393, 2013
[3] A. Griewank, and A. Walther, "Algorithm 799: Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation", ACM Transactions on Mathematical Software 26(1), 19--45, 2000

$$
\def\curl#1{\left\{#1\right\}}
\def\vek#1{\mathbf{#1}}
$$
lll-posed problems arise often in the context of scientific applications in which one cannot directly observe the object or quantity of interest. However, indirect observations or measurements can be made, and the observable data $y$ can be represented as the wanted observation $x$ being acted upon by an operator $\mathcal{A}$. Thus we want to solve the operator equation \begin{equation}\label{eqn.Txy} \mathcal{A} x = y, \end{equation} (1) often formulated in some Hilbert space $H$ with $\mathcal{A}:H\rightarrow H$ and $x,y\in H$. The difficulty then is that these problems are generally ill-posed, and thus $x$ does not depend continuously on the on the right-hand side. As $y$ is often derived from measurements, one has instead a perturbed $y^{\delta}$ such that ${y - y^{\delta}}_{H}<\delta$. Thus due to the ill-posedness, solving (1) with $y^{\delta}$ is not guaranteed to produce a meaningful solution. One such class of techniques to treat such problems are the Tikhonov-regularization methods. One seeks in reconstructing the solution to balance fidelity to the data against size of some functional evaluation of the reconstructed image (e.g., the norm of the reconstruction) to mitigate the effects of the ill-posedness. For some $\lambda>0$, we solve \begin{equation}\label{eqn.tikh} x_{\lambda} = \textrm{argmin}_{\widetilde{x}\in H}\left\lbrace{\left\|{b - A\widetilde{x}} \right\|_{H}^{2} + \lambda \left\|{\widetilde{x}}\right\|_{H}^{2}} \right\rbrace. \end{equation} In this talk, we discuss some new strategies for treating discretized versions of this problem. Here, we consider a discreditized, finite dimensional version of (1), \begin{equation}\label{eqn.Axb} Ax =  b \mbox{ with }  A\in \mathbb{R}^{n\times n}\mbox{ and } b\in\mathbb{R}^{n}, \end{equation} which inherits a discrete version of ill conditioning from [1]. We propose methods built on top of the Arnoldi-Tikhonov method of Lewis and Reichel, whereby one builds the Krylov subspace \begin{equation}
\mathcal{K}_{j}(\vek A,\vek w) = {\rm span\,}\curl{\vek w,\vek A\vek w,\vek A^{2}\vek w,\ldots,\vek A^{j-1}\vek w}\mbox{ where } \vek w\in\curl{\vek b,\vek A\vek b}
\end{equation}
and solves the discretized Tikhonov minimization problem projected onto that subspace. We propose to extend this strategy to setting of augmented Krylov subspace methods. Thus, we project onto a sum of subspaces of the form $\mathcal{U} + \mathcal{K}_{j}$ where $\mathcal{U}$ is a fixed subspace and $\mathcal{K}_{j}$ is a Krylov subspace. It turns out there are multiple ways to do this leading to different algorithms. We will explain how these different methods arise mathematically and demonstrate their effectiveness on a few example problems. Along the way, some new mathematical properties of the Arnoldi-Tikhonov method are also proven.

Since the legendary 1972 encounter of H. Montgomery and F. Dyson at tea time in Princeton, a statistical correspondence of the non-trivial zeros of the Riemann Zeta function with eigenvalues of high-dimensional random matrices has emerged. Surrounded by many deep conjectures, there is a striking analogyto the energy levels of a quantum billiard system with chaotic dynamics. Thanks 
to extensive calculation of Riemann zeros by A. Odlyzko, overwhelming numerical evidence has been found for the quantum analogy. The statistical accuracy provided by an enormous dataset of more than one billion zeros reveals distinctive finite size effects. Using the physical analogy, a precise prediction of these effects was recently accomplished through the numerical evaluation of operator determinants and their perturbation series (joint work with P. Forrester and A. Mays, Melbourne).
 

The present work concerns the approximation of the solution map associated to the parametric Helmholtz boundary value problem, i.e., the map which associates to each (real) wavenumber belonging to a given interval of interest the corresponding solution of the Helmholtz equation. We introduce a single-point Least Squares (LS) rational Padé-type approximation technique applicable to any meromorphic Hilbert space-valued univariate map, and we prove the uniform convergence of the Padé approximation error on any compact subset of the interval of interest that excludes any pole. We also present a simplified and more efficient version, named Fast LS-Padé, applicable to Helmholtz-type parametric equations with normal operators.

The LS-Padé techniques are then employed to approximate the frequency response map associated to various parametric time-harmonic wave problems, namely, a transmission/reflection problem, a scattering problem and a problem in high-frequency regime. In all cases we establish the meromorphy of the frequency response map. The Helmholtz equation with stochastic wavenumber is also considered. In particular, for Lipschitz functionals of the solution, and their corresponding probability measures, we establish weak convergence of the measure derived from the LS-Padé approximant to the true one. Two-dimensioanl numerical tests are performed, which confirm the effectiveness of the approximation method.As of the dates

 Joint work with: Francesca Bonizzoni and  Ilaria Perugia (Uni. Vienna), Davide Pradovera (EPFL)

Random matrices now play a role in many areas of theoretical, applied, and computational mathematics. Therefore, it is desirable to have tools for studying random matrices that are flexible, easy to use, and powerful. Over the last fifteen years, researchers have developed a remarkable family of results, called matrix concentration inequalities, that balance these criteria. This talk offers an invitation to the field of matrix concentration inequalities and their applications.

We first consider multilevel Monte Carlo and stochastic collocation methods for determining statistical information about an output of interest that depends on the solution of a PDE with inputs that depend on random parameters. In our context, these methods connect a hierarchy of spatial grids to the amount of sampling done for a given grid, resulting in dramatic acceleration in the convergence of approximations. We then consider multifidelity methods for the same purpose which feature a variety of models that have different fidelities. For example, we could have coarser grid discretizations, reduced-order models, simplified physics, surrogates such as interpolants, and, in principle, even experimental data. No assumptions are made about the fidelity of the models relative to the “truth” model of interest so that unlike multilevel methods, there is no a priori model hierarchy available. However, our approach can still greatly accelerate the convergence of approximations.

In the past few decades, power grids across the world have become dependent on markets that aim to efficiently match supply with demand at all times via a variety of pricing and auction mechanisms. These markets are based on models that capture interactions between producers, transmission and consumers. Energy producers typically maximize profits by optimally allocating and scheduling resources over time. A dynamic equilibrium aims to determine prices and dispatches that can be transmitted over the electricity grid to satisfy evolving consumer requirements for energy at different locations and times. Computation allows large scale practical implementations of socially optimal models to be solved as part of the market operation, and regulations can be imposed that aim to ensure competitive behaviour of market participants.

Questions remain that will be outlined in this presentation.

Firstly, the recent explosion in the use of renewable supply such as wind, solar and hydro has led to increased volatility in this system. We demonstrate how risk can impose significant costs on the system that are not modeled in the context of socially optimal power system markets and highlight the use of contracts to reduce or recover these costs. We also outline how battery storage can be used as an effective hedging instrument.

Secondly, how do we guarantee continued operation in rarely occuring situations and when failures occur and how do we price this robustness?

Thirdly, how do we guarantee appropriate participant behaviour? Specifically, is it possible for participants to develop strategies that move the system to operating points that are not socially optimal?

Fourthly, how do we ensure enough transmission (and generator) capacity in the long term, and how do we recover the costs of this enhanced infrastructure?
 

Advances in manufacturing technologies, most prominently in additive manufacturing or 3d printing, are making it possible to fabricate highly optimised products with increasing geometric and hierarchical complexity. This talk will introduce our ongoing work on design optimisation that combines CAD-compatible geometry representations, multiresolution geometry processing techniques and immersed finite elements with classical shape and topology calculus. As example applications,the shape optimisation of mechanical structures and electromechanical components, and the topology optimisation of lattice-skin structures will be discussed.

The development of reduced order models for complex applications, offering the promise for rapid and accurate evaluation of the output of complex models under parameterized variation, remains a very active research area. Applications are found in problems which require many evaluations, sampled over a potentially large parameter space, such as in optimization, control, uncertainty quantification and applications where near real-time response is needed.

However, many challenges remain to secure the flexibility, robustness, and efficiency needed for general large-scale applications, in particular for nonlinear and/or time-dependent problems.

After giving a brief general introduction to reduced order models, we discuss developments in two different directions. In the first part, we discuss recent developments of reduced methods that conserve chosen invariants for nonlinear time-dependent problems. We pay particular attention to the development of reduced models for Hamiltonian problems and propose a greedy approach to build the basis. As we shall demonstrate, attention to the construction of the basis must be paid not only to ensure accuracy but also to ensure stability of the reduced model. Time permitting, we shall also briefly discuss how to extend the approach to include more general dissipative problems through the notion of port-Hamiltonians, resulting in reduced models that remain stable even in the limit of vanishing viscosity and also touch on extensions to Euler and Navier-Stokes equations.

The second part of the talk discusses the combination of reduced order modeling for nonlinear problems with the use of neural networks to overcome known problems of on-line efficiency for general nonlinear problems. We discuss the general idea in which training of the neural network becomes part of the offline part and demonstrate its potential through a number of examples, including for the incompressible Navier-Stokes equations with geometric variations.

This work has been done with in collaboration with B.F. Afkram (EPFL, CH), N. Ripamonti EPFL, CH) and S. Ubbiali (USI, CH).

In this talk we will introduce and analyse a class of robust numerical methods for nonlocal possibly nonlinear diffusion and convection-diffusion equations. Diffusion and convection-diffusion models are popular in Physics, Chemistry, Engineering, and Economics, and in many models the diffusion is anomalous or nonlocal. This means that the underlying “particle" distributions are not Gaussian, but rather follow more general Levy distributions, distributions that need not have second moments and can satisfy (generalised) central limit theorems. We will focus on models with nonlinear possibly degenerate diffusions like fractional Porous Medium Equations, Fast Diffusion Equations, and Stefan (phase transition) Problems, with or without convection. The solutions of these problems can be very irregular and even possess shock discontinuities. The combination of nonlinear problems and irregular solutions makes these problems challenging to solve numerically.
The methods we will discuss are monotone finite difference quadrature methods that are robust in the sense that they “always” converge. By that we mean that under very weak assumptions, they converge to the correct generalised possibly discontinuous generalised solution. In some cases we can also obtain error estimates. The plan of the talk is: 1. to give a short introduction to the models, 2. explain the numerical methods, 3. give results and elements of the analysis for pure diffusion equations, and 4. give results and ideas of the analysis for convection-diffusion equations. 
 

A very common problem in Science is that we have some Data Observations and we are interested in either approximating the function underlying the data or computing some quantity of interest about this function.  This talk will discuss what are best algorithms for such tasks and how we can evaluate the performance of any such algorithm.
 

In this talk I will review recent results on the analysis of shock-capturing-type methods applied to convection-dominated problems. The method of choice is a variant of the Algebraic Flux-Correction (AFC) scheme. This scheme has received some attention over the last two decades due to its very satisfactory numerical performance. Despite this attention, until very recently there was no stability and convergence analysis for it. Thus, the purpose of the works reviewed in this talk was to bridge that gap. The first step towards the full analysis of the method is a rewriting of it as a nonlinear edge-based diffusion method. This writing makes it possible to present a unified analysis of the different variants of it. So, minimal assumptions on the components of the method are stated in such a way that the resulting scheme satisfies the Discrete Maximum Principle (DMP) and is convergence. One property that will be discussed in detail is the linearity preservation. This property has been linked to the good performance of methods of this kind. We will discuss in detail its role and the impact of it in the overall convergence of the method. Time permitting, some results on a posteriori error estimation will also be presented. 
This talk will gather contributions with A. Allendes (UTFSM, Chile), E. Burman (UCL, UK), V. John (WIAS, Berlin), F. Karakatsani (Chester, UK), P. Knobloch (Prague, Czech Republic), and 
R. Rankin (U. of Nottingham, China).

The development of reduced order models for complex applications, offering the promise for rapid and accurate evaluation of the output of complex models under parameterized variation, remains a very active research area. Applications are found in problems which require many evaluations, sampled over a potentially large parameter space, such as in optimization, control, uncertainty quantification and applications where near real-time response is needed.

However, many challenges remain to secure the flexibility, robustness, and efficiency needed for general large-scale applications, in particular for nonlinear and/or time-dependent problems.

After giving a brief general introduction to reduced order models, we discuss developments in two different directions. In the first part, we discuss recent developments of reduced methods that conserve chosen invariants for nonlinear time-dependent problems. We pay particular attention to the development of reduced models for Hamiltonian problems and propose a greedy approach to build the basis. As we shall demonstrate, attention to the construction of the basis must be paid not only to ensure accuracy but also to ensure stability of the reduced model. Time permitting, we shall also briefly discuss how to extend the approach to include more general dissipative problems through the notion of port-Hamiltonians, resulting in reduced models that remain stable even in the limit of vanishing viscosity and also touch on extensions to Euler and Navier-Stokes equations.

The second part of the talk discusses the combination of reduced order modeling for nonlinear problems with the use of neural networks to overcome known problems of on-line efficiency for general nonlinear problems. We discuss the general idea in which training of the neural network becomes part of the offline part and demonstrate its potential through a number of examples, including for the incompressible Navier-Stokes equations with geometric variations.

This work has been done with in collaboration with B.F. Afkram (EPFL, CH), N. Ripamonti EPFL, CH) and S. Ubbiali (USI, CH).