Past Computational Mathematics and Applications Seminar

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.

Function approximation, as a goal in itself or as an ingredient in scientific computing, typically relies on having a basis. However, in many cases of interest an obvious basis is not known or is not easily found. Even if it is, alternative representations may exist with much fewer degrees of freedom, perhaps by mimicking certain features of the solution into the “basis functions" such as known singularities or phases of oscillation. Unfortunately, such expert knowledge typically doesn’t match well with the mathematical properties of a basis: it leads instead to representations which are either incomplete or overcomplete. In turn, this makes a problem potentially unsolvable or ill-conditioned. We intend to show that overcomplete representations, in spite of inherent ill-conditioning, often work wonderfully well in numerical practice. We explore a theoretical foundation for this phenomenon, use it to devise ground rules for practitioners, and illustrate how the theory and its ramifications manifest themselves in a number of 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.

We use energy estimates to derive new bounds on the eigenvalues of a generic form of double saddle-point matrices, with and without regularization terms. Results related to inertia and algebraic multiplicity of eigenvalues are also presented. The analysis includes eigenvalue bounds for preconditioned matrices based on block-diagonal Schur complement-based preconditioners, and it is shown that in this case the eigenvalues are clustered within a few intervals bounded away from zero. The analytical observations are linked to a few multiphysics problems of interest. This is joint work with Susanne Bradley.

---

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.

Variational phase-field models of fracture have been at the center of a multidisciplinary effort involving a large community of mathematicians, mechanicians, engineers, and computational scientists over the last 25 years or so.

I will start with a modern interpretation of Griffith's classical criterion as a variational principle for a free discontinuity energy and will recall some of the milestones in its analysis. Then, I will introduce the phase-field approximation per se and describe its numerical implementation. I illustrate how phase-field models have led to major breakthroughs in the predictive simulation of fracture in complex situations.

I then will turn my attention to current issues, with a specific emphasis on crack nucleation in nominally brittle materials. I will recall the fundamental incompatibility between Griffith’s theory and nucleation criteria based on a stress yield surface: the strength vs. toughness paradox. I will then present several attempts at addressing this issue within the realm of phase-fracture and discuss their respective strengths and weaknesses. 

--

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.

Uncertainty Quantification through Markov Chain Monte Carlo (MCMC) can be prohibitively expensive for target probability densities with expensive likelihood functions, for instance when the evaluation it involves solving a Partial Differential Equation (PDE), as is the case in a wide range of engineering applications. Multilevel Delayed Acceptance (MLDA) with an Adaptive Error Model (AEM) is a novel approach, which alleviates this problem by exploiting a hierarchy of models, with increasing complexity and cost, and correcting the inexpensive models on-the-fly. The method has been integrated within the open-source probabilistic programming package PyMC3 and is available in the latest development version.

In this talk I will talk about the problems with the Multilevel Markov Chain Monte Carlo (Dodwell et al. 2015). In so we will prove detailed balance for Adaptive Multilevel Delayed Acceptance, as well as showing that multilevel variance reduction can be achieved without bias, not possible in the original MLMCMC framework.

I will talk about our implementation in the latest version of pymc3, and demonstrate how for classical inverse problem benchmarks the AMLDA sampler offers huge computational savings (> factor of 100 fold speed up).

Finally I will talk heuristically about new / future research, in which we seek to develop parallel strategies for this inherently sequential sampler, as well as point to interesting applied application areas in which the method is proving particular effective.

--

This talk will be in person.

Computing spectral properties of operators is fundamental in the sciences, with applications in quantum mechanics, signal processing, fluid mechanics, dynamical systems, etc. However, the infinite-dimensional problem is infamously difficult (common difficulties include spectral pollution and dealing with continuous spectra). This talk introduces classes of practical resolvent-based algorithms that rigorously compute a zoo of spectral properties of operators on Hilbert spaces. We also discuss how these methods form part of a broader programme on the foundations of computation. The focus will be computing spectra with error control and spectral measures, for general discrete and differential operators. Analogous to eigenvalues and eigenvectors, these objects “diagonalise” operators in infinite dimensions through the spectral theorem. The first is computed by an algorithm that approximates resolvent norms. The second is computed by building convolutions of appropriate rational functions with the measure via the resolvent operator (solving shifted linear systems). The final part of the talk provides purely data-driven algorithms that compute the spectral properties of Koopman operators, with convergence guarantees, from snapshot data. Koopman operators “linearise” nonlinear dynamical systems, the price being a reduction to an infinite-dimensional spectral problem (c.f. “Koopmania”, describing their surge in popularity). The talk will end with applications of these new methods in several thousand state-space dimensions.

