Past Numerical Analysis Group Internal Seminar

The Ciarlet definition of a finite element has been core to our understanding of the finite element method since its inception. It has proved particularly useful in structuring the implementation of finite element software. However, the definition does not encapsulate all the details required to uniquely implement an element, meaning each user of the definition (whether a researcher or software package) must make further mathematical assumptions to produce a working system. 

The talk presents a new definition built on Ciarlet’s that addresses these concerns. The novel definition forms the core of a new piece of software in development, FUSE, which allows the users to consider the choice of finite element as part of the data they are working with. This is a new implementation strategy among finite element software packages, and we will discuss some potential benefits of the development.

We present a new ensemble framework for boosting the performance of overparameterized unfolding networks solving the compressed sensing problem. We combine a state-of-the-art overparameterized unfolding network with a continuation technique, to warm-start a crucial quantity of the said network's architecture; we coin the resulting continued network C-DEC. Moreover, for training and evaluating C-DEC, we incorporate the log-cosh loss function, which enjoys both linear and quadratic behavior. Finally, we numerically assess C-DEC's performance on real-world images. Results showcase that the combination of continuation with the overparameterized unfolded architecture, trained and evaluated with the chosen loss function, yields smoother loss landscapes and improved reconstruction and generalization performance of C-DEC, consistently for all datasets.

Inversion problems, such as full waveform inversion (FWI), based on wave propagation, are computationally costly optimization processes used in many applications, ranging from seismic imaging to brain tomography. In most of these uses, high-order methods are required for both accuracy and computational efficiency. Within finite element methods (FEM), using high(er)-order can provide accuracy and the usage of flexible meshes. However, FEM are rarely employed in connection with unstructured simplicial meshes because of the computational cost and complexity of code implementation. They are used frequently with quadrilateral or hexahedral spectral finite elements, but the mesh adaptivity on those elements has not yet been fully explored. In this work, we address these challenges by developing software that leverages accurate higher-order mass-lumped simplicial elements with a mesh-adaption parameter, allowing us to take advantage of the computational efficiency of newer mass-lumped simplicial elements together with waveform-adapted meshes and the accuracy of higher-order function spaces. We also calculate these mesh-related parameters and develop software for high-order spectral element methods, allowing mesh flexibility. We will also discuss future developments. The open-source code was implemented using the Firedrake framework and the Unified Form Language (UFL), a mathematical-based domain specific language, allowing flexibility in a wide range of wave-based problems. 

In this talk, I will present low-regularity numerical methods for nonlinear dispersive PDEs, with unfiltered schemes analyzed in Sobolev spaces and filtered schemes in discrete Bourgain spaces, offering effective handling of low-regularity and even rough solutions. I will highlight the significance of exploring structure-preserving low-regularity schemes, as this is a crucial area for further research.

A number of algorithms are now available---including Halko-Martinsson-Tropp, interpolative decomposition, CUR, generalized Nystrom, and QR with column pivoting---for computing a low-rank approximation of matrices. Some methods come with extremely strong guarantees, while others may fail with nonnegligible probability. We present methods for efficiently estimating the error of the approximation for a specific instantiation of the methods. Such certificate allows us to execute "responsibly reckless" algorithms, wherein one tries a fast, but potentially unstable, algorithm, to obtain a potential solution; the quality of the solution is then assessed in a reliable fashion, and remedied if necessary. This is joint work with Gunnar Martinsson. 

Time permitting, I will ramble about other topics in Randomised NLA. 

An adaptive regularization algorithm for unconstrained nonconvex optimization is presented in
which the objective function is never evaluated, but only derivatives are used and without prior knowledge of Lipschitz constant.  This algorithm belongs to the class of adaptive regularization methods, for which optimal worst-case complexity results are known for the standard framework where the objective function is evaluated. It is shown in this paper that these excellent complexity bounds are also valid for the new algorithm. Theoretical analysis of both exact and stochastic cases are discussed and  new probabilistic conditions on tensor derivatives are proposed.  Initial experiments on large binary classification highlight the merits of our method.

High-order tensor methods employing local Taylor approximations have attracted considerable attention for convex and nonconvex optimisation. The pth-order adaptive regularisation (ARp) approach builds a local model comprising a pth-order Taylor expansion and a (p+1)th-order regularisation term, delivering optimal worst-case global and local convergence rates. However, for p≥2, subproblem minimisation can yield multiple local minima, and while a global minimiser is recommended for p=2, effectively identifying a suitable local minimum for p≥3 remains elusive.
This work extends interpolation-based updating strategies, originally proposed for p=2, to cases where p≥3, allowing the regularisation parameter to adapt in response to interpolation models. Additionally, it introduces a new prerejection mechanism to discard unfavourable subproblem minimisers before function evaluations, thus reducing computational costs for p≥3.
Numerical experiments, particularly on Chebyshev-Rosenbrock problems with p=3, indicate that the proper use of different minimisers can significantly improve practical performance, offering a promising direction for designing more efficient high-order methods.

