Block Krylov Subspace Methods Research Articles

Krylov subspace methods to accelerate kernel machines on graphs

In classical frameworks as the Euclidean space, positive definite kernels as well as their analytic properties are explicitly available and can be incorporated directly in kernel-based learning algorithms. This is different if the underlying domain is a discrete irregular graph. In this case, respective kernels have to be computed in a preliminary step in order to apply them inside a kernel machine. Typically, such a kernel is given as a matrix function of the graph Laplacian. Its direct calculation leads to a high computational burden if the size of the graph is very large. In this work, we investigate five different block Krylov subspace methods to obtain cheaper iterative approximations of these kernels. We will investigate convergence properties of these Krylov subspace methods and study to what extent these methods are able to preserve the symmetry and positive definiteness of the original kernels they are approximating. We will further discuss the computational complexity and the memory requirements of these methods, as well as possible implications for the kernel predictors in machine learning.

Adaptively restarted block Krylov subspace methods with low-synchronization skeletons

With the recent realization of exascale performance by Oak Ridge National Laboratory’s Frontier supercomputer, reducing communication in kernels like QR factorization has become even more imperative. Low-synchronization Gram-Schmidt methods, first introduced in Świrydowicz et al. (Numer. Lin. Alg. Appl. 28(2):e2343, 2020), have been shown to improve the scalability of the Arnoldi method in high-performance distributed computing. Block versions of low-synchronization Gram-Schmidt show further potential for speeding up algorithms, as column-batching allows for maximizing cache usage with matrix-matrix operations. In this work, low-synchronization block Gram-Schmidt variants from Carson et al. (Linear Algebra Appl. 638:150–195, 2022) are transformed into block Arnoldi variants for use in block full orthogonalization methods (BFOM) and block generalized minimal residual methods (BGMRES). An adaptive restarting heuristic is developed to handle instabilities that arise with the increasing condition number of the Krylov basis. The performance, accuracy, and stability of these methods are assessed via a flexible benchmarking tool written in MATLAB. The modularity of the tool additionally permits generalized block inner products, like the global inner product.

Cross-Interactive Residual Smoothing for Global and Block Lanczos-Type Solvers for Linear Systems with Multiple Right-Hand Sides

Global and block Krylov subspace methods are efficient iterative solvers for large sparse linear systems with multiple right-hand sides. However, global or block Lanczos-type solvers often exhibit large oscillations in the residual norms and may have a large residual gap relating to the loss of attainable accuracy of the approximations. Conventional residual smoothing schemes suppress these oscillations but cannot improve the attainable accuracy, whereas a recent residual smoothing scheme enables the improvement of the attainable accuracy for single right-hand side Lanczos-type solvers. The underlying concept of this scheme is that the primary and smoothed sequences of the approximations and residuals influence one another, thereby avoiding the severe propagation of rounding errors. In the present study, we extend this cross-interactive residual smoothing to the case of solving linear systems with multiple right-hand sides. The resulting smoothed methods can reduce the residual gap with a low additional cost compared to their original counterparts. We demonstrate the effectiveness of the proposed approach through rounding error analysis and numerical experiments.

Randomized block Krylov methods for approximating extreme eigenvalues

Randomized block Krylov subspace methods form a powerful class of algorithms for computing the extreme eigenvalues of a symmetric matrix or the extreme singular values of a general matrix. The purpose of this paper is to develop new theoretical bounds on the performance of randomized block Krylov subspace methods for these problems. For matrices with polynomial spectral decay, the randomized block Krylov method can obtain an accurate spectral norm estimate using only a constant number of steps (that depends on the decay rate and the accuracy). Furthermore, the analysis reveals that the behavior of the algorithm depends in a delicate way on the block size. Numerical evidence confirms these predictions.

Global symplectic Lanczos method with application to matrix exponential approximation

It is well-known that the symplectic Lanczos method is an efficient tool for computing a few eigenvalues of large and sparse Hamiltonian matrices. A variety of block Krylov subspace methods were introduced by Lopez and Simoncini to compute an approximation of $exp(M)V$ for a given large square Hamiltonian matrix $M$ and a tall and skinny matrix $V$ that preserves the geometric property of $V$. For the same purpose, in this paper, we have proposed a new method based on a global version of the symplectic Lanczos algorithm, called the global $J$-Lanczos method ($GJ$-Lanczos). To the best of our knowledge, this is probably the first adaptation of the symplectic Lanczos method in the global case. Numerical examples are given to illustrate the effectiveness of the proposed approach.

Randomized block Krylov subspace methods for trace and log-determinant estimators

We present randomized algorithms based on block Krylov subspace methods for estimating the trace and log-determinant of Hermitian positive semi-definite matrices. Using the properties of Chebyshev polynomials and Gaussian random matrix, we provide the error analysis of the proposed estimators and obtain the expectation and concentration error bounds. These bounds improve the corresponding ones given in the literature. Numerical experiments are presented to illustrate the performance of the algorithms and to test the error bounds.

An efficient approach for solving saddle point problems using block structure