The roots of a function represented by its Chebyshev expansion are known to be the eigenvalues of the so-called colleague matrix, which is a Hessenberg matrix that is the sum of a symmetric tridiagonal matrix and a rank 1 perturbation. The rootfinding problem is thus reformulated as an eigenproblem, making the computation of the eigenvalues of such matrices a subject of significant practical interest. To obtain the roots with the maximum possible accuracy, the eigensolver used must posess a somewhat subtle form of stability.

In this talk, I will discuss a recently constructed algorithm for the diagonalization of colleague matrices, satisfying the relevant stability requirements.  The scheme has CPU time requirements proportional to n^2, with n the dimensionality of the problem; the storage requirements are proportional to n. Furthermore, the actual CPU times (and storage requirements) of the procedure are quite acceptable, making it an approach of choice even for small-scale problems. I will illustrate the performance of the algorithm with several numerical examples.

--

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.

The past few years have been an exciting time for my work related to rational approximation.  This talk will present four developments:

1. AAA approximation (2016, with Nakatsukasa & Sète)
2. Root-exponential convergence and tapered exponential clustering (2020, with Nakatsukasa & Weideman)
3. Lightning (2017-2020, with Gopal & Brubeck)
4. Log-lightning (2020-21, with Nakatsukasa & Baddoo)

Two other topics will not be discussed:

X. AAA-Lawson approximation (2018, with Nakatsukasa)
Y. AAA-LS approximation (2021, with Costa)

Randomized NLA methods have recently gained popularity because of their easy implementation, computational efficiency, and numerical robustness. We propose a randomized version of a well-established FEAST eigenvalue algorithm that enables computing the eigenvalues of the Hermitian matrix pencil $(\textbf{A},\textbf{B})$ located in the given real interval $\mathcal{I} \subset [\lambda_{min}, \lambda_{max}]$. In this talk, we will present deterministic as well as probabilistic error analysis of the accuracy of approximate eigenpair and subspaces obtained using the randomized FEAST algorithm. First, we derive bounds for the canonical angles between the exact and the approximate eigenspaces corresponding to the eigenvalues contained in the interval $\mathcal{I}$. Then, we present bounds for the accuracy of the eigenvalues and the corresponding eigenvectors. This part of the analysis is independent of the particular distribution of an initial subspace, therefore we denote it as deterministic. In the case of the starting guess being a Gaussian random matrix, we provide more informative, probabilistic error bounds. Finally, we will illustrate numerically the effectiveness of all the proposed error bounds.

---

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 recent advances in sampling methods for positive semidefinite (PSD) matrix approximation. In particular, I will show how new techniques based on recursive leverage score sampling yield a surprising algorithmic result: we give a method for computing a near optimal k-rank approximation to any n x n PSD matrix in O(n * k^2) time. When k is not too large, our algorithm runs in sublinear time -- i.e. it does not need to read all entries of the matrix. This result illustrates the ability of randomized methods to exploit the structure of PSD matrices and go well beyond what is possible with traditional algorithmic techniques. I will discuss a number of current research directions and open questions, focused on applications of randomized methods to sublinear time algorithms for structured matrix problems.

--

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.

One of the great challenges of neural networks is to understand how they work.  For example: does a neuron encode a meaningful signal on its own?  Or is a neuron simply an undistinguished and arbitrary component of a feature vector space?  The tension between the neuron doctrine and the population coding hypothesis is one of the classical debates in neuroscience. It is a difficult debate to settle without an ability to monitor every individual neuron in the brain.

Within artificial neural networks we can examine every neuron. Beginning with the simple proposal that an individual neuron might represent one internal concept, we conduct studies relating deep network neurons to human-understandable concepts in a concrete, quantitative way: Which neurons? Which concepts? Are neurons more meaningful than an arbitrary feature basis? Do neurons play a causal role? We examine both simplified settings and state-of-the-art networks in which neurons learn how to represent meaningful objects within the data without explicit supervision.

