Past Computational Mathematics and Applications Seminar

Differentiable programming is the underpinning technology for the AI revolution. It allows neural networks to be programmed in very high level user code while still achieving very high performance for both the evaluation of the network and, crucially, its derivatives. The Firedrake project applies exactly the same concepts to the simulation of physical phenomena modelled with partial differential equations (PDEs). By exploiting the high level mathematical abstraction offered by the finite element method, users are able to write mathematical operators for the problem they wish to solve in Python. The high performance parallel implementations of these operators are then automatically generated, and composed with the PETSc solver framework to solve the resulting PDE. However, because the symbolic differential operators are available as code, it is possible to reason symbolically about them before the numerical evaluation. In particular, the operators can be differentiated with respect to their inputs, and the resulting derivative operators composed in forward or reverse order. This creates a differentiable programming paradigm congruent with (and compatible with) machine learning frameworks such as Pytorch and JAX. 

In this presentation, David Ham will present Firedrake in the context of differentiable programming, and show how this enables productivity, capability and performance to be combined in a unique way. I will also touch on the mechanism that enables Firedrake to be coupled with Pytorch and JAX.

Please note this talk will take place at Rutherford Appleton Laboratory, Harwell Campus, Didcot. 

In this talk I will review some recent results concerning the qualitative behaviour of symplectic integrators applied to Hamiltonian PDEs, such as the nonlinear wave equation or Schrödinger equations. 

Additionally, I will discuss the problem of numerical resonances, the existence of modified energy and the existence and stability of numerical solitons over long times. 

These are works with B. Grébert, D. Bambusi, G. Maierhofer and K. Schratz. 

Can one recover a matrix from only matrix-vector products? If so, how many are needed? We will consider the matrix recovery problem for the class of hierarchical rank-structured matrices. This problem arises in scientific machine learning, where one wishes to recover the solution operator of a PDE from only input-output pairs of forcing terms and solutions. Peeling algorithms are the canonical method for recovering a hierarchical matrix from matrix-vector products, however their recursive nature poses a potential stability issue which may deteriorate the overall quality of the approximation. Our work resolves the open question of the stability of peeling. We introduce a robust version of peeling and prove that it achieves low error with respect to the best possible hierarchical approximation to any matrix, allowing us to analyze the performance of the algorithm on general matrices, as opposed to exactly hierarchical ones. This analysis relies on theory for low-rank approximation, as well as the surprising result that the Generalized Nystrom method is more accurate than the randomized SVD algorithm in this setting. 

The Landau equation stands as one of the fundamental equations in kinetic theory and plays a key role in plasma physics. However, computing it presents significant challenges due to the complexity of the Landau operator,  the dimensionality, and the need to preserve the physical properties of the solution. In this presentation, I will introduce deep learning assisted particle methods aimed at addressing some of these challenges. These methods combine the benefits of traditional structure-preserving techniques with the approximation power of neural networks, aiming to handle high dimensional problems with minimal training. 

 Many optimal control, estimation and design problems can be formulated as so-called dynamic optimization problems, which are optimization problems with differential equations and other constraints. State-of-the-art methods based on collocation, which enforce the differential equations at only a finite set of points, can struggle to solve certain dynamic optimization problems, such as those with high-index differential algebraic equations, consistent overdetermined constraints or problems with singular arcs. We show how numerical methods based on integrating the differential equation residuals can be used to solve dynamic optimization problems where collocation methods fail. Furthermore, we show that integrated residual methods can be computationally more efficient than direct collocation.

This seminar takes place at RAL (Rutherford Appleton Lab). 

We introduce a novel approach to global optimization  via continuous-time dynamic programming and Hamilton-Jacobi-Bellman (HJB) PDEs. For non-convex, non-smooth objective functions,  we reformulate global optimization as an infinite horizon, optimal asymptotic stabilization control problem. The solution to the associated HJB PDE provides a value function which corresponds to a (quasi)convexification of the original objective.  Using the gradient of the value function, we obtain a  feedback law driving any initial guess towards the global optimizer without requiring derivatives of the original objective. We then demonstrate that this HJB control law can be integrated into other global optimization frameworks to improve its performance and robustness. 

