Past Computational Mathematics and Applications Seminar

We give a strongly polynomial algorithm for minimum cost generalized flow, and as a consequence, for all linear programs with at most two nonzero entries per row, or at most two nonzero entries per column. While strongly polynomial algorithms for the primal and dual feasibility problems have been known for a long time, various combinatorial approaches used for those problems did not seem to carry over to the minimum-cost variant.

Our approach is to show that the ‘subspace layered least squares’ interior point method, an earlier joint work with Allamigeon, Dadush, Loho and Natura requires only a strongly polynomial number of iterations for minimum cost generalized flow. We achieve this by bounding the straight line complexity, introduced in the same paper. The talk will give an overview of the interior point method as well as the combinatorial straight-line complexity analysis for the particular setting. This is joint work with Daniel Dadush, Zhuan Khye Koh, Bento Natura, and Neil Olver.

This talk explores recent advancements in stress and flux-based finite element methods. It focuses on addressing the limitations of traditional finite elements, in order to describe complex material behavior and engineer new metamaterials.

Stress and flux-based finite element methods are particularly useful in error estimation, laying the groundwork for adaptive refinement strategies. This concept builds upon the hypercircle theorem [1], which states that in a specific energy space, both the exact solution and any admissible stress field lie on a hypercircle. However, the construction of finite element spaces that satisfy admissible states for complex material behavior is not straightforward. It often requires a relaxation of specific properties, especially when dealing with non-symmetric stress tensors [2] or hyperelastic materials.

Alternatively, methods that directly approximate stresses can be employed, offering high accuracy of the stress fields and adherence to physical conservation laws. However, when approximating eigenvalues, this significant benefit for the solution's accuracy implies that the solution operator cannot be compact. To address this, the solution operator must be confined to a subset of the solution that excludes the stresses. Yet, due to compatibility conditions, the trial space for the other solution components typically does not yield the desired accuracy. The second part of this talk will therefore explore the Least-Squares method as a remedy to these challenges [3].

To conclude this talk, we will emphasize the integration of those methods within global solution strategies, with a particular focus on the challenges regarding model order reduction methods [4].

[1] W. Prager, J. Synge. Approximations in elasticity based on the concept of function space.

Quarterly of Applied Mathematics 5(3), 1947.

[2] FB, K. Bernhard, M. Moldenhauer, G. Starke. Weakly symmetric stress equilibration and a posteriori error estimation for linear elasticity, Numerical Methods for Partial Differential Equations 37(4), 2021.

[3] FB, D. Boffi. First order least-squares formulations for eigenvalue problems, IMA Journal of Numerical Analysis 42(2), 2023.

[4] FB, D. Boffi, A. Halim. A reduced order model for the finite element approximation of eigenvalue problems,Computer Methods in Applied Mechanics and Engineering 404, 2023.

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

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

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

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

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

This talk concerns the design and analysis of adaptive FEM-based solution strategies for partial differential equations (PDEs) with uncertain or parameter-dependent inputs. We present two conceptually different strategies: one is projection-based (stochastic Galerkin FEM) and the other is sampling-based (stochastic collocation FEM). These strategies have emerged and become popular as effective alternatives to Monte-Carlo sampling in the context of (forward) uncertainty quantification. Both stochastic Galerkin and stochastic collocation approximations are typically represented as finite (sparse) expansions in terms of a parametric polynomial basis with spatial coefficients residing in finite element spaces. The focus of the talk is on multilevel approaches where different spatial coefficients may reside in different finite element spaces and, therefore, the underlying spatial approximations are allowed to be refined independently from each other.

We start with a more familiar setting of projection-based methods, where exploiting the Galerkin orthogonality property and polynomial approximations in terms of an orthonormal basis facilitates the design and analysis of adaptive algorithms. We discuss a posteriori error estimation as well as the convergence and rate optimality properties of the generated adaptive multilevel Galerkin approximations for PDE problems with affine-parametric coefficients. We then show how these ideas of error estimation and multilevel adaptivity can be applied in a non-Galerkin setting of stochastic collocation FEM, in particular, for PDE problems with non-affine parameterization of random inputs and for problems with parameter-dependent local spatial features.

The talk is based on a series of joint papers with Dirk Praetorius (TU Vienna), Leonardo Rocchi (Birmingham), Michele Ruggeri (University of Strathclyde, Glasgow), David Silvester (Manchester), and Feng Xu (Manchester).

We develop an algorithm for parameter-free stochastic convex optimization (SCO) whose rate of convergence is only a double-logarithmic factor larger than the optimal rate for the corresponding known-parameter setting. In contrast, the best previously known rates for parameter-free SCO are based on online parameter-free regret bounds, which contain unavoidable excess logarithmic terms compared to their known-parameter counterparts. Our algorithm is conceptually simple, has high-probability guarantees, and is also partially adaptive to unknown gradient norms, smoothness, and strong convexity. At the heart of our results is a novel parameter-free certificate for the step size of stochastic gradient descent (SGD), and a time-uniform concentration result that assumes no a-priori bounds on SGD iterates.

