Numerical Solution of Double Saddle-Point Systems

Double saddle-point systems are drawing increasing attention in the past few years, due to the importance of multiphysics applications and recent advances in the solution of indefinite linear systems. In this talk I will describe some of the numerical properties of these linear systems. Eigenvalue bounds will be presented, and we will show what role Schur complements and nested Schur complements play in preconditioning. Some numerical experiments on the Stokes-Darcy equations will be described.

When is accurate and efficient expression evaluation and linear algebra possible?

When, and at what cost, can we solve linear algebra problems, or even just evaluate multivariate polynomials, with high relative accuracy? Suppose we only know that we can perform the 4 basic operations (+, -, *, /) with high relative accuracy, i.e. rnd(a op b) = (a op b)*(1+delta) with |delta| <= epsilon << 1. Which 2 of the following 3 problems can we solve with high relative accuracy?

  (1) Evaluate the Motzkin polynomial p(x,y,z) = z^6 + x^2*y^2*(x^2 + y^2 - 3*z^2)

  (2) Compute the eigenvalues of the Cauchy matrix C(i,j) = 1/(x(i) + x(j)), given 0 < x(1) < x(2) < ... , an example of which is the Hilbert matrix.

  (3) Evaluate x+y+z.

We will answer this question, and many related ones, and talk about the open problem of creating a decision procedure that would take a polynomial, and decide whether an accurate algorithm for it exists. Properties of the variety of the polynomial will play a central role. We will also discuss some properties of floating point arithmetic that will let us expand the list of basic operations above, and so expand the set of problems we can solve accurately.

Markov chains, graph Laplacians, and column subset selection

The column subset selection problem (CSSP) is an essential target of research in low-rank approximation. I will describe how our particular motivation in reduced order models of Markov chains has led to new theory and algorithms for this long-standing problem. Our problem setting, relevant to ubiquitous dynamical systems in chemistry, physics, and computer science, is nicely mappable to solving the CSSP for the (suitably regularized) inverse of the respective graph Laplacian. Using contour integral and probabilistic arguments, I will first show how the dynamical error of reversible Markov chain compression can be formulated and simply bounded in terms of Nyström approximation. Then I will describe how nuclear maximization is an especially attractive method for the CSSP in this setting, showing new matrix-free algorithms based on randomized diagonal estimation via approximate Cholesky factorization. I will close by demonstrating our algorithms' empirical performance and theoretical guarantees, detailing a novel combination of ideas from submodularity, probability theory, and determinantal point processes. Time permitting, I will outline our successes in extending our methods and error analysis to general CX, Nyström, and CUR approximation.

Rank-revealing factorizations: applications, algorithms, and theory

Rank-revealing matrix factorizations play a key role in applications ranging from structured low-rank approximation to computing localized basis functions in computational quantum chemistry. We will discuss some motivating applications to illustrate the breadth of desiderata for such factorizations. This will be followed by a brief overview of algorithms (where we will restrict our discussion to pivoted QR factorizations for simplicity) and highlights of historical and recent theoretical developments. We will close by highlighting some open problems in the area.

A randomized linear algebra perspective on geometric optimization methods in scientific machine learning

Subsampled natural gradient (SNG) and Gauss-Newton (SGN) methods have demonstrated impressive performance for parametric optimization problems in scientific machine learning, including neural network wavefunctions and physics-informed neural networks. However, their development has been hindered by a lack of theoretical understanding. To make progress on this problem, we study these algorithms on a simplified parametric optimization problem involving a linear model and a quadratic loss function. In this setting we show that both SNG and SGN can be interpreted as sketch-and-project methods and can thus enjoy rigorous convergence guarantees for arbitrary batch sizes. This interpretation also explains a recently proposed accelerated variant of SNG and clarifies how this accelerated variant can be extended to the setting of SGN. Finally, these explorations inspire a number of new questions about sketch-and-project algorithms, some of which have been resolved and many of which remain open.

 Use of Johnson–Lindenstrauss (JL) sketching operators in faster hierarchical matrix construction

We present an adaptive, partially matrix-free, hierarchical matrix construction framework using a broader class of Johnson–Lindenstrauss (JL) sketching operators. On the theoretical side, we extend the earlier concentration bounds to all JL sketching operators and examine this bound for specific classes of such operators including the original Gaussian sketching operators, subsampled randomized Hadamard transform (SRHT) and the sparse Johnson-Lindenstrauss transform (SJLT). We discuss the implementation details of applying SJLT and SRHT efficiently. Then we demonstrate experimentally that using SJLT or SRHT instead of Gaussian sketching operators leads to up to 2.5 times speedups of the serial HSS construction implementation in the STRUMPACK C++ library. Additionally, we discuss the implementation of a parallel distributed HSS construction that leverages Gaussian or SJLT sketching operators. We observe a performance improvement of up to 35 times when using SJLT sketching operators over Gaussian sketching operators. The generalized algorithm allows users to select their own JL sketching operators with theoretical lower bounds on the size of the operators which may lead to faster runtime with similar HSS construction accuracy.

Revisiting MR^3 for the Bidiagonal SVD

The Multiple Relatively Robust Representations (MRRR or MR^3) algorithm is optimal for the symmetric tridiagonal eigenvalue problem -- that is, it can compute k eigenpairs (with numerically orthogonal eigenvectors) in only O(nk) operations.  Accordingly, MR^3 yields the fastest algorithm for the symmetric eigenvalue problem when paired with standard tridiagonal reduction. Given that the bidiagonal SVD and symmetric tridiagonal eigenvalue problems are two sides of the same coin  -- and indeed the fastest algorithms for computing an arbitrary SVD begin with bidiagonal reduction -- it seems natural that MR^3 should imply an optimal algorithm for this problem as well. Nevertheless, it has yet to be successfully deployed in the SVD setting. This is in spite of the fact that one of the key insights of MR^3 is that representing a tridiagonal as a product of bidiagonals ensures high accuracy. In this talk we revisit this problem, tracing the origins of MR^3, the pitfalls of using it to compute the bidiagonal SVD, and the many efforts to make it work anyways.