Showing papers in "SIAM Journal on Matrix Analysis and Applications in 2000"
TL;DR: There is a strong analogy between several properties of the matrix and the higher-order tensor decomposition; uniqueness, link with the matrix eigenvalue decomposition, first-order perturbation effects, etc., are analyzed.
Abstract: We discuss a multilinear generalization of the singular value decomposition. There is a strong analogy between several properties of the matrix and the higher-order tensor decomposition; uniqueness, link with the matrix eigenvalue decomposition, first-order perturbation effects, etc., are analyzed. We investigate how tensor symmetries affect the decomposition and propose a multilinear generalization of the symmetric eigenvalue decomposition for pair-wise symmetric tensors.
4,101 citations
TL;DR: A multilinear generalization of the best rank-R approximation problem for matrices, namely, the approximation of a given higher-order tensor, in an optimal least-squares sense, by a tensor that has prespecified column rank value, rowRank value, etc.
Abstract: In this paper we discuss a multilinear generalization of the best rank-R approximation problem for matrices, namely, the approximation of a given higher-order tensor, in an optimal least-squares sense, by a tensor that has prespecified column rank value, row rank value, etc For matrices, the solution is conceptually obtained by truncation of the singular value decomposition (SVD); however, this approach does not have a straightforward multilinear counterpart We discuss higher-order generalizations of the power method and the orthogonal iteration method
1,638 citations
TL;DR: The problem of finding good preconditioners for the numerical solution of indefinite linear systems is considered and special emphasis is put on preconditionsers that have a 2 × 2 block structure and that incorporate the (1,2 and (2,1) blocks of the original matrix.
Abstract: The problem of finding good preconditioners for the numerical solution of indefinite linear systems is considered. Special emphasis is put on preconditioners that have a 2 × 2 block structure and that incorporate the (1,2) and (2,1) blocks of the original matrix. Results concerning the spectrum and form of the eigenvectors of the preconditioned matrix and its minimum polynomial are given. The consequences of these results are considered for a variety of Krylov subspace methods. Numerical experiments validate these conclusions.
365 citations
TL;DR: This work presents a common framework for efficient algorithms that regularize after this second projection rather than before it, and shows that determining regularization parameters based on the final projected problem rather than on the original discretization has firmer justification and often involves less computational expense.
Abstract: Numerical solution of ill-posed problems is often accomplished by discretization (projection onto a finite dimensional subspace) followed by regularization. If the discrete problem has high dimension, though, typically we compute an approximate solution by projecting the discrete problem onto an even smaller dimensional space, via iterative methods based on Krylov subspaces. In this work we present a common framework for efficient algorithms that regularize after this second projection rather than before it. We show that determining regularization parameters based on the final projected problem rather than on the original discretization has firmer justification and often involves less computational expense. We prove some results on the approximate equivalence of this approach to other forms of regularization, and we present numerical examples.
321 citations
TL;DR: A restarted variant of the Lanczos method for symmetric eigenvalue problems named the thick-restart Lanczo method is proposed, able to retain an arbitrary number of Ritz vectors from the previous iterations with a minimal restarting cost.
Abstract: In this paper, we propose a restarted variant of the Lanczos method for symmetric eigenvalue problems named the thick-restart Lanczos method. This new variant is able to retain an arbitrary number of Ritz vectors from the previous iterations with a minimal restarting cost. Since it restarts with Ritz vectors, it is simpler than similar methods, such as the implicitly restarted Lanczos method. We carefully examine the effects of the floating-point round-off errors on stability of the new algorithm and present an implementation of the partial reorthogonalization scheme that guarantees accurate Ritz values with a minimal amount of reorthogonalization. We also show a number of heuristics on deciding which Ritz pairs to save during restart in order to maximize the overall performance of the thick-restart Lanczos method.
301 citations
TL;DR: This work considers bipartite matching algorithms for computing permutations of a sparse matrix so that the diagonal of the permuted matrix has entries of large absolute value and considers scaling techniques to increase the relative values of the diagonal entries.
Abstract: We consider bipartite matching algorithms for computing permutations of a sparse matrix so that the diagonal of the permuted matrix has entries of large absolute value. We discuss various strategies for this and consider their implementation as computer codes. We also consider scaling techniques to further increase the relative values of the diagonal entries. Numerical experiments show the effect of the reorderings and the scaling on the solution of sparse equations by a direct method and by preconditioned iterative techniques.
280 citations
TL;DR: This paper provides an iterative algorithm to jointly approximately diagonalize K Hermitian positive definite matrices and calculates the matrix B which minimizes the criterionsum_{k=1}^K n_k, nk being positive numbers, which is a measure of the deviation from diagonality of the matrices BCkB*.
Abstract: This paper provides an iterative algorithm to jointly approximately diagonalize K Hermitian positive definite matrices ${\bf\Gamma}_1$, \dots, ${\bf\Gamma}_K$. Specifically, it calculates the matrix B which minimizes the criterion $\sum_{k=1}^K n_k [\log \det\diag(\B\C_k\B^*) - \log\det(\B\C_k\B^*)]$, nk being positive numbers, which is a measure of the deviation from diagonality of the matrices BCkB*$. The convergence of the algorithm is discussed and some numerical experiments are performed showing the good performance of the algorithm.
257 citations
TL;DR: A block generalization of the 1-norm power method underlying the estimator is derived and developed into a practical algorithm applicable to both real and complex matrices, and it is shown how it can be used to efficiently approximate 1- norm pseudospectra.
Abstract: The matrix 1-norm estimation algorithm used in LAPACK and various other software libraries and packages has proved to be a valuable tool. However, it has the limitations that it offers the user no control over the accuracy and reliability of the estimate and that it is based on level 2 BLAS operations. A block generalization of the 1-norm power method underlying the estimator is derived here and developed into a practical algorithm applicable to both real and complex matrices. The algorithm works with n × t matrices, where t is a parameter. For t=1 the original algorithm is recovered, but with two improvements (one for real matrices and one for complex matrices). The accuracy and reliability of the estimates generally increase with t and the computational kernels are level 3 BLAS operations for t > 1. The last t-1 columns of the starting matrix are randomly chosen, giving the algorithm a statistical flavor. As a by-product of our investigations we identify a matrix for which the 1-norm power method takes the maximum number of iterations. As an application of the new estimator we show how it can be used to efficiently approximate 1-norm pseudospectra.
164 citations
TL;DR: This work presents algorithms that use implicit restarting in order to retain information that deflates the smallest eigenvalues and thus improves the convergence of the generalized minimum residual method.
Abstract: The generalized minimum residual method (GMRES) is well known for solving large nonsymmetric systems of linear equations. It generally uses restarting, which slows the convergence. However, some information can be retained at the time of the restart and used in the next cycle. We present algorithms that use implicit restarting in order to retain this information. Approximate eigenvectors determined from the previous subspace are included in the new subspace. This deflates the smallest eigenvalues and thus improves the convergence. The subspace that contains the approximate eigenvectors is itself a Krylov subspace, but not with the usual starting vector. The implicitly restarted FOM algorithm includes standard Ritz vectors in the subspace. The eigenvalue portion of its calculations is equivalent to Sorensen's IRA algorithm. The implicitly restarted GMRES algorithm uses harmonic Ritz vectors. This algorithm also gives a new approach to computing interior eigenvalues.
155 citations
TL;DR: It is shown that the Lagrangian dual based on relaxing the constraints XXT=I and the seemingly redundant constraints XT X=I has a zero duality gap, which has natural applications to quadratic assignment and graph partitioning problems, as well as the problem of minimizing the weighted sum of the largest eigenvalues of a matrix.
Abstract: Quadratically constrained quadratic programs (QQPs) play an important modeling role for many diverse problems These problems are in general NP hard and numerically intractable Lagrangian relaxations often provide good approximate solutions to these hard problems Such relaxations are equivalent to semidefinite programming relaxations
For several special cases of QQP, eg, convex programs and trust region subproblems, the Lagrangian relaxation provides the exact optimal value, ie, there is a zero duality gap However, this is not true for the general QQP, or even the QQP with two convex constraints, but a nonconvex objective
In this paper we consider a certain QQP where the quadratic constraints correspond to the matrix orthogonality condition XXT=I For this problem we show that the Lagrangian dual based on relaxing the constraints XXT=I and the seemingly redundant constraints XT X=I has a zero duality gap This result has natural applications to quadratic assignment and graph partitioning problems, as well as the problem of minimizing the weighted sum of the largest eigenvalues of a matrix We also show that the technique of relaxing quadratic matrix constraints can be used to obtain a strengthened semidefinite relaxation for the max-cut problem
140 citations
TL;DR: A transformation-free form of this method that exploits incomplete Denman--Beavers square root iterations and aims for a specified accuracy (ignoring roundoff) is presented.
Abstract: The standard inverse scaling and squaring algorithm for computing the matrix logarithm begins by transforming the matrix to Schur triangular form in order to facilitate subsequent matrix square root and Pade approximation computations. A transformation-free form of this method that exploits incomplete Denman--Beavers square root iterations and aims for a specified accuracy (ignoring roundoff) is presented. The error introduced by using approximate square roots is accounted for by a novel splitting lemma for logarithms of matrix products. The number of square root stages and the degree of the final Pade approximation are chosen to minimize the computational work. This new method is attractive for high-performance computation since it uses only the basic building blocks of matrix multiplication, LU factorization and matrix inversion.
TL;DR: Newton's method and a class of basic fixed-point iterations can be used to find its minimal positive solution whenever it has a positive solution of any equation in this class of nonsymmetric algebraic Riccati equations.
Abstract: We consider the iterative solution of a class of nonsymmetric algebraic Riccati equations, which includes a class of algebraic Riccati equations arising in transport theory. For any equation in this class, Newton's method and a class of basic fixed-point iterations can be used to find its minimal positive solution whenever it has a positive solution. The properties of these iterative methods are studied and some practical issues are addressed. An algorithm is then proposed to find the minimal positive solution efficiently. Numerical results are also given.
TL;DR: The backward greedy algorithm is shown to be optimal for the subset selection problem in the sense that it selects the "correct" subset of columns from A if the perturbation of the data vector b is small enough.
Abstract: The following linear inverse problem is considered: Given a full column rank m × n data matrix A, and a length m observation vector b, find the best least-squares solution to A x = b with at most r < n nonzero components. The backward greedy algorithm computes a sparse solution to A x = b by removing greedily columns from A until r columns are left. A simple implementation based on a QR downdating scheme using Givens rotations is described. The backward greedy algorithm is shown to be optimal for the subset selection problem in the sense that it selects the "correct" subset of columns from A if the perturbation of the data vector b is small enough. The results generalize to any other norm of the residual.
TL;DR: A matrix with positive row sums and all its off-diagonal elements bounded above by their corresponding row means is called a B-matrix and it is proved that the class of B-Matrices is a subset of theclass of P-matrices.
Abstract: A matrix with positive row sums and all its off-diagonal elements bounded above by their corresponding row means is called a B-matrix. It is proved that the class of B-matrices is a subset of the class of P-matrices. Properties of B-matrices are used to localize the real eigenvalues of a real matrix and the real parts of all eigenvalues of a real matrix.
TL;DR: New theoretical results for the conjugate gradient method are given and a new version is proposed that generates a Krylov subspace which can be efficiently recycled thanks to orthogonal projections in subsequent systems.
Abstract: Many scientific applications require one to solve successively linear systems Ax=b with different right-hand sides $b$ and a symmetric positive definite matrix A. The conjugate gradient method applied to the first system generates a Krylov subspace which can be efficiently recycled thanks to orthogonal projections in subsequent systems. A modified conjugate gradient method is then applied with a specific initial guess and initial descent direction and a modified descent direction during the iterations. This paper gives new theoretical results for this method and proposes a new version. Numerical experiments show the efficacy of our method even for quite different right-hand sides.
TL;DR: A new set of algorithms for computation of matrix rational interpolants and one-sided matrix greatest common divisors, suitable for computation in exact arithmetic domains where growth of coefficients in intermediate computations is a central concern.
Abstract: We present a new set of algorithms for computation of matrix rational interpolants and one-sided matrix greatest common divisors. Examples of these interpolants include Pade approximants, Newton--Pade, Hermite--Pade, and simultaneous Pade approximants, and more generally M-Pade approximants along with their matrix generalizations. The algorithms are fast and compute all solutions to a given problem. Solutions for all (possibly singular) subproblems along offdiagonal paths in a solution table are also computed by stepping around singular blocks on a path corresponding to "closest" regular interpolation problems.
The algorithms are suitable for computation in exact arithmetic domains where growth of coefficients in intermediate computations is a central concern. This coefficient growth is avoided by using fraction-free methods. At the same time, the methods are fast in the sense that they are at least an order of magnitude faster than existing fraction-free methods for the corresponding problems. The methods make use of linear systems having a special striped Krylov structure.
TL;DR: These techniques are based on an analysis and manipulation of the underlying graph structure, using the framework developed in an earlier paper on rank-1 modifications, to develop sparse techniques for updating the factorization of a sparse symmetric positive definite matrix.
Abstract: Given a sparse symmetric positive definite matrix $\mathbf{AA}\tr$ and an associated sparse Cholesky factorization $\mathbf{LDL}\tr$ or $\mathbf{LL}\tr$, we develop sparse techniques for updating the factorization after either adding a collection of columns to A or deleting a collection of columns from A. Our techniques are based on an analysis and manipulation of the underlying graph structure, using the framework developed in an earlier paper on rank-1 modifications [T. A. Davis and W. W. Hager, SIAM J. Matrix Anal. Appl., 20 (1999), pp. 606--627]. Computationally, the multiple-rank update has better memory traffic and executes much faster than an equivalent series of rank-1 updates since the multiple-rank update makes one pass through L computing the new entries, while a series of rank-1 updates requires multiple passes through L.
TL;DR: It is proved that the first m eigenvalues of Pn-1Tn tend to r and the last m tend to R, for any fixed m, and the exact limit value of the condition number of the preconditioned matrices is computed.
Abstract: We study the asymptotic behavior of the eigenvalues of Hermitian n × n block Toeplitz matrices Tn, with k × k blocks, as n tends to infinity. No hypothesis is made concerning the structure of the blocks. Such matrices {Tn} are generated by the Fourier coefficients of a Hermitian matrix-valued function $f\in L^2$, and we study the distribution of their eigenvalues for large n, relating their behavior to some properties of f as a function; in particular, we show that the distribution of the eigenvalues converges to a limit $\mu_f$, and we explicitly compute $\mu_f$ in terms of f, showing that $\int F\, d\mu_f=1/k\int\tr F(f)$. Some consequences of this distribution and some localization results for the eigenvalues of Tn are discussed. We also study the eigenvalues of the preconditioned matrices {Pn-1Tn}, where the sequence {Pn} is generated by a positive definite matrix-valued function p. We show that the spectrum of any Pn-1Tn is contained in the interval [r,R], where r is the smallest and R the largest eigenvalue of p-1f. We also prove that the first m eigenvalues of Pn-1Tn tend to r and the last m tend to R, for any fixed m. Finally, the exact limit value of the condition number of the preconditioned matrices is computed.
TL;DR: It is shown that a similar technique works for star embeddings when the Laplacian has a zero Dirichlet boundary condition, though the related eigenvalues in this case are reciprocals of each other.
Abstract: Graph embeddings are useful in bounding the smallest nontrivial eigenvalues of Laplacian matrices from below. For an n × n Laplacian, these embedding methods can be characterized as follows: The lower bound is based on a clique embedding into the underlying graph of the Laplacian. An embedding can be represented by a matrix $\Gamma$; the best possible bound based on this embedding is $n/\lambda_{\max} (\Gamma^T \Gamma)$, where $\lambda_{\max}$ indicates the largest eigenvalue of the specified matrix. However, the best bounds produced by embedding techniques are not tight; they can be off by a factor proportional to log2 n for some Laplacians.
We show that this gap is a result of the representation of the embedding: By including edge directions in the embedding matrix representation $\Gamma$, it is possible to find an embedding such that $\Gamma^T \Gamma$ has eigenvalues that can be put into a one-to-one correspondence with the eigenvalues of the Laplacian. Specifically, if $\lambda$ is a nonzero eigenvalue of either matrix, then $n / \lambda$ is an eigenvalue of the other. Simple transformations map the corresponding eigenvectors to each other. The embedding that produces these correspondences has a simple description in electrical terms if the underlying graph of the Laplacian is viewed as a resistive circuit. We also show that a similar technique works for star embeddings when the Laplacian has a zero Dirichlet boundary condition, though the related eigenvalues in this case are reciprocals of each other. In the zero Dirichlet boundary case, the embedding matrix $\Gamma$ can be used to construct the inverse of the Laplacian. Finally, we connect our results with previous techniques for producing bounds and provide an example.
TL;DR: It is shown that in the two three- term recurrences the contributions of the local roundoff errors to the analyzed gaps may be dramatically amplified while propagating through the algorithm, explaining the well-known behavior of three-term-based versions of the biconjugate gradient method.
Abstract: It has been widely observed that Krylov space solvers based on two three-term recurrences can give significantly less accurate residuals than mathematically equivalent solvers implemented with three two-term recurrences. In this paper we attempt to clarify and justify this difference theoretically by analyzing the gaps between recursively and explicitly computed residuals.
It is shown that, in contrast with the two-term recurrences analyzed by Sleijpen, van der Vorst, and Fokkema [ Numer. Algorithms, 7 (1994), pp. 75--109] and Greenbaum [SIAM J. Matrix Anal. Appl., 18 (1997), pp. 535--551], in the two three-term recurrences the contributions of the local roundoff errors to the analyzed gaps may be dramatically amplified while propagating through the algorithm. This result explains, for example, the well-known behavior of three-term-based versions of the biconjugate gradient method, where large gaps between recursively and explicitly computed residuals are not uncommon. For the conjugate gradient method, however, such a devastating behavior---although possible---is not observed frequently in practical computations, and the difference between two-term and three-term implementations is usually moderate or small. This can also be explained by our results.
TL;DR: A fast algorithm for the basic deconvolution problem is developed due to the low displacement rank of the involved matrices and the sparsity of the generators and Monte-Carlo simulations indicate the superior statistical performance of the structured total least squares estimator compared to other estimators such as the ordinary total least square estimator.
Abstract: In this paper we develop a fast algorithm for the basic deconvolution problem. First we show that the kernel problem to be solved in the basic deconvolution problem is a so-called structured total least squares problem. Due to the low displacement rank of the involved matrices and the sparsity of the generators, we are able to develop a fast algorithm. We apply the new algorithm on a deconvolution problem arising in a medical application in renography. By means of this example, we show the increased computational performance of our algorithm as compared to other algorithms for solving this type of structured total least squares problem. In addition, Monte-Carlo simulations indicate the superior statistical performance of the structured total least squares estimator compared to other estimators such as the ordinary total least squares estimator.
TL;DR: It is shown that the optimal Ak and Bk are banded Toeplitz matrices, and an efficient algorithm for computing the approximation is provided.
Abstract: This paper considers the problem of finding n × n matrices Ak and Bk that minimize $||T - \sum A_k \otimes B_k||_F$, where $\otimes$ denotes Kronecker product and T is a banded n × n block Toeplitz matrix with banded n × n Toeplitz blocks. It is shown that the optimal Ak and Bk are banded Toeplitz matrices, and an efficient algorithm for computing the approximation is provided. An image restoration problem from the Hubble Space Telescope (HST) is used to illustrate the effectiveness of an approximate SVD preconditioner constructed from the Kronecker product decomposition.
TL;DR: In this paper, Bhatia and Kittaneh showed that the weak log-majorization relations hold for any complex number z. This sharpens some results due to R.Bhatia et al.
Abstract: Let Mn be the space of n × n complex matrices. For $A\in M_n,$ let $s(A)\equiv (s_1(A),\dotsc,\break s_n(A)),$ where $s_1(A)\ge\cdots\ge s_n(A)$ are the singular values of A. We prove that if $A,B\in M_n$ are positive semidefinite, then (i) $s_j(A-B)\le s_j(A\oplus B), j=1,2, . . . ,n, and (ii) the weak log-majorization relations $s(A-|z|B)\prec_{wlog} s(A+zB)\prec_{wlog} s(A+|z|B)$ hold for any complex number z. This sharpens some results due to R. Bhatia and F. Kittaneh.
TL;DR: A new class of SAI smoothers is proposed in this paper and their analytic smoothing factors for constant coefficient PDEs are presented, demonstrating the effectiveness of this new technique.
Abstract: Various forms of sparse approximate inverses (SAI) have been shown to be useful for preconditioning. Their potential usefulness in a parallel environment has motivated much interest in recent years. However, the capability of an approximate inverse in eliminating the local error has not yet been fully exploited in multigrid algorithms. A careful examination of the iteration matrices of these approximate inverses indicates their superiority in smoothing the high-frequency error in addition to their inherent parallelism. We propose a new class of SAI smoothers in this paper and present their analytic smoothing factors for constant coefficient PDEs. The following are several distinctive features that make this technique special:
By adjusting the quality of the approximate inverse, the smoothing factor can be improved accordingly. For hard problems, this is useful. In contrast to the ordering sensitivity of other smoothing techniques, this technique is ordering independent. In general, the sequential performance of many superior parallel algorithms is not very competitive. This technique is useful in both parallel and sequential computations.
Our theoretical and numerical results have demonstrated the effectiveness of this new technique.
TL;DR: This paper describes, in an asymptotic sense, the region containing those eigenvalues that are well approximated by the Ritz values, and it is characterized by an extremal problem in potential theory which was first considered by Rakhmanov.
Abstract: When discussing the convergence properties of the Lanczos iteration method for the real symmetric eigenvalue problem, Trefethen and Bau noted that the Lanczos method tends to find eigenvalues in regions that have too little charge when compared to an equilibrium distribution. In this paper a quantitative version of this rule of thumb is presented. We describe, in an asymptotic sense, the region containing those eigenvalues that are well approximated by the Ritz values. The region depends on the distribution of eigenvalues and on the ratio between the size of the matrix and the number of iterations, and it is characterized by an extremal problem in potential theory which was first considered by Rakhmanov. We give examples showing the connection with the equilibrium distribution.
TL;DR: This paper shows that it is an operator splitting preconditioner and clusters eigenvalues for the normal equation matrix thus ensuring a fast convergence of the conjugate gradient normal method.
Abstract: Preconditioning techniques for dense linear systems arising from singular boundary integral equations are described and analyzed. A particular class of approximate inverse based preconditioners related to the mesh neighbor methods is known to be efficient. This paper shows that it is an operator splitting preconditioner and clusters eigenvalues for the normal equation matrix thus ensuring a fast convergence of the conjugate gradient normal method. Clustering of the eigenvalues of the preconditioned matrix and fast convergence of the generalized minimal residual method are also observed. For the type of problems considered, we demonstrate a crucial connection between two essential features of eigenvalue clustering for a sparse preconditioner---approximate inversion for a small cluster radius and operator splitting for a small cluster size. Experimental results from several boundary integral equations are presented.
TL;DR: A class of preconditioners based on the augmented matrix approach considered earlier by two of the present authors are devoted, which is able to improve significantly an eigenvalue distribution of certain severely ill-conditioned algebraic systems by using properly chosen projection matrices, which correct the low-frequency components in the spectrum.
Abstract: The present work is devoted to a class of preconditioners based on the augmented matrix approach considered earlier by two of the present authors It presents some generalizations of the subspace-correction schemes studied earlier and gives a brief comparison of the developed technique with a somewhat similar "deflation" algorithm
The developed preconditioners are able to improve significantly an eigenvalue distribution of certain severely ill-conditioned algebraic systems by using properly chosen projection matrices, which correct the low-frequency components in the spectrum One of the main advantages of the proposed approach is the possibility of using inexact solvers within the projectors Another attractive feature of the developed method is that it can be easily combined with other preconditioners, for instance, those which correct the high-frequency eigenmodes
TL;DR: The partial fraction method is found to be the best method overall and the benefits it brings to a transformation-free form of the inverse scaling and squaring method recently proposed by Cheng, Higham, Kenney, and Laub are illustrated.
Abstract: The inverse scaling and squaring method for evaluating the logarithm of a matrix takes repeated square roots to bring the matrix close to the identity, computes a Pade approximant, and then scales back. We analyze several methods for evaluating the Pade approximant, including Horner's method (used in some existing codes), suitably customized versions of the Paterson--Stockmeyer method and Van Loan's variant, and methods based on continued fraction and partial fraction expansions. The computational cost, storage, and numerical accuracy of the methods are compared. We find the partial fraction method to be the best method overall and illustrate the benefits it brings to a transformation-free form of the inverse scaling and squaring method recently proposed by Cheng, Higham, Kenney, and Laub [ SIAM J. Matrix Anal. Appl., 22 (2001), pp. 1112--1125]. We comment briefly on how the analysis carries over to the matrix exponential.
TL;DR: In this paper, the Lanczos basis is used for the solution of symmetric indefinite linear systems, by solving a reduced system in one way or another, and it is shown that the method of solution may lead, under certain circumstances, to large additional errors, which are not corrected by continuing the iteration process.
Abstract: The three-term Lanczos process for a symmetric matrix leads to bases for Krylov subspaces of increasing dimension. The Lanczos basis, together with the recurrence coefficients, can be used for the solution of symmetric indefinite linear systems, by solving a reduced system in one way or another. This leads to well-known methods: MINRES (minimal residual), GMRES (generalized minimal residual), and SYMMLQ (symmetric LQ). We will discuss in what way and to what extent these approaches differ in their sensitivity to rounding errors.
In our analysis we will assume that the Lanczos basis is generated in exactly the same way for the different methods, and we will not consider the errors in the Lanczos process itself. We will show that the method of solution may lead, under certain circumstances, to large additional errors, which are not corrected by continuing the iteration process.
Our findings are supported and illustrated by numerical examples.
TL;DR: An implicitly restarted symplectic Lanczos method for the symplectic eigenvalue problem is presented and is used to compute some eigenvalues and eigenvectors of large and sparse symplectic operators.
Abstract: An implicitly restarted symplectic Lanczos method for the symplectic eigenvalue problem is presented. The Lanczos vectors are constructed to form a symplectic basis. The inherent numerical difficulties of the symplectic Lanczos method are addressed by inexpensive implicit restarts. The method is used to compute some eigenvalues and eigenvectors of large and sparse symplectic operators.