Past Numerical Analysis Group Internal Seminar

Choosing where to stop an iteration or how far to increase a model complexity parameter is a recurring problem in numerical computation and data analysis. Typical symptoms are diminishing returns, a noise-dominated floor, and overfitting---accordingly, many heuristics seek an elbow or plateau beyond which further effort is not worthwhile.  AAA rational approximation provides a sharp instance of this difficulty when constructing rational approximations from noisy data, where the error often decreases rapidly at first and then fluctuates in a noisy band.  Standard AAA has no mechanism to recognize this regime and may continue iterating until a preset degree cap is reached. We thus propose noiseChop, a noise-aware stopping rule designed to run online alongside AAA. The method is inspired by Chebfun's standardChop but is tailored to AAA by using quantities already available during the iteration---a monotone envelope of the $\infty$-norm nonlinear error and the linearized error from the Loewner least squares step.  
The method first detects evidence of stagnation and then selects an early cutoff degree that achieves good accuracy without chasing noise. Numerical tests illustrate robust behavior across several functions, sample sets, and noise levels. The method is soon to be available as an optional feature in Chebfun's AAA code.

Katherine Pearce is going to talk about: 'Randomized Algorithms for Tensor CUR Approximations in Attention Mechanisms'

Attention mechanisms are a central component of transformer models that capture contextual relationships between tokens in large language models. Although many of the underlying computations (e.g., query, key, and value embeddings in multi-head attention) are inherently multi-way, classical transformer models are built on matrix-based formulations. In this talk, we discuss several ways that tensorial structure can be imposed on and exploited in attention mechanisms of transformer models. We describe how tensor-based attention can capture higher-order contextual relationships among tokens. We then explore how randomized algorithms to compute tensor CUR decompositions may be used to accelerate computations in tensor-based attention and reduce storage requirements.

Lorenzo Lazzarino will talk about: 'Error estimations for randomized low-rank approximations'

Randomized algorithms in numerical linear algebra have proven to be effective in ameliorating issues of scalability when working with large matrices, efficiently producing accurate low-rank approximations. A key remaining challenge, however, is to efficiently assess the approximation accuracy of randomized methods without additional expensive matrix accesses.

In this talk, we discuss a posteriori error estimation strategies for randomized low-rank approximations, with a focus on estimators that can be constructed from the same data used to compute the approximation or without matrix global accesses. These can serve both as certification tools and as algorithmic building blocks, enabling adaptive approximations and informed trade-offs between accuracy and computational cost. As a motivation and a case study, we include a discussion on spectromicroscopy experiments.

India Marsden will talk about: 'Expanding the definition of a finite element: groups, complexes and software'

The finite element method is a flexible framework to discretise and solve partial differential equations which has been applied to many problems across science and engineering, for example weather modelling and battery design. A core feature of the success of the finite element method, the Ciarlet definition of the components of a finite element has been used for many years. The experience of these decades (and the subsequent implementations) has exposed several key deficiencies. In particular, Ciarlet’s definition is missing information about the global continuity of the mesh and how the degrees of freedom map to each other under the relative orientation of the mesh entities. This information is necessary to implement the finite element method, leaving scope for a new definition.

We propose a new definition to handle these issues and incorporate the constantly growing landscape of new elements. This new definition also aims to encapsulate more information about the elements, such as the symmetries, incorporating ideas from Group Theory. Through this work, we hope to produce a robust, thorough definition that allows processes such as implementation-independent serialisation of finite element data.

Alongside this new definition, we will discuss the new software FUSE, which provides a domain specific language for the definition and enables elements defined in this way to be used in high performance simulation using the finite element package Firedrake. 

Jing-Yuan Wang is going to talk about: 'A Runtime-Data-Driven Enhancement Preconditioner for PCG for a Sequence of SPD Linear Systems'
 