The widespread need to solve large-scale linear systems has sparked a growing interest in randomized techniques. One such class of techniques is known as iterative random sketching methods (e.g., Randomized Block Kaczmarz and Randomized Block Coordinate Descent). These methods "sketch" the linear system to generate iterative, easy-to-compute updates to a solution. By working with sketches, these methods can often enable more efficient memory operations, potentially leading to faster performance for large-scale problems. Unfortunately, tracking the progress of these methods still requires computing the full residual of the linear system, an operation that undermines the benefits of the solvers. In practice, this cost is mitigated by occasionally computing the full residual, typically after an epoch. However, this approach sacrifices real-time progress tracking, resulting in wasted computations. In this talk, we use statistical techniques to develop a progress estimation procedure that provides inexpensive, accurate real-time progress estimates at the cost of a small amount of uncertainty that we effectively control.

Multicomponent flows, i.e. mixtures, can be modeled effectively using the Onsager-Stefan-Maxwell (OSM) equations. The OSM equations can account for a wide variety of phenomena such as diffusive, convective, non-ideal mixing, thermal, pressure and electrochemical effects for steady and transient multicomponent flows. I will first introduce the general OSM framework and a finite element discretisation for multicomponent diffusion of ideal gasses. Then I will show two ways of preconditioning the multicomponent diffusion problem for various boundary conditions. Time permitting, I will also discuss how this can be extended to the non-ideal, thermal, and nonisobaric settings.

Tensor methods are methods for unconstrained continuous optimization that can incorporate derivative information of up to order p > 2 by computing a step based on the pth-order Taylor expansion at each iteration. The most important among them are regularization-based tensor methods which have been shown to have optimal worst-case iteration complexity of finding an approximate minimizer. Moreover, as one might expect, this worst-case complexity improves as p increases, highlighting the potential advantage of tensor methods. Still, the global complexity results only guarantee pessimistic sublinear rates, so it is natural to ask how local rates depend on the order of the Taylor expansion p. In the case of strongly convex functions and a fixed regularization parameter, the answer is given in a paper by Doikov and Nesterov from 2022: we get pth-order local convergence of function values and gradient norms. 
The value of the regularization parameter in their analysis depends on the Lipschitz constant of the pth derivative. Since this constant is not usually known in advance, adaptive regularization methods are more practical. We extend the local convergence results to locally strongly convex functions and fully adaptive methods. 
We discuss how for p > 2 it becomes crucial to select the "right" minimizer of the regularized local model in each iteration to ensure all iterations are eventually successful. Counterexamples show that in particular the global minimizer of the subproblem is not suitable in general. If the right minimizer is used, the pth-order local convergence is preserved, otherwise the rate stays superlinear but with an exponent arbitrarily close to one depending on the algorithm parameters.

The Parker conjecture, which explores whether magnetic fields in perfectly conducting plasmas can develop tangential discontinuities during magnetic relaxation, remains an open question in astrophysics. Helicity conservation provides a topological barrier against topologically nontrivial initial data relaxing to a trivial solution. Preserving this mechanism is therefore crucial for numerical simulation.  

This paper presents an energy- and helicity-preserving finite element discretization for the magneto-frictional system for investigating the Parker conjecture. The algorithm enjoys a discrete version of the topological mechanism and a discrete Arnold inequality. 
We will also discuss extensions to domains with nontrivial topology.

This is joint work with Prof Patrick Farrell, Dr Kaibo Hu, and Boris Andrews

In Bayesian inverse problems, it is common to consider several hyperparameters that define the prior and the noise model that must be estimated from the data. In particular, we are interested in linear inverse problems with additive Gaussian noise and Gaussian priors defined using Matern covariance models. In this case, we estimate the hyperparameters using the maximum a posteriori (MAP) estimate of the marginalized posterior distribution. 

However, this is a computationally intensive task since it involves computing log determinants.  To address this challenge, we consider a stochastic average approximation (SAA) of the objective function and use the preconditioned Lanczos method to compute efficient function evaluation approximations. 