Following this inquiry in computer vision leads us to insights about the computational structure of practical deep networks that enable several new applications, including semantic manipulation of objects in an image; understanding of the sparse logic of a classifier; and quick, selective editing of generalizable rules within a fully trained generative network.  It also presents an unanswered mathematical question: why is such disentanglement so pervasive?

In the talk, we challenge the notion that the internal calculations of a neural network must be hopelessly opaque. Instead, we propose to tear back the curtain and chart a path through the detailed structure of a deep network by which we can begin to understand its logic.

--

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 modified Newton methods for solving nonlinear systems of PDEs. These methods introduce additional variables before deriving the Newton linearization. These variables can then often be eliminated analytically before solving the Newton system, such that existing solvers can be adapted easily and the computational cost does not increase compared to a standard Newton method. The resulting algorithms yield favorable convergence properties. After illustrating the ideas on a simple example, I will show its application for the solution of incompressible Stokes flow problems with viscoplastic constitutive relation, where the additionally introduced variable is the stress tensor. These models are commonly used in earth science models. This is joint work with Johann Rudi (Argonne) and Melody Shih (NYU).

--

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.

An insightful exercise might be to ask what is the most important idea in linear algebra. Our first answer would not be eigenvalues or linearity, it would be “matrix factorizations.” We will discuss a blueprint to generate 53 inter-related matrix factorizations (times 2) most of which appear to be new. The underlying mathematics may be traced back to Cartan (1927), Harish-Chandra (1956), and Flensted-Jensen (1978) . We will discuss the interesting history. One anecdote is that Eugene Wigner (1968) discovered factorizations such as the SVD in passing in a way that was buried and only eight authors have referenced that work. Ironically Wigner referenced Sigurður Helgason (1962) but Wigner did not recognize his results in Helgason's book. This work also extends upon and completes open problems posed by Mackey² & Tisseur (2003/2005).

Classical results of Random Matrix Theory concern exact formulas from the Hermite, Laguerre, Jacobi, and Circular distributions. Following an insight from Freeman Dyson (1970), Zirnbauer (1996) and Duenez (2004/5) linked some of these classical ensembles to Cartan's theory of Symmetric Spaces. One troubling fact is that symmetric spaces alone do not cover all of the Jacobi ensembles. We present a completed theory based on the generalized Cartan distribution. Furthermore, we show how the matrix factorization obtained by the generalized Cartan decomposition G=K₁AK₂ plays a crucial role in sampling algorithms and the derivation of the joint probability density of A.

Joint work with Sungwoo Jeong

--

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.

In this talk, we present some effective methods for distributing points for approximating analytic functions with prescribed decay on a strip region including the real axis. Such functions appear when we use numerical methods with variable transformations. Typical examples of such methods are provided by single-exponential (SE) or double-exponential (DE) sinc formulas, in which variable transformations yield single- or double-exponential decay of functions on the real axis. It has been known that the formulas are nearly optimal on a Hardy space with a single- or double-exponential weight on the strip region, which is regarded as a space of transformed functions by the variable transformations.

Recently, we have proposed new approximation formulas that outperform the sinc formulas. For constructing them, we use an expression of the error norm (a.k.a. worst-case error) of an n-point interpolation operator in the weighted Hardy space. The expression is closely related to potential theory, and optimal points for interpolation correspond to an equilibrium measure of an energy functional with an external field. Since a discrete version of the energy becomes convex in the points under a mild condition about the weight, we can obtain good points with a standard optimization technique. Furthermore, with the aid of the formulation with the energy, we can find approximate distributions of the points theoretically.

[References]
- K. Tanaka, T. Okayama, M. Sugihara: Potential theoretic approach to design of accurate formulas for function approximation in symmetric weighted Hardy spaces, IMA Journal of Numerical Analysis Vol. 37 (2017), pp. 861-904.

- K. Tanaka, M. Sugihara: Design of accurate formulas for approximating functions in weighted Hardy spaces by discrete energy minimization, IMA Journal of Numerical Analysis Vol. 39 (2019), pp. 1957-1984.