In this work, we propose a runtime-data-driven enhancement preconditioner for improving the convergence of a preconditioned conjugate gradient method for solving a sequence of symmetric positive definite linear systems of equations. The methodology is designed for the situation where a subset of the systems has been solved and the convergence is considered too slow. In such a situation, data generated from the solved problems (residual vectors, intermediate solution vectors, approximate error vectors) are first analyzed by an unsupervised learning algorithm as a 3-step process: (1) dimension reduction; (2) classification of the slow features; (3) construction of projections to each of the feature subspaces. Based on the results of the analysis, one or more enhancement preconditioners are constructed using projection matrices corresponding to the features extracted from the slow convergence subspaces. The enhancement preconditioners are additively incorporated into the existing preconditioners and are employed to solve other systems in the sequence. The enhancement preconditioner can be further enhanced when necessary by repeating this process. Numerical experiments for time-dependent problems, including parabolic and hyperbolic equations, and stochastic elliptic equations demonstrate that the proposed approach improves the convergence considerably for other systems in the sequence when classical preconditioners are insufficient.

Daniel Cortild is going to talk about: 'Regularization Methods for Hierarchical Programming'

We consider hierarchical variational inequality problems, or more generally, variational inequalities defined over the set of zeros of a monotone operator. This framework includes convex optimization over equilibrium constraints and equilibrium selection problems. In a real Hilbert space setting, we combine a Tikhonov regularization and a proximal penalization to develop a flexible double-loop method for which we prove asymptotic convergence and provide rate statements in terms of gap functions. Our method is flexible, and effectively accommodates a large class of structured operator splitting formulations for which fixed-point encodings are available. 

Joint work with Meggie Marschner, and Mathias Staudigl (University of Mannheim)

Speaker Jung Eun Huh will talk about: 'Adaptive preconditioning for linear least-squares problems via iterative CUR'

Large-scale linear least-squares problems arise in many areas of computational science and data analysis, where efficiency and scalability are crucial. In this talk, we introduce a randomized preconditioning framework for iterative solvers based on low-rank approximations of small sketches of the original problem. The key idea is to iteratively construct low-rank preconditioners that reshape the singular value distribution in a favourable way. By tightly coupling the preconditioning and Krylov solving phases within an iterative CUR decomposition -- a low-rank approximation built from selected of columns and rows of the original matrix -- the proposed algorithm achieves faster and earlier convergence than existing methods. The algorithm performs particularly well on problems that are large in both dimensions, as well as on sparse and ill-conditioned systems. 

This is a joint work with Coralia Cartis and Yuji Nakatsukasa.

Charles Parker II will be talking about: 'Structure-preserving finite elements and the convergence of augmented Lagrangian methods'

Problems with physical constraints, such as the incompressibility constraint for mass conservation in fluids or Gauss's laws for electric and magnetic fields, result in generalized saddle point systems. So-called structure-preserving finite elements respect the constraints pointwise, resulting in more physically accurate solutions that are typically robust with respect to some problem parameters. However, constructing these finite elements may involve complicated spaces for the Lagrange multiplier variables. Augmented Lagrangian methods (ALMs) provide one process to compute the solution without the need for an explicit basis for the Lagrange multiplier space. In this talk, we present new convergence estimates for a standard ALM method, sometimes called the iterated penalty method, applied to structure-preserving discretizations of linear saddle point systems.

Speaker Francesco Hrobat will talk about; 'Lanczos with compression for symmetric matrix Lyapunov equations'

Large-scale symmetric matrix Lyapunov equations arise in control theory, model order reduction, and the discretization of PDEs. State-of-the-art algorithms, such as standard and rational Krylov methods, aim to approximate the solution with a low-rank matrix. However, the standard polynomial Krylov method (also referred to as the Lanczos method) often converges slowly and faces a memory bottleneck as the dimension of the Lanczos basis grows. Conversely, rational Krylov alternatives, while effective for low-rank approximations, require the solution of expensive shifted linear systems involving a large coefficient matrix.

In this talk, I will present a low-memory variant of the Lanczos algorithm for solving symmetric Lyapunov equations. Our approach leverages a polynomial Krylov subspace while employing rational subspaces associated with small matrices to compress the Lanczos basis. This method accesses the large coefficient matrix exclusively through matrix-vector products and maintains fixed storage requirements. The resulting low-rank solution has a rank that is independent of the dimension of the underlying polynomial Krylov subspace.

Roy Makhlouf will talk about: 'Random Embeddings for Global Optimization: Convergence Results Beyond Low Effective Dimension'
 