We can therefore compute the MAP estimate of the hyperparameters efficiently by building a preconditioner which can be updated cheaply for new values of the hyperparameters; and by leveraging numerical linear algebra tools to reuse information efficiently for computing approximations of the gradient evaluations.  We demonstrate the performance of our approach on inverse problems from tomography. 

In recent years, some progress has been made in the development of finite element complexes, particularly in the discretization of BGG complexes in two and three dimensions, including Hessian complexes, elasticity complexes, and divdiv complexes. In this talk, I will discuss distributional complexes in two and three dimensions. These complexes are simply constructed using geometric concepts such as vertices, edges, and faces, and they share the same cohomology as the complexes at the continuous level, which reflects that the discretization is structure preserving. The results can be regarded as a tensor generalization of the Whitney forms of the finite element exterior calculus. This talk is based on joint work with Snorre Christiansen (Oslo), Kaibo Hu (Edinburgh), and Qian Zhang (Michigan).

We introduce a new methodology based on the multirevolution idea for constructing integrators for stochastic differential equations in the situation where the fast oscillations themselves are driven by a Stratonovich noise. Applications include in particular highly-oscillatory Kubo oscillators and spatial discretizations of the nonlinear Schrödinger equation with fast white noise dispersion. We construct a method of weak order two with computational cost and accuracy both independent of the stiffness of the oscillations. A geometric modification that conserves exactly quadratic invariants is also presented. If time allows, we will discuss ongoing work on uniformly accurate methods for such systems. This is a joint work with Gilles Vilmart.

The backward error analysis is an important part of the perturbation theory and it is particularly useful for the study of the reliability of the numerical methods. We focus on the backward error for nonlinear eigenvalue problems. In this talk, the matrix-valued function is given as a linear combination of scalar functions multiplying matrix coefficients, and the perturbation is done on the coefficients. We provide theoretical results about the backward error of a set of approximate eigenpairs. Indeed, small backward errors for separate eigenpairs do not imply small backward errors for a set of approximate eigenpairs. In this talk, we provide inexpensive upper bounds, and a way to accurately compute the backward error by means of direct computations or through Riemannian optimization. We also discuss how the backward error can be determined when the matrix coefficients of the matrix-valued function have particular structures (such as symmetry, sparsity, or low-rank), and the perturbations are required to preserve them. For special cases (such as for symmetric coefficients), explicit and inexpensive formulas to compute the perturbed matrix coefficients are also given. This is a joint work with Leonardo Robol (University of Pisa).

Dispersive nonlinear partial differential equations can be used to describe a range of physical systems, from water waves to spin states in ferromagnetism. The numerical approximation of solutions with limited differentiability (low-regularity) is crucial for simulating fascinating phenomena arising in these systems including emerging structures in random wave fields and dynamics of domain wall states, but it poses a significant challenge to classical algorithms. Recent years have seen the development of tailored low-regularity integrators to address this challenge. Inherited from their description of physicals systems many such dispersive nonlinear equations possess a rich geometric structure, such as a Hamiltonian formulation and conservation laws. To ensure that numerical schemes lead to meaningful results, it is vital to preserve this structure in numerical approximations. This, however, results in an interesting dichotomy: the rich theory of existent structure-preserving algorithms is typically limited to classical integrators that cannot reliably treat low-regularity phenomena, while most prior designs of low-regularity integrators break geometric structure in the equation. In this talk, we will outline recent advances incorporating structure-preserving properties into low-regularity integrators. Starting from simple discussions on the nonlinear Schrödinger and the Korteweg–de Vries equation we will discuss the construction of such schemes for a general class of dispersive equations before demonstrating an application to the simulation of low-regularity vortex filaments. This is joint work with Yvonne Alama Bronsard, Valeria Banica, Yvain Bruned and Katharina Schratz.

The interplay between stochastic processes and optimal control has been extensively explored in the literature. With the recent surge in the use of diffusion models, stochastic processes have increasingly been applied to sample generation. This talk builds on the log transform, known as the Cole-Hopf transform in Brownian motion contexts, and extends it within a more abstract framework that includes a linear operator. Within this framework, we found that the well-known relationship between the Cole-Hopf transform and optimal transport is a particular instance where the linear operator acts as the infinitesimal generator of a stochastic process. We also introduce a novel scenario where the linear operator is the adjoint of the generator, linking to Bayesian inference under specific initial and terminal conditions. Leveraging this theoretical foundation, we develop a new algorithm, named the HJ-sampler, for Bayesian inference for the inverse problem of a stochastic differential equation with given terminal observations. The HJ-sampler involves two stages: solving viscous Hamilton-Jacobi (HJ) partial differential equations (PDEs) and sampling from the associated stochastic optimal control problem. Our proposed algorithm naturally allows for flexibility in selecting the numerical solver for viscous HJ PDEs. We introduce two variants of the solver: the Riccati-HJ-sampler, based on the Riccati method, and the SGM-HJ-sampler, which utilizes diffusion models. Numerical examples demonstrate the effectiveness of our proposed methods. This is an ongoing joint work with Zongren Zou, Jerome Darbon, and George Em Karniadakis.

