Lars Grasedyck, Daniel Kressner, Christine Tobler
During the last years, low-rank tensor approximation has been established as a new tool in scientific computing to address large-scale linear and multilinear algebra problems, which would be intractable by classical techniques. This survey attempts to give a literature overview of current developments in this area, with an emphasis on function-related tensors.
Noah Amsel, Yves Baumann, Paul Beckman et al. · berkeley
This document presents a series of open questions arising in matrix computations, i.e., the numerical solution of linear algebra problems. It is a result of working groups at the workshop Linear Systems and Eigenvalue Problems, which was organized at the Simons Institute for the Theory of Computing program on Complexity and Linear Algebra in Fall 2025. The complexity and numerical solution of linear algebra problems is a crosscutting area between theoretical computer science and numerical analysis. The value of the particular problem formulations here is that they were produced via discussions between researchers from both groups. The open questions are organized in five categories: iterative solvers for linear systems, eigenvalue computation, low-rank approximation, randomized sketching, and other areas including tensors, quantum systems, and matrix functions.
Daniel Kressner, Stefano Massei, Leonardo Robol
Linear matrix equations, such as the Sylvester and Lyapunov equations, play an important role in various applications, including the stability analysis and dimensionality reduction of linear dynamical control systems and the solution of partial differential equations. In this work, we present and analyze a new algorithm, based on tensorized Krylov subspaces, for quickly updating the solution of such a matrix equation when its coefficients undergo low-rank changes. We demonstrate how our algorithm can be utilized to accelerate the Newton method for solving continuous-time algebraic Riccati equations. Our algorithm also forms the basis of a new divide-and-conquer approach for linear matrix equations with coefficients that feature hierarchical low-rank structure, such as HODLR, HSS, and banded matrices. Numerical experiments demonstrate the advantages of divide-and-conquer over existing approaches, in terms of computational time and memory consumption.
Daniel Kressner, Robert Luce
The computation of the matrix exponential is a ubiquitous operation in numerical mathematics, and for a general, unstructured $n\times n$ matrix it can be computed in $\mathcal{O}(n^3)$ operations. An interesting problem arises if the input matrix is a Toeplitz matrix, for example as the result of discretizing integral equations with a time invariant kernel. In this case it is not obvious how to take advantage of the Toeplitz structure, as the exponential of a Toeplitz matrix is, in general, not a Toeplitz matrix itself. The main contribution of this work are fast algorithms for the computation of the Toeplitz matrix exponential. The algorithms have provable quadratic complexity if the spectrum is real, or sectorial, or more generally, if the imaginary parts of the rightmost eigenvalues do not vary too much. They may be efficient even outside these spectral constraints. They are based on the scaling and squaring framework, and their analysis connects classical results from rational approximation theory to matrices of low displacement rank. As an example, the developed methods are applied to Merton's jump-diffusion model for option pricing.
Minhong Chen, Daniel Kressner
Recursive blocked algorithms have proven to be highly efficient at the numerical solution of the Sylvester matrix equation and its generalizations. In this work, we show that these algorithms extend in a seamless fashion to higher-dimensional variants of generalized Sylvester matrix equations, as they arise from the discretization of PDEs with separable coefficients or the approximation of certain models in macroeconomics. By combining recursions with a mechanism for merging dimensions, an efficient algorithm is derived that outperforms existing approaches based on Sylvester solvers.
Daniel Kressner, Emre Mengi, Ivica Nakic et al.
We consider the distance from a (square or rectangular) matrix pencil to the nearest matrix pencil in 2-norm that has a set of specified eigenvalues. We derive a singular value optimization characterization for this problem and illustrate its usefulness for two applications. First, the characterization yields a singular value formula for determining the nearest pencil whose eigenvalues lie in a specified region in the complex plane. For instance, this enables the numerical computation of the nearest stable descriptor system in control theory. Second, the characterization partially solves the problem posed in [Boutry et al. 2005] regarding the distance from a general rectangular pencil to the nearest pencil with a complete set of eigenvalues. The involved singular value optimization problems are solved by means of BFGS and Lipschitz-based global optimization algorithms.
Alice Cortinovis, Daniel Kressner, Stefano Massei
The problem of finding a $k \times k$ submatrix of maximum volume of a matrix $A$ is of interest in a variety of applications. For example, it yields a quasi-best low-rank approximation constructed from the rows and columns of $A$. We show that such a submatrix can always be chosen to be a principal submatrix if $A$ is symmetric semidefinite or diagonally dominant. Then we analyze the low-rank approximation error returned by a greedy method for volume maximization, cross approximation with complete pivoting. Our bound for general matrices extends an existing result for symmetric semidefinite matrices and yields new error estimates for diagonally dominant matrices. In particular, for doubly diagonally dominant matrices the error is shown to remain within a modest factor of the best approximation error. We also illustrate how the application of our results to cross approximation for functions leads to new and better convergence results.
Daniel Kressner
Bivariate matrix functions provide a unified framework for various tasks in numerical linear algebra, including the solution of linear matrix equations and the application of the Fréchet derivative. In this work, we propose a novel tensorized Krylov subspace method for approximating such bivariate matrix functions and analyze its convergence. While this method is already known for some instances, our analysis appears to result in new convergence estimates and insights for all but one instance, Sylvester matrix equations.
Daniel Kressner, André Uschmajew
Low-rank tensor approximation techniques attempt to mitigate the overwhelming complexity of linear algebra tasks arising from high-dimensional applications. In this work, we study the low-rank approximability of solutions to linear systems and eigenvalue problems on Hilbert spaces. Although this question is central to the success of all existing solvers based on low-rank tensor techniques, very few of the results available so far allow to draw meaningful conclusions for higher dimensions. In this work, we develop a constructive framework to study low-rank approximability. One major assumption is that the involved linear operator admits a low-rank representation with respect to the chosen tensor format, a property that is known to hold in a number of applications. Simple conditions, which are shown to hold for a fairly general problem class, guarantee that our derived low-rank truncation error estimates do not deteriorate as the dimensionality increases.
Jonas Ballani, Daniel Kressner, Michael Peters
In this paper, we introduce and analyze a new low-rank multilevel strategy for the solution of random diffusion problems. Using a standard stochastic collocation scheme, we first approximate the infinite dimensional random problem by a deterministic parameter-dependent problem on a high-dimensional parameter domain. Given a hierarchy of finite element discretizations for the spatial approximation, we make use of a multilevel framework in which we consider the differences of the solution on two consecutive finite element levels in the collocation points. We then address the approximation of these high-dimensional differences by adaptive low-rank tensor techniques. This allows to equilibrate the error on all levels by exploiting analytic and algebraic properties of the solution at the same time. We arrive at an explicit representation in a low-rank tensor format of the approximate solution on the entire parameter domain, which can be used for, e.g., the direct and cheap computation of statistics. Numerical results are provided in order to illustrate the approach.
Zvonimir Bujanović, Lars Karlsson, Daniel Kressner
The QZ algorithm for computing eigenvalues and eigenvectors of a matrix pencil $A - λB$ requires that the matrices first be reduced to Hessenberg-triangular (HT) form. The current method of choice for HT reduction relies entirely on Givens rotations regrouped and accumulated into small dense matrices which are subsequently applied using matrix multiplication routines. A non-vanishing fraction of the total flop count must nevertheless still be performed as sequences of overlapping Givens rotations alternately applied from the left and from the right. The many data dependencies associated with this computational pattern leads to inefficient use of the processor and poor scalability. In this paper, we therefore introduce a fundamentally different approach that relies entirely on (large) Householder reflectors partially accumulated into block reflectors, by using (compact) WY representations. Even though the new algorithm requires more floating point operations than the state of the art algorithm, extensive experiments on both real and synthetic data indicate that it is still competitive, even in a sequential setting. The new algorithm is conjectured to have better parallel scalability, an idea which is partially supported by early small-scale experiments using multi-threaded BLAS. The design and evaluation of a parallel formulation is future work.
Nicola Guglielmi, Mutti-Ur Rehman, Daniel Kressner
A novel method for approximating structured singular values (also known as mu-values) is proposed and investigated. These quantities constitute an important tool in the stability analysis of uncertain linear control systems as well as in structured eigenvalue perturbation theory. Our approach consists of an inner-outer iteration. In the outer iteration, a Newton method is used to adjust the perturbation level. The inner iteration solves a gradient system associated with an optimization problem on the manifold induced by the structure. Numerical results and comparison with the well-known Matlab function mussv, implemented in the Matlab Control Toolbox, illustrate the behavior of the method.
Matthias Bolten, Karsten Kahl, Daniel Kressner et al.
Markov chains that describe interacting subsystems suffer, on the one hand, from state space explosion but lead, on the other hand, to highly structured matrices. In this work, we propose a novel tensor-based algorithm to address such tensor structured Markov chains. Our algorithm combines a tensorized multigrid method with AMEn, an optimization-based low-rank tensor solver, for addressing coarse grid problems. Numerical experiments demonstrate that this combination overcomes the limitations incurred when using each of the two methods individually. As a consequence, Markov chain models of unprecedented size from a variety of applications can be addressed.
Fabio Matti, Martin Skovgaard Andersen, Tianshi Chen et al.
We propose a novel Krylov subspace method for estimating the finite impulse response (FIR) of a one-dimensional linear time-invariant systems. The method approximates the system's FIR using a kernel-based formulation combined with hyperparameter selection based on maximum likelihood estimation (MLE), which requires repeated evaluation of two terms: The data fit $\boldsymbol{y}^{\top} (λ\boldsymbol{I} + \boldsymbol{A})^{-1} \boldsymbol{y}$ and the model complexity $\log(\det (λ\boldsymbol{I} + \boldsymbol{A}))$, where $\boldsymbol{A}$ is a certain positive semidefinite matrix that admits fast matrix--vector products and $λ> 0$ is a regularization parameter. Instead of approximating these two quantities separately, we jointly approximate them using a single augmented Krylov subspace for $\boldsymbol{A}$. One major benefit of augmentation is that we obtain accelerated convergence when approximating the data fit quadratic form, through implicit preconditioning. Thanks to the shift invariance of Krylov subspaces, the extracted approximations can be used to evaluate the MLE objective for many values of $λ$ at little additional cost. We derive error bounds for the approximations, reflecting the benefits of augmentation demonstrated through multiple numerical experiments.
Daniel Kressner, Jonas Latz, Stefano Massei et al.
Many techniques for data science and uncertainty quantification demand efficient tools to handle Gaussian random fields, which are defined in terms of their mean functions and covariance operators. Recently, parameterized Gaussian random fields have gained increased attention, due to their higher degree of flexibility. However, especially if the random field is parameterized through its covariance operator, classical random field discretization techniques fail or become inefficient. In this work we introduce and analyze a new and certified algorithm for the low-rank approximation of a parameterized family of covariance operators which represents an extension of the adaptive cross approximation method for symmetric positive definite matrices. The algorithm relies on an affine linear expansion of the covariance operator with respect to the parameters, which needs to be computed in a preprocessing step using, e.g., the empirical interpolation method. We discuss and test our new approach for isotropic covariance kernels, such as Matérn kernels. The numerical results demonstrate the advantages of our approach in terms of computational time and confirm that the proposed algorithm provides the basis of a fast sampling procedure for parameter dependent Gaussian random fields.
Daniel Kressner, Ana Susnjara
The efficient and accurate QR decomposition for matrices with hierarchical low-rank structures, such as HODLR and hierarchical matrices, has been challenging. Existing structure-exploiting algorithms are prone to numerical instability as they proceed indirectly, via Cholesky decompositions or a block Gram-Schmidt procedure. For a highly ill-conditioned matrix, such approaches either break down in finite-precision arithmetic or result in significant loss of orthogonality. Although these issues can sometimes be addressed by regularization and iterative refinement, it would be more desirable to have an algorithm that avoids these detours and is numerically robust to ill-conditioning. In this work, we propose such an algorithm for HODLR matrices. It achieves accuracy by utilizing Householder reflectors. It achieves efficiency by utilizing fast operations in the HODLR format in combination with compact WY representations and the recursive QR decomposition by Elmroth and Gustavson. Numerical experiments demonstrate that our newly proposed algorithm is robust to ill-conditioning and capable of achieving numerical orthogonality down to the level of roundoff error.
Bernhard Beckermann, Daniel Kressner, Marcel Schweitzer
We consider the task of updating a matrix function $f(A)$ when the matrix $A\in{\mathbb C}^{n \times n}$ is subject to a low-rank modification. In other words, we aim at approximating $f(A+D)-f(A)$ for a matrix $D$ of rank $k \ll n$. The approach proposed in this paper attains efficiency by projecting onto tensorized Krylov subspaces produced by matrix-vector multiplications with $A$ and $A^*$. We prove the approximations obtained from $m$ steps of the proposed methods are exact if $f$ is a polynomial of degree at most $m$ and use this as a basis for proving a variety of convergence results, in particular for the matrix exponential and for Markov functions. We illustrate the performance of our method by considering various examples from network analysis, where our approach can be used to cheaply update centrality and communicability measures.
Daniel Kressner, Robert Luce, Francesco Statti
We study the problem of computing the matrix exponential of a block triangular matrix in a peculiar way: Block column by block column, from left to right. The need for such an evaluation scheme arises naturally in the context of option pricing in polynomial diffusion models. In this setting a discretization process produces a sequence of nested block triangular matrices, and their exponentials are to be computed at each stage, until a dynamically evaluated criterion allows to stop. Our algorithm is based on scaling and squaring. By carefully reusing certain intermediate quantities from one step to the next, we can efficiently compute such a sequence of matrix exponentials.
Erna Begovic, Daniel Kressner
This paper is concerned with low multilinear rank approximations to antisymmetric tensors, that is, multivariate arrays for which the entries change sign when permuting pairs of indices. We show which ranks can be attained by an antisymmetric tensor and discuss the adaption of existing approximation algorithms to preserve antisymmetry, most notably a Jacobi algorithm. Particular attention is paid to the important special case when choosing the rank equal to the order of the tensor. It is shown that this case can be addressed with an unstructured rank-$1$ approximation. This allows for the straightforward application of the higher-order power method, for which we discuss effective initialization strategies.
Dorina Thanou, Xiaowen Dong, Daniel Kressner et al.
Effective information analysis generally boils down to properly identifying the structure or geometry of the data, which is often represented by a graph. In some applications, this structure may be partly determined by design constraints or pre-determined sensing arrangements, like in road transportation networks for example. In general though, the data structure is not readily available and becomes pretty difficult to define. In particular, the global smoothness assumptions, that most of the existing works adopt, are often too general and unable to properly capture localized properties of data. In this paper, we go beyond this classical data model and rather propose to represent information as a sparse combination of localized functions that live on a data structure represented by a graph. Based on this model, we focus on the problem of inferring the connectivity that best explains the data samples at different vertices of a graph that is a priori unknown. We concentrate on the case where the observed data is actually the sum of heat diffusion processes, which is a quite common model for data on networks or other irregular structures. We cast a new graph learning problem and solve it with an efficient nonconvex optimization algorithm. Experiments on both synthetic and real world data finally illustrate the benefits of the proposed graph learning framework and confirm that the data structure can be efficiently learned from data observations only. We believe that our algorithm will help solving key questions in diverse application domains such as social and biological network analysis where it is crucial to unveil proper geometry for data understanding and inference.
Daniel Kressner, Ana Susnjara
We consider the approximate computation of spectral projectors for symmetric banded matrices. While this problem has received considerable attention, especially in the context of linear scaling electronic structure methods, the presence of small relative spectral gaps challenges existing methods based on approximate sparsity. In this work, we show how a data-sparse approximation based on hierarchical matrices can be used to overcome this problem. We prove a priori bounds on the approximation error and propose a fast algo- rithm based on the QDWH algorithm, along the works by Nakatsukasa et al. Numerical experiments demonstrate that the performance of our algorithm is robust with respect to the spectral gap. A preliminary Matlab implementation becomes faster than eig already for matrix sizes of a few thousand.
Ana Susnjara, Nathanael Perraudin, Daniel Kressner et al.
Signal-processing on graphs has developed into a very active field of research during the last decade. In particular, the number of applications using frames constructed from graphs, like wavelets on graphs, has substantially increased. To attain scalability for large graphs, fast graph-signal filtering techniques are needed. In this contribution, we propose an accelerated algorithm based on the Lanczos method that adapts to the Laplacian spectrum without explicitly computing it. The result is an accurate, robust, scalable and efficient algorithm. Compared to existing methods based on Chebyshev polynomials, our solution achieves higher accuracy without increasing the overall complexity significantly. Furthermore, it is particularly well suited for graphs with large spectral gaps.
Daniel Kressner, Michael Steinlechner, Bart Vandereycken
The numerical solution of partial differential equations on high-dimensional domains gives rise to computationally challenging linear systems. When using standard discretization techniques, the size of the linear system grows exponentially with the number of dimensions, making the use of classic iterative solvers infeasible. During the last few years, low-rank tensor approaches have been developed that allow to mitigate this curse of dimensionality by exploiting the underlying structure of the linear operator. In this work, we focus on tensors represented in the Tucker and tensor train formats. We propose two preconditioned gradient methods on the corresponding low-rank tensor manifolds: A Riemannian version of the preconditioned Richardson method as well as an approximate Newton scheme based on the Riemannian Hessian. For the latter, considerable attention is given to the efficient solution of the resulting Newton equation. In numerical experiments, we compare the efficiency of our Riemannian algorithms with other established tensor-based approaches such as a truncated preconditioned Richardson method and the alternating linear scheme. The results show that our approximate Riemannian Newton scheme is significantly faster in cases when the application of the linear operator is expensive.
Petar Sirković, Daniel Kressner
This work is concerned with approximating the smallest eigenvalue of a parameter-dependent Hermitian matrix $A(μ)$ for many parameter values $μ\in \mathbb{R}^P$. The design of reliable and efficient algorithms for addressing this task is of importance in a variety of applications. Most notably, it plays a crucial role in estimating the error of reduced basis methods for parametrized partial differential equations. The current state-of-the-art approach, the so called Successive Constraint Method (SCM), addresses affine linear parameter dependencies by combining sampled Rayleigh quotients with linear programming techniques. In this work, we propose a subspace approach that additionally incorporates the sampled eigenvectors of $A(μ)$ and implicitly exploits their smoothness properties. Like SCM, our approach results in rigorous lower and upper bounds for the smallest eigenvalues on $D$. Theoretical and experimental evidence is given to demonstrate that our approach represents a significant improvement over SCM in the sense that the bounds are often much tighter, at negligible additional cost.