Timely optimization problems are high-dimensional, calling for dimensionality reduction techniques to solve them efficiently. The random embedding strategy, which optimizes the objective along a low-dimensional subspace of the search space, is arguably the simplest possible dimensionality reduction method. Recent works quantify the probability of success of this strategy to solve the original problem by lower bounding the probability of a random subspace to intersect the set of approximate global minimizers. These works showed that, when the objective has low effective dimension (i.e., is only varying along a low-dimensional subspace of the search space), random embeddings of sufficiently large dimension solve the original high-dimensional problem with probability one. In this work, we relax the low effective dimension assumption by considering objectives with anisotropic variability, namely, Lipschitz continuous functions whose Lipschitz constant is small (though nonzero) when the function is restricted to a high-dimensional subspace. Exploiting tools from stochastic geometry, we lower bound the probability for a random subspace to intersect the set of approximate global minimizers of these objectives, hence, the probability of random embeddings to succeed in solving (approximately) the original global optimization problem. Our findings offer deeper insights into the role of the dimension of the optimization problem in this probability of success.

The computation of accurate low-rank matrix approximations is central to improving the scalability of various techniques in machine learning, uncertainty quantification, and control. Traditionally, low-rank approximations are constructed using SVD-based approaches such as truncated SVD or RandomizedSVD. Although these SVD approaches---especially RandomizedSVD---have proven to be very computationally efficient, other low-rank approximation methods can offer even greater performance. One such approach is the CUR decomposition, which forms a low-rank approximation using direct row and column subsets of a matrix. Because CUR uses direct matrix subsets, it is also often better able to preserve native matrix structures like sparsity or non-negativity than SVD-based approaches and can facilitate data interpretation in many contexts. This paper introduces IterativeCUR, which draws on previous work in randomized numerical linear algebra to build a new algorithm that is highly competitive compared to prior work: (1) It is adaptive in the sense that it takes as an input parameter the desired tolerance, rather than an a priori guess of the numerical rank. (2) It typically runs significantly faster than both existing CUR algorithms and techniques such as RandomizedSVD, in particular when these methods are run in an adaptive rank mode. Its asymptotic complexity is  $\mathcal{O}(mn + (m+n)r^2 + r^3)$ for an $m\times n$ matrix of numerical rank $r$. (3) It relies on a single small sketch from the matrix that is successively downdated as the algorithm proceeds.

Speaker Edward Tansley will talk about: 'Low-rank functions in machine learning'

Functions that vary along a low-dimensional subspace of their input space, often called multi-index or low-rank functions, frequently arise in machine learning. Understanding how such structure emerges can provide insight into the learning dynamics of neural networks. One line of work that explores how networks learn low-rank data representations is the Neural Feature Ansatz (NFA), which states that after training, the Gram matrix of the first-layer weights of a deep network is proportional to some power of the average gradient outer product (AGOP) of the network with respect to its inputs. Existing results prove this relationship for 2-layer linear networks under balanced initialization. In this work, we extend these results to general L-layer linear networks and remove the assumption of balanced initialization for networks trained with weight decay.

Irina-Beatrice Haas will talk about; 'Sharp error bounds for approximate eigenvalues and singular values from subspace methods'
 

Subspace methods are commonly used for finding approximate eigenvalues and singular values of large-scale matrices. Once a subspace is found, the Rayleigh-Ritz method (for symmetric eigenvalue problems) and Petrov-Galerkin projection (for singular values) are the de facto method for extraction of eigenvalues and singular values. In this work we derive error bounds for approximate eigenvalues obtained via the Rayleigh-Ritz process. Our bounds are quadratic in the residual corresponding to each Ritz value while also being robust to clustered Ritz values, which is a key improvement over existing results. We apply these bounds to several methods for computing eigenvalues and singular values, including Krylov methods and randomized algorithms.

Dr Jari Fowkes will talk about; 'A Very Short Introduction to Ptychographic Image Reconstruction'

I will present a very short introduction to the mathematics behind the scientific imaging technique known as ptychography, starting with a brief overview of the physics model and the various simplifications required, before moving on to the main ptychography inverse problem and the three principal classes of optimization algorithms currently being used in practice. 

Casey Garner will talk about; 'General Matrix Optimization'

