Forward Error Analysis Research Articles

Accurate Evaluation of Polynomials in Legendre Basis

This paper presents a compensated algorithm for accurate evaluation of a polynomial in Legendre basis. Since the coefficients of the evaluated polynomial are fractions, we propose to store these coefficients in two floating point numbers, such as double-double format, to reduce the effect of the coefficients’ perturbation. The proposed algorithm is obtained by applying error-free transformation to improve the Clenshaw algorithm. It can yield a full working precision accuracy for the ill-conditioned polynomial evaluation. Forward error analysis and numerical experiments illustrate the accuracy and efficiency of the algorithm.

Robust explicit MPC design under finite precision arithmetic

We propose a design methodology for explicit Model Predictive Control (MPC) that guarantees hard constraint satisfaction in the presence of finite precision arithmetic errors. The implementation of complex digital control techniques, like MPC, is becoming increasingly adopted in embedded systems, where reduced precision computation techniques are embraced to achieve fast execution and low power consumption. However, in a low precision implementation, constraint satisfaction is not guaranteed if infinite precision is assumed during the algorithm design. To enforce constraint satisfaction under numerical errors, we use forward error analysis to compute an error bound on the output of the embedded controller. We treat this error as a state disturbance and use this to inform the design of a constraint-tightening robust controller. Benchmarks with a classical control problem, namely an inverted pendulum, show how it is possible to guarantee, by design, constraint satisfaction for embedded systems featuring low precision, fixed-point computations.

Control-Theoretic Forward Error Analysis of Iterative Numerical Algorithms

It has been known for at least five decades that control theory can be used to study iterative algorithms. However, little work can be found in the control systems literature on numerical algorithms, especially on the study of finite precision effects. In this technical note, we consider numerical iterative algorithms in finite precision as dynamical systems and study the effects of finite precision using control theory. By using the control tools of input-to-state stability and results from the study of quantization in control systems, we present new systematic ways to find bounds on the forward error for iterative algorithms. The advantages of the proposed schemes are shown by applying them to find bounds for the classical iterative methods for solving a system of linear equations.

Efficient computation of the matrix cosine

Trigonometric matrix functions play a fundamental role in second order differential equation systems. This work presents an algorithm for computing the cosine matrix function based on Taylor series and the cosine double angle formula. It uses a forward absolute error analysis providing sharper bounds than existing methods. The proposed algorithm had lower cost than state-of-the-art algorithms based on Hermite matrix polynomial series and Padé approximants with higher accuracy in the majority of test matrices.

Accurate evaluation of the [formula omitted]-th derivative of a polynomial and its application

This paper presents a compensated algorithm for the evaluation of the k-th derivative of a polynomial in power basis. The proposed algorithm makes it possible the direct evaluation without obtaining the k-th derivative expression of the polynomial itself, with a very accurate result to all but the most ill-conditioned evaluation. Forward error analysis and running error analysis are performed by an approach based on the data dependency graph. Numerical experiments illustrate the accuracy and efficiency of the algorithm.

Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm

A popular method for computing the matrix logarithm is the inverse scaling and squaring method, which essentially carries out the steps of the scaling and squaring method for the matrix exponential in reverse order. Here we make several improvements to the method, putting its development on a par with our recent version [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970--989] of the scaling and squaring method for the exponential. In particular, we introduce backward error analysis to replace the previous forward error analysis; obtain backward error bounds in terms of the quantities $\|A^p\|^{1/p}$, for several small integer $p$, instead of $\|A\|$; and use special techniques to compute the argument of the Padé approximant more accurately. We derive one algorithm that employs a Schur decomposition, and thereby works with triangular matrices, and another that requires only matrix multiplications and the solution of multiple right-hand side linear systems. Numerical experiments show the new algorithms to be generally faster and more accurate than their existing counterparts and suggest that the Schur-based method is the method of choice for computing the matrix logarithm.

Validated computation of certain hypergeometric functions

We present an efficient algorithm for the validated high-precision computation of real continued fractions, accurate to the last digit. The algorithm proceeds in two stages. In the first stage, computations are done in double precision. A forward error analysis and some heuristics are used to obtain an a priori error estimate. This estimate is used in the second stage to compute the fraction to the requested accuracy in high precision (adaptively incrementing the precision for reasons of efficiency). A running error analysis and techniques from interval arithmetic are used to validate the result. As an application, we use this algorithm to compute Gauss and confluent hypergeometric functions when one of the numerator parameters is a positive integer.

On the numerical stability of algorithmic differentiation

