Memory‐efficient Arnoldi algorithms for linearizations of matrix polynomials in Chebyshev basis

Novel memory‐efficient Arnoldi algorithms for solving matrix polynomial eigenvalue problems are presented. More specifically, we consider the case of matrix polynomials expressed in the Chebyshev basis, which is often numerically more appropriate than the standard monomial basis for a larger degree d. The standard way of solving polynomial eigenvalue problems proceeds by linearization, which increases the problem size by a factor d. Consequently, the memory requirements of Krylov subspace methods applied to the linearization grow by this factor. In this paper, we develop two variants of the Arnoldi method that build the Krylov subspace basis implicitly, in a way that only vectors of length equal to the size of the original problem need to be stored. The proposed variants are generalizations of the so‐called quadratic Arnoldi method and two‐level orthogonal Arnoldi procedure methods, which have been developed for the monomial case. We also show how the typical ingredients of a full implementation of the Arnoldi method, including shift‐and‐invert and restarting, can be incorporated. Numerical experiments are presented for matrix polynomials up to degree 30 arising from the interpolation of nonlinear eigenvalue problems, which stem from boundary element discretizations of PDE eigenvalue problems. Copyright © 2013 John Wiley & Sons, Ltd.     

Rational Krylov for Nonlinear Eigenproblems,an Iterative Projection Method

In recent papers Ruhe suggested a rational Krylov method for nonlinear eigenproblems knitting together a secant method for linearizing the nonlinear problem and the Krylov method for the linearized problem. In this note we point out that the method can be understood as an iterative projection method. Similarly to the Arnoldi method the search space is expanded by the direction from residual inverse iteration. Numerical methods demonstrate that the rational Krylov method can be accelerated considerably by replacing an inner iteration by an explicit solver of projected problems.     

Changing poles in the rational Lanczos method for the Hermitian eigenvalue problem

Applications such as the modal analysis of structures and acoustic cavities require a number of eigenvalues and eigenvectors of large‐scale Hermitian eigenvalue problems. The most popular method is probably the spectral transformation Lanczos method. An important disadvantage of this method is that a change of pole requires a complete restart. In this paper, we investigate the use of the rational Krylov method for this application. This method does not require a complete restart after a change of pole. It is shown that the change of pole can be considered as a change of Lanczos basis. The major conclusion of this paper is that the method is numerically stable when the poles are chosen in between clusters of the approximate eigenvalues. Copyright © 2001 John Wiley & Sons, Ltd.     

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.     

The implicit application of a rational filter in the RKS method

The implicitly restarted Arnoldi method implicitly applies a polynomial filter to the Arnoldi vectors by use of orthogonal transformations. In this paper, an implicit filtering by rational functions is proposed for the rational Krylov method. This filtering is performed in an efficient way. Two applications are considered. The first one is the filtering of unwanted eigenvalues using exact shifts. This approach is related to the use of exact shifts in the implicitly restarted Arnoldi method. Second, eigenvalue problems can have an infinite eigenvalue without physical relevance. This infinite eigenvalue can corrupt the eigensolution. An implicit filtering is proposed for avoiding such corruptions. The work of Gorik De Samblanx and Adhemar Bultheel was supported by the National Fund for Scientific Research (NFWO), project Lanczos, grant #2.0042.93 and by the Human Capital and Mobility project ROLLS of the European Community under contract ERBCHRXCT930416. The research by Karl Meerbergen was supported by the Belgian programme on Interuniversity Poles of Attraction (IUAP 17), initiated by the Belgian State—Prime Minister's Service—Federal Office for Scientific, Technical and Cultural Affairs and the project Iterative Methods in Scientific Computing, contract number HCM network CHRCCT93-0420, coordinated by CERFACS, Toulouse, France.     

A method for rational Chebyshev approximation of rational functions on the unit disk and on the unit interval

A method for construction of CF approximants in some cases of rational approximation of a rational function f on the unit disk and on the unit interval is presented. The inverted square root of the greatest positive eigenvalue and a corresponding eigenvector of an eigenvalue problem defined by the coefficients of f gives the solution.     

Structured eigenvalue problems for rational gauss quadrature

