On the loss of orthogonality in the second-order Arnoldi process

This note is devoted to the rounding error analysis of the second-order Arnoldi process for constructing an orthonormal basis of the second-order Krylov subspace. The effect of the rounding errors on the orthogonality of the derived vectors is given.     

Block Krylov Subspace Methods for Large Algebraic Riccati Equations

In the present paper, we present numerical methods for the computation of approximate solutions to large continuous-time and discrete-time algebraic Riccati equations. The proposed methods are projection methods onto block Krylov subspaces. We use the block Arnoldi process to construct an orthonormal basis of the corresponding block Krylov subspace and then extract low rank approximate solutions. We consider the sequential version of the block Arnoldi algorithm by incorporating a deflation technique which allows us to delete linearly and almost linearly dependent vectors in the block Krylov subspace sequences. We give some theoretical results and present numerical experiments for large problems.     

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.     

A model for understanding numerical stability

^** Email: bornemann{at}ma.tum.de We present a model of roundoff error analysis that combinessimplicity with predictive power. Though not considering allsources of roundoff within an algorithm, the model is relatedto a recursive roundoff error analysis and therefore is capableof correctly predicting stability or instability of an algorithm.By means of nontrivial examples, such as the componentwise backwardstability analysis of Gaussian elimination with a single iterativerefinement step, we demonstrate that the model even yields quantitativebackward error bounds that show all the known problem-dependentterms with the exception of dimension-dependent constants. Themodel can serve as a convenient tool for teaching or as a heuristicdevice to discover stability results before entering a furtherdetailed analysis.     

Roundoff Error Estimates of the Modified Gram–Schmidt Algorithm with Column Pivoting

In this paper we provide backward and forward roundoff error estimates of the modified Gram–Schmidt algorithm with column pivoting. It turns out that the row-wise growth factors of this algorithm are bounded. Furthermore, if the coefficient matrix is well conditioned, then with properly chosen tolerance , this algorithm can also correctly determine the numerical rank of the coefficient matrix. We also derive a forward roundoff error estimate of this algorithm for the solution of the least squares problem.     

Accurate solution of polynomial equations using Macaulay resultant matrices

We propose an algorithm for solving two polynomial equations in two variables. Our algorithm is based on the Macaulay resultant approach combined with new techniques, including randomization, to make the algorithm accurate in the presence of roundoff error. The ultimate computation is the solution of a generalized eigenvalue problem via the QZ method. We analyze the error due to roundoff of the method, showing that with high probability the roots are computed accurately, assuming that the input data (that is, the two polynomials) are well conditioned. Our analysis requires a novel combination of algebraic and numerical techniques.

     

A statistical model of roundoff error for varying length floating-point arithmetic

A model is presented which explains the behavior of the roundoff error in a result quantity when computing precision is varied. A set of hypotheses concerning this a posteriori model is tested in a matrix inversion algorithm. The characteristics of the algorithms where the error model is valid are discussed. As an application of the model, the usual estimation procedure for roundoff error consisting of comparing the results computed in two different precisions is analyzed statistically.     

Alternatives for parallel Krylov subspace basis computation

Numerical methods related to Krylov subspaces are widely used in large sparse numerical linear algebra. Vectors in these subspaces are manipulated via their representation onto orthonormal bases. Nowadays, on serial computers, the method of Arnoldi is considered as a reliable technique for constructing such bases. However, although easily parallelizable, this technique is not as scalable as expected for communications. In this work we examine alternative methods aimed at overcoming this drawback. Since they retrieve upon completion the same information as Arnoldi's algorithm does, they enable us to design a wide family of stable and scalable Krylov approximation methods for various parallel environments. We present timing results obtained from their implementation on two distributed-memory multiprocessor supercomputers: the Intel Paragon and the IBM Scalable POWERparallel SP2. © 1997 John Wiley & Sons, Ltd.     

Approximation of the linear combination of $\varphi$-functions using the block shift-and-invert Krylov subspace method

In this paper, we develop an algorithm in which the block shift-and-invert Krylov subspace method can be employed for approximating the linear combination of the matrix exponential and related exponential-type functions. Such evaluation plays a major role in a class of numerical methods known as exponential integrators. We derive a low-dimensional matrix exponential to approximate the objective function based on the block shift-and-invert Krylov subspace methods. We obtain the error expansion of the approximation, and show that the variants of its first term can be used as reliable a posteriori error estimates and correctors. Numerical experiments illustrate that the error estimates are efficient and the proposed algorithm is worthy of further study.     

IMinpert： An Incomplete Minimum Perturbation Algorithm for Large Unsymmetric Linear Systems