In contrast to integration, the differentiation of a function is an ill-conditioned process, if only an oracle is available for its pointwise evaluation. That is, unrelated small variations in the value of the composite function are allowed at nearly identical arguments. In contrast, we show here that, if the function is defined by an evaluation procedure as a composition of arithmetic operations and elementary functions, then automatic, or algorithmic differentiation is backward stable in the sense of Wilkinson. More specifically, the derivative values obtained are exact for a perturbation of the elementary components at the level of the machine precision. We also provide a forward error analysis for both the forward and reverse mode. The theoretical analysis is confirmed by numerical experiments.

Accurate evaluation of a polynomial in Chebyshev form

This paper presents a compensated algorithm to accurately evaluate a polynomial expressed in Chebyshev basis of the first and second kind with floating-point coefficients. The principle is to apply error-free transformations to improve the traditional Clenshaw algorithm. The new algorithm is as accurate as the Clenshaw algorithm performed in twice the working precision. Forward error analysis and numerical experiments illustrate the accuracy and properties of the proposed algorithm.

Accurate matrix exponential computation to solve coupled differential models in engineering

The matrix exponential plays a fundamental role in linear systems arising in engineering, mechanics and control theory. This work presents a new scaling-squaring algorithm for matrix exponential computation. It uses forward and backward error analysis with improved bounds for normal and nonnormal matrices. Applied to the Taylor method, it has presented a lower or similar cost compared to the state-of-the-art Padé algorithms with better accuracy results in the majority of test matrices, avoiding Padé’s denominator condition problems.

Gerschgorin’s theorem for generalized eigenvalue problems in the Euclidean metric

We present Gerschgorin-type eigenvalue inclusion sets applicable to generalized eigenvalue problems. Our sets are defined by circles in the complex plane in the standard Euclidean metric, and are easier to compute than known similar results. As one application we use our results to provide a forward error analysis for a computed eigenvalue of a diagonalizable pencil.

Accurate evaluation of a polynomial and its derivative in Bernstein form

In computer aided geometric design a polynomial is usually represented in Bernstein form. This paper presents two compensated algorithms to accurately evaluate a polynomial and its derivative in Bernstein form with floating point coefficients. The principle is to apply error-free transformation to improve traditional de Casteljau algorithm. Forward error analysis and numerical experiments illustrate the accuracy of our algorithms.

On direct elimination methods for solving the equality constrained least squares problem

Two closely related methods for solving the least squares problem with equality constraints (LSE) are considered. The first is the direct elimination (DE) method that is implemented using Modified Gram-Schmidt (MGS) for minimizing ‖Ax − c‖2 subject to the constraints ‖Bx − d‖2 = min. The second is the MGS for solving the following weighted least squares problem where we take the limit τ → ∞ analytically. The two methods are proved to be equivalent. We provide the backward and forward round-off error analysis for the DE method via MGS, and illustrate our results with numerical experiments.

Perturbation behavior of a multiple eigenvalue in generalized Hermitian eigenvalue problems

We show that a multiple eigenvalue has different sensitivities under perturbations in a generalized Hermitian eigenvalue problem. Our result provides a solution to a question raised by Stewart and Sun. We also show how this difference of sensitivities plays a role in the eigenvalue forward error analysis after the Rayleigh-Ritz process, for which we present an approach that provides tight bounds.

Certification of bounds on expressions involving rounded operators

Gappa is a tool designed to formally verify the correctness of numerical software and hardware. It uses interval arithmetic and forward error analysis to bound mathematical expressions that involve rounded as well as exact operators. It then generates a theorem and its proof for each verified enclosure. This proof can be automatically checked with a proof assistant, such as Coq or HOL Light. It relies on a large companion library of facts that we have developed. This Coq library provides theorems dealing with addition, multiplication, division, and square root, for both fixed- and floating-point arithmetics. Gappa uses multiple-precision dyadic fractions for the endpoints of intervals and performs forward error analysis on rounded operators when necessary. When asked, Gappa reports the best bounds it is able to reach for a given expression in a given context. This feature can be used to identify where the set of facts and automatic techniques implemented in Gappa becomes insufficient. Gappa handles seamlessly additional properties expressed as interval properties or rewriting rules in order to establish more intricate bounds. Recent work showed that Gappa is suited to discharge proof obligations generated for small pieces of software. They may be produced by third-party tools and the first applications of Gappa use proof obligations written by designers or obtained from traces of execution.

Fast matrix multiplication is stable

