New directions in applied linear algebra, numerical methods for PDEs, and applications


John Pearson, University of Edinburgh

9 - 11 April 2018

ICMS, 15 South College Street, Edinburgh, EH8 9AA

The accurate numerical solution of partial differential equations and optimization problems, including through the use of linear algebra techniques for solving the resulting matrix systems, presents vast challenges for researchers in applied mathematics and engineering. In this workshop, we consider the fast and efficient discretization and solution of PDEs arising from a wide variety of scientific applications, including optimization problems where PDEs act as constraints.

This workshop will consist of presentations from a number of leading researchers, aimed to encourage an exchange of ideas to further advance the state of the art in the numerical solution of PDEs, applied linear algebra, and computational optimization.

Invited speakers:

Dr Lehel Banjai, Heriot-Watt University
Prof. Luca Bergamaschi, University of Padua
Dr Silvia Gazzola, University of Bath
Dr Ben Goddard, University of Edinburgh
Prof. Jacek Gondzio, University of Edinburgh
Dr Stefan Güttel, University of Manchester
Dr Elias Jarlebring, KTH
Dr Robert Luce, École Polytechnique Fedérale de Lausanne
Dr John Pearson, University of Edinburgh
Dr Margherita Porcelli, Università degli Studi di Firenze
Dr Mariya Ptashnyk, Heriot-Watt University
Dr Tyrone Rees, Rutherford Appleton Laboratory
Dr Carola Schönlieb, University of Cambridge
Prof. Jennifer Scott, Rutherford Appleton Laboratory & University of Reading
Prof. Zdenek Strakos, Charles University in Prague
Prof. Andy Wathen, University of Oxford
Dr Alemseged Gabrehiwot Weldeyesus, University of Edinburgh
Prof. Walter Zulehner, Johannes Kepler University Linz
Dr Kostas Zygalakis, University of Edinburgh



If you would like to participate in this workshop please follow this link to register. It is expected that the workshop will start at lunchtime on Monday 9 April and finish at around 15.30 on Wednesday 11 April. Lunch and refreshments will be provided each day. Participants – other than invited speakers – are expected to make their own arrangements for accommodation in Edinburgh.


All talks will be held in the Newhaven Lecture Theatre. The lecture theatre has a built in computer, data projector, and visualiser/document camera. In addition, there are two blackboards. The projector and one board may be used simultaneously. We advise speakers that, where possible, you bring your talk on a memory stick/USB to put onto our computer in advance of your session - either at the start of the day or during coffee/lunch breaks. ICMS staff can help with this. It is possible for you to use your own laptops but it is then your own responsibility to ensure that resolutions are changed on the laptops if necessary (ICMS will not change the resolution on our equipment). If you use a Mac we expect you to bring your own adaptor.


If you are travelling from overseas you may require an entry visa. A European visa does not guarantee entry to the UK. Please use this link to the UK Visa Site to find out if you need a visa and if so how to apply for one.


Please note that it is your responsibility to have adequate travel insurance to cover medical and other emergencies that may occur on your trip.

A taxi directly from the airport will cost approximately 20.00 to 25.00 GBP to the city centre for a one-way journey. There is also a bus service direct from the airport to the city centre which will cost 4.50 GBP single or 7.50 GBP return - the Airlink 100. This is a frequent service (every 10 minutes during peak times) and will bring you close to Waverley Railway Station, only a short walk to the workshop venue.

Lothian Buses charge £1.70 for a single, £4.00 for a day ticket. Please note that the exact fare is required and no change is given.

If travelling by train, please note that Edinburgh has several railway stations - Waverley Railway Station being the main station and closest to the workshop venue at 15 South College Street. If you alight at Edinburgh Waverley, the workshop venue is an easy 10 minute walk over North and South Bridge. The second large railway station is called Haymarket and is at the West End of the city centre. Please be aware that there is also an Edinburgh Park railway station but this is at the west side of the city some way from the city centre. For ICMS, use Edinburgh Waverley.