- S. Hayakawa, K. Tanaka: Convergence analysis of approximation formulas for analytic functions via duality for potential energy minimization, arXiv:1906.03133.

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.

Computational science is facing several major challenges with rapidly changing highly complex heterogeneous computer architectures. To meet these challenges and yield fast and efficient performance, solvers need to be easily portable. Algebraic multigrid (AMG) methods have great potential to achieve good performance, since they have shown excellent numerical scalability for a variety of problems. However, their implementation on emerging computer architectures, which favor structure, presents new challenges. To face these difficulties, we have considered modularization of AMG, that is breaking AMG components into smaller kernels to improve portability as well as the development of new algorithms to replace components that are not suitable for GPUs. Another way to achieve performance on accelerators is to increase structure in algorithms. This talk will discuss new algorithmic developments, including a new class of interpolation operators that consists of simple matrix operations for unstructured AMG and efforts to develop a semi-structured AMG method.

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.

The bubble transform was a concept introduced by Richard Falk and me in a paper published in The Foundations of Computational Mathematics in 2016. From a simplicial mesh of a bounded domain in $R^n$ we constructed a map which decomposes scalar valued functions into a sum of local bubbles supported on appropriate macroelements.The construction is done without reference to any finite element space, but has the property that the standard continuous piecewise polynomial spaces are invariant. Furthermore, the transform is bounded in $L^2$ and $H^1$, and as a consequence we obtained a new tool for the understanding of finite element spaces of arbitrary polynomial order. The purpose of this talk is to review the previous results, and to discuss how to generalize the construction to differential forms such that the corresponding properties hold. In particular, the generalized transform will be defined such that it commutes with the exterior derivative.

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.

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.

We develop a trust-region method for minimizing the sum of a smooth term f and a nonsmooth term h, both of which can be nonconvex. Each iteration of our method minimizes a possibly nonconvex model of f+h in a trust region. The model coincides with f+h in value and subdifferential at the center. We establish global convergence to a first-order stationary point when f satisfies a smoothness condition that holds, in particular, when it has Lipschitz-continuous gradient, and h is proper and lower semi-continuous. The model of h is required to be proper, lower-semi-continuous and prox-bounded. Under these weak assumptions, we establish a worst-case O(1/ε^2) iteration complexity bound that matches the best known complexity bound of standard trust-region methods for smooth optimization. We detail a special instance in which we use a limited-memory quasi-Newton model of f and compute a step with the proximal gradient method, resulting in a practical proximal quasi-Newton method. We describe our Julia implementations and report numerical results on inverse problems from sparse optimization and signal processing. Our trust-region algorithm exhibits promising performance and compares favorably with linesearch proximal quasi-Newton methods based on convex models.

This is joint work with Aleksandr Aravkin and Robert Baraldi.

-

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.

Dissipative Hamiltonian systems are an important class of dynamical systems that arise in all areas of science and engineering. They are a special case of port-Hamiltonian control systems. When the system is linearized arround a stationary solution one gets a linear dissipative Hamiltonian typically differential-algebraic system. Despite the fact that the system looks unstructured at first sight, it has remarkable properties.  Stability and passivity are automatic, spectral structures for purely imaginary eigenvalues, eigenvalues at infinity, and even singular blocks in the Kronecker canonical form are very restricted and furthermore the structure leads to fast and efficient iterative solution methods for asociated linear systems. When port-Hamiltonian systems are subject to (structured) perturbations, then it is important to determine the minimal allowed perturbations so that these properties are not preserved. The computation of these structured distances to instability, non-passivity, or non-regularity, is typically a very hard non-convex optimization problem. However, in the context of dissipative Hamiltonian systems, the computation becomes much easier and can even be implemented efficiently for large scale problems in combination with model reduction techniques. We will discuss these distances and the computational methods and illustrate the results via an industrial problem.

 

--

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.

--

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.