Additionally, we present theoretical and numerical results for a dynamic step size schedule for SGD based on a variant of this idea. On a broad range of vision and language transfer learning tasks our methods performance is close to that of SGD with tuned learning rate. Also, a per-layer variant of our algorithm approaches the performance of tuned ADAM.

This talk is based on papers with Yair Carmon and Maor Ivgi.

Geometric evolutions can arise as simple models or fundamental building blocks in various applications with moving boundaries and time-dependent domains, such as grain boundaries in materials or deforming cell boundaries. Mesh-based methods require adaptation and smoothing, particularly in the case of strong deformations. We consider finite element schemes based on classical approaches for geometric evolution equations but augmented with the gradient of the Dirichlet energy or a variant of it, which is known to produce a tangential mesh movement beneficial for the mesh quality. We focus on the one-dimensional case, where convergence of semi-discrete schemes can be proved, and discuss two cases. For networks forming triple junctions, it is desirable to keep the impact of any additional, mesh smoothing terms on the geometric evolution as small as possible, which can be achieved with a perturbation approach. Regarding the elastic flow of curves, the Dirichlet energy can serve as a replacement of the usual penalty in terms of the length functional in that, modulo rescaling, it yields the same minimisers in the long run.

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

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

Optimization problems subject to PDE constraints constitute a mathematical tool that can be applied to a wide range of scientific processes, including fluid flow control, medical imaging, option pricing, biological and chemical processes, and electromagnetic inverse problems, to name a few. These problems involve minimizing some function arising from a particular physical objective, while at the same time obeying a system of PDEs which describe the process. It is necessary to obtain accurate solutions to such problems within a reasonable CPU time, in particular for time-dependent problems, for which the “all-at-once” solution can lead to extremely large linear systems.

In this talk we consider iterative methods, in particular Krylov subspace methods, to solve such systems, accelerated by fast and robust preconditioning strategies. In particular, we will survey several new developments, including block preconditioners for fluid flow control problems, a circulant preconditioning framework for solving certain optimization problems constrained by fractional differential equations, and multiple saddle-point preconditioners for block tridiagonal linear systems. We will illustrate the benefit of using these new approaches through a range of numerical experiments.

This talk is based on work with Santolo Leveque (Scuola Normale Superiore, Pisa), Spyros Pougkakiotis (Yale University), Jacek Gondzio (University of Edinburgh), and Andreas Potschka (TU Clausthal).

The need to solve linear systems of equations is ubiquitous in scientific computing. Powerful methods for preconditioning such systems have been developed in cases where we can exploit knowledge of the origin of the linear system; a recent example from the solution of systems from PDEs is the Gen-EO domain decomposition method which works well, but requires a non-trival amount of knowledge of the underlying problem to implement.  

In this talk I will present a new spectral coarse space that can be constructed in a fully-algebraic way, in contrast to most existing spectral coarse spaces, and will give a theoretical convergence result for Hermitian positive definite diagonally dominant matrices. Numerical experiments and comparisons against state-of-the-art preconditioners in the multigrid community show that the resulting two-level Schwarz preconditioner is efficient, especially for non-self-adjoint operators. Furthermore, in this case, our proposed preconditioner outperforms state-of-the-art preconditioners.

This is joint work with Hussam Al Daas, Pierre Jolivet and Jennifer Scott.

Randomized least squares is a promising method but not yet widely used in practice. We show an example of its use for finding low-rank canonical polyadic (CP) tensor decompositions for large sparse tensors. This involves solving a sequence of overdetermined least problems with special (Khatri-Rao product) structure.

In this work, we present an application of randomized algorithms to fitting the CP decomposition of sparse tensors, solving a significantly smaller sampled least squares problem at each iteration with probabilistic guarantees on the approximation errors. We perform sketching through leverage score sampling, crucially relying on the fact that the problem structure enable efficient sampling from overestimates of the leverage scores with much less work. We discuss what it took to make the algorithm practical, including general-purpose improvements.

Numerical results on real-world large-scale tensors show the method is faster than competing methods without sacrificing accuracy.

*This is joint work with Brett Larsen, Stanford University.

This work considers weighted and preconditioned GMRES. The objective is to provide a way of choosing the preconditioner and the inner product, also called weight, that ensure fast convergence. The main focus of the article is on Hermitian preconditioning (even for non-Hermitian problems).

It is indeed proposed to choose a Hermitian preconditioner H, and to apply GMRES in the inner product induced by H. If moreover, the problem matrix A is positive definite, then a new convergence bound is proved that depends only on how well H preconditions the Hermitian part of A, and on a measure of how non-Hermitian A is. In particular, if a scalable preconditioner is known for the Hermitian part of A, then the proposed method is also scalable. I will also illustrate this result numerically.

I will reflect on my time as Professor of Numerical Analysis.

