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.     

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.     

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.     

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.

     

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.

     

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.     

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.     

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.     

On some topological properties of numerical algorithms

A result quantity in a numerical algorithm is considered as a function of the input data, roundoff and truncation errors. In order to investigate this functional relationship using the methods of mathematical analysis a structural model of the numerical algorithm calledR-automaton is introduced. It is shown that the functional dependence defined by anR-automaton is a continuous rational function in a neighborhood of any data point except in a point set, the Lebesgue measure of which is zero. An effective general-purpose algorithm is presented to compute the derivative of any result quantity with respect to the individual roundoff and truncation errors. Some ways of generalizing theR-automation model without losing the results achieved are finally suggested.     

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.     

Interval arithmetic for guaranteed bounds in linear programming

In this note, we show how interval arithmetic can be used to give a solution to the linear programming problem which is guaranteed to be on the safe side of the true solution, where roundoff error is taken into account.This research was supported in part by the National Research Council of Canada.     

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.     

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.     

Multiplicative Error Analysis of Matrix Transformation Algorithms

The author's recently introduced relative error measure forvectors is applied to the error analysis of algorithms whichproceed by successive transformation of a matrix. Instead ofmodelling the roundoff errors at each stage by A: = T(A)+E onemodels them by A: =e^E T(A) where E is a small linear transformation.This can simplify analyses considerably. Applications to theparallel Jacobi method for eigenvalues, and to Gaussian elimination,are given.     

New algorithms for solving odes by pseudospectral method

To compute derivatives using matrix vector multiplication method, new algorithms were introduced in [1, 2]. By these algorithms, we reduced roundoff error in computing derivatives using Chebyshev collocation methods (CCM). In this paper, some applications of these algorithms are presented.     

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.     

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.     

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.

     

Implementing Taylor models arithmetic with floating-point arithmetic

The implementation of Taylor models arithmetic may use floating-point arithmetic to benefit from the speed of the floatingpoint implementation. The issue is then to take into account the roundoff errors. Here, we assume that the floating-point arithmetic is compliant with the IEEE-754 standard. We show how to get tight bounds of the roundoff errors, and more generally how to get high accuracy for the coefficients as well as for the bounds on the roundoff errors     

An Algorithm for image removals and decompositions without inverse matrices

Partial Differential Equation (PDE) based methods in image processing have been actively studied in the past few years. One of the effective methods is the method based on a total variation introduced by Rudin, Oshera and Fatemi (ROF) [L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268]. This method is a well known edge preserving model and an useful tool for image removals and decompositions. Unfortunately, this method has a nonlinear term in the equation which may yield an inaccurate numerical solution. To overcome the nonlinearity, a fixed point iteration method has been widely used. The nonlinear system based on the total variation is induced from the ROF model and the fixed point iteration method to solve the ROF model is introduced by Dobson and Vogel [D.C. Dobson, C.R. Vogel, Convergence of an iterative method for total variation denoising, SIAM J. Numer. Anal. 34 (5) (1997) 1779–1791]. However, some methods had to compute inverse matrices which led to roundoff error. To address this problem, we developed an efficient method for solving the ROF model. We make a sequence like Richardson’s method by using a fixed point iteration to evade the nonlinear equation. This approach does not require the computation of inverse matrices. The main idea is to make a direction vector for reducing the error at each iteration step. In other words, we make the next iteration to reduce the error from the computed error and the direction vector.