The restarted shift‐and‐invert Krylov method for matrix functions

In this paper, the numerical evaluation of matrix functions expressed in partial fraction form is addressed. The shift‐and‐invert Krylov method is analyzed, with special attention to error estimates. Such estimates give insights into the selection of the shift parameter and lead to a simple and effective restart procedure. Applications to the class of Mittag–Leffler functions are presented. Copyright © 2012 John Wiley & Sons, Ltd.     

Krylov eigenvalue strategy using the FEAST algorithm with inexact system solves

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

Krylov subspace methods for eigenvalues with special properties and their analysis for normal matrices

In this paper we propose a general approach by which eigenvalues with a special property of a given matrix A can be obtained. In this approach we first determine a scalar function ψ: C → C whose modulus is maximized by the eigenvalues that have the special property. Next, we compute the generalized power iterations uinj + 1 = ψ(A)u_j, j = 0, 1,…, where u₀ is an arbitrary initial vector. Finally, we apply known Krylov subspace methods, such as the Arnoldi and Lanczos methods, to the vector u_n for some sufficiently large n. We can also apply the simultaneous iteration method to the subspace span{x⁽ⁿ⁾₁,…,x⁽ⁿ⁾_k} with some sufficiently large n, where x^(j+1)_m = ψ(A)x^(j)_m, j = 0, 1,…, m = 1,…, k. In all cases the resulting Ritz pairs are approximations to the eigenpairs of A with the special property. We provide a rather thorough convergence analysis of the approach involving all three methods as n → ∞ for the case in which A is a normal matrix. We also discuss the connections and similarities of our approach with the existing methods and approaches in the literature.     

Convergence conditions for a restarted GMRES method augmented with eigenspaces

We consider the GMRES(m,k) method for the solution of linear systems Ax=b, i.e. the restarted GMRES with restart m where to the standard Krylov subspace of dimension m the other subspace of dimension k is added, resulting in an augmented Krylov subspace. This additional subspace approximates usually an A‐invariant subspace. The eigenspaces associated with the eigenvalues closest to zero are commonly used, as those are thought to hinder convergence the most. The behaviour of residual bounds is described for various situations which can arise during the GMRES(m,k) process. The obtained estimates for the norm of the residual vector suggest sufficient conditions for convergence of GMRES(m,k) and illustrate that these augmentation techniques can remove stagnation of GMRES(m) in many cases. All estimates are independent of the choice of an initial approximation. Conclusions and remarks assessing numerically the quality of proposed bounds conclude the paper. Copyright © 2004 John Wiley & Sons, Ltd.     

Globally and superlinearly convergent inexact Newton–Krylov algorithms for solving nonsmooth equations

This paper presents some variants of the inexact Newton method for solving systems of nonlinear equations defined by locally Lipschitzian functions. These methods use variants of Newton's iteration in association with Krylov subspace methods for solving the Jacobian linear systems. Global convergence of the proposed algorithms is established under a nonmonotonic backtracking strategy. The local convergence based on the assumptions of semismoothness and BD‐regularity at the solution is characterized, and the way to choose an inexact forcing sequence that preserves the rapid convergence of the proposed methods is also indicated. Numerical examples are given to show the practical viability of these approaches. Copyright © 2009 John Wiley & Sons, Ltd.     

Inexact Perturbed Newton Methods and Applications to a Class of Krylov Solvers

Inexact Newton methods are variant of the Newton method in which each step satisfies only approximately the linear system (Ref. 1). The local convergence theory given by the authors of Ref. 1 and most of the results based on it consider the error terms as being provided only by the fact that the linear systems are not solved exactly. The few existing results for the general case (when some perturbed linear systems are considered, which in turn are not solved exactly) do not offer explicit formulas in terms of the perturbations and residuals. We extend this local convergence theory to the general case, characterizing the rate of convergence in terms of the perturbations and residuals.The Newton iterations are then analyzed when, at each step, an approximate solution of the linear system is determined by the following Krylov solvers based on backward error minimization properties: GMRES, GMBACK, MINPERT. We obtain results concerning the following topics: monotone properties of the errors in these Newton–Krylov iterates when the initial guess is taken 0 in the Krylov algorithms; control of the convergence orders of the Newton–Krylov iterations by the magnitude of the backward errors of the approximate steps; similarities of the asymptotical behavior of GMRES and MINPERT when used in a converging Newton method. At the end of the paper, the theoretical results are verified on some numerical examples.     

A COCR method for solving complex symmetric linear systems

The Conjugate Orthogonal Conjugate Gradient (COCG) method has been recognized as an attractive Lanczos-type Krylov subspace method for solving complex symmetric linear systems; however, it sometimes shows irregular convergence behavior in practical applications. In the present paper, we propose a Conjugate A$A$-Orthogonal Conjugate Residual (COCR) method, which can be regarded as an extension of the Conjugate Residual (CR) method. Numerical examples show that COCR often gives smoother convergence behavior than COCG.     

On preconditioned Krylov subspace methods for discrete convection–diffusion problems

