Deflated block Krylov subspace methods for large scale eigenvalue problems

We discuss a class of deflated block Krylov subspace methods for solving large scale matrix eigenvalue problems. The efficiency of an Arnoldi-type method is examined in computing partial or closely clustered eigenvalues of large matrices. As an improvement, we also propose a refined variant of the Arnoldi-type method. Comparisons show that the refined variant can further improve the Arnoldi-type method and both methods exhibit very regular convergence behavior.     

A block inverse-free preconditioned Krylov subspace method for symmetric generalized eigenvalue problems

The inverse-free preconditioned Krylov subspace method of Golub and Ye [G.H. Golub, Q. Ye, An inverse free preconditioned Krylov subspace method for symmetric generalized eigenvalue problems, SIAM J. Sci. Comp. 24 (2002) 312-334] is an efficient algorithm for computing a few extreme eigenvalues of the symmetric generalized eigenvalue problem. In this paper, we first present an analysis of the preconditioning strategy based on incomplete factorizations. We then extend the method by developing a block generalization for computing multiple or severely clustered eigenvalues and develop a robust black-box implementation. Numerical examples are given to illustrate the analysis and the efficiency of the block algorithm.     

On computing PageRank via lumping the Google matrix

Computing Google’s PageRank via lumping the Google matrix was recently analyzed in [I.C.F. Ipsen, T.M. Selee, PageRank computation, with special attention to dangling nodes, SIAM J. Matrix Anal. Appl. 29 (2007) 1281–1296]. It was shown that all of the dangling nodes can be lumped into a single node and the PageRank could be obtained by applying the power method to the reduced matrix. Furthermore, the stochastic reduced matrix had the same nonzero eigenvalues as the full Google matrix and the power method applied to the reduced matrix had the same convergence rate as that of the power method applied to the full matrix. Therefore, a large amount of operations could be saved for computing the full PageRank vector.     

A modified SSOR iterative method for augmented systems

Golub, Wu and Yuan [G.H. Golub, X. Wu, J.Y. Yuan, SOR-like methods for augmented systems, BIT 41 (2001) 71–85] have presented the SOR-like algorithm to solve augmented systems. In this paper, we present the modified symmetric successive overrelaxation (MSSOR) method for solving augmented systems, which is based on Darvishi and Hessari’s work above. We derive its convergence under suitable restrictions on the iteration parameter, determine its optimal iteration parameter and the corresponding optimal convergence factor under certain conditions. Finally, we apply the MSSOR method to solve augmented systems.     

A method based on Rayleigh quotient gradient flow for extreme and interior eigenvalue problems

Recently, a continuous method has been proposed by Golub and Liao as an alternative way to solve the minimum and interior eigenvalue problems. According to their numerical results, their method seems promising. This article is an extension along this line. In this article, firstly, we convert an eigenvalue problem to an equivalent constrained optimization problem. Secondly, using the Karush-Kuhn-Tucker conditions of this equivalent optimization problem, we obtain a variant of the Rayleigh quotient gradient flow, which is formulated by a system of differential-algebraic equations. Thirdly, based on the Rayleigh quotient gradient flow, we give a practical numerical method for the minimum and interior eigenvalue problems. Finally, we also give some numerical experiments of our method, the Golub and Liao method, and EIGS (a Matlab implementation for computing eigenvalues using restarted Arnoldi’s method) for some typical eigenvalue problems. Our numerical experiments indicate that our method seems promising for most test problems.     

On generalized successive overrelaxation methods for augmented linear systems