Fourth-order elliptic problems arise in a variety of applications from thin plates to phase separation to liquid crystals. A conforming Galerkin discretization requires a finite dimensional subspace of H2, which in turn means that conforming finite element subspaces are C1-continuous. In contrast to standard H1-conforming C0-elements, C1-elements, particularly those of high order, are less understood from a theoretical perspective and are not implemented in many existing finite element codes. In this talk, we address the implementation of the elements. In particular, we present algorithms that compute C1-finite element approximations to fourth-order elliptic problems and which only require elements with at most C0-continuity. The algorithms are suitable for use in almost all standard finite element packages. Iterative methods and preconditioners for the subproblems in the algorithm will also be presented.

We present a spectral method that converges exponentially for a variety of fractional integral equations on a closed interval. The method uses an orthogonal fractional polynomial basis that is obtained from an appropriate change of variable in classical Jacobi polynomials. For a problem arising from time-fractional heat and wave equations, we elaborate the complexities of three spectral methods, among which our method is the most performant due to its superior stability. We present algorithms for building the fractional integral operators, which are applied to the orthogonal fractional polynomial basis as matrices. 

The main component of flow machines is the rotor; however, there may also be stationary parts surrounding the rotor, which are the diffuser blades. In order to consider these two parts simultaneously, the most intuitive approach is to perform a transient flow simulation; however, the computational cost is relatively high. Therefore, one possible approach is the Multiple Reference Frame (MRF) approach, which considers two directly coupled zones: one for the rotating reference frame (for the rotor blades) and one for the stationary reference frame (for the diffuser blades). When taking into account topology optimisation, some changes are required in order to take both rotating and stationary parts simultaneously in the design, which also leads to changes in the composition of the multi-objective function. Therefore, the topology optimisation method is formulated for MRF while also proposing this new multi-objective function. An integer variable-based optimisation algorithm is considered, with some adjustments for the MRF case. Some numerical examples are presented.

Our research explores the application of reinforcement learning (RL) strategies to solve complex combinatorial research problems, specifically the Job-shop Scheduling Problem (JSP) and the Stochastic Vehicle Routing Problem with Time Windows (SVRP). For JSP, we utilize Curriculum Learning (CL) to enhance the performance of dispatching policies. This approach addresses the significant optimality gap in existing end-to-end solutions by structuring the training process into a sequence of increasingly complex tasks, thus facilitating the handling of larger, more intricate instances. Our study introduces a size-agnostic model and a novel strategy, the Reinforced Adaptive Staircase Curriculum Learning (RASCL), which dynamically adjusts difficulty levels during training, focusing on the most challenging instances. Experimental results on Taillard and Demirkol datasets show that our approach reduces the average optimality gap to 10.46% and 18.85%, respectively.

For SVRP, we propose an end-to-end framework employing an attention-based neural network trained through RL to minimize routing costs while addressing uncertain travel costs and demands, alongside specific customer delivery time windows. This model outperforms the state-of-the-art Ant-Colony Optimization algorithm by achieving a 1.73% reduction in travel costs and demonstrates robustness across diverse environmental settings, making it a valuable baseline for future research. Both studies mark advancements in the application of machine learning techniques to operational research.

We consider the problem of approximating singular values of a matrix when provided with approximations to the leading singular vectors. In particular, we focus on the Generalized Nystrom (GN) method, a commonly used low-rank approximation, and its error in extracting singular values. Like other approaches, the GN approximation can be interpreted as a perturbation of the original matrix. Up to orthogonal transformations, this perturbation has a peculiar structure that we wish to exploit. Thus, we use the Jordan-Wieldant Theorem and similarity transformations to generalize a matrix perturbation theory result on eigenvalues of a perturbed Hermitian matrix. Finally, combining the above,  we can derive a bound on the GN singular values approximation error. We conclude by performing preliminary numerical examples. The aim is to heuristically study the sharpness of the bound, to give intuitions on how the analysis can be used to compare different approaches, and to provide ideas on how to make the bound computable in practice.