Since our early days in mathematics, we have been aware of two important characteristics of a matrix, namely, its coordinates and its spectrum. We have also witnessed the growth of matrix optimization models from matrix completion to semidefinite programming; however, only recently has the question of solving matrix optimization problems with general spectral and coordinate constraints been studied. In this talk, we shall discuss recent work done to study these general matrix optimization models and how they relate to topics such as Riemannian optimization, approximation theory, and more.

In this talk, we employ a level-set method to define complex computational domains and propose a ghost nodal finite element strategy tailored for two distinct applications. In the first part, we introduce a model for a Poisson-Nernst-Planck system that accounts for the correlated motion of positive and negative ions through Coulomb interactions. For very short Debye lengths, one can adopt the so called Quasi-Neutral limit which drastically simplifies the system, reducing it to a diffusion equation for a single carriers with effective diffusion coefficient. This approach, while simplifying the mathematical model, can limit the scope of numerical simulations, as it may not capture the full range of behaviors near the Quasi-Neutral limit. Our goal is therefore to design an Asymptotic Preserving (AP)  to handle both regimes: the full system when the Debye length is small but non-negligible, and the Quasi-Neutral regime as the Debye length approaches zero. In the second part, we study the formation of biological transportation networks governed by a nonlinear elliptic equation for the pressure coupled with a reaction-diffusion parabolic equation for the conductivity tensor. We compute numerical solutions using the proposed ghost nodal finite element method, which shows that the network becomes highly intricate and its branches extend over large portions of the domain.

An inexact framework for high-order adaptive regularization methods is presented, in which approximations may be used for the pth-order tensor, based on lower-order derivatives. Between each recalculation of the pth-order derivative approximation, a high-order secant equation can be used to update the pth-order tensor as proposed in (Welzel 2022) or the approximation can be kept constant in a lazy manner. When refreshing the pth-order tensor approximation after m steps, an exact evaluation of the tensor or a finite difference approximation can be used with an explicit discretization stepsize. For all the newly adaptive regularization variants, we retrieve standard complexity bound to reach a second-order stationary point.  Discussions on the number of oracle calls for each introduced variant are also provided. When p = 2, we obtain a second-order method that uses quasi-Newton approximations with optimal number of iterations bound. 

I will present a parametric finite element formulation for structure-preserving numerical methods. The approach introduces two scalar Lagrange multipliers and evolution equations for surface energy and volume, ensuring that the resulting schemes maintain the underlying geometric and physical structures. To illustrate the method, I will discuss two applications: surface diffusion and two-phase Stokes flow. By combining piecewise linear finite elements in space with structure-preserving second-order time discretizations, we obtain fully discrete schemes of high temporal accuracy. Numerical experiments confirm that the proposed methods achieve the expected accuracy while preserving surface energy and volume.

In continuum kinetic models of quasineutral plasmas, binary collisions between particles are represented by the bilinear Fokker-Planck collision operator. In full-F kinetic models, which solve for the entire particle probability distribution function, it is important to correctly capture this operator, which pushes the system towards thermodynamic equilibrium. We show a multi-species, conservative, finite element implementation of this operator, using the continuum Galerkin representation, in the Julia programming language. A Jacobian-free-Newton-Krylov solver is used to implement a backward-Euler time advance. We present several example problems that demonstrate the performance of the implementation, and we speculate on future applications.

In the numerical simulation of incompressible flows and elastic materials, it is often desirable to design discretisation schemes that preserve key structural properties of the underlying physical model. In particular, the conservation of angular momentum plays a critical role in accurately capturing rotational effects, and is closely tied to the symmetry of the stress tensor. Classical formulations such as the Stokes equations or linear elasticity can exhibit significant discrepancies when this symmetry is weakly enforced or violated at the discrete level.

This work focuses on mixed finite element methods that impose the symmetry of the stress tensor strongly, thereby ensuring exact conservation of angular momentum in the absence of body torques and couple stresses. We systematically study the effect of this constraint in both incompressible Stokes flow and linear elasticity, including anisotropic settings inspired by liquid crystal polymer networks. Through a series of benchmark problems—ranging from rigid body motions to transversely isotropic materials—we demonstrate the advantages of angular-momentum-preserving discretisations, and contrast their performance with classical elements.

Our findings reveal that strong symmetry enforcement not only leads to more robust a priori error estimates and pressure-independent velocity approximations, but also more reliable physical predictions in scenarios where angular momentum conservation is critical.