For the augmented system of linear equations, Golub, Wu and Yuan recently studied an SOR-like method (BIT 41(2001)71–85). By further accelerating it with another parameter, in this paper we present a generalized SOR (GSOR) method for the augmented linear system. We prove its convergence under suitable restrictions on the iteration parameters, and determine its optimal iteration parameters and the corresponding optimal convergence factor. Theoretical analyses show that the GSOR method has faster asymptotic convergence rate than the SOR-like method. Also numerical results show that the GSOR method is more effective than the SOR-like method when they are applied to solve the augmented linear system. This GSOR method is further generalized to obtain a framework of the relaxed splitting iterative methods for solving both symmetric and nonsymmetric augmented linear systems by using the techniques of vector extrapolation, matrix relaxation and inexact iteration. Besides, we also demonstrate a complete version about the convergence theory of the SOR-like method. Subsidized by The Special Funds For Major State Basic Research Projects (No. G1999032803) and The National Natural Science Foundation (No. 10471146), P.R. China     

The computation of bounds for the norm of the error in the conjugate gradient algorithm

In this paper we consider computing estimates of the norm of the error in the conjugate gradient (CG) algorithm. Formulas were given in a paper by Golub and Meurant (1997). Here, we first prove that these expressions are indeed upper and lower bounds for the A-norm of the error. Moreover, starting from these formulas, we investigate the computation of the l ₂-norm of the error. Finally, we define an adaptive algorithm where the approximations of the extreme eigenvalues that are needed to obtain upper bounds are computed when running CG leading to an improvement of the upper bounds for the norm of the error. Numerical experiments show the effectiveness of this algorithm. This revised version was published online in August 2006 with corrections to the Cover Date.     

A subspace preconditioning algorithm for eigenvector/eigenvalue computation

We consider the problem of computing a modest number of the smallest eigenvalues along with orthogonal bases for the corresponding eigenspaces of a symmetric positive definite operatorA defined on a finite dimensional real Hilbert spaceV. In our applications, the dimension ofV is large and the cost of invertingA is prohibitive. In this paper, we shall develop an effective parallelizable technique for computing these eigenvalues and eigenvectors utilizing subspace iteration and preconditioning forA. Estimates will be provided which show that the preconditioned method converges linearly when used with a uniform preconditioner under the assumption that the approximating subspace is close enough to the span of desired eigenvectors.     

A global harmonic Arnoldi method for large non-Hermitian eigenproblems with an application to multiple eigenvalue problems

The global Arnoldi method can be used to compute exterior eigenpairs of a large non-Hermitian matrix A, but it does not work well for interior eigenvalue problems. Based on the global Arnoldi process that generates an F-orthonormal basis of a matrix Krylov subspace, we propose a global harmonic Arnoldi method for computing certain harmonic F-Ritz pairs that are used to approximate some interior eigenpairs. We propose computing the F-Rayleigh quotients of the large non-Hermitian matrix with respect to harmonic F-Ritz vectors and taking them as new approximate eigenvalues. They are better and more reliable than the harmonic F-Ritz values. The global harmonic Arnoldi method inherits convergence properties of the harmonic Arnoldi method applied to a larger matrix whose distinct eigenvalues are the same as those of the original given matrix. Some properties of the harmonic F-Ritz vectors are presented. As an application, assuming that A is diagonalizable, we show that the global harmonic Arnoldi method is able to solve multiple eigenvalue problems both in theory and in practice. To be practical, we develop an implicitly restarted global harmonic Arnoldi algorithm with certain harmonic F-shifts suggested. In particular, this algorithm can be adaptively used to solve multiple eigenvalue problems. Numerical experiments show that the algorithm is efficient for the eigenproblem and is reliable for quite ill-conditioned multiple eigenproblems.     

On nonsymmetric saddle point matrices that allow conjugate gradient iterations

Linear systems in saddle point form are usually highly indefinite,which often slows down iterative solvers such as Krylov subspace methods. It has been noted by several authors that negating the second block row of a symmetric indefinite saddle point matrix leads to a nonsymmetric matrix ${{\mathcal A}}Linear systems in saddle point form are usually highly indefinite,which often slows down iterative solvers such as Krylov subspace methods. It has been noted by several authors that negating the second block row of a symmetric indefinite saddle point matrix leads to a nonsymmetric matrix whose spectrum is entirely contained in the right half plane. In this paper we study conditions so that is diagonalizable with a real and positive spectrum. These conditions are based on necessary and sufficient conditions for positive definiteness of a certain bilinear form,with respect to which is symmetric. In case the latter conditions are satisfied, there exists a well defined conjugate gradient (CG) method for solving linear systems with . We give an efficient implementation of this method, discuss practical issues such as error bounds, and present numerical experiments. In memory of Gene Golub (1932–2007), our wonderful friend and colleague, who had a great interest in the conjugate gradient method and the numerical solution of saddle point problems. The work of J?rg Liesen was supported by the Emmy Noether-Program and the Heisenberg-Program of the Deutsche Forschungsgemeinschaft.     

