Improved Backward Error Bounds for LU and Cholesky Factorizations

Abstract

Assuming standard floating-point arithmetic (in base β, precision p) and barring underflow and overflow, classical rounding error analysis of the LU or Cholesky factorization of an n × n matrix A provides backward error bounds of the form |∆A| γn| L|| U | or |∆A| γ n+1 | R T || R|. Here, L, U , and R denote the computed factors, and γn is the usual fraction nu/(1−nu) = nu+O(u 2) with u the unit roundoff. Similarly, when solving an n × n triangular system T x = b by substitution, the computed solution x satisfies (T + ∆T) x = b with |∆T | γn|T |. All these error bounds contain quadratic terms in u and limit n to satisfy either nu < 1 or (n + 1)u < 1. We show in this paper that the constants γn and γ n+1 can be replaced by nu and (n + 1)u, respectively, and that the restrictions on n can be removed. To get these new bounds the main ingredient is a general framework for bounding expressions of the form |ρ − s|, where s is the exact sum of a floating-point number and n − 1 real numbers, and where ρ is a real number approximating the computed sum s. By instantiating this framework with suitable values of ρ, we obtain improved versions of the well-known Lemma 8.4 from [N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002] (used for analyzing triangular system solving and LU factorization) and of its Cholesky variant. All our results hold for rounding to nearest with any tie-breaking strategy and no matter what the order of summation.

Improved error bounds for floating-point products and Horner’s scheme

“…without restriction on the integer k and, in [9], similar improvements have been obtained for the residuals of the computed LU and Cholesky factors as well as for triangular system solutions. A similar result was recently shown by Graillat, Lefèvre, and Muller [1] for binary arithmetic: This improves the classical Wilkinson-type estimate | r − x k+1 | γ k |x k+1 |.…”

“…A similar result was recently shown by Graillat, Lefèvre, and Muller [1] for binary arithmetic: This improves the classical Wilkinson-type estimate | r − x k+1 | γ k |x k+1 |. They also note that for k ≈ u −1 the relative error on r can indeed be larger than ku, thus suggesting that in the case of integer powers, the price to be paid for the refined constant ku is a necessary restriction on the range of k. This is in contrast with bounds like (1.2) and the results in [4,9], where restrictions on k can be avoided.…”

Error Bounds for Computer Arithmetics

“…(12) Finally, it was shown in [36] that the concept of linearizing bounds by replacing γ k by ku is also true for some standard numerical linear algebra algorithms. If for some A P F mˆn with m ě n Gaussian elimination runs to completion, then the computed factorsL andÛ satisfy (comparison and absolute value to be understood entrywise)…”

“…By means of explicit examples (cf. [24]) it is easy to see that some restriction on n is mandatory for (36). But although it was common belief that this is the worst case, it could not be proved.…”

Exploiting Structure in Floating-Point Arithmetic

“…mul. nu + O(u 2 ) nu [24] Euclidean norm ( n 2 + 1)u + O(u 2 ) ( n 2 + 1)u [25] T x = b, A = LU nu + O(u 2 ) nu [39] A = R T R (n + 1)u + O(u 2 ) (n + 1)u [39] x [38] poly. eval.…”