This talk concerns applications of differential geometry in numerical optimization. They arise when the optimization problem can be formulated as finding an optimum of a real-valued cost function defined on a smooth nonlinear search space. Oftentimes, the search space is a "matrix manifold", in the sense that its points admit natural representations in the form of matrices. In most cases, the matrix manifold structure is due either to the presence of certain nonlinear constraints (such as orthogonality or rank constraints), or to invariance properties in the cost function that need to be factored out in order to obtain a nondegenerate optimization problem. Manifolds that come up in practical applications include the rotation group SO(3) (generation of rigid body motions from sample points), the set of fixed-rank matrices (low-rank models, e.g., in collaborative filtering), the set of 3x3 symmetric positive-definite matrices (interpolation of diffusion tensors), and the shape manifold (morphing).

In the recent years, the practical importance of optimization problems on manifolds has stimulated the development of geometric optimization algorithms that exploit the differential structure of the manifold search space. In this talk, we give an overview of geometric optimization algorithms and their applications, with an emphasis on the underlying geometric concepts and on the numerical efficiency of the algorithm implementations.

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.

Data scientists are often faced with the challenge of understanding a high dimensional data set organized as a table. These tables may have columns of different (sometimes, non-numeric) types, and often have many missing entries. In this talk, we discuss how to use low rank models to analyze these big messy data sets. Low rank models perform well --- indeed, suspiciously well — across a wide range of data science applications, including applications in social science, medicine, and machine learning. In this talk, we introduce the mathematics of low rank models, demonstrate a few surprising applications of low rank models in data science, and present a simple mathematical explanation for their effectiveness.

--

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.

Conventional computational methods often create a dilemma for fluid-structure interaction problems. Typically, solids are simulated using a Lagrangian approach with grid that moves with the material, whereas fluids are simulated using an Eulerian approach with a fixed spatial grid, requiring some type of interfacial coupling between the two different perspectives. Here, a fully Eulerian method for simulating structures immersed in a fluid will be presented. By introducing a reference map variable to model finite-deformation constitutive relations in the structures on the same grid as the fluid, the interfacial coupling problem is highly simplified. The method is particularly well suited for simulating soft, highly-deformable materials and many-body contact problems, and several examples will be presented.

This is joint work with Ken Kamrin (MIT).

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.

Computer-based simulation of partial differential equations (PDEs) involves approximating the unknowns and relies on suitable description of geometrical entities such as the computational domain and its properties. The Finite Element Method (FEM) is by large the most popular technique for the computer-based simulation of PDEs and hinges on the assumption that discretized domain and unknown fields are both represented by piecewise polynomials, on tetrahedral or hexahedral partitions. In reality, the simulation of PDEs is a brick within a workflow where, at the beginning, the geometrical entities are created, described and manipulated with a geometry processor, often through Computer-Aided Design systems (CAD), and then used for the simulation of the mechanical behaviour of the designed object. This workflow is often repeated many times as part of a shape optimisation loop. Within this loop, the use of FEM on CAD geometries (which are mainly represented through their boundaries) calls then for (re-) meshing and re-interpolation techniques that often require human intervention and result in inaccurate solutions and lack of robustness of the whole process. In my talk, I will present the mathematical counterpart of this problem, I will discuss the mismatch in the mathematical representations of geometries and PDEs unknowns and introduce a promising framework where geometric objects and PDEs unknowns are represented in a compatible way. Within this framework, the challenges to be addressed in order to construct robust PDE solvers are many and I will discuss some of them. Mathematical results will besupported by numerical validation.

Composite materials make up over 50% of recent aircraft constructions. They are manufactured from very thin fibrous layers  (~10^-4 m) and even  thinner resin interfaces (~10^-5 m). To achieve the required strength, a particular layup sequence of orientations of the anisotropic fibrous layers is used. During manufacturing, small localised defects in the form of misaligned fibrous layers can occur in composite materials, adding an additional level of complexity. After FE discretisation the model exhibits multiple scales and large spatial variations in model parameters. Thus the resultant linear system of equations can be very ill-conditioned and extremely large. The limitations of commercially available modelling tools for solving these problems has led us to the implementation of a robust and scalable preconditioner called GenEO for parallel Krylov solvers. I will discuss using the GenEO coarse space as an effective multiscale model for the fine-scale displacement and stress fields. For the coarse space construction, GenEO computes generalised eigenvectors of the local stiffness matrices on the overlapping subdomains and builds an approximate coarse space by combining the smallest energy eigenvectors on each subdomain via a partition of unity.

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.