We perform forward error analysis for a large class of recursive matrix multiplication algorithms in the spirit of Bini and Lotti [Numer. Math. 36:63–72, 1980]. As a consequence of our analysis, we show that the exponent of matrix multiplication (the optimal running time) can be achieved by numerically stable algorithms. We also show that new group-theoretic algorithms proposed in Cohn and Umans [Foundations of Computer Science, 44th Annual IEEE Symposium, pp. 438–449, 2003] and Cohn et al. [Foundations of Computer Science, 46th Annual IEEE Symposium, pp. 379–388, 2005] are all included in the class of algorithms to which our analysis applies, and are therefore numerically stable. We perform detailed error analysis for three specific fast group-theoretic algorithms.

Evaluation algorithms for multivariate polynomials in Bernstein–Bézier form

The evaluation of multivariate polynomials of n variables in Bernstein–Bézier form is considered. A forward error analysis for the corresponding de Casteljau algorithm and the VS algorithm is performed. We also include algorithms that simultaneously evaluate the polynomial and provide “a posteriori” error bounds, without increasing significantly the computational cost. The sharpness of our running error bounds is shown in the case of trivariate polynomials.

Evaluation of the derivative of a polynomial in Bernstein form

In geometric modelling and computer-aided design a polynomial is usually represented in Bernstein form. Forward error analysis for the computation of the derivative of a polynomial in Bernstein form is performed. The corresponding running error is also analyzed. We have compared our results with those corresponding to a polynomial represented in power basis and using Legendre and Chebyshev forms. Since the conversions between the Bernstein and the power bases and between Bernstein and Chebyshev polynomials were already studied in the literature, for the sake of completeness the conversion between Bernstein and Chebyshev polynomials of the first kind is analyzed in an appendix.

The Orthogonal Rayleigh Quotient Iteration (ORQI) method

This paper presents a new method for computing all the eigenvectors of a real n× n symmetric band matrix T. The algorithm computes an orthogonal matrix Q=[ q 1,…, q n] and a diagonal matrix Λ=diag{ λ 1,…, λ n } such that TQ= QΛ. The basic ideas are rather simple. Assume that q 1,…, q k−1 and λ 1,…, λ k−1 have already been computed. Then q k is obtained via the Rayleigh Quotient Iteration (RQI) method. Starting from an arbitrary vector u 0 the RQI method generates a sequence of vectors u ℓ , ℓ=1,2,… , and a sequence of scalars ρ ℓ, ℓ=0,1,2,… The theory tells us that these two sequences converge (almost always) to an eigenpair (ρ *, u *) . The appeal of the RQI method comes from the observation that the final rate of convergence is cubic. Furthermore, if the starting point is forced to satisfy q T i u 0=0 for i=1,…, k−1, as our method does, then all the coming vectors, u ℓ,ℓ=1,2,… , and their limit point, u * , should stay orthogonal to q 1,…, q k−1 . In practice orthogonality is lost because of rounding errors. This difficulty is resolved by successive orthogonalization of u ℓ against q 1,…, q k−1 . The key for effective implementation of the algorithm is to use a selective orthogonalization scheme in which u ℓ is orthogonalized only against “close” eigenvectors. That is, u ℓ is orthogonalized against q i only if | ρ ℓ− λ i |⩽ γ where γ is a small threshold value, e.g., γ=∥ T∥ ∞/1000. An essential feature of the proposed orthogonalization scheme is the use of reorthogonalization. The ORQI method is supported by forward and backward error analysis. Preliminary experiments on medium-size problems ( n⩽1000) are quite encouraging. The average number of iterations per eigenvector was less than 13, while the overall number of flops required for orthogonalizations is often below n 3/2.

Quasiconvex relaxations based on interval arithmetic

Interval analysis provides a tool for (i) forward error analysis, (ii) estimating and controlling rounding and approximation errors automatically, and (iii) proving existence and uniqueness of solutions. In this context the terms self-validating methods, inclusion methods orverification methods are in use. In this paper, we present a new self-validating method for solving global constrained optimization problems. This method is based on the construction of quasiconvex lower bound and quasiconcave upper bound functions of a given function, the latter defined by an arithmetical expression. No further assumptions about the nonlinearities of the given function are necessary. These lower and upper bound functions are rigorous by using the tools of interval arithmetic. In its easiest form they are constructed by taking appropriate linear and/or quadratical estimators which yield quasiconvex/quasiconcave bound functions. We show how these bound functions can be used to define rigorous quasiconvex relaxations for constrained global optimization problems and nonlinear systems. These relaxations can be incorporated in a branch and bound framework yielding a self-validating method.