Tensor decomposition express a tensor as a linear combination of elementary tensors. They have applications in chemometrics, computer science, machine learning, psychometrics, and signal processing. Their uniqueness properties render them suitable for data analysis tasks in which the elementary tensors are the quantities of interest. However, in applications, the idealized mathematical model is corrupted by measurement errors. For a robust interpretation of the data, it is therefore imperative to quantify how sensitive these elementary tensors are to perturbations of the whole tensor. I will give an overview of recent results on the condition number of tensor decompositions, established with my collaborators C. Beltran, P. Breiding, and N. Dewaele.

Data-driven reduced-order modeling aims at constructing models describing the underlying dynamics of unknown systems from measurements. This has become an increasingly preeminent discipline in the last few years. It is an essential tool in situations when explicit models in the form of state space formulations are not available, yet abundant input/output data are, motivating the need for data-driven modeling. Depending on the underlying physics, dynamical systems can inherit differential structures leading to specific physical interpretations. In this work, we concentrate on systems that are described by differential equations and possess linear dynamics. Extensions to more complicated, nonlinear dynamics are also possible and will be briefly covered here if time permits.

The methods developed in our study use rational approximation based on Loewner matrices. Starting with the approach by Antoulas and Anderson in '86, and moving forward to the one by Mayo and Antoulas in '07, the Loewner framework (LF) has become an established methodology in the model reduction and reduced-order modeling community. It is a data-driven approach in the sense that what is needed to compute the reduced models is solely data, i.e., samples of the system's transfer function. As opposed to conventional intrusive methods that require an actual large-scale model to reduce (described by many differential equations), the LF only needs measurements in compressed format. In the former category of approaches, we mention balanced truncation (BT), arguably one of the most prevalent model reduction methods. Introduced in the early 80s, this method constructs reduced-order models (ROMs) by using balancing and truncating steps (with respect to classical system theory concepts such as controllability and observability). We show that BT can be reinterpreted as a data-driven approach, by using again the Loewner matrix as a central ingredient. By making use of quadrature approximations of certain system theoretical quantities (infinite Gramian matrices), a novel method called QuadBT (quadrature-based BT) is introduced by G., Gugercin, and Beattie in '22. We show parallels with the LF and, if time permits, certain recent extensions of QuadBT. Finally, all theoretical considerations are validated on various numerical test cases.

On supercomputers that exist today, achieving even close to the peak performance is incredibly difficult if not impossible for many applications. Techniques designed to improve the performance of matrix computations - making computations less expensive by reorganizing an algorithm, making intentional approximations, and using lower precision - all introduce what we can generally call ``inexactness''. The questions to ask are then:

1. With all these various sources of inexactness involved, does a given algorithm still get close enough to the right answer?
2. Given a user constraint on required accuracy, how can we best exploit and balance different types of inexactness to improve performance?

Studying the combination of different sources of inexactness can thus reveal not only limitations, but also new opportunities for developing algorithms for matrix computations that are both fast and provably accurate. We present few recent results toward this goal, icluding mixed precision randomized decompositions and mixed precision sparse approximate inverse preconditioners.

Linear poroelasticity models have important applications in biology and geophysics. In particular, the well-known Biot consolidation model describes the coupled interaction between the linear response of a porous elastic medium saturated with fluid and a diffusive fluid flow within it, assuming small deformations. This is the starting point for modeling human organs in computational medicine and for modeling the mechanics of permeable
rock in geophysics. Finite element methods for Biot’s consolidation model have been widely studied over the past four decades.
In the first part of the talk, we discuss a posteriori error estimators for locking-free mixed finite element approximation of Biot’s consolidation model. The simplest of these is a conventional residual-based estimator. We establish bounds relating the estimated and true errors, and show that these are independent of the physical parameters. The other two estimators require the solution of local problems. These local problem estimators are also shown to be reliable, efficient and robust. Numerical results are presented that
validate the theoretical estimates, and illustrate the effectiveness of the estimators in guiding adaptive solution algorithms.
In the second part of talk, we discuss a novel locking-free stochastic Galerkin mixed finite element method for the Biot consolidation model with uncertain Young’s modulus and hydraulic conductivity field. After introducing a five-field mixed variational formulation of the standard Biot consolidation model, we discuss stochastic Galerkin mixed finite element approximation, focusing on the issue of well-posedness and efficient linear algebra for the discretized system. We introduce a new preconditioner for use with MINRES and
establish eigenvalue bounds. Finally, we present specific numerical examples to illustrate the efficiency of our numerical solution approach.

Finally, we discuss some remarks related to non-conforming approximation of Biot’s consolidation model.