Monday 9 April


Registration with Lunch in the Chapterhouse, Level 1




Jacek Gondzio (University of Edinburgh)
Iterative solver for linear systems arising in interior point methods for semidefinite


Elias Jarlebring (Royal Institute of Technology, Sweden)
Large-scale time-periodic time-delay systems and Broyden's method for nonlinear eigenvalue problems


Tea/Coffee in the Chapterhouse, Level 1


Stefan Güttel (University of Manchester)
Compressing variable-coefficient exterior Helmholtz problems via RKFIT


Kostas Zygalakis (University of Edinburgh)
Explicit stabilised  Runge-Kutta methods for optimization


Margherita Porcelli (Università degli Studi di Firenze)
Preconditioned Newton-like methods for PDE-constrained optimization with sparsity constraints


Informal wine reception in the Chapterhouse, Level 1


Tuesday 10 April


Silvia Gazzola (University of Bath)
Recent advances in iterative regularization methods


Luca Bergamaschi  (University of Padua)
Low rank update of preconditioners for sequences of linear systems


Tea/Coffee in the Chapterhouse, Level 1

11:15- 11:50

Alemseged Weldeyesus (University of Edinburgh)
Interior point methods for large-scale problems arising from structural optimization


Carola-Bibiane Schönlieb (Cambridge University)
Using multiple preconditioners in a Krylov subspace method


Lunch & discussion provided in the Chapterhouse, Level 1


Walter Zulehner (Johannes Kepler University Linz)
On a new mixed formulation of Kirchhoff and Reissner-Mindlin plates


Robert Luce (École polytechnique fédérale de Lausanne & Gurobi Optimization)
A fast Cholesky factorization for Toeplitz-plus-Hankel matrices


Tea/Coffee in the Chapterhouse, Level 1


Lehel Banjai (Heriot-Watt University)
Tensor Finite Element Methods for the Fractional Laplacian 


John Pearson (University of Edinburgh)
Linear algebra for the solution of PDE-constrained optimization problems


Workshop dinner 19:00 at  Blonde 71-75 St. Leonard's St, Edinburgh, EH8 9QR


Wednesday 11 April


Jennifer Scott (Rutherford Appleton Laboratory & Reading)
The challenge of large-scale linear least-squares problems


Ben Goddard (University of Edinburgh)
Pseudospectral methods for integro-PDEs


Tea/Coffee in the Chapterhouse, Level 1


Mariya Ptashnyk (Heriot-Watt University)
Multiscale modelling and simulations of plant tissue biomechanics


Andy Wathen (Oxford University)
Preconditioning for nonsymmetric problems


Lunch & discussion provided in the Chapterhouse, Level 1


Zdenek Strakos (Charles University Prague)
Decomposition into subspaces and operator preconditioning


Tyrone Rees (Rutherford Appleton Laboratory)
Using multiple preconditioners in a Krylov subspace method


Closing remarks



Lehel Banjai
Tensor Finite Element Methods for the Fractional Laplacian
In this talk we discuss the analysis and fast implementation of numerical algorithms for the solution of fractional linear elliptic problems on a bounded domain with a homogeneous Dirichlet boundary condition. We focus on the spectral fractional Laplacian and the corresponding Caffarelli-Silvestre extension. We make use of tensor product finite element spaces with $hp$ finite elements in the extended direction and piecewise linear or $hp$ finite elements in the bounded domain and prove convergence estimates. We describe fast methods for solving the resulting linear system resulting in algorithms of almost optimal complexity. Numerical results confirm the theoretical predictions.
This is joint work with J. Melenk, R. Nochetto, E. Otárola, A. Salgado, and Ch. Schwab

Luca Bergamaschi
Low rank update of preconditioners for sequences of linear systems
PDF of abstract

