Session III.1 - Numerical Linear Algebra
Talks
Monday, June 19, 14:00 ~ 14:30
Universal Matrix Sparsifiers and Fast Deterministic Algorithms for Singular Value Approximation
Cameron Musco
University of Massachusetts Amherst, USA - This email address is being protected from spambots. You need JavaScript enabled to view it.
Given a matrix $A \in \mathbb{R}^{n \times n}$ which is normalized so that its entries are bounded in magnitude by $1$, it is well-known that randomly sparsifying $A$ to be supported on a subset of just $\tilde{O}(n/\epsilon^2)$ of its entries, yields an approximation $A_S$ with $\| A - A_S\|_2 \le \epsilon n$ with high probability, where $\|\cdot\|_2$ is the spectral norm. We show that for positive semidefinite (PSD) matrices, no randomness is needed at all in this statement. Namely, there exists a fixed subset of $\tilde{O}(n/\epsilon^2)$ entries that acts as a $universal\, sparsifier$: $\| A - A_S\|_2 \le \epsilon n$ holds simultaneously for every bounded entry PSD matrix. One can view this result as a significant extension of a spectral graph expander nearly matching the Ramanujan bound.
We leverage the existence of such universal sparsifiers to give the first $deterministic\, algorithms$ for several central linear algebraic problems, including singular value and singular vector approximation and positive semidefiniteness testing, that run in faster than matrix multiplication time on worst-case input matrices. This partially addresses a significant gap between randomized and deterministic algorithms for fast linear algebraic computation. We also extend our results to non-PSD matrices and to deterministic algorithms that approximate $A$ with a non-sparse matrix. We prove lower bounds showing that many of our results are nearly optimal.
Joint work with Rajarshi Bhattacharjee (University of Massachusetts Amherst), Gregory Dexter (Purdue University), Archan Ray (University of Massachusetts Amherst) and David P. Woodruff (Carnegie Mellon University).
Monday, June 19, 14:30 ~ 15:00
XTrace: Making the most of every sample in stochastic trace estimation
Joel Tropp
Caltech, USA - This email address is being protected from spambots. You need JavaScript enabled to view it.
The implicit trace estimation problem asks for an approximation of the trace of a square matrix, accessed via matrix--vector products (matvecs). This talk introduces new randomized algorithms, \textsc{XTrace} and \textsc{XNysTrace}, for the trace estimation problem by exploiting both variance reduction and the exchangeability principle. For a fixed budget of matvecs, numerical experiments show that the new methods can achieve errors that are orders of magnitude smaller than existing algorithms. A theoretical analysis confirms the benefits by offering a precise description of the performance of these algorithms as a function of the spectrum of the input matrix.
Available from arXiv:2301.07825.
Joint work with Ethan N. Epperly (Caltech) and Robert J. Webber (Caltech).
Monday, June 19, 15:00 ~ 15:30
Solving the Parametric Eigenvalue Problem by Taylor Series and Chebyshev Expansion
Thomas Mach
University of Potsdam, Germany - This email address is being protected from spambots. You need JavaScript enabled to view it.
We discuss two approaches to solving the parametric (or stochastic) eigenvalue problem. One of them uses a Taylor expansion and the other a Chebyshev expansion. The parametric eigenvalue problem assumes that the matrix $A$ depends on a parameter $\mu$, where $\mu$ might be a random variable. Consequently, the eigenvalues and eigenvectors are also functions of $\mu$. We compute a Taylor approximation of these functions about $\mu_{0}$ by iteratively computing the Taylor coefficients. The complexity of this approach is $O(n^{3})$ for all eigenpairs, if the derivatives of $A(\mu)$ at $\mu_{0}$ are given. The Chebyshev expansion works similarly. We first find an initial approximation iteratively which we then refine with Newton's method. This second method is more expensive but provides a good approximation over the whole interval of the expansion instead around a single point. We present numerical experiments confirming the complexity and demonstrating that the approaches are capable of tracking eigenvalues at intersection points. Further experiments shed light on the limitations of the Taylor expansion approach with respect to the distance from the expansion point $\mu_{0}$.
This presentation is based on a paper available on arXiv, https://arxiv.org/abs/2302.03661, that has been submitted for publication.
Joint work with Melina A. Freitag (University of Potsdam, Germany).
Monday, June 19, 15:30 ~ 16:00
Perfect Shifted QR for Rank Structured Pencils
Raf Vandebril
KU Leuven, Belgium - This email address is being protected from spambots. You need JavaScript enabled to view it.
It is known that executing a perfect shifted QR step via the implicit QR algorithm performs poor in terms of stability, and accuracy. Typically several steps are required before deflation actually takes place. This behavior can be remedied by determining the similarity transformation via the associated eigenvector. Similar techniques can be deduced for the QZ algorithm and the rational QZ algorithm. In this talk we generalize this even further and present an approach for executing a perfect shifted QR step for general rank structured pencils. We prove that the rank structures of the matrices involved in the pencil are preserved, and moreover we examine the preservation of spectral properties of subblocks in the pencil.
Joint work with Nicola Mastronardi, Marc Van Barel, Raf Vandebril, and Paul Van Dooren.
Monday, June 19, 16:30 ~ 17:00
Speeding up Krylov subspace methods for matrix functions via randomization
Alice Cortinovis
Stanford University, United States - This email address is being protected from spambots. You need JavaScript enabled to view it.
In this talk we consider the computation of the action of a matrix function f(A), such as the matrix exponential or the matrix square root, on a vector b. For a general matrix A, this can be done by computing the compression of A onto a suitable Krylov subspace. Such compression is usually computed by forming an orthonormal basis of the Krylov subspace using the Arnoldi method. In this talk, we propose to compute (non-orthonormal) bases in a faster way and to use a fast randomized algorithm for least-squares problems to compute the compression of A onto the Krylov subspace. We present some numerical examples which show that our algorithms can be faster than the standard Arnoldi method while achieving comparable accuracy.
Joint work with Daniel Kressner (EPFL, Switzerland) and Yuji Nakatsukasa (University of Oxford, UK).
Monday, June 19, 17:00 ~ 17:30
Gradient-Type Subspace Iteration Methods for Symmetric Eigen-Problem
Yousef Saad
University of Minnesota, USA - This email address is being protected from spambots. You need JavaScript enabled to view it.
Computing invariant subspaces is at the core of many applications, from machine learning to signal processing, and control theory, to name just a few examples. Often one wishes to compute the subspace associated with eigenvalues located at one end of the spectrum, i.e., either the largest or the smallest eigenvalues. It is also quite common that the data at hand undergoes frequent changes and that one is required to keep updating or tracking the target invariant subspace. For such problems, the subspace iteration algorithm is ideally suited. The talk will explores variants of the subspace iteration algorithm for computing approximate invariant subspaces. The talk will start by reviewing the subspace iteration approach and then discuss a few variants that exploit gradient-type techniques combined with a Grassmann manifold viewpoint. A gradient descent/ascent method as well as a conjugate gradient technique will be described and analyzed. Numerical experiments will illustrate the behavior of the algorithm on a few examples.
Joint work with Bart Vandereycken (University of Geneva) and Foivos Alimisis (University of Geneva).
Monday, June 19, 17:30 ~ 18:00
The limit empirical spectral distribution of random matrix polynomials
Vanni Noferini
Aalto University, Finland - This email address is being protected from spambots. You need JavaScript enabled to view it.
Three famous classic results concern the distributions of the roots of a random polynomial and the eigenvalues of a random matrix or pencil. Under relatively mild assumptions on the distribution of the coefficients, the former is known to converge to the uniform distribution on the unit circle when the degree $k$ approaches infinity. Under similarly unrestrictive assumptions on the distributions of the entries, the latter is known to converge to the uniform distribution on the unit disk (when the entries have mean $0$ and variance $1/n$) when the size $n$ approaches infinity. Several mathematicians have also independently derived the distribution of the generalised eigenvalues of a random pencil: in this case, the distribution is uniform on the Riemann sphere. However, until the present work nothing was to my knowledge known about the eigenvalues of a general random matrix polynomial, which could be thought as an intermediate case between a random matrix or pencil ($k=1$) and a random polynomial ($n=1$).
In this talk I plan to first give some gentle introduction, thought for non-experts on random variables and random matrices, to the known results mentioned above. I will then move on to describe recent new results that we obtained about the limit spectral distributions of a random matrix polynomial, both in the regime $k \rightarrow \infty$ and in the case $n \rightarrow \infty$. After discussing the (easier) nonmonic case, I will also comment on what changes if the random matrix polynomial is assumed to be monic, i.e., having all random coefficients except the leading one which is taken to be the identity matrix.
The main tools from our results come both from random matrix theory and from deterministic matrix theory. In particular, we exploit (1) the replacement principle, proven by Tao, Vu and Krishnapur, (2) the logarithmic potential approach as proposed by Girko, and (3) classic perturbation theory results.
Joint work with Giovanni Barbarino.
Tuesday, June 20, 14:00 ~ 14:30
On the Effectiveness of Single Vector Krylov Methods for Low-Rank Approximation
Christopher Musco
New York University - This email address is being protected from spambots. You need JavaScript enabled to view it.
Krylov subspace methods are a ubiquitous tool for computing near-optimal rank $k$ approximations of large matrices. While ``large block'' Krylov methods with block size at least $k$ give the best known theoretical guarantees, block size one (a single vector) or a small constant is often preferred in practice. Despite their popularity, we lack theoretical bounds on the performance of such ``small block'' Krylov methods for low-rank approximation.
In this talk I will talk about a recent result that addresses this gap between theory and practice by proving that small block Krylov methods essentially match all known low-rank approximation guarantees for large block methods. Via a black-box reduction we show, for example, that the standard single vector Krylov method run for $t$ iterations obtains the same spectral norm and Frobenius norm error bounds as a Krylov method with block size $\ell \geq k$ run for $O(t/\ell)$ iterations, up to a logarithmic dependence on the smallest gap between sequential singular values. That is, for a given number of matrix-vector products, single vector methods are essentially as effective as the best choice of large block size.
By combining our result with tail-bounds on eigenvalue gaps in random matrices, we prove that the dependence on the smallest singular value gap can be eliminated if the input matrix is perturbed by a small random matrix.
Joint work with Raphael Meyer (New York University) and Cameron Musco (University of Massachusetts Amherst).
Tuesday, June 20, 14:30 ~ 15:00
Condition numbers related to low-rank matrix approximation
Nick Vannieuwenhoven
KU Leuven, Belgium - This email address is being protected from spambots. You need JavaScript enabled to view it.
I will present the Rice condition number of the best rank-$r$ matrix approximation problem, i.e., of computing the rank-$r$ truncated singular value decomposition. This condition number measures, in the Frobenius norm, how much an infinitesimal perturbation to an arbitrary input matrix is amplified in the perturbation of the best rank-$r$ approximation. I also present a Rice condition number of the problem of taking a matrix to a selected left or right singular subspace, or, equivalently, to a projector thereon. This condition number measures the input error in the Frobenius norm and the output error in the chordal, Grassmann, and Procrustes distances on the Grassmannian manifold of linear subspaces, rigorizing and extending prior, classic results in numerical linear algebra. We are able to derive explicit closed-form formulas of these condition numbers in terms of the singular values of the input matrix for both of these classic numerical linear algebra problems through an interesting voyage in Riemannian geometry.
Joint work with Paul Breiding (University of Osnabruck, Germany) in the first part.
Tuesday, June 20, 15:00 ~ 15:30
A nested divide-and-conquer method for tensor Sylvester equations with positive definite hierarchically semiseparable coefficients
Leonardo Robol
University of Pisa, Italy - This email address is being protected from spambots. You need JavaScript enabled to view it.
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(I\otimes \dots\otimes I \otimes A_1+\dots + A_d\otimes I \otimes\dots \otimes I\right)x=b, $$ where the matrices $A_t\in\mathbb R^{n\times n}$ are symmetric 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) + \log(\kappa)^2 + \log(\kappa) \log(\epsilon^{-1})))$ where $\kappa$ is the condition number of the linear system, and $\epsilon$ the target accuracy. 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.
Joint work with Stefano Massei (University of Pisa).
Tuesday, June 20, 15:30 ~ 16:00
Structured matrix recovery from matrix-vector products
Alex Townsend
Cornell University, United States - This email address is being protected from spambots. You need JavaScript enabled to view it.
Can one recover a structured matrix efficiently from only matrix-vector products? If so, how many are needed? In this talk, we will describe algorithms to recover structured matrices, such as tridiagonal, Toeplitz-like, and hierarchical low-rank, from matrix-vector products. In particular, we derive an algorithm for recovering an unknown $N\times N$ hierarchical low-rank matrix from only $\mathcal{O}((k+p)\log N)$ matrix-vector products that is stable with high probability, where $k$ is the rank of the off-diagonal blocks, and $p$ is a small oversampling parameter. We do this by carefully constructing randomized input vectors for our matrix-vector products that use the hierarchical structure of the matrix. While existing algorithms for hierarchical matrix recovery use a recursive “peeling" procedure based on elimination, our approach uses a recursive projection procedure,
Joint work with Diana Halikias (Cornell University).
Tuesday, June 20, 16:30 ~ 17:00
Randomization techniques for solving linear systems of equations and eigenvalue problems
Laura Grigori
EPFL, Switzerland - This email address is being protected from spambots. You need JavaScript enabled to view it.
In this talk we discuss recent progress in using randomization for solving linear systems of equations and eigenvalue problems. We present a randomized version of the Gram-Schmidt process for orthogonalizing a set of vectors and its usage in the Arnoldi iteration. This leads to introducing new Krylov subspace methods for solving large scale linear systems of equations and eigenvalue problems. The new methods retain the numerical stability of classic Krylov methods while reducing communication and being more efficient on modern massively parallel computers.
Joint work with O. Balabanov (Alpines, Inria Paris and LJLL, Sorbonne University), M. Beaupere (Alpines, Inria Paris and LJLL, Sorbonne University), V. Lederer (Alpines, Inria Paris and LJLL, Sorbonne University) and E. Timsit (Alpines, Inria Paris and LJLL, Sorbonne University).
Tuesday, June 20, 17:00 ~ 17:30
A new backward error analysis framework for GMRES
Bastien Vieublé
The University of Manchester, United Kingdom - This email address is being protected from spambots. You need JavaScript enabled to view it.
GMRES is a well-known iterative solver for the solution of linear systems $Ax=b$. Its combination with approximate computing techniques such as low precision or randomization has recently brought much attention. This is because these techniques can offer significant gains in speed and energy and enlarge the scale of problems we are able to solve. GMRES can employ approximate computing in many different ways: at the preconditioner level, at the restart level, at the orthogonalization level, etc. Moreover, considering the wide variety of preconditioners that can be used (e.g., ILU, polynomial, Jacobi, iterative solver, etc.) and the different orthogonalization algorithms and their variants (e.g., Householder, MGS, CGS2, etc.), the combinatory of possibilities is gigantic. Generally, existing backward error analyses of GMRES are specialized to one of these combinations and cannot be extended to others. For this reason, in this talk, we present a general framework for backward error analysis of GMRES that aims at gathering all these combinations inside a common and coherent analysis, in addition to deriving specialized bounds on the backward and forward errors of a given GMRES algorithm. We finally show how to apply these tools on a new mixed precision preconditioned restarted GMRES algorithm that uses six independent precision parameters and applies the preconditioner in a lower precision than the matrix-vector product with $A$. We motivate this new mixed precision scheme by numerical experiments on real-life SuiteSparse matrices.
Joint work with Alfredo Buttari (IRIT, France), Nicholas J. Higham (The University of Manchester, UK) and Théo Mary (LIP6, France).
Tuesday, June 20, 17:30 ~ 18:00
Construction of Lobatto and Kronrod rules from orthogonal Laurent polynomials
Carl Jagels
Hanover College, United States - This email address is being protected from spambots. You need JavaScript enabled to view it.
The construction of Gaussian rules can be viewed as a spectral decomposition of the tridiagonal Jacobi matrix generated by the Lanczos process. Radau and Lobatto rules follow from a modification of this matrix. Analogous rules exist for Laurent polynomials, polynomials that contain reciprocal powers. The analog of the tridiagonal recursion matrix is a pentadiagonal matrix. This talk discusses augmentations of the pentadiagonal matrix that yield Lobatto and Kronrod rules.
Tuesday, June 20, 18:00 ~ 18:30
Are randomized NLA algorithms numerically stable?
Yuji Nakatsukasa
University of Oxford, United Kingdom - This email address is being protected from spambots. You need JavaScript enabled to view it.
Randomized algorithms in numerical linear algebra can provide efficient means for, for example, computing low-rank approximations, and solving least-squares problems, linear systems, and eigenvalue problems. An aspect that often gets overlooked is: are these methods numerically (backward) stable? In this talk I will highlight that the situation is quite subtle: for example, sketch-to-precondition methods for least-squares are not stable, and methods for low-rank approximation may require a careful implementation to be stable. Part of the talk is based on joint work with Maike Meier (University of Oxford, UK), Alex Townsend (Cornell University, USA), Marcus Webb (University of Manchester, UK).
Wednesday, June 21, 14:00 ~ 14:30
On the optimal value for the conditioning of the eigenvector problem
Carlos Beltrán
Universidad de Cantabria, Spain - This email address is being protected from spambots. You need JavaScript enabled to view it.
The condition number for eigenvalue--eigenvector pair computations is a well--studied quantity measuring how much the solution changes when the input is slightly perturbed. The eigenvalues of symmetric matrices are optimally conditioned, but, which is a perfectly conditioned matrix w.r.t. eigenvector computations? In this talk we answer this nontrivial question with exact first order asymptotic.
The author is partially supported by Grant PID2020-113887GB-I00 funded by MCIN/ AEI /10.13039/501100011033
Joint work with Joint work with Laurent Bétermin (Université Claude Bernard Lyon 1), Peter Grabner (Technische Universität Graz) and Stefan Steinerberger (University of Washington).
Wednesday, June 21, 14:30 ~ 15:00
Avoiding discretization issues for nonlinear eigenvalue problems
Matthew Colbrook
University of Cambridge, UK - This email address is being protected from spambots. You need JavaScript enabled to view it.
The first step when solving an infinite-dimensional eigenvalue problem is often to discretize it. We show that one must be extremely careful when discretizing nonlinear eigenvalue problems. Using examples, we show that discretization can: (1) introduce spurious eigenvalues, (2) entirely miss spectra, and (3) bring in severe ill-conditioning. While there are many eigensolvers for solving matrix nonlinear eigenvalue problems, we propose a solver for general holomorphic infinite-dimensional nonlinear eigenvalue problems that avoids discretization issues, which we prove is stable and converges. Moreover, we provide an algorithm that computes the problem's pseudospectra with explicit error control, allowing verification of computed spectra. The algorithm and numerical examples are publicly available in infNEP, which is a software package written in MATLAB. Finally, we show how this work fits into the bigger picture of recent progress in the foundations of infinite-dimensional computations.
Joint work with Alex Townsend (Cornell University).
Wednesday, June 21, 15:00 ~ 15:30
Strongly minimal self-conjugate linearizations for polynomial and rational matrices
Froilán Dopico
Universidad Carlos III de Madrid, Spain - This email address is being protected from spambots. You need JavaScript enabled to view it.
We present the concept of strongly minimal linearization and prove that we can always construct strongly minimal linearizations of an arbitrary rational matrix from its Laurent expansion around the point at infinity, which happens to be the case for polynomial matrices expressed in the monomial basis. If the rational matrix has a particular self-conjugate structure we show how to construct strongly minimal linearizations that preserve it. The structures that are considered are the Hermitian and skew-Hermitian rational matrices with respect to the real line, and the para-Hermitian and para-skew-Hermitian matrices with respect to the imaginary axis. We pay special attention to the construction of strongly minimal linearizations for the particular case of structured polynomial matrices. The proposed constructions lead to efficient numerical algorithms for constructing strongly minimal linearizations. The fact that they are valid for any rational matrix is an improvement on any other previous approach for constructing other classes of structure preserving linearizations, which are not valid for any structured rational or polynomial matrix. The use of the recent concept of strongly minimal linearization is the key for getting such generality.
Joint work with María C. Quintana (Aalto University, Finland) and Paul Van Dooren (Université catholique de Louvain, Belgium).
Wednesday, June 21, 15:30 ~ 16:00
Randomized matrix-free quadrature
Tyler Chen
New York University, USA - This email address is being protected from spambots. You need JavaScript enabled to view it.
We discuss randomized matrix-free quadrature algorithms for spectrum and spectral sum approximation. The algorithms studied are characterized by the use of a Krylov subspace method to approximate independent and identically distributed samples of $\mathbf{v}^*f(\mathbf{A})\mathbf{v}$, where $\mathbf{v}$ is an isotropic random vector, $\mathbf{A}$ is a Hermitian matrix, and $f(\mathbf{A})$ is a matrix function. This class of algorithms includes the kernel polynomial method and stochastic Lanczos quadrature, two widely used methods for approximating spectra and spectral sums. We will provide a unified framework for understanding these algorithms, and provide examples which shed light on the commonalities and tradeoffs between them.
Joint work with Thomas Trogdon (University of Washington) and Shashanka Ubaru (IBM Watson).
Wednesday, June 21, 16:30 ~ 17:00
Role extraction for digraphs via neighbourhood pattern similarity
Giovanni Barbarino
aalto university, Finland - This email address is being protected from spambots. You need JavaScript enabled to view it.
The analysis of large graphs frequently assumes that there is an underlying structure in the graph that allows us to represent it in a simpler manner. A typical example of this is the detection of communities, which groups together nodes with similar connectivity properties. Yet, many graph structures cannot be modelled using communities: for example, arrowhead and tree graph structures, which appear in overlapping communities, human protein-protein interaction networks, and food and web networks. These more general types of network structures can be modelled as role structures, and the process of finding them is called the role extraction problem or block modelling. In this presentation, we show that a particular spectral method, using the so-called Neighbourhood Pattern Similarity matrices, allows us to recover asymptotically the block structure of a general directed graph with a stochastic block model (SBM) structure. In fact, the increasing gap between the dominant and dominated eigenvalues of the similarity matrix guarantees the asymptotic correct identification of the number of different roles and enables us to compute an asymptotic estimation of the total misclassification error. This is in particular one of the few existing results of correctness for the spectral clustering algorithm applied to directed graph with a SBM structure.
Joint work with Vanni Noferini (Aalto University, Finland) and Paul Van Dooren (Université catholique de Louvain, Belgium).
Wednesday, June 21, 17:00 ~ 17:30
What part of a numerical problem is ill-conditioned?
Nick Dewaele
KU Leuven, Belgium - This email address is being protected from spambots. You need JavaScript enabled to view it.
Many numerical problems with input $x$ and output $y$ can be formulated as an system of equations $F(x, y) = 0$ where the goal is to solve for $y$. The condition number measures the change of $y$ for small perturbations to $x$. From this numerical problem, one can derive a (typically underdetermined) subproblem by omitting any number of constraints from $F$. We propose a condition number for underdetermined systems that relates the condition number of a numerical problem to those of its subproblems. We illustrate the use of our technique by computing the condition of two numerical linear algebra problems that do not have a condition number in the classical sense: the decomposition of a low-rank matrix into unstructured factors and Tucker decomposition.
Joint work with Nick Vannieuwenhoven (KU Leuven, Belgium).
Wednesday, June 21, 17:30 ~ 18:00
Randomized Contour Integral Methods for Eigenvalue Problems
Agnieszka Miedlar
Virginia Tech, United States of America - This email address is being protected from spambots. You need JavaScript enabled to view it.
Randomized NLA methods have recently gained popularity because of their easy implementation, computational efficiency, and numerical robustness. In this talk, we present the analysis of a randomized version of a well-established FEAST algorithm that enables computing the eigenvalues of the Hermitian matrix pencil $(A, B)$ located in the given real interval $\mathcal{I} \subset [\lambda_{min}, \lambda_{max}]$. First, we establish new structural as well as probabilistic error analysis of the accuracy of approximate eigenpairs and subspaces obtained using the randomized FEAST algorithm, i.e., bounds for the canonical angles between the exact and the approximate eigenspaces, and for the accuracy of the eigenvalues and the corresponding eigenvectors. Since this part of the analysis is independent of the particular distribution of an initial subspace, we denote it as structural. In the case of the starting guess being a Gaussian random matrix, we provide more informative, probabilistic error bounds. Our modified algorithm allows to improve the accuracy of orthogonalization when $B$ is ill-conditioned, efficiently apply the rational filter by using MPGMRES-Sh [Bakhos, Kitanidis, Ladenheim, Saibaba and Szyld, 2016] method to accelerate solving shifted linear systems and estimate the eigenvalue counts in a given interval. Finally, we illustrate numerically the effectiveness of presented error bounds and proposed algorithmic modifications.
Posters
Sparse and symmetry-preserving compression of tensor trains arising in quantum chemistry
Siwar Badreddine
Sorbonne university, France - This email address is being protected from spambots. You need JavaScript enabled to view it.
The quantum chemical density matrix renormalization group (QC-DMRG) is regarded as a powerful method for the resolution of the ground state energy of the many-body electronic Schrödinger equation. The construction of the quantum chemical Hamiltonian operator in matrix product or tensor train (TT) form is at the core of the QC-DMRG algorithm. Although the ranks of the exact TT representation grow quadratically with the system size $d$, $O(d^2)$, it can be reduced using the TT-SVD algorithm. However, this might destroy the sparsity as well as the symmetries of the original operator if it is not carefully done.
This study shows that the TT format of the Hamiltonians in QC-DMRG has an inherent structure that can be exploited to get a sparse and symmetry-preserving scheme. It is well known that the quantum chemical Hamiltonian has a variety of symmetries that can be exploited to design a more efficient TT representa- tion. This not only results in reducing the computational cost but also preserves the properties of the original operator that are crucial for an accurate simulation. Therefore, we demonstrate that enforcing abelian symmetries such as the particle number results in a sparse block structured form in the TT cores of the TT representation of the Hamiltonian operator, leading to a symmetry-preserving, more structured, and organized representation, allowing for faster computation and manipulation.
In conclusion, this study highlights the potential for exploiting the new structure of the TT-cores of the TT representation of the Hamiltonian operator which can be used to reduce the cost of operations involved in the QC-DMRG algorithm.
Joint work with Eric Cances (Ecole des Ponts ParisTech and Inria), Mi-Song Dupuy ( Laboratoire Jacques-Louis LIONS in Sorbonne Université (Paris, France)) and Laura Grigori (Laboratoire Jacques-Louis LIONS in Sorbonne Université (Paris, France) and Inria).
Randomized Joint Diagonalization of Symmetric Matrices
Haoze He
École Polytechnique Fédérale de Lausanne (EPFL), Switzerland - This email address is being protected from spambots. You need JavaScript enabled to view it.
Given a family of nearly commuting symmetric matrices, we consider the task of computing an orthogonal matrix that nearly diagonalizes every matrix in the family. In this work, we propose and analyze randomized joint diagonalization (RJD) for performing this task. RJD applies a standard eigenvalue solver to random linear combinations of the matrices. Unlike existing optimization-based methods, RJD is simple to implement and leverages existing high-quality linear algebra software packages. Our main novel contribution is to prove robust recovery: Given a family that is $\epsilon$-close to a commuting family, RJD jointly diagonalizes this family, with high probability, up to an error of norm $\mathcal{O}(\epsilon)$. No other existing method is known to enjoy such a universal robust recovery guarantee. We also discuss how the algorithm can be further improved by deflation techniques and demonstrate its state-of-the-art performance by numerical experiments with synthetic and real-world data.
Joint work with Daniel Kressner (EPFL).
Fast and direct inverse nonequispaced Fourier transforms
Melanie Kircheis
Chemnitz University of Technology, Germany - This email address is being protected from spambots. You need JavaScript enabled to view it.
The well-known discrete Fourier transform (DFT) can easily be generalized to arbitrary nodes in the spatial domain. The fast procedure for this generalization is referred to as nonequispaced fast Fourier transform (NFFT). Various applications such as MRI, solution of PDEs, etc., are interested in the inverse problem, i.,e., computing Fourier coefficients from given nonequispaced data. In contrast to iterative procedures, where multiple iteration steps are needed for computing a solution, we focus especially on so-called direct inversion methods. We review density compensation techniques and introduce a new scheme that leads to an exact reconstruction for trigonometric polynomials. In addition, we consider a matrix optimization approach using Frobenius norm minimization to obtain an inverse NFFT.
Joint work with Daniel Potts (Chemnitz University of Technology, Germany).
Fast computations with arrowhead and diagonal-plus-rank-k matrices over associative fields
Ivan Slapničar
University of Split, FESB, Croatia - This email address is being protected from spambots. You need JavaScript enabled to view it.
We consider arrowhead and diagonal-plus-rank-$k$ matrices (DPRk) in $\mathbb{F}^{n\times n}$ where $\mathbb{F}\in\{\mathbb{R},\mathbb{C},\mathbb{H}\}$, and $\mathbb{H}$ is a non-commutative field of quaternions. Such matrices arise in many applications and computations with such matrices are parts of important linear algebra algorithms. We give formulas for matrix-vector multiplications, determinants, and inverses for both types of matrices. The formulas are unified in the sense that the same formulas hold in both, commutative and noncommutative associative algebras. Each formula requires $O(n)$ arithmetic operations. Most of the formulas hold for block matrices, as well. Further, we derived efficient $O(n^2)$ eigensolvers for quaternionic arrowhead and DPRk matrices. The eigensolvers use a version of Wielandt deflation technique. All algorithms are elegantly implemented using the polymorphism feature of the modern high-level computing language Julia. The code is available at https://github.com/ivanslapnicar/MANAA. Our results complement and extend the existing results for commutative fields.
Joint work with Nevena Jakovčević Stor and Thaniporn Chaysri.
Achieving scalability with a Hermitian preconditioner for a class of non-Hermitian matrices
Nicole Spillane
CNRS, Ecole Polytechnique, France - This email address is being protected from spambots. You need JavaScript enabled to view it.
This work considers the convergence of GMRES for non-singular problems. Preconditioning and weighted norms within GMRES are considered.
The main focus is on Hermitian preconditioning (even for non-Hermitian problems). It is proposed to choose a Hermitian preconditioner H and to apply GMRES in the inner product induced by H. If moreover, the problem matrix A is positive definite, then a new convergence bound is proved that depends only on how well H preconditions the Hermitian part of A, and on how non-Hermitian A is.
In particular, if a scalable preconditioner is known for the Hermitian part of A, then the proposed method is also scalable. This result is illustrated numerically.
Reference: Nicole Spillane. Hermitian Preconditioning for a class of Non-Hermitian Linear Systems. 2023. https://hal.science/hal-04028590v1/document
Probabilistic bounds on best rank-one approximation ratio
Josue Tonelli-Cueto
The University of Texas at San Antonio, United States - This email address is being protected from spambots. You need JavaScript enabled to view it.
In this poster, we provide new upper and lower bounds on the minimum possible ratio of the spectral and Frobenius norms of a (partially) symmetric tensor. In the particular case of general tensors our result recovers a known upper bound. For symmetric tensors our upper bound unveils that the ratio of norms has the same order of magnitude as the trivial lower bound $1/n^{\frac{d-1}{2}}$, when the order of a tensor $d$ is fixed and the dimension of the underlying vector space $n$ tends to infinity. However, when $n$ is fixed and $d$ tends to infinity, our lower bound is better than $1/n^{\frac{d-1}{2}}$.
Joint work with Khazhgali Kozhasov (Technische Universität Braunschweig, Germany).