This paper focuses on saddle point problems with a 2-by-2 block coefficient matrix. When the number of columns in the upper-right block and the number of rows in the lower-left block of the coefficient matrix is large, the convergence behavior of Krylov subspace methods for the saddle point problems tends to be poor even if the upper-left block is a well-conditioned matrix. In this paper, an efficient approach for solving the saddle point problems using block structure of the problems is proposed. The most time-consuming part of our proposed approach is the solution of a linear system with multiple right-hand sides. To solve the linear system with multiple right-hand sides efficiently, we propose to apply Block Krylov subspace methods to this linear system. Numerical experiments show that the proposed approach with Block Krylov subspace methods can solve the saddle point problems more efficiently than the conventional approach in terms of the number of iterations and the computation time.

Compress‐and‐restart block Krylov subspace methods for Sylvester matrix equations

SummaryBlock 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.

The tensor t‐function: A definition for functions of third‐order tensors

SummaryA definition for functions of multidimensional arrays is presented. The definition is valid for third‐order tensors in the tensor t‐product formalism, which regards third‐order tensors as block circulant matrices. The tensor function definition is shown to have similar properties as standard matrix function definitions in fundamental scenarios. To demonstrate the definition's potential in applications, the notion of network communicability is generalized to third‐order tensors and computed for a small‐scale example via block Krylov subspace methods for matrix functions. A complexity analysis for these methods in the context of tensors is also provided.

Development of the Block BiCGGR2 method for linear systems with multiple right-hand sides

Block Krylov subspace methods are efficient iterative methods for solving linearsystems with multiple right-hand sides. The Block BiCGGR method is one of the Block Krylov subspace methods, and can generate highly accurate approximatesolutions. However, when the number of right-hand sides is large, this method becomes numerically unstable. In this paper, the Block BiCGGR2 method is proposed. This method is a numerically stable and highly accurate method than the Block BiCGGR method and the Block BiCGSTAB method. Numerical results illustrate the robustness and the accuracy of the approximate solutions of the Block BiCGGR2 method.

Preconditioned global Krylov subspace methods for solving saddle point problems with multiple right-hand sides

In the present paper, we propose a preconditioned global approach as a new strategy to solve linear systems with several right-hand sides coming from saddle point problems. The preconditioner is obtained by replacing a (2,2)-block in the original saddle-point matrix A by another well-chosen block. We apply the global GMRES method to solve this new problem with several right-hand sides and give some convergence results. Moreover, we analyze the eigenvalue distribution and the eigenvectors of the proposed preconditioner when the first block is positive definite. We also compare different preconditioned global Krylov subspace algorithms (CG, MINRES, FGMRES, GMRES) with preconditioned block (CG, GMRES) algorithms. Numerical results show that our preconditioned global GMRES method is competitive with other preconditioned global Krylov subspace and preconditioned block Krylov subspace methods for solving saddle point problems with several right-hand sides.

Investigation of Block Krylov Subspace Methods with Dummy Right-Hand Sides for Solving Linear Systems Obtained By Extended Element-Free Galerkin Method

Investigation of Block Krylov Subspace Methods with Dummy Right-Hand Sides for Solving Linear Systems Obtained By Extended Element-Free Galerkin Method

Weighted and flexible versions of block CMRH method for solving nonsymmetric linear systems with multiple right-hand sides

Block Krylov subspace methods are the most popular algorithms for solving large nonsymmetric linear systems with multiple right-hand sides. One of them is the block CMRH method. This method generates a (non orthogonal) basis of the Krylov subspace through the block Hessenberg process. To accelerate the convergence of the block CMRH method, we will introduce two new methods. First, we present the block CMRH method with weighting strategy. In this method, the block CMRH method uses a different product at each restart. Second, we introduce a flexible version of the block CMRH algorithm that allows varying preconditioning at every step of the algorithm. Numerical experiments illustrate the benefits of the presented methods.

A spectrally preconditioned and initially deflated variant of the restarted block GMRES method for solving multiple right-hand sides linear systems

The solution of large linear systems with multiple right-hand sides given simultaneously is required in many large-scale scientific and engineering applications modelled by either partial differential or boundary integral equations. Block Krylov subspace methods are attractive to use for this problem class as they can overcome the memory bottleneck of direct methods and they perform block matrix-vector products, achieving high computational efficiency on modern cache-based computer architectures. In this paper we introduce variants of the block GMRES method that combine initial deflation and eigenvalue recycling strategies to remedy some typical convergence problems of block Krylov solvers. The new class of block iterative solvers has the ability to handle the approximate linear dependence of the block of right-hand sides and exploits approximate invariant subspaces recycled over the iterations to mitigate the bad effects that small eigenvalues can have on the convergence, by adapting an existing preconditioner. We illustrate the numerical behavior of the spectrally updated and initially deflated block GMRES method on a set of linear systems arising from the discretization of the Dirac equation and of boundary integral equations in electromagnetics scattering.

Krylov eigenvalue strategy using the FEAST algorithm with inexact system solves