Silvia Gazzola
Recent advances in iterative regularization methods
In this talk we will survey some classical and well-established Krylov methods that can be employed as regularization methods for the solution of large-scale ill-posed problems, with a particular focus on imaging applications. We will then introduce a variety of new approaches, based on the idea of introducing flexible prior-conditioning within Krylov subspaces, with the goal of incorporating nonnegativity or sparsity constraints into the solution subspace. This is a joint work with James Nagy, Julianne Chung, and Malena Sabate Landman

Ben Goddard
Pseudospectral methods for integro-PDEs
Motivated by applications in statistical mechanics and fluid dynamics, we propose a novel, efficient pseudospectral collocation scheme for non-local, non-linear, integro-PDEs. In particular, we compute the non-local, integral terms in real space with the help of a specialised Gauss quadrature. Due to the exponential accuracy of the quadrature, and a convenient choice of collocation points near interfaces, we can use grids with a significantly lower number of nodes than most other reported methods.  In addition, the method can be directly applied to both (non-periodic) bounded and unbounded domains.  We demonstrate the effectiveness of the scheme by applying it to examples from soft matter and complex fluids.
Joint work with Serafim Kalliadasis, Andreas Nold, Nikos Savva and Peter Yatsyshin.

Jacek Gondizo
Iterative solver for linear systems arising in interior point methods for semidefinite programming
Interior point methods for Semidefinite Programming (SDP) face a difficult linear algebra subproblem. We propose an iterative scheme to solve the Newton equation system arising in SDP. It relies on a new preconditioner which exploits well the sparsity of matrices. Theoretical insights into the method will be provided. The computational results for the method applied to large maximum cut and matrix completion problems will be reported. 
This is a joint work with Stefania Bellavia and Margherita Porcelli. 

Stefan Guetell
Compressing variable-coefficient exterior Helmholtz problems via RKFIT
The efficient discretization of Helmholtz problems on unbounded domains is a challenging task, in particular, when the wave medium is non-homogeneous. We present a new numerical approach for compressing finite difference discretizations of such problems, thereby giving rise to efficient perfectly matched layers for non-homogeneous media. Our approach is based on the solution of a non-linear rational least squares problem using the RKFIT method. We show how the solution of this least squares problem can be converted into an accurate finite difference grid within a rational Krylov framework. This is joint work with Vladimir Druskin (Schlumberger-Doll Research, Boston) and Leonid Knizhnerman (Central Geophysical Expedition, Moscow).

Elias Jarlebring
Large-scale time-periodic time-delay systems and Broyden's method for nonlinear eigenvalue problems
Broyden's method is a general iterative method commonly used for nonlinear systems of equations, when very little information is available about the problem. We develop an approach based on Broyden's method for nonlinear eigenvalue problems and adapt it to the analysis of a partial differential equation coupled with a particular type of time-delay systems. The time-delay system is large-scale and periodic, i.e., $x'(t)=A(t)x(t)+B(t)x(t-\tau)$, where $A(t),B(t)\in\mathbb{R}^{n\times n}$ are periodic and large, i.e., $n\gg 1$ stemming from the discretization of the PDE. We adapt our approach by using that the characteristic equation of the time-periodic system is a nonlinear eigenvalue problem where the matrix-vector action can be computed by solving an ordinary differential equation. Accurate solution of this differential equation is computationally expensive. The flexibility of our Broyden approach allows us to naturally handle the trade-off between accuracy and computation time for this problem. Further improvements of the algorithm is possible by exploiting the structure of the Jacobian matrix and allows us to incorporate it into the algorithm to improve convergence. The algorithm exhibits local superlinear convergence for simple eigenvalues, and we characterize the convergence. A specific problem in machine tool milling, coupled with a the discretization of a partial-differential equation modeling vibrations in the workpiece, is used to illustrate the applicability of the approach.

Robert Luce
A fast Cholesky factorization for Toeplitz-plus-Hankel matrices
Motivated by a problem in seismic reflection tomography, we consider the problem of computing the Cholesky factorization of a symmetric positive definite Toeplitz-plus-Hankel (TpH) matrix.  We give a simple and elementary derivation of a generalized displacement operator for such matrices, which yields a displacement rank of four.  Hence the generalized Schur algorithm can be used to compute the Cholesky factor with quadratic operation complexity instead of a cubic one.  We present an efficient implementation of this algorithm, as well as an extension to the block TpH case.

