Analysis of a Classical Matrix Preconditioning Algorithm

We study a classical iterative algorithm for balancing matrices in the $L_\infty$ norm via a scaling transformation. This algorithm, which goes back to Osborne and Parlett \& Reinsch in the 1960s, is implemented as a standard preconditioner in many numerical linear algebra packages. Surprisingly, despite its widespread use over several decades, no bounds were known on its rate of convergence. In this paper we prove that, for any irreducible $n\times n$ (real or complex) input matrix~$A$, a natural variant of the algorithm converges in $O(n^3\log(n\rho/\varepsilon))$ elementary balancing operations, where $\rho$ measures the initial imbalance of~$A$ and $\varepsilon$ is the target imbalance of the output matrix. (The imbalance of~$A$ is $\max_i |\log(a_i^{\text{out}}/a_i^{\text{in}})|$, where $a_i^{\text{out}},a_i^{\text{in}}$ are the maximum entries in magnitude in the $i$th row and column respectively.) This bound is tight up to the $\log n$ factor. A balancing operation scales the $i$th row and column so that their maximum entries are equal, and requires $O(m/n)$ arithmetic operations on average, where $m$ is the number of non-zero elements in~$A$. Thus the running time of the iterative algorithm is $\tilde{O}(n^2m)$. This is the first time bound of any kind on any variant of the Osborne-Parlett-Reinsch algorithm. We also prove a conjecture of Chen that characterizes those matrices for which the limit of the balancing process is independent of the order in which balancing operations are performed.Comment: The previous version (1) (see also STOC'15) handled UB ("unique balance") input matrices. In this version (2) we extend the work to handle all input matrice

Low-rank approximation in the Frobenius norm by column and row subset selection

A CUR approximation of a matrix A is a particular type of low-rank approximation A approx CUR, where C and R consist of columns and rows of A, respectively. One way to obtain such an approximation is to apply column subset selection to A and AT . In this work, we describe a numerically robust and much faster variant of the column subset selection algorithm proposed by Deshpande and Rademacher, which guarantees an error close to the best approximation error in the Frobenius norm. For cross approximation, in which U is required to be the inverse of a submatrix of A described by the intersection of C and R, we obtain a new algorithm with an error bound that stays within a factor k + 1 of the best rank-k approximation error in the Frobenius norm. To the best of our knowledge, this is the first deterministic polynomial-time algorithm for which this factor is bounded by a polynomial in k. Our derivation and analysis of the algorithm is based on derandomizing a recent existence result by Zamarashkin and Osinsky. To illustrate the versatility of our new column subset selection algorithm, an extension to low multilinear rank approximations of tensors is provided as well. © 2020 Society for Industrial and Applied Mathematics

On randomized trace estimates for indefinite matrices with an application to determinants

Randomized trace estimation is a popular and well-studied technique that approximates the trace of a large-scale matrix B by computing the average of xTBx for many samples of a random vector X. Often, B is symmetric positive definite (SPD) but a number of applications give rise to indefinite B. Most notably, this is the case for log-determinant estimation, a task that features prominently in statistical learning, for instance in maximum likelihood estimation for Gaussian process regression. The analysis of randomized trace estimates, including tail bounds, has mostly focused on the SPD case. In this work, we derive new tail bounds for randomized trace estimates applied to indefinite B with Rademacher or Gaussian random vectors. These bounds significantly improve existing results for indefinite B, reducing the number of required samples by a factor n or even more, where n is the size of B. Even for an SPD matrix, our work improves an existing result by Roosta-Khorasani and Ascher (Found Comput Math, 15(5):1187–1212, 2015) for Rademacher vectors. This work also analyzes the combination of randomized trace estimates with the Lanczos method for approximating the trace of f(B). Particular attention is paid to the matrix logarithm, which is needed for log-determinant estimation. We improve and extend an existing result, to not only cover Rademacher but also Gaussian random vectors

Hierarchical adaptive low-rank format with applications to discretized partial differential equations

A novel framework for hierarchical low-rank matrices is proposed that combines an adaptive hierarchical partitioning of the matrix with low-rank approximation. One typical application is the approximation of discretized functions on rectangular domains; the flexibility of the format makes it possible to deal with functions that feature singularities in small, localized regions. To deal with time evolution and relocation of singularities, the partitioning can be dynamically adjusted based on features of the underlying data. Our format can be leveraged to efficiently solve linear systems with Kronecker product structure, as they arise from discretized partial differential equations (PDEs). For this purpose, these linear systems are rephrased as linear matrix equations and a recursive solver is derived from low-rank updates of such equations. We demonstrate the effectiveness of our framework for stationary and time-dependent, linear and nonlinear PDEs, including the Burgers' and Allen-Cahn equations

A block Newton method for nonlinear eigenvalue problems

We consider matrix eigenvalue problems that are nonlinear in the eigenvalue parameter. One of the most fundamental differences from the linear case is that distinct eigenvalues may have linearly dependent eigenvectors or even share the same eigenvector. This has been a severe hindrance in the development of general numerical schemes for computing several eigenvalues of a nonlinear eigenvalue problem, either simultaneously or subsequently. The purpose of this work is to show that the concept of invariant pairs offers a way of representing eigenvalues and eigenvectors that is insensitive to this phenomenon. To demonstrate the use of this concept in the development of numerical methods, we have developed a novel block Newton method for computing such invariant pairs. Algorithmic aspects of this method are considered and a few academic examples demonstrate its viability. © Springer-Verlag 2009

Low-rank tensor approximation for Chebyshev interpolation in parametric option pricing

Treating high dimensionality is one of the main challenges in the development of computational methods for solving problems arising in finance, where tasks such as pricing, calibration, and risk assessment need to be performed accurately and in real-time. Among the growing literature addressing this problem, Gass et al. [Finance Stoch., 22 (2018), pp. 701--731] propose a complexity reduction technique for parametric option pricing based on Chebyshev interpolation. As the number of parameters increases, however, this method is affected by the curse of dimensionality. In this article, we extend this approach to treat high-dimensional problems: Additionally, exploiting low-rank structures allows us to consider parameter spaces of high dimensions. The core of our method is to express the tensorized interpolation in the tensor train format and to develop an efficient way, based on tensor completion, to approximate the interpolation coefficients. We apply the new method to two model problems: American option pricing in the Heston model and European basket option pricing in the multidimensional Black--Scholes model. In these examples, we treat parameter spaces of dimensions up to 25. The numerical results confirm the low-rank structure of these problems and the effectiveness of our method compared to advanced techniques

Perturbation bounds for isotropic invariant subspaces of skew-Hamiltonian matrices

We investigate the behavior of isotropic invariant subspaces of skew-Hamiltonian matrices under structured perturbations. It is shown that finding a nearby subspace is equivalent to solving a certain quadratic matrix equation. This connection is used to derive meaningful error bounds and condition numbers that can be used to judge the quality of invariant subspaces computed by strongly backward stable eigensolvers

Compress-and-Restart Block Krylov Subspace Methods for Sylvester Matrix Equations

Block Krylov subspace methods (KSMs) comprise building blocks in many state-of-the-art solvers for large-scale matrix equations as they arise, for example, from the discretization of partial differential equations. While extended and rational block Krylov subspace methods provide a major reduction in iteration counts over polynomial block KSMs, they also require reliable solvers for the coefficient matrices, and these solvers are often iterative methods themselves. It is not hard to devise scenarios in which the available memory, and consequently the dimension of the Krylov subspace, is limited. In such scenarios for linear systems and eigenvalue problems, restarting is a well-explored technique for mitigating memory constraints. In this work, such restarting techniques are applied to polynomial KSMs for matrix equations with a compression step to control the growing rank of the residual. An error analysis is also performed, leading to heuristics for dynamically adjusting the basis size in each restart cycle. A panel of numerical experiments demonstrates the effectiveness of the new method with respect to extended block KSMs

A periodic Krylov-Schur algorithm for large matrix products

Stewart's recently introduced Krylov-Schur algorithm is a modification of the implicitly restarted Arnoldi algorithm which employs reordered Schur decompositions to perform restarts and deflations in a numerically reliable manner. This paper describes a variant of the Krylov-Schur algorithm suitable for addressing eigenvalue problems associated with products of large and sparse matrices. It performs restarts and deflations via reordered periodic Schur decompositions and, by taking the product structure into account, it is capable of achieving qualitatively better approximations to eigenvalues of small magnitude