References:
1. A. Khan, D. J. Silvester, Robust a posteriori error estimation for mixed finite
element approximation of linear poroelsticity, IMA Journal of Numerical Analysis, Oxford University Press, 41 (3), 2021, 2000-2025.
2. A. Khan, C. E. Powell, Parameter-robust stochastic Galerkin approxination for linear poroelasticity with uncertain inputs, SIAM Journal on Scientific Computing (SISC), 43 (4), 2021, B855-B883.
3. A. Khan, P. Zanotti, A nonsymmetric approach and a quasi-optimal and robust discretization for the Biot’s model. Mathematics of Computation, 91 (335), 2022, 1143-1170.
4. V. Anaya, A. Khan, D. Mora, R. Ruiz-Baier, Robust a posteriori error analysis for rotation-based formulations of the elasticity/poroelasticity coupling, SIAM Journal
on Scientific Computing (SISC), 2022.

We propose to use stochastic Riemannian coordinate descent on the Stiefel manifold for deep neural network training. The algorithm rotates successively two columns of the matrix, an operation that can be efficiently implemented as a multiplication by a Givens matrix. In the case when the coordinate is selected uniformly at random at each iteration, we prove the convergence of the proposed algorithm under standard assumptions on the loss function, stepsize and minibatch noise. Experiments on benchmark deep neural network training problems are presented to demonstrate the effectiveness of the proposed algorithm.

Correlation operators are used in data assimilation algorithms
to weight the contribution of prior and observation information.
Efficient implementation of these operators is therefore crucial for
operational implementations. Diffusion-based correlation operators are popular in ocean data assimilation, but can require a large number of serial matrix-vector products. An all-at-once formulation removes this requirement, and offers the opportunity to exploit modern computer architectures. High quality preconditioners for the all-at-once approach are well-known, but impossible to apply in practice for the
high-dimensional problems that occur in oceanography. In this talk we
consider a nested preconditioning approach which retains many of the
beneficial properties of the ideal analytic preconditioner while
remaining affordable in terms of memory and computational resource.

Polynomial splines over simplicial meshes in R^n (triangulations in 2D, tetrahedral meshes in 3D, and so on) sometimes have extra orders of smoothness at a vertex. This property is known as supersmoothness, and plays a role both in the construction of macroelements and in the finite element method.
Supersmoothness depends both on the number of simplices that meet at the vertex and their geometric configuration.

In this talk we review what is known about supersmoothness of polynomial splines and then discuss the more general setting of splines whose individual pieces are any infinitely smooth functions.

This is joint work with Kaibo Hu.

Skeletal joints are often modelled as two adjacent layers of poroviscoelastic cartilage that are permitted to slide past each other.  The talk will begin by outlining a mathematical model that may be used, focusing on two unusual features of the model: (i) the solid component of the poroviscoelastic body has a charged surface that ionises the fluid within the pores, generating a swelling pressure; and (ii) appropriate conditions are required at the interface between the two adjacent layers of cartilage.  The remainder of the talk will then address various theoretical and practical issues in computing a finite element solution of the governing equations.

In this talk, we first review the de Rham complex and the finite element exterior calculus, a cohomological framework for structure-preserving discretisation of PDEs. From de Rham complexes, we derive other complexes with applications in elasticity, geometry and general relativity. The derivation, inspired by the Bernstein-Gelfand-Gelfand (BGG) construction, also provides a general machinery to establish results for tensor-valued problems (e.g., elasticity) from de Rham complexes (e.g., electromagnetism and fluid mechanics). We discuss some applications and progress in this direction, including mechanics models and the construction of bounded homotopy operators (Poincaré integrals) and finite elements.

The neutron transport equation is a linear Boltzmann-type PDE that models radiative transfer processes, and fission nuclear reactions. The computation of the largest eigenvalue of this Boltzmann operator is crucial in nuclear safety studies but it has classically been formulated only at a discretized level, so the predictive capabilities of such computations are fairly limited. In this talk, I will give an overview of the modeling for this equation, as well as recent analysis that leads to an infinite dimensional formulation of the eigenvalue problem. We leverage this point of view to build a numerical scheme that comes with a rigorous, a posteriori estimation of the error between the exact, infinite-dimensional solution, and the computed one.

I will start with two very brief surveys.  First is a class of problems, namely variational inequalities (VIs), which generalize PDE problems, and second is a class of solver algorithms, namely full approximation storage (FAS) nonlinear multigrid for PDEs.  Motivation for applying FAS to VIs is demonstrated in the standard mathematical model for glacier surface evolution, a very general VI problem relevant to climate modeling.  (Residuals for this nonlinear and non-local VI problem are computed by solving a Stokes model.)  Some existing nonlinear multilevel VI schemes, based on global (Newton) linearization would seem to be less suited to such general VI problems.  From this context I will sketch some work-in-progress toward the scalable solutions of nonlinear and nonlocal VIs by an FAS-type multilevel method.

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

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

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

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

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

Learning dynamical models from data plays a vital role in engineering design, optimization, and predictions. Building models describing the dynamics of complex processes (e.g., weather dynamics, reactive flows, brain/neural activity, etc.) using empirical knowledge or first principles is frequently onerous or infeasible. Therefore, system identification has evolved as a scientific discipline for this task since the 1960ies. Due to the obvious similarity of approximating unknown functions by artificial neural networks, system identification was an early adopter of machine learning methods. In the first part of the talk, we will review the development in this area until now.

