Schedule for: 23w5108 - Perspectives on Matrix Computations: Theoretical Computer Science Meets Numerical Analysis

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.

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 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!

Solving systems of non-autonomous ordinary differential equations (ODE) is a crucial and often challenging problem. Recently a new approach was introduced based on a generalization of the Volterra composition. The method transforms the ODE system into the solution of a structured linear system. Numerical examples illustrate the method's efficacy and its properties.

A common task in NLA is low-rank approximation, which can be done by approximating the top eigenvectors of a matrix. Block Krylov Iteration gives the best known theoretical guarantees and empirical performance for this task. Though, it's conceptually unclear how we should choose the block size. The best theoretical guarantees in the prior work requires using large block sizes, at least block size k when finding a rank-k approximation -- and possibly much larger. However, in practice, a small constant block size independent of k typically performs much better than such large block sizes. When the block size is exactly 1, we call this "Single Vector Krylov". We resolve this unclarity and this theory-practice gap by proving that, when counting matrix-vector products instead of iteration complexity, Single Vector Krylov matches the rate of convergence of the best possible choice of block size, up to just a mild logarithmic dependence on eigenvalue gaps. Moreover, for a wide family of matrices, Single Vector Krylov converges faster, improving from $O( k \cdot log(\frac1\varepsilon) )$ MatVecs to $O( k+log(\frac1\varepsilon) )$ MatVecs. The core of our proof is the observation that the Krylov Subpsace generated by Single Vector Krylov initialized from random vector is exactly equal to the Krylov Subspace generated by Block Krylov Iteration initialized to a structured starting block. We conclude with implications of this theory for simplifying algorithms, smoothed analysis, and experimental validation.

Linear systems with a tensor product structure arise naturally when considering the discretization of Laplace type differential equations or, more generally, multidimensional operators with separable coefficients. In this work, we focus on the numerical solution of linear systems of the form $$ \left(A_1\otimes I\otimes\dots\otimes I+\dots + I\otimes\dots I\otimes A_d\right)x=b, $$ where the matrices $A_t\in\mathbb R^{n\times n}$ are Hermitian positive definite and belong to the class of hierarchically semiseparable matrices. We propose and analyze a nested divide-and-conquer scheme, based on the technology of low-rank updates, that attains the quasi-optimal computational cost $\mathcal O(n^d\log(n))$. Our theoretical analysis highlights the role of inexactness in the nested calls of our algorithm and provides worst case estimates for the amplification of the residual norm. The performances are validated on 2D and 3D case studies.

Over the past two decades, areas of algorithms traditionally considered combinatorial have benefitted enormously from the study of provably worst-case efficient numerical solvers for structured linear systems. This talk aims to give a picture of the current state-of-the-art for solving linear systems under various structural assumptions, resource goals, and representations of numbers.

We study several variants of decomposing a symmetric matrix into a sum of a low-rank positive semidefinite matrix and a diagonal matrix. Such decompositions have applications in factor analysis and they have been studied for many decades. On the one hand, we prove that when the rank of the positive semidefinite matrix in the decomposition is bounded above by an absolute constant, the problem can be solved in polynomial time. On the other hand, we prove that, in general, these problems as well as their certain approximation versions are all NP-hard. Finally, we prove that many of these low-rank decomposition problems are complete in the first-order theory of the reals; i.e., given any system of polynomial equations, we can write down a low-rank decomposition problem in polynomial time so that the original system has a solution iff our corresponding decomposition problem has a feasible solution of certain (lowest) rank.

In this talk, we study the growth factor in Gaussian elimination under complete pivoting and show a number of new results, including improved lower bounds, estimates for growth in floating point arithmetic, and bounds for growth for matrices with restricted entries.

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.

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

A framework is introduced for fitting trigonometric rational functions to data in a way that is robust to various forms of corruption, including additive Gaussian noise, perturbed sampling grids, and missing data. Our approach combines a variant of Prony’s method with a modified version of the AAA algorithm. Using representations in both frequency and time space, a collection of algorithms is described for adaptively computing with trigonometric rationals. This includes procedures for algebraic operations, differentiation, convolution, and more.

