A simple derivation of the Hansen-Bliek-Rohn-Ning-Kearfott

enclosure for linear interval equations

Arnold Neumaier

Institut für Mathematik Universität Wien Strudlhofgasse 4 A-1090 Wien Austria email: neum@cma.univie.ac.at WWW: http://solon.cma.univie.ac.at/˜ neum/

Abstract. Recently, Ning & Kearfott derived a formula for the interval enclosure of the solution set of linear systems of equations with uncertain data ranging in intervals, in the case when the coefficient matrix is an Hmatrix. The enclosure is optimal when the midpoint matrix is diagonal, and when the midpoint is the identity, it reduces to the optimal method for enclosing preconditioned systems found by Hansen and Bliek and simplified by Rohn. An elementary proof of this formula is given using only simple properties of H-matrices and Schur complements. The new proof gives additional insight into why the theorem is true. It is also shown how to preserve rigor in the enclosure when finite precision arithmetic is used.

revised version, November 1998

Keywords: linear interval equations, optimal enclosure, H-matrix, Schur complement, finite precision arithmetic, directed rounding

1991 MSC Classification: 65G10

1 Introduction

An interesting and unexpected result by Hansen [2] gave an explicit formula for the the interval hull of the solution set of a linear system of interval equations in the special case that the midpoint matrix is the identity. Since such matrices arise by preconditioning more general linear interval equations, this result is of considerable interest. The result was independently found by Bliek [1] in Chapter 4.4 of his unpublished Ph.D. thesis, but both Hansen’s and Bliek’s proofs were not completely rigorous. Rohn [5] gave an impeccable proof of the result and simplified the formula, reducing the amount of work needed to compute the interval hull. Recently, Ning & Kearfott [4] generalized Rohn’s result to a formula giving an enclosure of the solution set of linear interval equations whose coefficient matrix is an H-matrix, and showed that the enclosure is optimal when the midpoint matrix is diagonal. In this paper, we give a more direct and elementary proof of this formula, together with implementation details in finite precision arithmetic.

In the following, concepts and notation are as in Neumaier [3], except that interval quantities are bold face. In particular, inequalities, absolute values and operations inf, sup, max, min are interpreted componentwise, ⊓ ⊔S = [inf S, sup S] is the interval hull of a bounded set S ⊂ ℝ^{n}, and

Σ(A, b) = {x ∈ ℝ^{n} | Ax = b for some A ∈ A, b ∈ b}

denotes the solution set of a system of linear interval equations with coefficient matrix A and right hand side b. A square interval matrix A ∈ 𝕀ℝ^{n×n} is called an H-matrix iff there is a vector v > 0 with ⟨A⟩v > 0. Here ⟨A⟩ is the comparison matrix with entries

⟨A⟩_{ii} = ⟨A_{ii}⟩ = min{|α| | α ∈ A_{ii}},
⟨A⟩_{ik} = −|A_{ik}| = − max{|α| | α ∈ A_{ik}} for k ≠ i.

Among the many properties of H-matrices discussed in Neumaier [3, Section 3.6-3.7], we need the fact that every matrix A contained in an interval Hmatrix A is nonsingular, the comparison matrix ⟨A⟩ is an H-matrix, too (even an M-Matrix), and

|A^{−1}| ≤ ⟨A⟩^{−1}.  (1)

We write A^{(i)} for the submatrix of a square matrix A obtained by dropping the ith row A_{i:} and the ith column A_{:i}, and b^{(i)} for the subvector of a vector b obtained by dropping the ith component b_{i}. Similarly, A_{i:}^{(i)} and A_{:i}^{(i)} denote the ith row and ith column of A with the entry A_{ii} dropped.

2 The Ning-Kearfott Theorem

Ning & Kearfott [4] proved the following theorem by reducing it to Rohn’s version of Hansen’s result.

Theorem 2.1 Let A ∈ 𝕀ℝ^{n×n} be an H-matrix, b ∈ 𝕀ℝ^{n} a right hand side,

u = ⟨A⟩^{−1}|b|, d_{i} = (⟨A⟩^{−1})_{ii}  (2)

and

α_{i} = ⟨A_{ii}⟩ − 1/d_{i}, β_{i} = u_{i}/d_{i} − |b_{i}|.  (3)

Then ⊓ ⊔Σ(A, b) is contained in the vector x with components

x_{i} = \frac{b_{i} + [−β_{i}, β_{i}]}{A_{ii} + [−α_{i}, α_{i}]}.  (4)

Moreover, if the midpoint of A is diagonal then ⊓ ⊔Σ(A, b) = x.

The new proof is based on two simple observations.

Lemma 2.2 If A^{(i)} is nonsingular and Ax = b then

(A_{ii} − A_{i:}^{(i)}(A^{(i)})^{−1}A_{:i}^{(i)}) x_{i} = b_{i} − A_{i:}^{(i)}(A^{(i)})^{−1}b^{(i)}.  (5)

Proof. We separate the ith variable of x and the ith equation of Ax = b and find