Deflation is a technique to remove a solution to a problem so that other solutions to this problem can subsequently be found. The most prominent instance is deflation we see in eigenvalue solvers, but recent interest has been in deflation of rootfinding problems from nonlinear PDEs with many isolated solutions (spearheaded by Farrell and collaborators). 

In this talk I’ll show you recent results on deflation techniques for optimisation algorithms with many local minima, focusing on the Gauss—Newton algorithm for nonlinear least squares problems.  I will demonstrate advantages of these techniques instead of the more obvious approach of applying deflated Newton’s method to the first order optimality conditions and present some proofs that these algorithms will avoid the deflated solutions. Along the way we will see an interesting generalisation of Woodbury’s formula to least squares problems, something that should be more well known in Numerical Linear Algebra (joint work with Güttel, Nakatsukasa and Bloor Riley).

Main preprint: https://arxiv.org/abs/2409.14438. 

WoodburyLS preprint: https://arxiv.org/abs/2406.15120

There is a mystery at the heart of operator learning: how can one recover a non-self-adjoint operator from data without probing the adjoint? Current practical approaches suggest that one can accurately recover an operator while only using data generated by the forward action of the operator without access to the adjoint. However, naively, it seems essential to sample the action of the adjoint for learning time-dependent PDEs. 

In this talk, we will first explore connections with low-rank matrix recovery problems in numerical linear algebra. Then, we will show that one can approximate a family of non-self-adjoint infinite-dimensional compact operators via projection onto a Fourier basis without querying the adjoint.

We present an exponential-integrator-based multi-index Monte Carlo (MIMC) method for the weak approximation of mild solutions to semilinear stochastic partial differential equations (SPDEs). Theoretical results on multi-index coupled solutions of the SPDE are provided, demonstrating their stability and the satisfaction of multiplicative error estimates. Leveraging this theory, we develop a tractable MIMC algorithm. Numerical experiments illustrate that MIMC outperforms alternative approaches, such as multilevel Monte Carlo, particularly in low-regularity settings.

We are now able to solve many partial differential equation problems that were well beyond reach when I started in academia. Some of this success is due to computer hardware but much is due to algorithmic advances. 

I will give a personal perspective of the development of computational methodology in this area over my career thus far. 

Large Language Models (LLMs) have demonstrated impressive achievements. However, recent research has shown that their deeper layers often contribute minimally, with effectiveness diminishing as layer depth increases. This pattern presents significant opportunities for model compression. 

In the first part of this seminar, we will explore how this phenomenon can be harnessed to improve the efficiency of LLM compression and parameter-efficient fine-tuning. Despite these opportunities, the underutilization of deeper layers leads to inefficiencies, wasting resources that could be better used to enhance model performance. 

The second part of the talk will address the root cause of this ineffectiveness in deeper layers and propose a solution. We identify the issue as stemming from the prevalent use of Pre-Layer Normalization (Pre-LN) and introduce Mix-Layer Normalization (Mix-LN) with combined Pre-LN and Post-LN as a new approach to mitigate this training deficiency.

Kinetic and mean-field equations are central to understanding complex systems across fields such as classical physics, engineering, and the socio-economic sciences. Efficiently solving these equations remains a significant challenge due to their high dimensionality and the need to preserve key structural properties of the models. 

In this talk, we will focus on recent advancements in deterministic numerical methods, which provide an alternative to particle-based approaches (such as Monte Carlo or particle-in-cell methods) by avoiding stochastic fluctuations and offering higher accuracy. We will discuss strategies for designing these methods to reduce computational complexity while preserving fundamental physical properties and maintaining efficiency in stiff regimes. 
Special attention will be given to the role of these methods in addressing multi-scale problems in rarefied gas dynamics and plasma physics. Time permitting, we will also touch on emerging techniques for uncertainty quantification in these systems.