These insights advocate for the broader adoption of structure-preserving methods in computational continuum mechanics, especially in applications sensitive to rotational invariants.

The typical energy estimate for the Navier-Stokes equations provides a bound for the gradient of the velocity; energy-stable numerical methods that preserve this estimate preserve this bound. However, the bound scales with the Reynolds number (Re) causing solutions to be numerically unstable (i.e. exhibit spurious oscillations) on under-resolved meshes. The dissipation of enstrophy on the other hand provides, in the transient 2D case, a bound for the gradient that is independent of Re.

We propose a finite-element integrator for the Navier-Stokes equations that preserves the evolution of both the energy and enstrophy, implying gradient bounds that are, in the 2D case, independent of Re. Our scheme is a mixed velocity-vorticity discretisation, making use of a discrete Stokes complex. While we introduce an auxiliary vorticity in the discretisation, the energy- and enstrophy-stability results both hold on the primal variable, the velocity; our scheme thus exhibits greater numerical stability at large Re than traditional methods.

We conclude with a demonstration of numerical results, and a discussion of the existence and uniqueness of solutions.

The nematic Helmholtz-Korteweg equation is a fourth-order scalar PDE modelling time-harmonic acoustic waves in nematic Korteweg fluids, such as nematic liquid crystals. Conforming discretizations typically require C1-conforming elements, for example the Argyris element, whose implementation is notoriously challenging - especially in three dimensions - and often demands a high polynomial degree. 
In this talk, we consider an alternative non-conforming C0-hybrid interior penalty method that is both stable and convergent for any polynomial degree greater than two. Classical C0-interior penalty methods employ an H1-conforming subspace and treat the non-conformity with respect to H2 with discontinuous Galerkin techniques. Building on this idea, we use hybridization techniques to improve the computational efficiency of the discretization. We provide a brief overview of the numerical analysis and show numerical examples, demonstrating the method's ability to capture anisotropic propagation of sound in two and three dimensions. 

High-order tensor methods for solving both convex and nonconvex optimization problems have recently generated significant research interest, due in part to the natural way in which higher derivatives can be incorporated into adaptive regularization frameworks, leading to algorithms with optimal global rates of convergence and local rates that are faster than Newton's method. On each iteration, to find the next solution approximation, these methods require the unconstrained local minimization of a (potentially nonconvex) multivariate polynomial of degree higher than two, constructed using third-order (or higher) derivative information, and regularized by an appropriate power of the change in the iterates. Developing efficient techniques for the solution of such subproblems is currently, an ongoing topic of research,  and this talk addresses this question for the case of the third-order tensor subproblem. In particular, we propose the CQR algorithmic framework, for minimizing a nonconvex Cubic multivariate polynomial with  Quartic Regularisation, by sequentially minimizing a sequence of local quadratic models that also incorporate both simple cubic and quartic terms.

The role of the cubic term is to crudely approximate local tensor information, while the quartic one provides model regularization and controls progress. We provide necessary and sufficient optimality conditions that fully characterise the global minimizers of these cubic-quartic models. We then turn these conditions into secular equations that can be solved using nonlinear eigenvalue techniques. We show, using our optimality characterisations, that a CQR algorithmic variant has the optimal-order evaluation complexity of $O(\epsilon^{-3/2})$ when applied to minimizing our quartically-regularised cubic subproblem, which can be further improved in special cases.  We propose practical CQR variants that judiciously use local tensor information to construct the local cubic-quartic models. We test these variants numerically and observe them to be competitive with ARC and other subproblem solvers on typical instances and even superior on ill-conditioned subproblems with special structure.

Spectro-microscopy is an experimental technique with great potential to science challenges such as the observation of changes over time in energy materials or environmental samples and investigations of the chemical state in biological samples. However, its application is often limited by factors like long acquisition times and radiation damage. We present two measurement strategies that significantly reduce experiment times and applied radiation doses. These strategies involve acquiring only a small subset of all possible measurements and then completing the full data matrix from the sampled measurements. The methods are data-driven, utilizing spectral and spatial importance subsampling distributions to select the most informative measurements. Specifically, we use data-driven leverage scores and adaptive randomized pivoting techniques. We explore raster importance sampling combined with the LoopASD completion algorithm, as well as CUR-based sampling where the CUR approximation also serves as the completion method. Additionally, we propose ideas to make the CUR-based approach adaptive. As a result, capturing as little as 4–6% of the measurements is sufficient to recover the same information as a conventional full scan.