SummaryThe FEAST eigenvalue algorithm is a subspace iteration algorithm that uses contour integration to obtain the eigenvectors of a matrix for the eigenvalues that are located in any user‐defined region in the complex plane. By computing small numbers of eigenvalues in specific regions of the complex plane, FEAST is able to naturally parallelize the solution of eigenvalue problems by solving for multiple eigenpairs simultaneously. The traditional FEAST algorithm is implemented by directly solving collections of shifted linear systems of equations; in this paper, we describe a variation of the FEAST algorithm that uses iterative Krylov subspace algorithms for solving the shifted linear systems inexactly. We show that this iterative FEAST algorithm (which we call IFEAST) is mathematically equivalent to a block Krylov subspace method for solving eigenvalue problems. By using Krylov subspaces indirectly through solving shifted linear systems, rather than directly using them in projecting the eigenvalue problem, it becomes possible to use IFEAST to solve eigenvalue problems using very large dimension Krylov subspaces without ever having to store a basis for those subspaces. IFEAST thus combines the flexibility and power of Krylov methods, requiring only matrix–vector multiplication for solving eigenvalue problems, with the natural parallelism of the traditional FEAST algorithm. We discuss the relationship between IFEAST and more traditional Krylov methods and provide numerical examples illustrating its behavior.

Block Krylov subspace methods for functions of matrices

A variety of block Krylov subspace methods have been successfully developed for linear systems and matrix equations. The application of block Krylov methods to compute matrix functions is, however, less established, despite the growing prevalence of matrix functions in scientific computing. Of particular importance is the evaluation of a matrix function on not just one but multiple vectors. The main contribution of this paper is a class of efficient block Krylov subspace methods tailored precisely to this task. With the full orthogonalization method (FOM) for linear systems forming the backbone of our theory, the resulting methods are referred to as B(FOM)2: block FOM for functions of matrices. Many other important results are obtained in the process of developing these new methods. Matrix-valued inner products are used to construct a general framework for block Krylov subspaces that encompasses already established results in the literature. Convergence bounds for B(FOM)2 are proven for Stieltjes functions applied to a class of matrices which are self-adjoint and positive definite with respect to the matrix-valued inner product. A detailed algorithm for B(FOM)2 with restarts is developed, whose efficiency is based on a recursive expression for the error, which is also used to update the solution. Numerical experiments demonstrate the power and versatility of this new class of methods for a variety of matrix-valued inner products, functions, and matrices.

A block Krylov subspace implementation of the time-parallel Paraexp method and its extension for nonlinear partial differential equations

A parallel time integration method for nonlinear partial differential equations is proposed. It is based on a new implementation of the Paraexp method for linear partial differential equations (PDEs) employing a block Krylov subspace method. For nonlinear PDEs the algorithm is based on our Paraexp implementation within a waveform relaxation. The initial value problem is solved iteratively on a complete time interval. Nonlinear terms are treated as source terms, provided by the solution from the previous iteration. At each iteration, the problem is decoupled into independent subproblems by the principle of superposition. The decoupled subproblems are solved fast by exponential integration, based on a block Krylov method. The new time integration is demonstrated for the one-dimensional advection–diffusion equation and the viscous Burgers equation. Numerical experiments confirm excellent parallel scaling for the linear advection–diffusion problem, and good scaling in case the nonlinear Burgers equation is simulated.

Block Krylov subspace methods for approximating the linear combination of [formula omitted]-functions arising in exponential integrators

In this paper we present numerical methods for computing the matrix functions arising in exponential integrators. The matrix functions are the linear combination of the form ∑j=0pφj(A)bj, where φj is related to exponential function. By reducing the corresponding systems of linear ODEs with polynomial inhomogeneity based on block Krylov subspace methods, two effective exponentials of low dimension matrix are derived to approximate the objective functions. Some useful a posteriori error estimates are established. The algorithms can be combined with time-stepping strategy to satisfy the accuracy requirements under a moderate dimension of the block Krylov subspace. Several numerical experiments are presented to demonstrate the efficiency of the methods.

Stochastic estimation of level density in nuclear shell-model calculations

An estimation method of the nuclear level density stochastically based on nuclear shell-model calculations is introduced. In order to count the number of the eigen-values of the shell-model Hamiltonian matrix, we perform the contour integral of the matrix element of a resolvent. The shifted block Krylov subspace method enables us its efficient computation. Utilizing this method, the contamination of center-of-mass motion is clearly removed.

On short recurrence Krylov type methods for linear systems with many right-hand sides

Block and global Krylov subspace methods have been proposed as methods adapted to the situation where one iteratively solves systems with the same matrix and several right hand sides. These methods are advantageous, since they allow to cast the major part of the arithmetic in terms of matrix-block vector products, and since, in the block case, they take their iterates from a potentially richer subspace. In this paper we consider the most established Krylov subspace methods which rely on short recurrences, i.e. BiCG, QMR and BiCGStab. We propose modifications of their block variants which increase numerical stability, thus at least partly curing a problem previously observed by several authors. Moreover, we develop modifications of the “global” variants which almost halve the number of matrix–vector multiplications. We present a discussion as well as numerical evidence which both indicate that the additional work present in the block methods can be substantial, and that the new “economic” versions of the “global” BiCG and QMR method can be considered as good alternatives to the BiCGStab variants.