AI tools like ChatGPT, Microsoft Copilot, GitHub Copilot, Claude and even older AI-enabled tools like Grammarly and MS Word, are becoming an everyday part of our research environment.  This last-minute opening up of a seminar slot due to the unfortunate illness of our intended speaker (who will hopefully re-schedule for next term) gives us an opportunity to discuss what this means for us as researchers; what are good helpful uses of AI, and are there uses of AI which we might view as inappropriate?  Please come ready to participate with examples of things which you have done yourselves with AI tools.

In this seminar I will begin by giving an overview of some problems in stochastic simulation and uncertainty quantification. I will then outline the Multilevel Monte Carlo for situations in which accurate simulations are very costly, but it is possible to perform much cheaper, less accurate simulations.  Inspired by the multigrid method, it is possible to use a combination of these to achieve the desired overall accuracy at a much lower cost.

Quantum computers are becoming a reality and current generations of machines are already well beyond the 50-qubit frontier. However, hardware imperfections still overwhelm these devices and it is generally believed the fault-tolerant, error-corrected systems will not be within reach in the near term: a single logical qubit needs to be encoded into potentially thousands of physical qubits which is prohibitive.

Due to limited resources, in the near term, hybrid quantum-classical protocols are the most promising candidates for achieving early quantum advantage and these need to resort to quantum error mitigation techniques. I will explain the basic concepts and introduce hybrid quantum-classical protocols are the most promising candidates for achieving early quantum advantage. These have the potential to solve real-world problems---including optimisation or ground-state search---but they suffer from a large number of circuit repetitions required to extract information from the quantum state. I will finally identify the most likely areas where quantum computers may deliver a true advantage in the near term.

Bálint Koczor

Associate Professor in Quantum Information Theory

Mathematical Institute, University of Oxford

webpage

Following the 2012 breakthrough in deep learning for classification and visions problems, the last decade has seen tremendous raise of interest in machine learning in a wider mathematical research community from foundational research through field specific analysis to applications. 

As data is at the core of any inverse problem, it was a natural direction for the field to investigate how machine learning could aid various aspects of inversion yielding numerous approaches from somewhat ad-hoc but very effective like learned unrolled methods to provably convergent learned regularisers with everything in between. In this talk I will review some on these developments through a lens of the research of our group.   

In this informal talk I will review some theoretical and practical aspects related to the finite element approximation of eigenvalue problems arising from PDEs.
The review will cover elliptic eigenvalue problems and eigenvalue problems in mixed form, with particular emphasis on the Maxwell eigenvalue problem.
Other topics can be discussed depending on the interests of the audience, including adaptive schemes, approximation of parametric problems, reduced order models.
 

This seminar will introduce 2D-MoSub, a derivative-free optimization method based on the subspace method and quadratic models, specifically tackling large-scale derivative-free problems. 2D-MoSub combines 2-dimensional quadratic interpolation models and trust-region techniques to update the points and explore the 2-dimensional subspace iteratively. Its framework includes constructing the interpolation set, building the quadratic interpolation model, performing trust-region trial steps, and updating the trust-region radius and subspace. Computation details and theoretical properties will be discussed. Numerical results demonstrate the advantage of 2D-MoSub.

Short Bio:
Pengcheng Xie, PhD (Chinese Academy of Sciences), is joining Lawrence Berkeley National Laboratory as a postdoctoral scholar specializing in mathematical optimization and numerical analysis. He has developed optimization methods, including 2D-MoSub and SUSD-TR. Pengcheng has published in major journals and presented at ISMP 2024 (upcoming), ICIAM 2023, and CSIAM 2022. He received the Hua Loo-keng scholarship in 2019 and the CAS-AMSS Presidential scholarship in 2023.
 

The classical finite element method uses piecewise-polynomial function spaces satisfying continuity and boundary conditions. Hybrid finite element methods, by contrast, drop these continuity and boundary conditions from the function spaces and instead enforce them weakly using Lagrange multipliers. The hybrid approach has several numerical and implementational advantages, which have been studied over the last few decades.