We will discuss recent theoretical results about the global rapid convergence of the QR algorithm, which in particular imply the existence of a randomized diagonalization algorithm that runs in $\tilde{O}(n^3)$ operations on any arbitrary input. This is joint work with Jess Banks and Nikhil Srivastava.

TBA

This talk explores the limitations and the successes of understanding algorithms via their performance on random matrices. For example, the conjugate gradient algorithm produces highly concentrated output from random input, indicating robustness. On the other hand, the power iteration produces much less concentrated output. This can be explained via the dependence of the algorithms on modified moments versus classical moments for a matrix-dependent measure. Yet, random input can give potentially misleading conclusions as the well-known instabilities of the Lanczos iteration are often missing in random situations. This is joint work with Percy Deift, Tyler Chen, Elliot Paquette and Xiucai Ding.

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.

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

Join for 1hr long Tour of the Iconic Banff Centre. Jim Olver will meet you outside the Professional Development Centre. Dress for the weather.

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.

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.

In recent years, stochastic trace estimation has proven to be a powerful algorithmic tool in numerical linear algebra, and theoretical computer science, machine learning, and computational quantum physics. However, despite the popularity of such estimators, the literature on the topic is highly fragmented. Indeed, while NLA and TCS typically attribute the estimators to several papers appearing around 1990, algorithms based on the same ideas date back to at least the 1970s. Moreover, stochastic trace estimation is, in some moral sense, a special instance of quantum typicality, an idea dating back at least to work of Schr\"odinger and von Neumann in the late 1920s and which regained attention in the mid 2000s. This talk uses trace estimation/typicality as a case study in research at the intersection of nearby fields. We discuss provide some history into the topic, as well as speculate as to whether this level of fragmentation is an inherent consequence of modern research, or whether it can been mitigated.

This work is concerned with computing low-rank approximations of a matrix function $f(A)$ for a large symmetric positive semi-definite matrix $A$, a task that arises in, e.g., statistical learning and inverse problems. The application of popular randomized methods, such as the randomized singular value decomposition or the Nyström approximation, to $f(A)$ requires multiplying $f(A)$ with a few random vectors. A significant disadvantage of such an approach, matrix-vector products with $f(A)$ are considerably more expensive than matrix-vector products with $A$, even when carried out only approximately via, e.g., the Lanczos method. In this work, we present and analyze funNyström, a simple and inexpensive method that constructs a low-rank approximation of $f(A)$ directly from a Nyström approximation of $A$, completely bypassing the need for matrix-vector products with $f(A)$. It is sensible to use funNyström whenever $f$ is monotone and satisfies $f(0) = 0$. Under the stronger assumption that $f$ is operator monotone, which includes the matrix square root $A^{1/2}$ and the matrix logarithm $\log(I+A)$, we derive probabilistic bounds for the error in the Frobenius, nuclear, and operator norms. These bounds confirm the numerical observation that funNyström tends to return an approximation that compares well with the best low-rank approximation of $f(A)$. Furthermore, compared to existing methods, funNyström requires significantly fewer matrix-vector products with $A$ to obtain a low-rank approximation of $f(A)$, without sacrificing accuracy or reliability. Our method is also of interest when estimating quantities associated with $f(A)$, such as the trace or the diagonal entries of $f(A)$. In particular, we propose and analyze funNyström++, a combination of funNyström with the recently developed Hutch++ method for trace estimation.

We describe a Lanczos or Arnoldi-based algorithm for approximating the product of a rational function of a matrix $A$ with a given vector $b$. This algorithm is optimal over the Krylov subspace, in a norm induced by the denominator of the rational function, and the approximation can be computed using information from a slightly larger Krylov subspace. It requires storage of just a few additional vectors, proportional to the degree of the denominator of the rational function. It can be used to construct near-optimal approximations for other matrix functions such as the matrix sign function.

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.