We study nonstationary iterative methods for solving preconditioned systems arising from discretizations of the convection–diffusion equation. The preconditioners arise from Gauss–Seidel methods applied to the original system. It is shown that the performance of the iterative solvers is affected by the relationship of the ordering of the underlying grid and the direction of the fow associated with the differential operator. Specifically, only those orderings that follow the fow give fast iterative solvers. © 1997 John Wiley & Sons, Inc. Numer Methods Partial Differential Eq 13 :321–330     

A new investigation of the extended Krylov subspace method for matrix function evaluations

For large square matrices A and functions f, the numerical approximation of the action of f(A) to a vector v has received considerable attention in the last two decades. In this paper we investigate theextended Krylov subspace method, a technique that was recently proposed to approximate f(A)v for A symmetric. We provide a new theoretical analysis of the method, which improves the original result for A symmetric, and gives a new estimate for A nonsymmetric. Numerical experiments confirm that the new error estimates correctly capture the linear asymptotic convergence rate of the approximation. By using recent algorithmic improvements, we also show that the method is computationally competitive with respect to other enhancement techniques. Copyright © 2009 John Wiley & Sons, Ltd.     

A parallel balance scheme for banded linear systems

A parallel algorithm is proposed for the solution of narrow banded non‐symmetric linear systems. The linear system is partitioned into blocks of rows with a small number of unknowns common to multiple blocks. Our technique yields a reduced system defined only on these common unknowns which can then be solved by a direct or iterative method. A projection based extension to this approach is also proposed for computing the reduced system implicitly, which gives rise to an inner–outer iteration method. In addition, the product of a vector with the reduced system matrix can be computed efficiently on a multiprocessor by concurrent projections onto subspaces of block rows. Scalable implementations of the algorithm can be devized for hierarchical parallel architectures by exploiting the two‐level parallelism inherent in the method. Our experiments indicate that the proposed algorithm is a robust and competitive alternative to existing methods, particularly for difficult problems with strong indefinite symmetric part. Copyright © 2001 John Wiley & Sons, Ltd.     

A generalization of parameterized inexact Uzawa method for generalized saddle point problems

Recently, a class of parameterized inexact Uzawa methods has been proposed for generalized saddle point problems by Bai and Wang [Z.-Z. Bai, Z.-Q. Wang, On parameterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl. 428 (2008) 2900–2932], and a generalization of the inexact parameterized Uzawa method has been studied for augmented linear systems by Chen and Jiang [F. Chen, Y.-L. Jiang, A generalization of the inexact parameterized Uzawa methods for saddle point problems, Appl. Math. Comput. (2008)]. This paper is concerned about a generalization of the parameterized inexact Uzawa method for solving the generalized saddle point problems with nonzero (2, 2) blocks. Some new iterative methods are presented and their convergence are studied in depth. By choosing different parameter matrices, we derive a series of existing and new iterative methods, including the preconditioned Uzawa method, the inexact Uzawa method, the SOR-like method, the GSOR method, the GIAOR method, the PIU method, the APIU method and so on. Numerical experiments are used to demonstrate the feasibility and effectiveness of the generalized parameterized inexact Uzawa methods.     

A parallel inexact newton method for stochastic programs with recourse

A parallel inexact Newton method with a line search is proposed for two-stage quadratic stochastic programs with recourse. A lattice rule is used for the numerical evaluation of multi-dimensional integrals, and a parallel iterative method is used to solve the quadratic programming subproblems. Although the objective only has a locally Lipschitz gradient, global convergence and local superlinear convergence of the method are established. Furthermore, the method provides an error estimate which does not require much extra computation. The performance of the method is illustrated on a CM5 parallel computer.This work was supported by the Australian Research Council and the numerical experiments were done on the Sydney Regional Centre for Parallel Computing CM5.     

A descent method with inexact linear search for nonsmooth equilibrium problems

A descent method with respect to the gap function is formulated and justified for the nonsmooth equilibrium problem. It uses the procedure of inexact linear search of the Armijo type. The proposed method converges under the same assumptions as the methods with exact linear search.     

A Power–Arnoldi algorithm for computing PageRank

The PageRank algorithm plays an important role in modern search engine technology. It involves using the classical power method to compute the principle eigenvector of the Google matrix representing the web link graph. However, when the largest eigenvalue is not well separated from the second one, the power method may perform poorly. This happens when the damping factor is sufficiently close to 1. Therefore, it is worth developing new techniques that are more sophisticated than the power method. The approach presented here, called Power–Arnoldi, is based on a periodic combination of the power method with the thick restarted Arnoldi algorithm. The justification for this new approach is presented. Numerical tests illustrate the efficiency and convergence behaviour of the new algorithm. Copyright © 2007 John Wiley & Sons, Ltd.     

Convergence estimates of nonrestarted and restarted block‐Lanczos methods