In this talk, we show how the hybrid perspective has yielded new insights—and new methods—in structure-preserving numerical PDEs. These include multisymplectic methods for Hamiltonian PDEs, charge-conserving methods for the Maxwell and Yang-Mills equations, and hybrid methods in finite element exterior calculus.

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.

Bilevel optimization has been part of machine learning for over 4 decades now, although perhaps not always in an obvious way. The interconnection between the two topics started appearing more clearly in publications since about 20 years now, and in the last 10 years, the number of machine learning applications of bilevel optimization has literally exploded. This rise of bilevel optimization in machine learning has been highly positive, as it has come with many innovations in the theoretical and numerical perspectives in understanding and solving the problem, especially with the rebirth of the implicit function approach, which seemed to have been abandoned at some point.
Overall, machine learning has set the bar very high for the whole field of bilevel optimization with regards to the development of numerical methods and the associated convergence analysis theory, as well as the introduction of efficient tools to speed up components such as derivative calculations among other things. However, it remains unclear how the techniques from the machine learning—based bilevel optimization literature can be extended to other applications of bilevel programming. 
For instance, many machine learning loss functions and the special problem structures enable the fulfillment of some qualification conditions that will fail for multiple other applications of bilevel optimization. In this talk, we will provide an overview of machine learning applications of bilevel optimization while giving a flavour of corresponding solution algorithms and their limitations. 
Furthermore, we will discuss possible paths for algorithms that can tackle more complicated machine learning applications of bilevel optimization, while also highlighting lessons that can be learned for more general bilevel programs.

In this talk, we consider the problem of approximating failure regions. More specifically, given a costly computational model with random parameters and a failure condition, our objective is to determine the parameter region in which the failure condition is likely to not be satisfied. In mathematical terms, this problem can be cast as approximating the level set of a probability density function. We solve this problem by dividing it into two: 1) The design of an efficient Monte Carlo strategy for probability estimation. 2) The construction of an efficient algorithm for level-set approximation. Following this structure, this talk is comprised of two parts:

In the first part, we present a new multi-output multilevel best linear unbiased estimator (MLBLUE) for approximating expectations. The advantage of this estimator is in its convenience and optimality: Given any set of computational models with known covariance structure, MLBLUE automatically constructs a provenly optimal estimator for any (finite) number of quantities of interest. Nevertheless, the optimality of MLBLUE is tied to its optimal set-up, which requires the solution of a nonlinear optimization problem. We show how the latter can be reformulated as a semi-definite program and thus be solved reliably and efficiently.

In the second part, we construct an adaptive level-set approximation algorithm for smooth functions corrupted by noise in $\mathbb{R}^d$. This algorithm only requires point value data and is thus compatible with Monte Carlo estimators. The algorithm is comprised of a criterion for level-set adaptivity combined with an a posteriori error estimator. Under suitable assumptions, we can prove that our algorithm will correctly capture the target level set at the same cost complexity of uniformly approximating a $(d-1)$-dimensional function.

In this talk, we discuss continuous in time dynamics for the problem of approaching the set of zeros of a single-valued monotone and continuous operator V . Such problems are motivated by minimax convexconcave and, in particular, by convex optimization problems with linear constraints. The central role is played by a second-order dynamical system that combines a vanishing damping term with the time derivative of V along the trajectory, which can be seen as an analogous of the Hessian-driven damping in case the operator is originating from a potential. We show that these methods exhibit fast convergence rates for kV (z(t))k as t ! +1, where z( ) denotes the generated trajectory, and for the restricted gap function, and that z( ) converges to a zero of the operator V . For the corresponding implicit and explicit discrete time models with Nesterov’s momentum, we prove that they share the asymptotic features of the continuous dynamics.

Extensions to variational inequalities and fixed-point problems are also addressed. The theoretical results are illustrated by numerical experiments on bilinear games and the training of generative adversarial networks.