Schedule for: 23w5004 - New Directions in Applied Linear Algebra

A buffet dinner is served daily between 5:30pm and 7:30pm in Vistas Dining Room, top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

A brief introduction to BIRS with important logistical information, technology instruction, and opportunity for participants to ask questions.

Linear programming is a central problem in computer science and applied mathematics with numerous applications across a wide range of domains, including machine learning and data science. Interior point methods (IPMs) are a common approach to solving linear programs with strong theoretical guarantees and solid empirical performance. The time complexity of these methods is dominated by the cost of solving a linear system of equations at each iteration. In common applications of linear programming, particularly in data science and scientific computing, the size of this linear system can become prohibitively large, requiring the use of iterative solvers which provide an approximate solution to the linear system. Approximately solving the linear system at each iteration of an IPM invalidates common analyses of IPMs and the theoretical guarantees they provide. In this talk we will discuss how randomized linear algebra can be used to design and analyze theoretically and practically efficient IPMs when using approximate linear solvers.

Double saddle-point systems are drawing increasing attention in the past few years, due to the importance of associated applications and the challenge in developing efficient numerical solvers. In this talk we describe some of their numerical properties and make some distinctions between this family and the more well-established class of classical 2-by-2 block-structured saddle-point systems. We derive eigenvalue bounds, expressed in terms of extremal eigenvalues and singular values of block sub-matrices. The analysis includes bounds for preconditioned matrices based on block diagonal preconditioners using Schur complements, and it is shown that in this case the eigenvalues are clustered within a few intervals bounded away from zero. Analysis for approximations of Schur complements is included.

I will present practical methods to apply block $\alpha$-circulant preconditioners approximately to an all-at-once system from ocean data assimilation. We apply the preconditioner via a nested iterative method within an outer Chebyshev semi-iteration. We present theory on the rate of convergence and use this to design an adapted approach which ensures each sub-problem is solved to the same tolerance. A panel of numerical experiments reveal that using an appropriate choice of $\alpha$ our new approaches are competitive in terms of outer iterations and matrix-vector products compared to current choices of preconditioner.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Meet in the PDC front desk for a guided tour of The Banff Centre campus.

Meet in foyer of TCPL to participate in the BIRS group photo. The photograph will be taken outdoors, so dress appropriately for the weather. Please don't be late, or you might not be in the official group photo!

Low rank structure is pervasive in real-world datasets. This talk shows how to accelerate the solution of fundamental computational problems, including eigenvalue decomposition, linear system solves, composite convex optimization, and stochastic optimization (including deep learning), by exploiting this low rank structure. We present a simple method based on randomized numerical linear algebra for efficiently computing approximate top eigendecompositions, which can be used to replace large matrices (such as Hessians and constraint matrices) with low rank surrogates that are faster to apply and invert. The resulting solvers for linear systems (NystromPCG), composite convex optimization (NysADMM), and deep learning (SketchySGD) demonstrate strong theoretical and numerical support, outperforming state-of-the-art methods in terms of speed and robustness to hyperparameters.

Anderson acceleration with window size $m$ is a nonlinear convergence acceleration mechanism for fixed-point iterative methods that is widely used in scientific computing, optimization and machine learning. For many applications $AA(m)$ dramatically improves the convergence speed, both when iteration counts are small, and asymptotically for large iteration counts. Nevertheless, there are still no known results yet that can bound or quantify the improvement in asymptotic convergence speed provided by windowed $AA(m)$. In this talk I will give an overview of what is known about the asymptotic convergence speed of windowed $AA(m)$ and highlight some recent results and open problems.

Model order reduction techniques are effective for solving parametrized models involving PDEs, including models of incompressible flow, where the constraint is the incompressibility constraint, and in optimal control, where the constraints themselves are PDEs. However, reduced models may fail to be inf-sup stable. We present a new approach for generating reduced bases in this scenario, using a so-called stacked reduced basis, which avoids some of the difficulties associated with inf-sup stability. We show that this approach is effective although in some circumstances it also requires stabilization, which can be done using either classic methods of penalization or through Petrov-Galerkin methods. Computational tests are presented for models based on PDE-constrained optimization and incompressible flow.

