Verified Error Bounds for Sparse Systems Part I: The Splitting of a Matrix into Two Factors

Abstract. Verification methods provide mathematically correct error bounds for the solution of a numerical problem. That includes the proof of solvability of the problem and often uniqueness of the solution within the computed bounds. There are many verification methods for standard problems in numerical analysis, including linear and nonlinear systems of equations, matrix decompositions, eigenproblems, local and global optimization, and ordinary and partial differential equations. Many of those are included in INTLAB, the MATLAB/Octave toolbox for reliable computing. Despite several efforts, the verified solution of general sparse linear systems was an open problem. There are satisfactory algorithms for systems with symmetric positive definite input matrix. To that end error bounds for the solution of [Formula: see text] with general matrix [Formula: see text] could be computed using [Formula: see text], but that reduces the applicability in double precision to matrices with condition number up to [Formula: see text]. We give in this note an algorithm to compute entrywise error bounds for the solution of general real or complex sparse systems with condition number up to the limit [Formula: see text]. Our algorithm splits into three subalgorithms for symmetric positive definite, symmetric indefinite, and general input matrix [Formula: see text]. It is based on a mathematically correct lower bound on the smallest singular value [Formula: see text]. A key point is a factorization [Formula: see text] such that [Formula: see text] and [Formula: see text] have identical sets of singular values with the smallest one close to [Formula: see text]. A mathematically correct lower bound on [Formula: see text] is then computed using [Formula: see text]. Numerical evidence suggests that bounds for the solution of a linear system are computed for condition numbers up to [Formula: see text], and that often the bounds for all entries are close to maximally accurate, i.e., the bounds differ by few bits. We present a new a priori error bound based on Perron–Frobenius theory on the residual of a floating-point Cholesky decomposition which improves on existing ones by some two orders of magnitude. Our three subalgorithms benefit from this new bound. Based on the results in this Part I, an alternative approach will be presented in Part II of this note. Those methods are somewhat simpler, but often slower. However, they seem to be more robust, i.e., produce verified inclusions where the methods of Part I fail. Both approaches for square linear systems will be used in Part II of this note to compute verified error bounds for the solution of least squares problems and for underdetermined linear systems. Moreover, inclusions of the solution of general real or complex systems of nonlinear equations with sparse Jacobi matrix are computed by transforming the problem into a linear system with point matrix and interval right-hand side.