The block‐Lanczos method serves to compute a moderate number of eigenvalues and the corresponding invariant subspace of a symmetric matrix. In this paper, the convergence behavior of nonrestarted and restarted versions of the block‐Lanczos method is analyzed. For the nonrestarted version, we improve an estimate by Saad by means of a change of the auxiliary vector so that the new estimate is much more accurate in the case of clustered or multiple eigenvalues. For the restarted version, an estimate by Knyazev is generalized by extending our previous results on block steepest descent iterations and single‐vector restarted Krylov subspace iterations. The new estimates can also be reformulated and applied to invert‐block‐Lanczos methods for solving generalized matrix eigenvalue problems.     

A unified rational Krylov method for elliptic and parabolic fractional diffusion problems

We present a unified framework to efficiently approximate solutions to fractional diffusion problems of stationary and parabolic type. After discretization, we can take the point of view that the solution is obtained by a matrix-vector product of the form $f^{\mathit{\tau }}(L)\mathbf{b}$, where $L$ is the discretization matrix of the spatial operator, $\mathbf{b}$ a prescribed vector, and $f^{\mathit{\tau }}$ a parametric function, such as a fractional power or the Mittag-Leffler function. In the abstract framework of Stieltjes and complete Bernstein functions, to which the functions we are interested in belong to, we apply a rational Krylov method and prove uniform convergence when using poles based on Zolotarëv's minimal deviation problem. The latter are particularly suited for fractional diffusion as they allow for an efficient query of the map $\mathit{\tau }\mapsto f^{\mathit{\tau }}(L)\mathbf{b}$ and do not degenerate as the fractional parameters approach zero. We also present a variety of both novel and existing pole selection strategies for which we develop a computable error certificate. Our numerical experiments comprise a detailed parameter study of space-time fractional diffusion problems and compare the performance of the poles with the ones predicted by our certificate.     

Tuned preconditioners for inexact two‐sided inverse and Rayleigh quotient iteration

Convergence results are provided for inexact two‐sided inverse and Rayleigh quotient iteration, which extend the previously established results to the generalized non‐Hermitian eigenproblem and inexact solves with a decreasing solve tolerance. Moreover, the simultaneous solution of the forward and adjoint problem arising in two‐sided methods is considered, and the successful tuning strategy for preconditioners is extended to two‐sided methods, creating a novel way of preconditioning two‐sided algorithms. Furthermore, it is shown that inexact two‐sided Rayleigh quotient iteration and the inexact two‐sided Jacobi‐Davidson method (without subspace expansion) applied to the generalized preconditioned eigenvalue problem are equivalent when a certain number of steps of a Petrov–Galerkin–Krylov method is used and when this specific tuning strategy is applied. Copyright © 2014 John Wiley & Sons, Ltd.     

Restarted Full Orthogonalization Method for Shifted Linear Systems

Restarted GMRES is known to be inefficient for solving shifted systems when the shifts are handled simultaneously. Variants have been proposed to enhance its performance. We show that another restarted method, restarted Full Orthogonalization Method (FOM), can effectively be employed. The total number of iterations of restarted FOM applied to all shifted systems simultaneously is the same as that obtained by applying restarted FOM to the shifted system with slowest convergence rate, while the computational cost grows only sub-linearly with the number of shifts. Numerical experiments are reported.     

A multigrid-based preconditioned Krylov subspace method for the Helmholtz equation with PML

In this paper, we generalize the complex shifted Laplacian preconditioner to the complex shifted Laplacian-PML preconditioner for the Helmholtz equation with perfectly matched layer (Helmholtz-PML equation). The Helmholtz-PML equation is discretized by an optimal 9-point difference scheme, and the preconditioned linear system is solved by the Krylov subspace method, especially by the biconjugate gradient stabilized method (Bi-CGSTAB). The spectral analysis of the linear system is given, and a new matrix-based interpolation operator is proposed for the multigrid method, which is used to approximately invert the preconditioner. The numerical experiments are presented to illustrate the efficiency of the preconditioned Bi-CGSTAB method with the multigrid based on the new interpolation operator, also, numerical results are given for comparing the performance of the new interpolation operator with that of classic bilinear interpolation operator and the one suggested in Erlangga et al. (2006) [10].     

Preconditioner updates for solving sequences of linear systems in matrix‐free environment

We present two new ways of preconditioning sequences of nonsymmetric linear systems in the special case where the implementation is matrix free. Both approaches are fully algebraic, they are based on the general updates of incomplete LU decompositions recently introduced in (SIAM J. Sci. Comput. 2007; 29 (5):1918–1941), and they may be directly embedded into nonlinear algebraic solvers. The first of the approaches uses a new model of partial matrix estimation to compute the updates. The second approach exploits separability of function components to apply the updated factorized preconditioner via function evaluations with the discretized operator. Experiments with matrix‐free implementations of test problems show that both new techniques offer useful, robust and black‐box solution strategies. In addition, they show that the new techniques are often more efficient in matrix‐free environment than either recomputing the preconditioner from scratch for every linear system of the sequence or than freezing the preconditioner throughout the whole sequence. Copyright © 2010 John Wiley & Sons, Ltd.