This is joint work with Kayla D. Davie (University of Maryland, College Park).

A buffet dinner is served daily between 5:30pm and 7:30pm in Vistas Dining Room, top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

In this talk I will outline recent trends in uncertainty quantification (UQ) with PDE-based models and connect the trends to linear algebra. The topics in UQ are: (1) advanced random field models such as hierarchical random fields, (2) dynamical systems with random coefficients and their sensitivity analysis, (3) reliability analysis and rare events. The associated linear algebra topics are: matrix factorization, matrix eigenvalue problems, numerical continuation, and dimension reduction.

We study the iterative solution of linear systems of equations arising from stochastic Galerkin finite element discretizations of symmetric saddle point problems with stochastic data. As preconditioners we consider block matrices with building blocks that have Kronecker product structure. We derive bounds for the eigenvalues of the preconditioned system matrix and show how the existence requirements of a Bramble-Pasciak-type conjugate gradient method can be met. Based on numerical experiments for Stokes flow with random viscosity, we compare a MINRES solver with block diagonal preconditioner to a conjugate gradient method with block triangular preconditioner.

Parameters in mathematical models for physical processes are often impossible to determine fully or accurately, and are hence subject to uncertainty. In this talk, we employ the multilevel Monte Carlo (MLMC) method to compute expected values of quantities of interest related to partial differential equations with random coefficients. We make use of the circulant embedding method for sampling from the coefficient, and devise a smoothing technique to further improve the computational complexity of the MLMC estimator.

Random feature methods have been successful in various machine learning tasks, are easy to compute, and come with theoretical accuracy bounds. They serve as an alternative approach to standard neural networks since they can represent similar function spaces without a costly training phase. However, for accuracy, random feature methods require more measurements than trainable parameters, limiting their use for data-scarce applications or problems in scientific machine learning. In this talk, we will introduce the sparse random feature expansion to obtain parsimonious random feature models. Specifically, we leverage ideas from compressive sensing to generate random feature expansions with theoretical guarantees even in the data-scarce setting. We also present a sparse random mode decomposition to extract intrinsic modes from challenging time-series data and propose a new method to predict the infectious using sparse random feature models and time-delay equations. Applications of our methods on identifying important variables in high-dimensional settings as well as on decomposing music pieces, visualizing black-hole mergers, and epidemiologic forecasting will be addressed.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

A buffet dinner is served daily between 5:30pm and 7:30pm in Vistas Dining Room, top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

A variety of problems in modern applications of optimization require a selection of a 'sparse' solution, a vector with preferably few nonzero entries. Such problems may originate from very different applications in computational statistics, signal or image processing or compressed sensing, finance, machine learning and discrete optimal transport, to mention just a few. Sparse approximation problems are often solved with dedicated and highly specialised first-order methods of optimization. In this talk I will argue that these problems may be very efficiently solved by interior point methods. The key to their success is a design of suitable preconditioners.

Optimal Transport (OT) is a branch of optimization theory that has prompted a number of successful applications in different fields in the last few years. However, the wide-spread application of OT ideas in science and engineering is still fragmentary, due to the lack of accurate and efficient numerical solvers. We will explore some recent advances and open issues in the construction of efficient preconditioners involved in the solution of PDE-based formulations of OT.

