Past Computational Mathematics and Applications Seminar

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.