The connection between Gauss quadrature rules and the algebraic eigenvalue problem for a Jacobi matrix was first exploited in the now classical paper by Golub and Welsch (Math. Comput. 23(106), 221–230, 1969). From then on many computational problems arising in the construction of (polynomial) Gauss quadrature formulas have been reduced to solving direct and inverse eigenvalue problems for symmetric tridiagonals. Over the last few years (rational) generalizations of the classical Gauss quadrature formulas have been studied, i.e., formulas integrating exactly in spaces of rational functions. This paper wants to illustrate that stable and efficient procedures based on structured numerical linear algebra techniques can also be devised for the solution of the eigenvalue problems arising in the field of rational Gauss quadrature.     

Rational Krylov sequence methods for eigenvalue computation

Algorithms to solve large sparse eigenvalue problems are considered. A new class of algorithms which is based on rational functions of the matrix is described. The Lanczos method, the Arnoldi method, the spectral transformation Lanczos method, and Rayleigh quotient iteration all are special cases, but there are also new algorithms which correspond to rational functions with several poles. In the simplest case a basis of a rational Krylov subspace is found in which the matrix eigenvalue problem is formulated as a linear matrix pencil with a pair of Hessenberg matrices.     

Efficient computation of the exponential operator for large,sparse, symmetric matrices

In this paper we compare Krylov subspace methods with Chebyshev series expansion for approximating the matrix exponential operator on large, sparse, symmetric matrices. Experimental results upon negative‐definite matrices with very large size, arising from (2D and 3D) FE and FD spatial discretization of linear parabolic PDEs, demonstrate that the Chebyshev method can be an effective alternative to Krylov techniques, especially when memory bounds do not allow the storage of all Ritz vectors. We also discuss the sensitivity of Chebyshev convergence to extreme eigenvalue approximation, as well as the reliability of various a priori and a posteriori error estimates for both methods. Copyright © 2000 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.     

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.     

Compact linearization for binary quadratic problems

We show that a well-known linearization technique initially proposed for quadratic assignment problems can be generalized to a broader class of quadratic 0–1 mixed-integer problems subject to assignment constraints. The resulting linearized formulation is more compact and tighter than that obtained with a more usual linearization technique. We discuss the application of the compact linearization to three classes of problems in the literature, among which the graph partitioning problem.      

A rank‐exploiting infinite Arnoldi algorithm for nonlinear eigenvalue problems

We consider the nonlinear eigenvalue problem M(λ)x = 0, where M(λ) is a large parameter‐dependent matrix. In several applications, M(λ) has a structure where the higher‐order terms of its Taylor expansion have a particular low‐rank structure. We propose a new Arnoldi‐based algorithm that can exploit this structure. More precisely, the proposed algorithm is equivalent to Arnoldi's method applied to an operator whose reciprocal eigenvalues are solutions to the nonlinear eigenvalue problem. The iterates in the algorithm are functions represented in a particular structured vector‐valued polynomial basis similar to the construction in the infinite Arnoldi method [Jarlebring, Michiels, and Meerbergen, Numer. Math., 122 (2012), pp. 169–195]. In this paper, the low‐rank structure is exploited by applying an additional operator and by using a more compact representation of the functions. This reduces the computational cost associated with orthogonalization, as well as the required memory resources. The structure exploitation also provides a natural way in carrying out implicit restarting and locking without the need to impose structure in every restart. The efficiency and properties of the algorithm are illustrated with two large‐scale problems. Copyright © 2016 John Wiley & Sons, Ltd.     

Superlinear convergence of the rational Arnoldi method for the approximation of matrix functions

A superlinear convergence bound for rational Arnoldi approximations to functions of matrices is derived. This bound generalizes the well-known superlinear convergence bound for the conjugate gradient method to more general functions with finite singularities and to rational Krylov spaces. A constrained equilibrium problem from potential theory is used to characterize a max-min quotient of a nodal rational function underlying the rational Arnoldi approximation, where an additional external field is required for taking into account the poles of the rational Krylov space. The resulting convergence bound is illustrated at several numerical examples, in particular, the convergence of the extended Krylov method for the matrix square root.     

On IGMRES: An incomplete generalized minimal residual method for large unsymmetric linear systems