We present a scalable method for large-scale elliptic PDE and bound constrained optimization. Such problems are a means to learn unknown aspects of PDE-based models and are particularly challenging as their optimality systems are coupled PDEs and complementarity conditions. We utilize an interior-point Newton method with a filter-line search strategy and solve the Newton linear systems with a block Gauss-Seidel preconditioned Krylov-subspace method. We show that, in certain problem regimes, the number of Krylov-subspace iterations is independent of discretization and ill-conditioning from the bound constraints. We conclude with scaling results on an example problem using the MFEM and hypre packages.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Probabilistic linear solvers are iterative solvers for linear systems that return probability measures. These measures are designed to capture uncertainty in the solution due to incomplete computation. In this talk we review the probabilistic linear solver BayesCG, and discuss ongoing research into applications of the solver. Specifically, we introduce a new application of the solver to the problem of multiple related linear systems, and discuss its role as a means of providing both warm starting and subspace recycling.

A buffet dinner is served daily between 5:30pm and 7:30pm in Vistas Dining Room, top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

We provide a computational framework for approximating a class of structured matrices; here, the term structure is very general. Our matrix-to-tensor invertible mapping allows us to pose the matrix-approximation problem as a tensor-approximation problem. Mapping the tensor approximant back to matrix space, we obtain a structured matrix approximation that can be expressed as a sum of structured Kronecker products or written in block low-rank form. We illustrate the ability of our method to uncover block-low-rank format from applications such as system identification and show we can uncover sum of structured Kronecker products structure on several matrices from the SuiteSparse collection.

The blind source separation problem is a classical signal processing problem. Mathematically it translates into the factorization of $X$, an observed signal matrix, into a mixing and a source signal matrix, i.e., $X=MS^T$. In [Domanov, 2016], the authors propose a criteria checking list to guarantee the decomposition generic uniqueness for a deterministic method class. From this and [De Lathauwer, 2006; Domanov, 2014] results, we design an algebraic algorithm to retrieve the unique factorization of $X$. Moreover, this algorithm finds application in the canonical polyadic tensor decomposition. The structure of the objects involved in the algorithm allows an efficient implementation.

This is joint work with Ignat Domanov and Lieven De Lathauwer.

Randomized algorithms in numerical linear algebra typically generate random vectors from standard distributions, such as Gaussian. However, in certain applications it may be advantageous to run computations with vectors compatible with the underlying structure of the problem. In this talk we discuss algorithms that use rank-one vectors, i.e., Kronecker products of two random vectors, which may allow for faster operations with matrices represented as short sums of Kronecker products. Possible applications include trace estimation and large-scale eigenvalue computation; we provide theoretical and numerical evidence that the use of rank-one instead of unstructured random vectors still leads to good estimates.

Sketching can be seen as a random dimensionality reduction able to preserving the main features of the original problem with probabilistic confidence. Sketching is one of the most promising tools to boost numerical computations although its understanding is still limited. In this talk we present the main features of sketching and how this can be combined with Krylov methods for the solution of large-scale linear systems. However, thanks to the novel sketched Arnoldi relation we will illustrate, the results discussed in this talk can be extended to the numerical evaluation of matrix functions and the solution of matrix equations.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Stellarators are a class of non-axisymmetric toroidal magnetic confinement systems for fusion relevant plasmas. Confinement depends on nested magnetic flux surfaces that act as transport barriers. Flow along magnetic field lines has a Hamiltonian structure, and we generically expect such invariant surfaces to exist in stellarators by KAM theory. However, many devices exhibit not only nested flux surfaces, but also magnetic islands and chaotic regions. Physicists often explore this visually via Poincaré plots that show the structure of the monodromy map for magnetic field line flow around a device. However, to optimize devices we want fast and automatic methods to identify such structure. In this talk, we describe several such analysis methods, with an emphasis on new methods with roots in linear algebra, vector sequence extrapolation, and system identification. We illustrate our method with a variety of examples.

