Computational Mathematics and Applications Seminar

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).

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.

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.

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.

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.

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.

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.

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

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.