A preconditioned and extrapolation-accelerated GMRES method for PageRank

In this paper, we propose and analyze GMRES-type methods for the PageRank computation. However, GMRES may converge very slowly or sometimes even diverge or break down when the damping factor is close to 1 and the dimension of the search subspace is low. We propose two strategies: preconditioning and vector extrapolation accelerating, to improve the convergence rate of the GMRES method. Theoretical analysis demonstrate the efficiency of the proposed strategies and numerical experiments show that the performance of the proposed methods is very much better than that of the traditional methods for PageRank problems.     

A block Arnoldi-type contour integral spectral projection method for solving generalized eigenvalue problems

For generalized eigenvalue problems, we consider computing all eigenvalues located in a certain region and their corresponding eigenvectors. Recently, contour integral spectral projection methods have been proposed for solving such problems. In this study, from the analysis of the relationship between the contour integral spectral projection and the Krylov subspace, we conclude that the Rayleigh–Ritz-type of the contour integral spectral projection method is mathematically equivalent to the Arnoldi method with the projected vectors obtained from the contour integration. By this Arnoldi-based interpretation, we then propose a block Arnoldi-type contour integral spectral projection method for solving the eigenvalue problem.     

Matrices, moments and quadrature II; How to compute the norm of the error in iterative methods

In this paper, we study the numerical computation of the errors in linear systems when using iterative methods. This is done by using methods to obtain bounds or approximations of quadratic formsu ^T A ⁻¹ u whereA is a symmetric positive definite matrix andu is a given vector. Numerical examples are given for the Gauss-Seidel algorithm. Moreover, we show that using a formula for theA-norm of the error from Dahlquist, Golub and Nash [1978] very good bounds of the error can be computed almost for free during the iterations of the conjugate gradient method leading to a reliable stopping criterion. The work of the first author was partially supported by NSF Grant CCR-950539.     

Implicitly restarted projection algorithm for solving optimization problems

An algorithm for solving the problem of minimizing a quadratic function subject to ellipsoidal constraints is introduced. This algorithm is based on the impHcitly restarted Lanczos method to construct a basis for the Krylov subspace in conjunction with a model trust region strategy to choose the step. The trial step is computed on the small dimensional subspace that lies inside the trust region.

One of the main advantages of this algorithm is the way that the Krylov subspace is terminated. We introduce a terminationcondition that allows the gradient to be decreased on that subspace.

A convergence theory for this algorithm is presented. It is shown that this algorithm is globally convergent and it shouldcope quite well with large scale minimization problems. This theory is sufficiently general that it holds for any algorithm that projects the problem on a lower dimensional subspace.     

Modified block preconditioners for the discretized time-harmonic Maxwell equations in mixed form

In this paper, building on the previous work by Greif and Schötzau [Preconditioners for the discretized time-harmonic Maxwell equations in mixed form, Numer. Linear Algebra Appl. 14 (2007) 281–297] and Benzi and Olshanskii [An augmented lagrangian-based approach to the Oseen problem, SIAM J. Sci. Comput. 28 (2006) 2095–2113], we present the improved preconditioning techniques for the iterative solution of the saddle point linear systems, which arise from the finite element discretization of the mixed formulation of the time-harmonic Maxwell equations. The modified block diagonal and triangular preconditioners considered are based on augmentation with using the symmetric nonsingular weighted matrix. We discuss the spectral properties of the preconditioned matrix in detail and generalize the results of the above-mentioned paper by Greif and Schötzau. Numerical experiments are given to demonstrate the efficiency of the presented preconditioners.     