For complex systems, identifying the full dynamics using system identification may still lead to high-dimensional models. For engineering tasks like optimization and control synthesis as well as in the context of digital twins, such learned models might still be computationally too challenging in the aforementioned multi-query scenarios. Therefore, it is desirable to identify compact approximate models from the available data. In the second part of this talk, we will therefore exploit that the dynamics of high-fidelity models often evolve in lowdimensional manifolds. We will discuss approaches learning representations of these lowdimensional manifolds using several ideas, including the lifting principle and autoencoders. In particular, we will focus on learning state-space representations that can be used in classical tools for computational engineering. Several numerical examples will illustrate the performance and limitations of the suggested approaches.

Tensors, also known as multiway arrays, have become ubiquitous as representations for operators or as convenient schemes for storing data. Yet, when it comes to compressing these objects or analyzing the data stored in them, the tendency is to ``flatten” or ``matricize” the data and employ traditional linear algebraic tools, ignoring higher dimensional correlations/structure that could have been exploited. Impediments to the development of equivalent tensor-based approaches stem from the fact that familiar concepts, such as rank and orthogonal decomposition, have no straightforward analogues and/or lead to intractable computational problems for tensors of order three and higher.

In this talk, we will review some of the common tensor decompositions and discuss their theoretical and practical limitations. We then discuss a family of tensor algebras based on a new definition of tensor-tensor products. Unlike other tensor approaches, the framework we derive based around this tensor-tensor product allows us to generalize in a very elegant way all classical algorithms from linear algebra. Furthermore, under our framework, tensors can be decomposed in a natural (e.g. ‘matrix-mimetic’) way with provable approximation properties and with provable benefits over traditional matrix approximation. In addition to several examples from recent literature illustrating the advantages of our tensor-tensor product framework in practice, we highlight interesting open questions and directions for future research.

TBA

Surface water waves of significant interest, such as tsunamis and solitary waves, are nonlinear and dispersive waves. Unluckily, the equations derived from first principles that describe the propagation of surface water waves, known as Euler's equations, are immensely hard to study. For this reason, several approximate systems have been proposed as mathematical alternatives. We show that among the numerous simplified systems of PDEs of water wave theory there is only one that is provably well-posed (in Hadamard’s sense) in bounded domains with slip-wall boundary conditions. We also show that the particular well-posed system obeys most of the physical laws that acceptable water wave equations must obey, and it is consistent with the Euler equations. For the numerical solution of our system we rely on a Galerkin/finite element method based on Nitsche's method for which we have proved its convergence. Validation with laboratory data is also presented.

Direct Multisearch (DMS) is a well-known multiobjective derivative-free optimization class of methods, with competitive computational implementations that are often successfully used for benchmark of new algorithms and in practical applications. As a directional direct search method, its structure is organized in a search step and a poll step, being the latter responsible for its convergence. A first implementation of DMS was released in 2010. Since then, the algorithmic class has continued to be analyzed from the theoretical point of view and new improvements have been proposed for the numerical implementation. Worst-case-complexity bounds have been derived, a search step based on polynomial models has been defined, and parallelization strategies have successfully improved the numerical performance of the code, which has also shown to be competitive for multiobjective derivative-based problems. In this talk we will survey the algorithmic structure of this class of optimization methods, the main theoretical properties associated to it and report numerical experiments that validate its numerical competitiveness.

Combining the classical theory of optimal transport with modern operator splitting techniques, I will present a new numerical method for nonlinear, nonlocal partial differential equations, arising in models of porous media,materials science, and biological swarming. Using the JKO scheme, along with the Benamou-Brenier dynamical characterization of the Wasserstein distance, we reduce computing the solution of these evolutionary PDEs to solving a sequence of fully discrete minimization problems, with strictly convex objective function and linear constraint. We compute the minimizer of these fully discrete problems by applying a recent, provably convergent primal dual splitting scheme for three operators. By leveraging the PDE’s underlying variational structure, ourmethod overcomes traditional stability issues arising from the strong nonlinearity and degeneracy, and it is also naturally positivity preserving and entropy decreasing. Furthermore, by transforming the traditional linear equality constraint, as has appeared in previous work, into a linear inequality constraint, our method converges in fewer iterations without sacrificing any accuracy. We prove that minimizers of the fully discrete problem converge to minimizers of the continuum JKO problem as the discretization is refined, and in the process, we recover convergence results for existing numerical methods for computing Wasserstein geodesics. Simulations of nonlinear PDEs and Wasserstein geodesics in one and two dimensions that illustrate the key properties of our numerical method will be shown.