Reliable and efficient iterative solvers for models of steady incompressible flow emerged in the early 1990s. Strategies based on block preconditioning of the underlying matrix operators using (algebraic or geometric) multigrid components have proved to be the key to realising mesh independent convergence (and optimal complexity) without the need for tuning parameters, particularly in the context of classical mixed finite element approximation. The focus of this contribution is on efficient solver strategies in cases where (an inf-sup) stable Taylor-Hood mixed approximation is augmented by a piecewise constant pressure in order to guarantee local conservation of mass. The augmentation leads to over-specification of the pressure solution requiring a redesign of the established solver technology. This enrichment process causes over-specification of the pressure, which complicates the design and implementation of efficient solvers for the resulting linear systems. We first describe the impact of this choice of pressure space on the matrices involved. Next, we show how to recover effective solvers for Stokes problems, using a preconditioner based on the singular pressure mass matrix, and for Oseen systems arising from linearising the Navier-Stokes equations, by using a two-stage pressure convection-diffusion strategy. This is joint work with Jennifer Pestana.

A buffet dinner is served daily between 5:30pm and 7:30pm in Vistas Dining Room, top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

The first part of the talk will provide a general introduction to low rank tensor formats for the data sparse representation of higher order tensors and multivariate functions (CP, Tucker, TT, HT). The second half of the talk is geared towards more recent results on multiparametric multigrid methods and challenges in approximation with regard to special distance measures (e.g. KL) and constraints on the set of admissible tensors (e.g. positivity).

We consider the problem of approximating the von Neumann entropy of a large and sparse symmetric positive semidefinite matrix $A$, defined as $\text{trace}(f(A))$ where $f(z) = -z \log(z)$, which is an important task in several fields, such as quantum information theory and network science. We consider the use of polynomial and rational Krylov subspace algorithms [3] within both randomized trace estimators [2] and probing methods based on graph colorings [1]. We develop error bounds and heuristics that are employed in the implementation of the algorithms. The performance of the methods is illustrated with several numerical experiments. The activities mentioned herein are performed in the framework of the project MUR_PNRR_PE_NQSTI_BELTRAM funded by the Italian Ministry for University and Research (MUR). A preprint is available on arXiv [4]. This is a joint work with Michele Benzi and Michele Rinelli.

References:

[1] A. Frommer, C. Schimmel, and M. Schweitzer, Analysis of probing techniques for sparse approximation and trace estimation of decaying matrix functions, SIAM Journal on Matrix Analysis and Applications, 2021, 42:3, 1290--1318.
[2] R. A. Meyer, C. Musco, C. Musco, and D. P. Woodruff, Hutch++: Optimal Stochastic Trace Estimation, Symposium on Simplicity in Algorithms (SOSA), 2021, 142--155.
[3] S. Güttel, Rational Krylov approximation of matrix functions: numerical methods and optimal pole selection, GAMM-Mitteilungen, 2013, 36: 8-31.
[4] M. Benzi, M. Rinelli, I. Simunec, Computation of the von Neumann entropy of large matrices via trace estimators and rational Krylov methods, arXiv preprint arXiv:2212.09642 (2022).

5-day workshop participants are welcome to use BIRS facilities (TCPL ) until 3 pm on Friday, although participants are still required to checkout of the guest rooms by 11AM.

We consider the solution of linear systems arising from \emph{parameter robust} discretizations of singularly perturbed convection-diffusion problems, whose accuracy is proven to be independent of the perturbation parameter. There are many studies of such methods, but few of the robust solution of resulting linear systems. I'll discuss challenges for both direct and iterative solvers for such systems, and propose a multigrid-based preconditioning strategy tuned to the matrix structure induced by using layer-adapted meshes. Supporting theoretical and numerical results will be presented.

This is joint work with Anh Thài Nhan (SCU) and Scott MacLachlan (MUN): SIMAX 43 (2022).

In this talk, we show how a recently introduced three-level variant of Parareal can be interpreted as a multigrid reduction in time (MGRIT) method for a particular choice of restriction/prolongation operators and relaxation scheme. Specifically, we show that when suitably initialized, the MGRIT V-cycle with F-relaxation and injection as the prolongation operator generates iterates that are identical to three-level Parareal. This equivalence allows one to analyze the convergence of three-level Parareal from two different perspectives, using either ODE or multigrid techniques.

This is joint work with Stephanie Friedhoff and Martin J. Gander