Rutherford Appleton Laboratory, nr Didcot

TBA; hosted at RAL. 

Joint work Marta Betcke and Bolin Pan

In photoacoustic tomography (PAT) with flat sensor, we routinely encounter two types of limited data. The first is due to using a finite sensor and is especially perceptible if the region of interest is large relatively to the sensor or located farther away from the sensor. In this talk we focus on the second type caused by a varying sensitivity of the sensor to the incoming wavefront direction which can be modelled as binary i.e. by a cone of sensitivity. Such visibility conditions result, in Fourier domain, in a restriction of the data to a bowtie, akin to the one corresponding to the range of the forward operator but further narrowed according to the angle of sensitivity. 

We show how we can separate the visible and invisible wavefront directions in PAT image and data using a directional frame like Curvelets, and how such decomposition allows for decoupling of the reconstruction involving application of expensive forward/adjoint solvers from the training problem. We present fast and stable approximate Fourier domain forward and adjoint operators for reconstruction of the visible coefficients for such limited angle problem and a tailored UNet matching both the multi-scale Curvelet decomposition and the partition into the visible/invisible directions for learning the invisible coefficients from a training set of similar data.

We present a general framework for preconditioning Hermitian positive definite linear systems based on the Bregman log determinant divergence. This divergence provides a measure of discrepancy between a preconditioner and a target matrix, giving rise to

the study of preconditioners given as the sum of a Hermitian positive definite matrix plus a low-rank correction. We describe under which conditions the preconditioner minimises the $\ell^2$ condition number of the preconditioned matrix, and obtain the low-rank 

correction via a truncated singular value decomposition (TSVD). Numerical results from variational data assimilation (4D-VAR) support our theoretical results.

We also apply the framework to approximate factorisation preconditioners with a low-rank correction (e.g. incomplete Cholesky plus low-rank). In such cases, the approximate factorisation error is typically indefinite, and the low-rank correction described by the Bregman divergence is generally different from one obtained as a TSVD. We compare these two truncations in terms of convergence of the preconditioned conjugate gradient method (PCG), and show numerous examples where PCG converges to a small tolerance using the proposed preconditioner, whereas PCG using a TSVD-based preconditioner fails. We also consider matrices arising from interior point methods for linear programming that do not admit such an incomplete factorisation by default, and present a robust incomplete Cholesky preconditioner based on the proposed methodology.

The talk is based on papers with Martin S. Andersen (DTU).

We use the term shape optimization when we want to find a minimizer of an objective function that assigns real values to shapes of domains. Solving shape optimization problems can be quite challenging, especially when the objective function is constrained to a PDE, in the sense that evaluating the objective function for a given domain shape requires first solving a boundary value problem stated on that domain. The main challenge here is that shape optimization methods must employ numerical methods capable of solving a boundary value problem on a domain that changes after each iteration of the optimization algorithm.

The first part of this talk will provide a gentle introduction to shape optimization. The second part of this talk will highlight how the finite element framework leads to automated numerical shape optimization methods, as realized in the open-source library fireshape. The talk will conclude with a brief overview of some academic and industrial applications of shape optimization.

Whilst we all enjoy travelling to exciting and far-off locations, the current climate crisis is making flights less and less attractive. But is there anything we can do about this? By plotting courses that make best use of atmospheric data to minimise aircraft fuel burn, airlines can not only save money on fuel, but also reduce emissions, whilst not significantly increasing flight times. In each case the route between London Heathrow Airport and John F Kennedy Airport in New York is considered.  Atmospheric data is taken from a re-analysis dataset based on daily averages from 1st December, 2019 to 29th February, 2020.

Initially Pontryagin’s minimum principle is used to find time minimal routes between the airports and these are compared with flight times along the organised track structure routes prepared by the air navigation service providers NATS and NAV CANADA for each day.  Efficiency of tracks is measured using air distance, revealing that potential savings of between 0.7% and 16.4% can be made depending on the track flown. This amounts to a reduction of 6.7 million kg of CO2 across the whole winter period considered.

In a second formulation, fixed time flights are considered, thus reducing landing delays.  Here a direct method involving a reduced gradient approach is applied to find fuel minimal flight routes either by controlling just heading angle or both heading angle and airspeed. By comparing fuel burn for each of these scenarios, the importance of airspeed in the control formulation is established.  

Finally dynamic programming is applied to the problem to minimise fuel use and the resulting flight routes are compared with those actually flown by 9 different models of aircraft during the winter of 2019 to 2020. Results show that savings of 4.6% can be made flying east and 3.9% flying west, amounting to 16.6 million kg of CO2 savings in total.

Thus large reductions in fuel consumption and emissions are possible immediately, by planning time or fuel minimal trajectories, without waiting decades for incremental improvements in fuel-efficiency through technological advances.
 

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.

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.

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.

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.

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.