Solving sparse linear systems is omnipresent in scientific computing. Direct approaches based on matrix factorization are very robust, and since they can be used as a black-box, it is easy for other software to use them. However, the memory requirement of direct approaches scales poorly with the problem size, and the algorithms underpinning sparse direct solvers software are poorly suited to parallel computation. Multilevel Domain decomposition (MDD) methods are among the most efficient iterative methods for solving sparse linear systems. One of the main technical difficulties in using efficient MDD methods (and most other efficient preconditioners) is that they require information from the underlying problem which prohibits them from being used as a black-box. This was the motivation to develop the widely used algebraic multigrid for example. I will present a series of recently developed robust and fully algebraic MDD methods, i.e., that can be constructed given only the coefficient matrix and guarantee a priori prescribed convergence rate. The series consists of preconditioners for sparse least-squares problems, sparse SPD matrices, general sparse matrices, and saddle-point systems. Numerical experiments illustrate the effectiveness, wide applicability, scalability of the proposed preconditioners. A comparison of each one against state-of-the-art preconditioners is also presented.

Physics-informed neural networks (PINNs) [2] are a solution method for solving boundary value problems based on differential equations (PDEs). The key idea of PINNs is to incorporate the residual of the PDE as well as boundary conditions into the loss function of the neural network. This provides a simple and mesh-free approach for solving problems relating to PDEs. However, a key limitation of PINNs is their lack of accuracy and efficiency when solving problems with larger domains and more complex, multi-scale solutions. 

In a more recent approach, Finite Basis Physics-Informed Neural Networks (FBPINNs) [1], the authors use ideas from domain decomposition to accelerate the learning process of PINNs and improve their accuracy in this setting. In this talk, we show how Schwarz-like additive, multiplicative, and hybrid iteration methods for training FBPINNs can be developed. Furthermore, we will present numerical experiments on the influence on convergence and accuracy of these different variants. 

This is joint work with Alexander Heinlein (Delft) and Benjamin Moseley (Oxford).

References 
1.    [1]  B. Moseley, A. Markham, and T. Nissen-Meyer. Finite basis physics- informed neural networks (FBPINNs): a scalable domain decomposition approach for solving differential equations. arXiv:2107.07871, 2021. 
2.    [2]  M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019.

For the first time, a method has become available for fast computation of near-best rational approximations on arbitrary sets in the real line or complex plane: the AAA algorithm (Nakatsukasa-Sète-T. 2018).  After a brief presentation of the algorithm this talk will focus on twenty demonstrations of the kinds of things we can do, all across applied mathematics, with a black-box rational approximation tool.
 

Hilbert complexes are chains of spaces linked by operators, with properties that are crucial to establishing the well-posedness of certain systems of partial differential equations. Designing stable numerical schemes for such systems, without resorting to nonphysical stabilisation processes, requires reproducing the complex properties at the discrete level. Finite-element complexes have been extensively developed since the late 2000's, in particular by Arnold, Falk, Winther and collaborators. These are however limited to certain types of meshes (mostly, tetrahedral and hexahedral meshes), which limits options for, e.g., local mesh refinement.

In this talk we will introduce the Discrete De Rham complex, a discrete version of one of the most popular complexes of differential operators (involving the gradient, curl and divergence), that can be applied on meshes consisting of generic polytopes. We will use a simple magnetostatic model to motivate the need for (continuous and discrete) complexes, then give a presentation of the lowest-order version of the complex and sketch its links with the CW cochain complex on the mesh. We will then briefly explain how this lowest-order version is naturally extended to an arbitrary-order version, and briefly present the associated properties (Poincaré inequalities, primal and adjoint consistency, commutation properties, etc.) that enable the analysis of schemes based on this complex.

The design and analysis of finite element methods for high Reynolds flow remains a challenging task, not least because of the difficulties associated with turbulence. In this talk we will first revisit some theoretical results on interior penalty methods using equal order interpolation for smooth solutions of the Navier-Stokes’ equations at high Reynolds number and show some recent computational results for turbulent flows.

Then we will focus on so called pressure robust methods, i.e. methods where the smoothness of the pressure does not affect the upper bound of error estimates for the velocity of the Stokes’ system. We will discuss how convection can be stabilized for such methods in the high Reynolds regime and, for the lowest order case, show an interesting connection to turbulence modelling.

The spread of a matrix is defined as the diameter of its spectrum. In this talk, we consider the problem of maximizing the spread of a symmetric non-negative matrix with bounded entries and discuss a number of recent results. This optimization problem is closely related to a pair of conjectures in spectral graph theory made by Gregory, Kirkland, and Hershkowitz in 2001, which were recently resolved by Breen, Riasanovsky, Tait, and Urschel. This talk will give a light overview of the approach used in this work, with a strong focus on ideas, many of which can be abstracted to more general matrix optimization problems.

We consider the Bayesian inverse problem of inferring the initial condition of a linear dynamical system from noisy output measurements taken after the initial time. In practical applications, the large dimension of the dynamical system state poses a computational obstacle to computing the exact posterior distribution. Balanced truncation is a system-theoretic method for model reduction which obtains an efficient reduced-dimension dynamical system by projecting the system operators onto state directions which simultaneously maximize energies defined by reachability and observability Gramians. We show that in our inference setting, the prior covariance and Fisher information matrices can be naturally interpreted as reachability and observability Gramians, respectively. We use these connections to propose a balancing approach to model reduction for the inference setting. The resulting reduced model then inherits stability properties and error bounds from system theory, and yields an optimal posterior covariance approximation. 