The truncated version of the generalized minimal residual method (GMRES), the incomplete generalized minimal residual method (IGMRES), is studied. It is based on an incomplete orthogonalization of the Krylov vectors in question, and gives an approximate or quasi-minimum residual solution over the Krylov subspace. A convergence analysis of this method is given, showing that in the non-restarted version IGMRES can behave like GMRES once the basis vectors of Krylov subspace generated by the incomplete orthogonalization are strongly linearly independent. Meanwhile, some relationships between the residual norms for IOM and IGMRES are established. Numerical experiments are reported to show convergence behavior of IGMRES and of its restarted version IGMRES(m). Project supported by the China State Key Basic Researches, the National Natural Science Foundation of China (Grant No. 19571014), the Doctoral Program (97014113), the Foundation of Returning Scholars of China and the Natural Science Foundation of Liaoning Province.     

Adaptive computation of smallest eigenvalues of self‐adjoint elliptic partial differential equations

We consider a new adaptive finite element (AFEM) algorithm for self‐adjoint elliptic PDE eigenvalue problems. In contrast to other approaches we incorporate the inexact solutions of the resulting finite‐dimensional algebraic eigenvalue problems into the adaptation process. In this way we can balance the costs of the adaptive refinement of the mesh with the costs for the iterative eigenvalue method. We present error estimates that incorporate the discretization errors, approximation errors in the eigenvalue solver and roundoff errors, and use these for the adaptation process. We show that it is also possible to restrict to very few iterations of a Krylov subspace solver for the eigenvalue problem on coarse meshes. Several examples are presented to show that this new approach achieves much better complexity than the previous AFEM approaches which assume that the algebraic eigenvalue problem is solved to full accuracy. Copyright © 2010 John Wiley & Sons, Ltd.     

Applications of realizations (aka linearizations) to free probability

We show how the combination of new “linearization” ideas in free probability theory with the powerful “realization” machinery – developed over the last 50 years in fields including systems engineering and automata theory – allows solving the problem of determining the eigenvalue distribution (or even the Brown measure, in the non-selfadjoint case) of noncommutative rational functions of random matrices when their size tends to infinity. Along the way we extend evaluations of noncommutative rational expressions from matrices to stably finite algebras, e.g. type II₁ von Neumann algebras, with a precise control of the domains of the rational expressions.The paper provides sufficient background information, with the intention that it should be accessible both to functional analysts and to algebraists.     

A rational Arnoldi process with applications

The rational Arnoldi process is a popular method for the computation of a few eigenvalues of a large non‐Hermitian matrix and for the approximation of matrix functions. The method is particularly attractive when the rational functions that determine the process have only few distinct poles , because then few factorizations of matrices of the form A ? z_jI have to be computed. We discuss recursion relations for orthogonal bases of rational Krylov subspaces determined by rational functions with few distinct poles. These recursion formulas yield a new implementation of the rational Arnoldi process. Applications of the rational Arnoldi process to the approximation of matrix functions as well as to the computation of eigenvalues and pseudospectra of A are described. The new implementation is compared to several available implementations.     

The Rational Krylov Algorithm for Nonlinear Matrix Eigenvalue Problems

It is shown how the rational Krylov algorithm can be applied to a matrix eigenvalue problem that is nonlinear in the eigenvalue parameter. Bibliography: 6 titles.     

Subspace methods for three‐parameter eigenvalue problems

We propose subspace methods for three‐parameter eigenvalue problems. Such problems arise when separation of variables is applied to separable boundary value problems; a particular example is the Helmholtz equation in ellipsoidal and paraboloidal coordinates. While several subspace methods for two‐parameter eigenvalue problems exist, their extensions to a three‐parameter setting seem challenging. An inherent difficulty is that, while for two‐parameter eigenvalue problems, we can exploit a relation to Sylvester equations to obtain a fast Arnoldi‐type method, such a relation does not seem to exist when there are three or more parameters. Instead, we introduce a subspace iteration method with projections onto generalized Krylov subspaces that are constructed from scratch at every iteration using certain Ritz vectors as the initial vectors. Another possibility is a Jacobi–Davidson‐type method for three or more parameters, which we generalize from its two‐parameter counterpart. For both approaches, we introduce a selection criterion for deflation that is based on the angles between left and right eigenvectors. The Jacobi–Davidson approach is devised to locate eigenvalues close to a prescribed target; yet, it often also performs well when eigenvalues are sought based on the proximity of one of the components to a prescribed target. The subspace iteration method is devised specifically for the latter task. The proposed approaches are suitable especially for problems where the computation of several eigenvalues is required with high accuracy. MATLAB implementations of both methods have been made available in the package MultiParEig (see http://www.mathworks.com/matlabcentral/fileexchange/47844-multipareig ).