A^{(i)}x^{(i)} + A_{:i}^{(i)}x_{i} = b^{(i)},  (6)
A_{i:}^{(i)}x^{(i)} + A_{ii}x_{i} = b_{i}.  (7)

Solving (6) for x^{(i)} gives x^{(i)} = (A^{(i)})^{−1}(b^{(i)} − A_{:i}^{(i)}x_{i}), and insertion into (7) and simplification gives (5). _{□} Lemma 2.3 Let A be an H-matrix. Then A^{(i)} is an H-matrix, too.

Proof. By definition, A is an H-matrix iff there is a vector v > 0 with ⟨A⟩v > 0. Dropping the ith row and separating the contribution of v_{i} gives ⟨A^{(i)}⟩v^{(i)} > |A_{:i}^{(i)}|v_{i} ≥ 0. Therefore ⟨A^{(i)}⟩v^{(i)} > 0, and A^{(i)} is an H-matrix. _{□}

Proof of Theorem 2.1. By Lemma 2.3, we may apply Lemma 2.2 to an arbitrary A ∈ A and b ∈ b, and find for the solution x of Ax = b the relation (5), i.e.,

(A_{ii} − α̃_{i})x_{i} = b_{i} − β̃_{i},  (8)

where

α̃_{i} = A_{i:}^{(i)}(A^{(i)})^{−1}A_{:i}^{(i)}, β̃_{i} = A_{i:}^{(i)}(A^{(i)})^{−1}b^{(i)}.  (9)

By (1) and the rules for absolute values we find

|α̃_{i}| ≤ |A_{i:}^{(i)}|⟨A^{(i)}⟩^{−1}|A_{:i}^{(i)}| =: α_{i},  (10)
|β̃_{i}| ≤ |A_{i:}^{(i)}|⟨A^{(i)}⟩^{−1}|b^{(i)}| =: β_{i}.  (11)

We first verify that these definitions for α_{i} and β_{i} agree with the expressions given in the theorem. If we apply Lemma 2.2 to the equation ⟨A⟩c = I_{:i} in place of Ax = b, where I_{:i} is the ith column of the identity matrix, we find

(⟨A⟩_{ii} − α_{i})c_{i} = 1.  (12)

Since c_{i} = (⟨A⟩^{−1})_{ii} = d_{i}, this shows that α_{i} has the value claimed in (3). And if we apply Lemma 2.2 to the equation ⟨A⟩u = |b| in place of Ax = b, we find

(⟨A⟩_{ii} − α_{i})u_{i} = |b_{i}| + β_{i}.  (13)

(The plus sign is due to a minus sign in the comparison matrix!) By solving (13) for β_{i} we see (using (12)) that β_{i} also has the values claimed in (3). To prove (4), we first note that relation (12) implies ⟨A⟩_{ii} > α_{i} since c_{i} = (⟨A⟩^{−1})_{ii} > 0. Therefore, the denominator in (4) cannot contain zero, and (4) follows from (8) – (11) since

x_{i} = \frac{b_{i} − β̃_{i}}{A_{ii} − α̃_{i}} ∈ \frac{b_{i} + [−β_{i}, β_{i}]}{A_{ii} + [−α_{i}, α_{i}]} = x_{i}.  (14)

For the proof of optimality we first note that we may multiply the rows of the linear systems by factors ±1 to make all mid b_{k} ≥ 0, and then multiply the columns of A and the corresponding variables by factors ±1 to make all mid A_{kk} ≥ 0. This choice of signs guarantees

sup b = |b|,  (15)
⟨A⟩_{kk} = \underline{A}_{kk} > 0 for all k,  (16)

since the diagonal elements of H-matrices are nonzero. If mid A is diagonal, A contains the matrix A := ⟨A⟩ and the matrices obtained from A by making the entries of the ith row and/or column nonnegative without changing their absolute value. This gives four choices for A ∈ A with

A^{(i)} = ⟨A^{(i)}⟩, A_{i:}^{(i)} = ±|A_{i:}^{(i)}|, A_{:i}^{(i)} = ±|A_{:i}^{(i)}|

for independent choices of the signs, and this remains valid no matter how we change A_{ii} ∈ A_{ii}. If we now choose b ∈ b such that b^{(i)} = |b^{(i)}| which is possible by (15), we see from (9)–(11) that α̃_{i} = ±α_{i}, β̃_{i} = ±β_{i} with arbitrary and independent choices of the sign. Since we may still give A_{ii}, b_{i} any values in A_{ii}, b_{i}, the extremes in (14) are attained for suitable A ∈ A, b ∈ b. _{□}

Remark. It is easy to see that the Ning-Kearfott enclosure is always at least as good as that provided by the Krawczyk inverse (Theorem 3.7.8 in [3]). However, if mid A is not diagonal, one often gets a further improvement by refining the Ning-Kearfott enclosure by a few Gauss-Seidel iterations. Since the latter only costs O(n^{2}) operations, it should be applied generically in applications where it is not apriori known that (apart from rounding errors) mid A is diagonal.

3 Finite Precision Calculations