There is a wide range of problems in continuum mechanics that involve transformation fronts, which are non-stationary interfaces between two different phases in a phase-transforming or a chemically-transforming material. From the mathematical point of view, the considered problems are represented by systems of non-linear PDEs with discontinuities across non-stationary interfaces, kinetics of which depend on the solution of the PDEs. Such problems have a significant industrial relevance – an example of a transformation front is the localised stress-affected chemical reaction in Li-ion batteries with Si-based anodes. Since the kinetics of the transformation fronts depends on the continuum fields, the transformation front propagation can be decelerated and even blocked by the mechanical stresses. This talk will focus on three topics: (1) the stability of the transformation fronts in the vicinity of the equilibrium position for the chemo-mechanical problem, (2) a fictitious-domain finite-element method (CutFEM) for solving non-linear PDEs with transformation fronts and (3) an applied problem of Si lithiation.

Motivated by the problem of reconstructing the electron density of a molecule from pulsed X-ray diffraction images (about 10e+9 per reconstruction), we develop a framework for analyzing the convergence to invariant measures of random fixed point iterations built from mappings that, while expansive, nevertheless possess attractive fixed points.  Building on techniques that we have established for determining rates of convergence of numerical methods for inconsistent nonconvex
feasibility, we lift the relevant regularities to the setting of probability spaces to arrive at a convergence analysis for noncontractive Markov operators.  This approach has many other applications, for instance the analysis of distributed randomized algorithms.
We illustrate the approach on the problem of solving linear systems with finite precision arithmetic.

 

That the linear systems arising upon the discretization of elliptic PDEs can be solved efficiently is well-known, and iterative solvers that often attain linear complexity (multigrid, Krylov methods, etc) have proven very successful. Interestingly, it has recently been demonstrated that it is often possible to directly compute an approximate inverse to the coefficient matrix in linear (or close to linear) time. The talk will argue that such direct solvers have several compelling qualities, including improved stability and robustness, the ability to solve certain problems that have remained intractable to iterative methods, and dramatic improvements in speed in certain environments.

After a general introduction to the field, particular attention will be paid to a set of recently developed randomized algorithms that construct data sparse representations of large dense matrices that arise in scientific computations. These algorithms are entirely black box, and interact with the linear operator to be compressed only via the matrix-vector multiplication.

In space dimension 2 we present a finite element complex for the deformation operator acting on vectorfields and the linearized curvature operator acting on symmetric 2 by 2 matrices. We also present the tools that were used in the construction, namely the BGG diagram chase and the framework of finite element systems. For this general framework we can prove a de Rham theorem on cohomology groups in the flat case and a Bianchi identity in the case with curvature.

Tensor structured linear operators play an important role in matrix equations and low-rank modelling. Motivated by this we consider the problem of approximating a matrix by a sum of Kronecker products. It is known that an optimal approximation in Frobenius norm can be obtained from the singular value decomposition of a rearranged matrix, but when the goal is to approximate the matrix as a linear map, an operator norm would be a more appropriate error measure. We present an alternating optimization approach for the corresponding approximation problem in spectral norm that is based on semidefinite programming, and report on its practical performance for small examples.
This is joint work with Venkat Chandrasekaran and Mareike Dressler.

The recent, rapid advances in modern biology and data science have opened up a whole range of challenging mathematical problems. In this talk I will discuss a class of interacting particle models with anisotropic repulsive-attractive interaction forces. These models are motivated by the simulation of fingerprint databases, which are required in forensic science and biometric applications. In existing models, the forces are isotropic and particle models lead to non-local aggregation PDEs with radially symmetric potentials. The central novelty in the models I consider is an anisotropy induced by an underlying tensor field. This innovation does not only lead to the ability to describe real-world phenomena more accurately, but also renders their analysis significantly harder compared to their isotropic counterparts. I will discuss the role of anisotropic interaction in these models, present a stability analysis of line patterns, and show numerical results for the simulation of fingerprints. I will also outline how very similar models can be used in data classification, where it is desirable to assign labels to points in a point cloud, given that a certain number of points is already correctly labeled.