Low-rank approximation of parameter-dependent matrices A(t) is an important task in the computational sciences, with applications in areas such as dynamical systems and the compression of series of images. In this talk, we introduce AdaCUR, an efficient randomised algorithm for computing low-rank approximations of parameter-dependent matrices using the CUR decomposition. The key idea of our approach is the ability to reuse column and row indices for nearby parameter values, improving efficiency. The resulting algorithm is rank-adaptive, provides error control, and has complexity that compares favourably with existing methods. This is joint work with Yuji Nakatsukasa.

Shape optimization is a rich field at the intersection of analysis, optimization, and engineering. It seeks to determine the optimal geometry of structures to minimize performance objectives, subject to physical constraints—often modeled by Partial Differential Equations (PDEs). Traditional approaches commonly assume that these constraints admit a unique solution for each candidate shape, implying a single-valued shape-to-solution map. However, many real-world structures exhibit multistability, where multiple stable configurations exist under identical physical conditions.

This research departs from the single-solution paradigm by investigating shape optimization in the presence of multiple solutions to the same PDE constraints. Focusing on a neo-Hookean hyperelastic model, we formulate an optimization problem aimed at controlling the energy jump between distinct solutions. Drawing on bifurcation theory, we develop a theoretical framework that interprets these solutions as continuous branches parameterized by shape variations. Building on this foundation, we implement a numerical optimization strategy and present numerical results that demonstrate the effectiveness of our approach.

Many applications in electromagnetism, magnetohydrodynamics, and pour media flow are well-posed in spaces from the 3D de Rham complex involving $H^1$, $H(curl)$, $H(div)$, and $L^2$. Discretizing these spaces with the usual conforming finite element spaces typically leads to discrete problems that are both structure-preserving and uniformly stable with respect to the mesh size and polynomial degree. Robust preconditioners/solvers usually require the inversion of subproblems or auxiliary problems on vertex, edge, or face patches of elements. For high-order discretizations, the cost of inverting these patch problems scales like $\mathcal{O}(p^9)$ and is thus prohibitively expensive. We propose a new set of basis functions for each of the spaces in the discrete de Rham complex that reduce the cost of the patch problems to $\mathcal{O}(p^6)$ complexity. By taking advantage of additional properties of the new basis, we propose further computationally cheaper variants of existing preconditioners. Various numerical examples demonstrate the performance of the solvers.

A beautiful example of a nonlinear eigenvalue problem is the determination of complex eigenvalues for wave scattering. This talk will show how nicely this can be done by applying AAA rational approximation to a scalarized resolvent sampled at a few real frequencies.  Even for a domain as elementary as a circle with a gap in it, such computations do not seem to have been done before. This is joint work with Oscar Bruno and Manuel Santana at Caltech.

Multicomponent fluids are mixtures of distinct chemical species (i.e. components) that interact through complex physical processes such as cross-diffusion and chemical reactions. Additional physical phenomena often must be accounted for when modelling these fluids; examples include momentum transport, thermality and (for charged species) electrical effects. Despite the ubiquity of chemical mixtures in nature and engineering, multicomponent fluids have received almost no attention from the finite element community, with many important applications remaining out of reach from numerical methods currently available in the literature. This is in spite of the fact that, in engineering applications, these fluids often reside in complicated spatial regions -- a situation where finite elements are extremely useful! In this talk, we present a novel class of high-order finite element methods for simulating cross-diffusion and momentum transport (i.e. convection) in multicomponent fluids. Our model can also incorporate local electroneutrality when the species carry electrical charge, making the numerical methods particularly desirable for simulating liquid electrolytes in electrochemical applications. We discuss challenges that arise when discretising the partial differential equations of multicomponent flow, as well as some salient theoretical properties of our numerical schemes. Finally, we present numerical simulations involving (i) the microfluidic non-ideal mixing of hydrocarbons and (ii) the transient evolution of a lithium-ion battery electrolyte in a Hull cell electrode.

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.