In order that the enclosure given by the Ning-Kearfott Theorem 2.1 remains rigorous in finite precision arithmetic, care is needed in the calculation of the quantities defined there. In particular, we need a rigorous upper bound B for ⟨A⟩^{−1}, in order to obtain rigorous upper bounds for α_{i} and β_{i} that render the subsequent interval calculations rigorous.

All platforms where interval arithmetic can be implemented have (hardware or emulated software) routines that set the rounding mode for subsequent operations to either upward rounding or downward rounding. For example, the Matlab toolbox Intlab of Rump [6] provides for this task two routines

SetRoundDown and SetRoundUp.

Because ⟨A⟩^{−1} is nonnegative by (1), a tight upper bound B for ⟨A⟩^{−1} can be constructed from an approximation B̃ to ⟨A⟩^{−1} and vectors v, w satisfying

I − ⟨A⟩B̃ ≤ ⟨A⟩vw^{T}  (17)

by taking

B = B̃ + vw^{T} .

(Proof: multiply by ⟨A⟩^{−1} ≥ 0.) Since A is an H-matrix, there is a vector v > 0 such that u = ⟨A⟩v > 0, and for such a vector, we can satisfy (17) with the vector w with tiny components

w_{k} = max_{i} \frac{−R_{ik}}{u_{i}} ,

where

R = ⟨A⟩B̃ − I

is the residual in the computation of the approximate inverse. To find an explicit such v, we may assume that a positive vector e is available whose components are proportional to natural scales for the components of x. Then v := B̃e ≈ ⟨A⟩^{−1}e > 0 if we indeed have an H-matrix, and then u = ⟨A⟩v ≈ e is positive if the approximation is good enough. Since conversely, the condition u = ⟨A⟩v > 0 implies that A is an H-matrix only when v ≥ 0, we take in the program the absolute value of the computed v to ensure the H-matrix property. In the absence of information about the relative magnitudes of the components of x we may simply choose the all-one vector for e, and this is the choice used in the program below. If we do the above calculations in finite precision arithmetic, we may calculate B̃ and v with arbitrary rounding, but u and R must underestimate the exact quantities, while w, B and the quantities calculated in (2) and (3) must overestimate the exact quantities. This is easily achieved by a proper setting of the rounding modes.

By looking at the operations involved, it is now easy to verify that the Matlab program

n = dim(A);                                  % dimension of A
dA = diag(A);    % save diagonal entries before overwriting A
A = compmat(A);        % A now contains the comparison matrix
B = inv(A);        % approximate inverse of comparison matrix
v = abs(B*ones(n,1));            % enforce nonnegativity of v
SetRoundDown;
u = A*v;           % is approximately equal to all-one vector
if ~all(min(u)>0),                    % check positivity of u
% A is numerically not an H-matrix
x = midrad(0,inf+zeros(n,1));
else
dAc = diag(A); % save diagonal entries before overwriting A
A = A*B-eye(n);        % A now contains the residual matrix
SetRoundUp;
w = zeros(1,n);
for i=1:n,
  w = max(w,(-A(i,:))/u(i));
end;
B = B+v*w;               % rigorous upper bound for exact B
u = B*abs(b);
d = diag(B);
alpha = dAc+(-1)./d;              % ensures upward rounding
beta = u./d-abs(b);
x = (b+midrad(0,beta))./(dA+midrad(0,alpha));
end;

(to be used together with the Intlab toolbox of Rump [6]) produces the desired rigorous bounds and hence a rigorous enclosure. It should not be difficult to rewrite this to equivalent code in other programming environments.

To calculate the Ning-Kearfott enclosure in the above way we need 2n^{3} + O(n^{2}) operations (counting additions and multiplications separately) for the approximate inversion and as many operations for the calculation of the residual, plus O(n^{2}) operations for the remaining calculations, giving a total of 4n^{3} + O(n^{2}) operations. Using a running error analysis for the inversion process, slightly suboptimal bounds for the residual can be obtained with O(n^{2}) operations, thus halving the work in high dimensions; see, e.g., Spellucci & Krier [7].

References

[1] C. Bliek, Computer methods for design automation, Ph.D. Thesis, Dept. of Ocean Engineering, Massachusetts Institute of Technology, 1992.

[2] E. Hansen, Bounding the solution of interval linear equations, SIAM J. Numer. Anal. 29 (1992), 1493-1503.

[3] A. Neumaier, Interval methods for systems of equations, Cambridge Univ. Press, Cambridge 1990.

[4] S. Ning and R. B. Kearfott, A comparison of some methods for solving linear interval equations, SIAM J. Numer. Anal. 34 (1997), 1289-1305.

[5] J. Rohn, Cheap and tight bounds: The recent result by E. Hansen can be made more efficient, Interval Computations 4 (1993), 13-21.

[6] S.M. Rump, INTLAB – INTerval LABoratory, Reliable Computing, to appear.

[7] P. Spellucci and N. Krier, Untersuchungen der Grenzgenauigkeit von Algorithmen zur Auflösung linearer Gleichungssysteme mit Fehlererfassung, pp. 288-297 in: Interval Mathematics (K. Nickel, ed.), Lecture Notes in Computer Science 29, Springer, Berlin 1975.