John Pearson
Linear algebra for the solution of PDE-constrained optimization problems
In this talk we examine the numerical solution of optimization problems constrained by systems of partial differential equations. Specifically, we build iterative methods for solving the matrix systems obtained upon discretization, accelerated by suitably constructed preconditioners. We survey applications of such solvers to a number of scientific processes, such as fluid flow control, pattern formation, and imaging problems.

Mariya Ptashnyk
Multiscale Modelling and Simulations of Plant Tissue Biomechanics
Analysis of interactions between mechanical properties and chemical processes, which influence the elasticity and extensibility of plant tissues, is important for better understanding of plant growth and development. Plant tissues are composed of cells surrounded by cell walls and connected to each other by a cross-linked pectin network of middle lamella. The main feature of plant cells are their walls, which must be strong to resist high internal hydrostatic pressure (turgor pressure) and flexible to permit growth. It is supposed that calcium-pectin cross-linking chemistry is one of the main regulators of plant cell wall elasticity and extension.
In the microscopic model for plant cell wall and tissue biomechanics we will consider the influence of the microscopic structure and chemical processes on the mechanical properties of plant tissues. The interplay between the mechanics and the chemistry will be defined by assuming that the elastic properties of the plant cell walls depend on the chemical processes and chemical reactions depend on mechanical stresses within the cell walls. To analyse the macroscopic behaviour of plant tissues, the macroscopic models for plant cell wall and tissue biomechanics will be derived using homogenization techniques. Numerical solutions for the macroscopic models will demonstrate the patterns in the interactions between mechanical stresses and chemical processes.

Margherita Porcelli
Preconditioned Newton-like methods for PDE-constrained optimization with sparsity constraints
PDE-constrained optimization aims at finding optimal setups for partial differential equations so that relevant quantities are minimized. Including sparsity promoting terms in the formulation of such problems results in more practically relevant computed controls but adds more challenges to the numerical solution of these problems. The needed L^1-terms as well as additional inclusion of box control constraints require the use of specialized second-order methods. We discuss the use of both semismooth Newton methods and interior point methods and propose preconditioners for the solution of the arising linear algebra phase. Preliminary numerical experiments illustrate the behaviour of the presented methods.

