Roundoff error analysis of algorithms based on Krylov subspace methods

We study the roundoff error propagation in an algorithm which computes the orthonormal basis of a Krylov subspace with Householder orthonormal matrices. Moreover, we analyze special implementations of the classical GMRES algorithm, and of the Full Orthogonalization Method. These techniques approximate the solution of a large sparse linear system of equations on a sequence of Krylov subspaces of small dimension. The roundoff error analyses show upper bounds for the error affecting the computed approximated solutions.This work was carried out with the financial contribution of the Human Capital and Mobility Programme of the European Union grant ERB4050PL921378.     

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.     

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.     

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.     

Adaptive solution of elliptic PDE-eigenvalue problems.

In this paper we introduce a new adaptive algorithm (AFEMLA) for elliptic PDE-eigenvalue problems. In contrast to other approaches the algebraic eigenvalue problem does not have to be solved to full accuracy. We incorporate the iterative solution of the resulting finite dimensional algebraic eigenvalue problems in the adaptation process in order to balance the cost 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. (© 2009 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

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.     

Optimization of functions whose values are subject to small errors

In this paper we consider the minimization of a function whose values can only be obtained with an error. For the case when the error has certain statistical properties this problem has been investigated by Kiefer and Wolfowitz (1) and Kushner (2, 3). Kushner has shown that a certain class of algorithms converge to a stationary point with probability one. Here a different approach is used. The error is assumed to have an upper bound and it is shown that a stationary point can be obtained to within a certain accuracy, dependent on the magnitude of the error. Our results are related to works concerning roundoff errors for one dimensional optimization (4) and solution of nonlinear equations (5). The algorithm we use can be regarded as an extension of the methods used in (6), (8) and (9).Supported by: Institutet för tillämpad matematik (Sweden) and the National Science Foundation (USA) under grant MPS 72-04787-a02.     

Truncation and roundoff errors in three-point approximations of first and second derivatives

We use linear combinations of Taylor expansions to develop three-point finite difference expressions for the first and second derivative of a function at a given node. We derive analytical expressions for the truncation and roundoff errors associated with these finite difference formulae. Using these error expressions, we find optimal values for the stepsize and the distribution of the three points, relative to the given node. The latter are obtained assuming that the three points are equispaced. For the first derivative approximation, the distribution of the points relative to the given node is not symmetrical, while it is so for the second derivative approximation. We illustrate these results with a numerical example in which we compute upper bounds on the roundoff error.     

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.     

Superlinear Convergence of a Stabilized SQP Method to a Degenerate Solution

We describe a slight modification of the well-known sequential quadratic programming method for nonlinear programming that attains superlinear convergence to a primal-dual solution even when the Jacobian of the active constraints is rank deficient at the solution. We show that rapid convergence occurs even in the presence of the roundoff errors that are introduced when the algorithm is implemented in floating-point arithmetic.     

A NUMERICALLY STABLE BLOCK MODIFIED GRAM-SCHMIDT ALGORITHM FOR SOLVING STIFF WEIGHTED LEAST SQUARES PROBLEMS

Recently, Wei in proved that perturbed stiff weighted pseudoinverses and stiff weighted least squares problems are stable, if and only if the original and perturbed coefficient matrices A and A^- satisfy several row rank preservation conditions. According to these conditions, in this paper we show that in general, ordinary modified Gram-Schmidt with column pivoting is not numerically stable for solving the stiff weighted least squares problem. We then propose a row block modified Gram-Schmidt algorithm with column pivoting, and show that with appropriately chosen tolerance, this algorithm can correctly determine the numerical ranks of these row partitioned sub-matrices, and the computed QR factor R^- contains small roundoff error which is row stable. Several numerical experiments are also provided to compare the results of the ordinary Modified Gram-Schmidt algorithm with column pivoting and the row block Modified Gram-Schmidt algorithm with column pivoting.     

Numerical investigation into the dependence of the Allen–Cahn equation on the free energy

In this paper, we study a direct parallel-in-time (PinT) algorithm for first- and second-order time-dependent differential equations. We use a second-order boundary value method as the time integrator. Instead of solving the corresponding all-at-once system iteratively, we diagonalize the time discretization matrix B, which yields a direct parallel implementation across all time steps. A crucial issue of this methodology is how the condition number (denoted by Cond₂(V )) of the eigenvector matrix V of B behaves as n grows, where n is the number of time steps. A large condition number leads to large roundoff error in the diagonalization procedure, which could seriously pollute the numerical accuracy. Based on a novel connection between the characteristic equation and the Chebyshev polynomials, we present explicit formulas for V and V^− 1, by which we prove that Cond\(_{2}(V)=\mathcal {O}(n^{2})\). This implies that the diagonalization process is well-conditioned and the roundoff error only increases moderately as n grows, and thus, compared to other direct PinT algorithms, a much larger n can be used to yield satisfactory parallelism. A fast structure-exploiting algorithm is also designed for computing the spectral diagonalization of B. Numerical results on parallel machine are given to support our findings, where over 60 times speedup is achieved with 256 cores.

     

An algorithm for computing a Padé approximant with minimal degree denominator

In this paper, a new definition of a reduced Padé approximant and an algorithm for its computation are proposed. Our approach is based on the investigation of the kernel structure of the Toeplitz matrix. It is shown that the reduced Padé approximant always has nice properties which the classical Padé approximant possesses only in the normal case. The new algorithm allows us to avoid the appearance of Froissart doublets induced by computer roundoff in the non-normal Padé table.     

Spectral collocation method and Darvishi’s preconditionings to solve the generalized Burgers–Huxley equation

In this paper, a numerical solution of the generalized Burgers–Huxley equation is presented. This is the application of spectral collocation method. To reduce roundoff error in this method we use Darvishi’s preconditionings. The numerical results obtained by this method have been compared with the exact solution. It can be seen that they are in a good agreement with each other, because errors are very small and figures of exact and numerical solutions are very similar.     

Almost optimal convergence of the point vortex method for vortex sheets using numerical filtering

Standard numerical methods for the Birkhoff-Rott equation for a vortex sheet are unstable due to the amplification of roundoff error by the Kelvin-Helmholtz instability. A nonlinear filtering method was used by Krasny to eliminate this spurious growth of round-off error and accurately compute the Birkhoff-Rott solution essentially up to the time it becomes singular. In this paper convergence is proved for the discretized Birkhoff-Rott equation with Krasny filtering and simulated roundoff error. The convergence is proved for a time almost up to the singularity time of the continuous solution. The proof is in an analytic function class and uses a discrete form of the abstract Cauchy-Kowalewski theorem. In order for the proof to work almost up to the singularity time, the linear and nonlinear parts of the equation, as well as the effects of Krasny filtering, are precisely estimated. The technique of proof applies directly to other ill-posed problems such as Rayleigh-Taylor unstable interfaces in incompressible, inviscid, and irrotational fluids, as well as to Saffman-Taylor unstable interfaces in Hele-Shaw cells.

     

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.     

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.     

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.