Computational Mathematics and Applications Seminar

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.

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.

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.

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)

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.

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.

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.