This paper gives the truncated version of the Minpert method:the incomplete minimum perturbation algorithm(IMinpert).It is based on an incomplete orthogonal- ization of the Krylov vectors in question,and gives a quasi-minimum backward error solution over the Krylov subspace.In order to make the practical implementation of IMinpert easy and convenient,we give another approximate version of the IMinpert method:A-IMinpert.Theoretical properties of the latter algorithm are discussed.Nu- merical experiments are reported to show the proposed method is effective in practice and is competitive with the Minpert algorithm.     

A modified fixed-point iterative algorithm for image restoration using fourth-order PDE model

In this paper, we propose a modified fixed point iterative algorithm to solve the fourth-order PDE model for image restoration problem. Compared with the standard fixed point algorithm, the proposed algorithm needn?t to compute inverse matrices so that it can speed up the convergence and reduce the roundoff error. Furthermore, we prove the convergence of the proposed algorithm and give some experimental results to illustrate its effectiveness by comparing with the standard fixed point algorithm, the time marching algorithm and the split Bregman algorithm.     

Improved seed methods for symmetric positive definite linear equations with multiple right‐hand sides

We consider symmetric positive definite systems of linear equations with multiple right‐hand sides. The seed conjugate gradient (CG) method solves one right‐hand side with the CG method and simultaneously projects over the Krylov subspace thus developed for the other right‐hand sides. Then the next system is solved and used to seed the remaining ones. Rounding error in the CG method limits how much the seeding can improve convergence. We propose three changes to the seed CG method: only the first right‐hand side is used for seeding, this system is solved past convergence, and the roundoff error is controlled with some reorthogonalization. We will show that results are actually better with only one seeding, even in the case of related right‐hand sides. Controlling rounding error gives the potential for rapid convergence for the second and subsequent right‐hand sides. Polynomial preconditioning can help reduce storage needed for reorthogonalization. The new seed methods are applied to examples including matrices from quantum chromodynamics. Copyright © 2013 John Wiley & Sons, Ltd.     

Numerical stability of nonequispaced fast Fourier transforms

This paper presents some new results on numerical stability for multivariate fast Fourier transform of nonequispaced data (NFFT). In contrast to fast Fourier transform (of equispaced data), the NFFT is an approximate algorithm. In a worst case study, we show that both approximation error and roundoff error have a strong influence on the numerical stability of NFFT. Numerical tests confirm the theoretical estimates of numerical stability.     

On solution of large systems of linear equations with block-Toeplitz banded matrices

A novel numerical algorithm is proposed for solving systems of linear equations with block-Toeplitz narrow-banded matrices. Its arithmetical cost is about double of that of the block cyclic reduction. The backward roundoff error stability of the method guarantees its reliability for the matrices that are not symmetric positive definite or diagonally dominant.     

A harmonic restarted Arnoldi algorithm for calculating eigenvalues and determining multiplicity

A restarted Arnoldi algorithm is given that computes eigenvalues and eigenvectors. It is related to implicitly restarted Arnoldi, but has a simpler restarting approach. Harmonic and regular Rayleigh-Ritz versions are possible.For multiple eigenvalues, an approach is proposed that first computes eigenvalues with the new harmonic restarted Arnoldi algorithm, then uses random restarts to determine multiplicity. This avoids the need for a block method or for relying on roundoff error to produce the multiple copies.     

An Arnoldi based algorithm for large algebraic Riccati equations

In the present work, we present a numerical method for the computation of approximate solutions for large continuous-time algebraic Riccati equations. The proposed method is a method of projection onto a matrix Krylov subspace. We use a matrix Arnoldi process to construct an orthonormal basis. We give some theoretical results and numerical experiments for large problems.     

Roundoff errors in floating-point summation

All possible schemes for the calculation of the sum ofn addends by means ofn–1 floating-point additions is considered and in case of positive addends it is shown how to use the Huffman algorithm to choose a scheme to obtain the least upper bound on the accumulated roundoff error in the result.     

Quasi-Minimal Residual Eigenpairs

We consider the extraction of approximate eigenpairs from a given Krylov subspace that are optimal, or at least quasi–optimal with respect to the normwise backward error. An algorithm based on Newton–Grassmannian optimization and a simple example are given; a useful graphical method of comparison with other extraction methods is briefly sketched. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Error reduction for higher derivatives of Chebyshev collocation method using preconditioning and domain decomposition

A new preconditioning method is investigated to reduce the roundoff error in computing derivatives using Chebyshev collocation methods(CCM). Using this preconditioning causes ratio of roundoff error of preconditioning method and CCM becomes small when N gets large. Also for accuracy enhancement of differentiation we use a domain decomposition approach. Error analysis shows that for this domain decomposition method error reduces proportional to the length of subintervals. Numerical results show that using domain decomposition and preconditioning simultaneously, gives super accurate approximate values for first derivative of the function and good approximate values for moderately high derivatives.     

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.