Tyrone Rees
Using multiple preconditioners in a Krylov subspace method
When solving large linear systems of equations from PDEs, we often have no option but to use an iterative method.  It is generally accepted that iterative methods need to be coupled with a suitable preconditioner in order to be effective.  However, it is rarely the case that there is a clear choice of `right' preconditioner for a given problem, and we are forced to make a choice. In this talk I will describe methods for using more than one preconditioner at each iteration of a Krylov subspace method.  In addition to the potential for decreasing the number of iterations, the preconditioners can be applied in parallel, allowing a further decrease in the time required to solve the linear system.
This is joint work with Chen Greif, Daniel Szyld, and Scott MacLachlan

Carola Schönlieb
Bi-level optimisation with non-smooth lower-level problems
I will speak about PDE constrained optimisation for parameter optimisation for variational models in imaging.

Jennifer Scott
The challenge of large-scale linear least-squares problems
In recent years, a variety of preconditioners have been proposed for use in solving large-scale linear least-squares problems. These include simple diagonal preconditioning, preconditioners based on  incomplete factorizations and stationary inner iterations used with Krylov subspace methods.  In this talk, we briefly review preconditioners for which software has been made available and present a numerical evaluation of them using performance profiles and a large set of problems arising from practical applications. We highlight problems that are particularly challenging and look at how we might tackle them.
This is work with Nick Gould and Miroslav Tuma.

Zdenek Strakos
Decomposition into subspaces and operator preconditioning
We will consider linear equations in the abstract infinite-dimensional Hilbert space setting with bounded, coercive and self-adjoint operators, which can represent, e.g., boundary value problems formulated via partial differential equations. Efficient numerical solution procedures often incorporate transformation of the original problem using preconditioning. Motivated, in particular, by the works of several authors published in the early 90's, this text will present an abstract formulation of operator preconditioning based on the idea of decomposition of a Hilbert space into a finite number of (infinite-dimensional) subspaces with focusing on the concepts of norm equivalence and spectral equivalence of infinite-dimensional operators.  The goal is to describe in a concise way the common principles behind various computational techniques using infinite-dimensional function spaces.

Andy Wathen
Preconditioning for nonsymmetric problems
The convergence of Krylov subspace methods for self-adjoint(symmetric) matrix linear systems is well known to depend only on eigenvalues of coefficient matrix and leads to the robust theoretical underpinning of many effective preconditioning approaches. Thus preconditioned Conjugate Gradients for symmetric definite problems and preconditioned MINRES for symmetric indefinite problems are reliably and widely used. We will make some comments on non-symmetric problems where preconditioning approaches must necessarily be heuristic in general. We consider, in particular, a class of Toeplitz-related matrices for which some progress has been made towards reliable preconditioning with convergence guarantees.
This is joint work with Sean Hon.

Alemseged Gebrehiwot Weldeyesus
Interior point methods for large-scale problems arising from structural optimization
Interior point methods for optimization are well-known for obtaining solutions within modest number of iterations. However, every iteration requires solving linear systems and this imposes new requirements when solving large-scale problems. We study novel techniques that can improve the efficiency of interior point methods for solving two classes of large-scale structural optimization problems, the so-called truss layout optimization and free material optimization.
In truss layout optimization, we determine optimal cross-sectional areas of the bars to find the lightest truss structure that can sustain a given set of loads. We consider problems that can be formulated as a linear program, and propose an interior point method coupled with several techniques contributing to its efficiency. These include column generation, exploiting algebraic structures, and the use of iterative methods to solve the resulting linear systems.
In Free material optimization, the design variable is the material stiffness tensor which is allowed to vary freely at each point of the design domain. The problems belong to semidefinite program. We describe a special purpose interior method where its implementation utilizes the special property that the many matrix inequality constraints involve matrices of small sizes.

Walter Zulehner
On a new mixed formulation of Kirchhoff and Reissner-Mindlin plates
In 3D linear elasticity it is well-known that stress tensor fields satisfying the homogeneous equilibrium equation can be expressed in terms of the Beltrami stress functions, provided the domain is topologically simple. There is a similar result for the tensor field of bending and twisting moments in plate models, if the mid-surface of the plate is simply connected. While the stress tensor field in 3D linear elasticity can be written as a second-order differential operator applied to the Beltrami stress functions, the tensor field of bending and twisting moments in plate models is a first-order differential operator applied to some 2D vector field.
We will show how this result can be used to reformulate the Kirchhoff and the Reissner-Mindlin plate models as well-posed second-order systems. For the Kirchhoff model, the second-order system consists of three (consecutively to solve) standard second-order elliptic problems. For the Reissner-Mindlin model we obtain a decomposition into two second-order elliptic problems and one saddle point problem. These subproblems can be treated by standard discretization methods and solution techniques.

Kostas Zygalakis
Explicit stabilised Runge-Kutta methods for optimization
Explicit stabilized Runge-Kutta (RK) methods are explicit one-step methods with extended stability domains along the negative real axis. These methods are intended for large systems of ordinary differential equations originating mainly from semidiscretization in space of parabolic or hyperbolic-parabolic equations. In this talk we will explore their applicability to optimization of strongly convex functions. In particular, for quadratic problems it is possible to show rigorously that for suitable choice of algorithmic parameters the convergence rate matches the one of the conjugate gradient. In the general case, numerical investigations indicate that the convergence rate remains the same and the corresponding algorithm is able to outperform state of the art optimization algorithms such as the Nesterov's accelerated method. Time permitting we might also discuss extension of these ideas to the stochastic context.