Solving large-scale Bayesian inverse problems governed by complex models suffers from the twin difficulties of the high dimensionality of the uncertain parameters and computationally expensive forward models. In this talk, we focus on 1. reducing the computational cost when solving these problems (via joint parameter and state dimension reduction) and 2. accounting for the error due to using a reduced order forward model (via Bayesian Approximation Error (BAE)).  To reduce the parameter dimension, we exploit the underlying problem structure (e.g., local sensitivity of the data to parameters, the smoothing properties of the forward model, the fact that the data contain limited information about the (infinite-dimensional) parameter field, and the covariance structure of the prior) and identify a likelihood-informed parameter subspace that shows where the change from prior to posterior is most significant. For the state dimension reduction, we employ a proper orthogonal decomposition (POD) combined with the discrete empirical interpolation method (DEIM) to approximate the nonlinear term in the forward model. We illustrate our approach with a model ice sheet inverse problem governed by the nonlinear Stokes equation for which the basal sliding coefficient field (a parameter that appears in a Robin boundary condition at the base of the geometry) is inferred from the surface ice flow velocity. The results show the potential to make the exploration of the full posterior distribution of the parameter or subsequent predictions more tractable.

This is joint work with Ki-Tae Kim (UC Merced), Benjamin Peherstorfer (NYU) and Tiangang Cui (Monash University).

Polynomial preconditioning for linear solvers is well-known but not frequently used in current scientific applications.  Furthermore, polynomial preconditioning has long been touted as well-suited for parallel computing; does this claim still hold in our new world of GPU-dominated machines?  We give details of the GMRES polynomial preconditioner and discuss its simple implementation, its properties such as eigenvalue remapping, and choices such as the starting vector and added roots.  We compare polynomial preconditioned GMRES to related methods such as FGMRES and full GMRES without restarting. We conclude with initial evaluations of the polynomial preconditioner for parallel and GPU computing, further discussing how polynomial preconditioning can become useful to real-word applications.

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

TBA

Let $A$ be an $n$ by $n$ matrix or a bounded linear operator on a complex Hilbert space $(H, \langle \cdot , \cdot \rangle , \| \cdot \|)$. A closed set $\Omega \subset \mathbb{C}$ is a $K$-spectral set for $A$ if the spectrum of $A$ is contained in $\Omega$ and if, for all rational functions $f$ bounded in $\Omega$, the following inequality holds:
\[\| f(A) \| \leq K \| f \|_{\Omega} ,\]
where $\| \cdot \|$ on the left denotes the norm in $H$ and $\| \cdot \|_{\Omega}$ on the right denotes the $\infty$-norm on $\Omega$. A simple way to obtain a $K$ value for a given set $\Omega$ is to use the Cauchy integral formula and replace the norm of the integral by the integral of the resolvent norm:
\[f(A) = \frac{1}{2 \pi i} \int_{\partial \Omega} ( \zeta I - A )^{-1}
f( \zeta )\,d \zeta \Rightarrow
\| f(A) \| \leq \frac{1}{2 \pi} \left( \int_{\partial \Omega}
\| ( \zeta I - A )^{-1} \|~| d \zeta | \right) \| f \|_{\Omega} .\]
Thus one can always take
\[K = \frac{1}{2 \pi} \int_{\partial \Omega} \| ( \zeta I - A )^{-1} \| | d \zeta | .\]
In M. Crouzeix and A. Greenbaum, Spectral sets: numerical range and beyond, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 1087-1101, different bounds on $K$ were derived.  I will show how these compare to that from the Cauchy integral formula for a variety of applications.  In case $A$ is a matrix and $\Omega$ is simply connected, we can numerically compute what we believe to be the optimal value for $K$ (and, at least, is a lower bound on $K$).  I will show how these values compare with the proven bounds as well.

(joint with  Michel Crouzeix and Natalie Wellen)
 

---

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

I will discuss new work on practically popular algorithms, including the kernel polynomial method (KPM) and moment matching method, for approximating the spectral density (eigenvalue distribution) of an n x n symmetric matrix A. We will see that natural variants of these algorithms achieve strong worst-case approximation guarantees: they can approximate any spectral density to epsilon accuracy in the Wasserstein-1 distance with roughly O(1/epsilon) matrix-vector multiplications with A. Moreover, we will show that the methods are robust to *in accuracy* in these matrix-vector multiplications, which allows them to be combined with any approximation multiplication algorithm. As an application, we develop a randomized sublinear time algorithm for approximating the spectral density of a normalized graph adjacency or Laplacian matrices. The talk will cover the main tools used in our work, which include random importance sampling methods and stability results for computing orthogonal polynomials via three-term recurrence relations.

---

A link for this talk will be sent to our mailing list a day or two in advance.  If you are not on the list and wish to be sent a link, please contact @email.

Climate- and weather prediction centres such as the Met Office rely on efficient numerical methods for simulating large scale atmospheric flow. One computational bottleneck in many models is the repeated solution of a large sparse system of linear equations. Preconditioning this system is particularly challenging for state-of-the-art discretisations, such as (mimetic) finite elements or Discontinuous Galerkin (DG) methods. In this talk I will present recent work on developing efficient multigrid preconditioners for practically relevant modelling codes. As reported in a REF2021 Industrial Impact Case Study, multigrid has already led to runtime savings of around 10%-15% for operational global forecasts with the Unified Model. Multigrid also shows superior performance in the Met Office next-generation LFRic model, which is based on a non-trivial finite element discretisation.