Estimating the largest singular values of large sparse matrices via modified moments

We describe a procedure for determining a few of the largest singular values of a large sparse matrix. The method by Golub and Kent which uses the method of modified moments for estimating the eigenvalues of operators used in iterative methods for the solution of linear systems of equations is appropriately modified in order to generate a sequence of bidiagonal matrices whose singular values approximate those of the original sparse matrix. A simple Lanczos recursion is proposed for determining the corresponding left and right singular vectors. The potential asynchronous computation of the bidiagonal matrices using modified moments with the iterations of an adapted Chebyshev semi-iterative (CSI) method is an attractive feature for parallel computers. Comparisons in efficiency and accuracy with an appropriate Lanczos algorithm (with selective re-orthogonalization) are presented on large sparse (rectangular) matrices arising from applications such as information retrieval and seismic reflection tomography. This procedure is essentially motivated by the theory of moments and Gauss quadrature.This author's work was supported by the National Science Foundation under grants NSF CCR-8717492 and CCR-910000N (NCSA), the U.S. Department of Energy under grant DOE DE-FG02-85ER25001, and the Air Force Office of Scientific Research under grant AFOSR-90-0044 while at the University of Illinois at Urbana-Champaign Center for Supercomputing Research and Development.This author's work was supported by the U.S. Army Research Office under grant DAAL03-90-G-0105, and the National Science Foundation under grant NSF DCR-8412314.     

Numerical Taylor expansions of invariant manifolds in large dynamical systems

Summary. In this paper we develop a numerical method for computing higher order local approximations of invariant manifolds, such as stable, unstable or center manifolds near steady states of a dynamical system. The underlying system is assumed to be large in the sense that a large sparse Jacobian at the equilibrium occurs, for which only a linear (black box) solver and a low dimensional invariant subspace is available, but for which methods like the QR–Algorithm are considered to be too expensive. Our method is based on an analysis of the multilinear Sylvester equations for the higher derivatives which can be solved under certain nonresonance conditions. These conditions are weaker than the standard gap conditions on the spectrum which guarantee the existence of the invariant manifold. The final algorithm requires the solution of several large linear systems with a bordered Jacobian. To these systems we apply a block elimination method recently developed by Govaerts and Pryce [12, 14]. Received March 12, 1996 / Revised version reveiced August 8, 1997     

On the construction of restricted-denominator exponential W-methods

In this paper we consider the practical construction of exponential W-methods for the solution of large stiff nonlinear initial value problems, based on the restricted-denominator rational approach for the computation of the functions of matrices required. This approach is employed together with the Krylov subspace method based on the Arnoldi algorithm. Two integrators are constructed and tested on some classical stiff equations arising from the semidiscretization of parabolic problems.     

求解PageRank问题的重启GMRES修正的多分裂迭代法

PageRank算法已经成为网络搜索引擎的核心技术。针对PageRank问题导出的线性方程组,首先将Krylov子空间方法中的重启GMRES(generalized minimal residual)方法与多分裂迭代(multi-splitting iteration,MSI)方法相结合,提出了一种重启GMRES修正的多分裂迭代法;然后,给出了该算法的详细计算流程和收敛性分析;最后,通过数值实验验证了该算法的有效性。     

Error analysis of a pairwise summation algorithm to compute the sample variance

Summary We give an error analysis of an algorithm for computing the sample variance due to Chan, Golub, and LeVeque [The American Statistician 7 (1983), pp. 242–247]. It is shown that this algorithm is numerically stable. The algorithm computes the sample variance (and the sample mean) using just one pass through the sample data. It is amenable to pairwise summation and thus requires onlyO(logn) parallel steps.Research supported by the Air Force Office of Scientific Research under grant no. AFOSR-88-0161 and by the Office of Naval Research under grant no. N00024-85-C-6041