A COMPARISON OF SOME METHODS FOR SOLVING LINEAR INTERVAL EQUATIONS^{∗}

S. NING AND R. B. KEARFOTT^{†}

Abstract. Certain cases in which the interval hull of a system of linear interval equations can be computed inexpensively are outlined. We extend a proposed technique of Hansen and Rohn with a formula that bounds the solution set of a system of equations whose coefficient matrix A = [\underline{A}, \overline{A}] is an H-matrix; when A is centered about a diagonal matrix, these bounds are the smallest possible (i.e. the bounds are then the solution hull). Hansen’s scheme also computes the solution hull when the linear interval system Ax = b = [\underline{b}, \overline{b}] is such that A is inverse positive and \underline{b} = −\overline{b} ≠ 0. Earlier results of others also imply that, when A is an M-matrix and b ≥ 0, b ≤ 0 or 0 ∈ b, interval Gaussian elimination (IGA) computes the hull. We also give a method of computing the solution hull inexpensively in many instances when A is inverse positive, given an outer approximation such as that obtained from IGA. Examples are used to compare these schemes under various conditions.

Key words. numerical linear algebra, interval computations, inverse-positive matrix, H-matrix, interval hull, bounding solution sets

AMS subject classifications. 15A06, 15A23, 65F05, 65G10, 65F35

1. Introduction. We assume familiarity with the elementary properties of interval arithmetic; for an introduction, see [1, 5, 11, 14, 15]. Also, we will not attempt to survey results on computation of the hull of the solution set of linear interval systems; the reader may wish to consult [14, Ch. 6] or Rohn’s recent review [17]. We will, however, introduce our particular notation and some background facts used in this article.

We will use boldface to denote intervals. We will use underscores to denote lower bounds of intervals and overscores to denote upper bounds of intervals. Subscripts will be used to denote components of vectors. For example, we may have:

x = (x_{1}, x_{2}, …, x_{n})^{T} ,

where x_{i} = [\underline{x}_{i}, \overline{x}_{i}]. In general, vectors will be denoted by lower case letters, while matrices will be denoted by upper case. For example, A = \{[\underline{A}_{i,j}, \overline{A}_{i,j}]\}_{i,j=1}^{n} will denote an n by n interval matrix. We use 𝕀ℝ^{n} to denote the set of all n dimensional interval vectors and 𝕀ℝ^{n×n} to denote the set of all n × n interval matrices. The quantity ⟨x⟩ will denote the mignitude of an interval, defined by ⟨x⟩ = min {|x| : x ∈ x}. Similarly, the magnitude of an interval is defined as |x| = max {|x| : x ∈ x}.

The radius of an interval x will be denoted by rad(x) = (\overline{x} − \underline{x})/2, and the radius of an interval vector x or matrix A, denoted rad(x) or rad(A), will be defined componentwise.

If v and w are n-dimensional interval vectors, then v > w means that the inequality is valid componentwise, i.e. \underline{v}_{i} > \overline{w}_{i} for i between 1 and n. Similar interpretation is given to v ≥ w, v < w and v ≤ w. Inequalities of interval matrices M > N, etc. are defined similarly.

If A ∈ 𝕀ℝ^{n×n} is a regular interval matrix, we define a matrix inverse of A by

A^{−1} = □{A^{−1} : A ∈ A}. (1)

∗ This work was supported in part by National Science Foundation grant CCR-9203730.

† Department of Mathematics, University of Southwestern Louisiana, U.S.L. Box 43026, Lafayette, LA 70504 USA. email: sxn0893@usl.edu, rbk@usl.edu

That is, A^{−1} is the smallest interval matrix that contains the set {A^{−1} : A ∈ A}. We call an interval matrix A inverse positive if A is regular and A^{−1} ≥ 0. An n × n interval matrix A is called an M-matrix iff A_{ij} ≤ 0 for all i ≠ j, and

Au > 0 for some positive vector u ∈ ℝ^{n}.

We use ⟨A⟩ to denote the comparison matrix of A, with entries

⟨A⟩_{ii} = ⟨A_{ii}⟩, ⟨A⟩_{ij} = −|A_{ij}| for i ≠ j.

We call A an H-matrix iff its comparison matrix ⟨A⟩ is an M-matrix.

Now we state [14, Proposition 3.6.6]:

Proposition 1.1. (Kuttler [10]) Let A = [\underline{A}, \overline{A}] ∈ 𝕀ℝ^{n×n}. If \underline{A} and \overline{A} are regular and \underline{A}^{−1} ≥ 0, \overline{A}^{−1} ≥ 0, then A is regular and

A^{−1} = [\overline{A}^{−1}, \underline{A}^{−1}] ≥ 0.

It is known that every M-matrix is inverse positive. Thus, Proposition 1.1 can be applied to M-matrices.

Suppose that A ∈ 𝕀ℝ^{n×n} is a regular interval matrix and b ∈ 𝕀ℝ^{n}. We refer to the set

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

as the solution set of the interval linear system

Ax = b. (3)

In Neumaier [14] and elsewhere, it is shown that the hull A^{H}b of (2), that is, the smallest interval vector that contains Σ(A, b), obeys

A^{H}b = □{x_{D} : |D| = I, inf(D(Ax_{D} − b)) = 0}. (4)

We will speak of the quantity A^{H}b as the hull of the interval linear system (3), or the solution hull of the interval linear system. In general, computation of the hull (4) requires solution of the 2^{n} equations

inf(D(Ax_{D} − b)) = 0, (5)

where D ranges over all D with |D| = I. In many special cases, however, not all 2^{n} vectors x_{D}, (|D| = I) are needed. We give a method to compute the exact hull of Equation (3) when A is centered about a diagonal matrix. We also consider how to compute the hull of Equation (3) when A is inverse positive. Our results extend and amplify results of Neumaier, Hansen and Rohn.

We state [14, Theorem 3.6.7] as follows.

Theorem 1.2. (Beeck [3]) Let A ∈ 𝕀ℝ^{n×n} be inverse positive. Then

A^{H}b = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = [\underline{x}, \overline{x}]

where A^{(1)}, A^{(2)} ∈ A are defined by

A_{ik}^{(1)} = \overline{A}_{ik} if \underline{x}_{k} ≥ 0 and A_{ik}^{(1)} = \underline{A}_{ik} otherwise,

A_{ik}^{(2)} = \overline{A}_{ik} if \overline{x}_{k} ≤ 0 and A_{ik}^{(2)} = \underline{A}_{ik} otherwise.

In particular,

A^{H}b = A^{−1}b = \{ \begin{matrix}[\overline{A}^{−1}\underline{b}, \underline{A}^{−1}\overline{b}] & if b ≥ 0, \\ [\underline{A}^{−1}\underline{b}, \underline{A}^{−1}\overline{b}] & if 0 ∈ b, \\ [\underline{A}^{−1}\underline{b}, \overline{A}^{−1}\overline{b}] & if b ≤ 0. \end{matrix}

In [6], E. R. Hansen proposed a scheme to estimate the hull of a preconditioned linear interval system. Suppose that

A = Ǎ + [−1, 1]rad(A) ∈ 𝕀ℝ^{n×n}

and b = [\underline{b}, \overline{b}] ∈ 𝕀ℝ^{n}. Consider the equation

Ax = b. (6)

Multiplying on the left by Ǎ^{−1}, we obtain the preconditioned equation

Mx = r, (7)

where r = Ǎ^{−1}b and M = Ǎ^{−1}A = I + [−1, 1]|Ǎ^{−1}| rad(A) = I + [−1, 1]Q is centered about an identity matrix I, where Q = |Ǎ^{−1}|rad(A) ≥ 0. As proposed in [6], the scheme required solution of 2n systems of linear equations sharing the same coefficient matrix. In [16], Rohn proposed and improvement that requires only a single floating point matrix inversion^{1}.

To unify the presentation, we present versions of Hansen’s and Rohn’s formulas with modified hypotheses. In particular, in Hansen’s original theorem, it was assumed that M was regular; this was assured by assuming M = [\underline{M}, \overline{M}] = [I − Q, I + Q] was diagonally dominant. In Rohn’s theorem, it was assumed that the spectral radius ρ(Q) < 1. In contrast, in both results, we assume that \underline{M}^{−1} ≥ 0. The condition \underline{M}^{−1} ≥ 0 implies ρ(Q) < 1, which in turn implies M is regular: Since Q ≥ 0, the Perron–Frobenius theorem ([14, Theorem 3.2.2]) implies that there exists a real nonnegative eigenvector x ≠ 0 with Qx = ρ(Q)x. Since \underline{M}x = (I −Q)x = (1−ρ(Q))x, x = (1 − ρ(Q))\underline{M}^{−1}x. From this, since x ≥ 0, x ≠ 0 and \underline{M}^{−1}x ≥ 0, 1 − ρ(Q) > 0, so ρ(Q) < 1. Combining this with [14, Proposition 4.1.1] then implies M is strongly regular, and a fortiori, regular.

Thus restated, Hansen’s result is

Theorem 1.3. (Hansen [6]) Suppose M = [\underline{M}, \overline{M}] of Equation (7) satisfies

\underline{M}^{−1} ≥ 0. (8)

Let

(9) s^{(i)} = \{ \begin{matrix}\overline{r}_{i} & for j = i \\ max\{−\underline{r}_{j}, \overline{r}_{j}\} & for j ≠ i; j = 1, …, n \end{matrix}

(10) t^{(i)} = \{ \begin{matrix}\underline{r}_{i} & for j = i \\ min\{\underline{r}_{j}, −\overline{r}_{j}\} & for j ≠ i; j = 1, …, n \end{matrix}

and

c_{i} = \frac{1}{2(\underline{M}^{−1})_{ii} − 1}. (11)

1 with appropriate rounding-out to maintain rigor

Then the hull of Equation (7) is

M^{H}r = [\underline{x}, \overline{x}] (12)

with

\underline{x}_{i} = \{ \begin{matrix}c_{i}e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} ≥ 0 \\ e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} < 0 \end{matrix} (13)

\overline{x}_{i} = \{ \begin{matrix}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} ≥ 0 \\ c_{i}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} < 0 \end{matrix} (14)

for i = 1, 2, …, n, where e_{i}^{T} is the unit vector whose i-th coordinate is 1 and all of whose other coordinates are 0.

Rohn’s result, with our hypothesis, is

Theorem 1.4. (Rohn [16]) As in Theorem 1.3, assume \underline{M} is inverse positive. Then the hull of Equation (7) is M^{H}r = [\underline{x}, \overline{x}], where

\underline{x}_{i} = min\{x̃_{i}, c_{i}x̃_{i}\},

\overline{x}_{i} = max\{x̃_{i}, c_{i}x̃_{i}\},

where

x̃_{i} = −x_{i}^{∗} + (\underline{M}^{−1})_{ii}(ř + |ř|)_{i}
x̃_{i} = x_{i}^{∗} + (\underline{M}^{−1})_{ii}(ř − |ř|)_{i}
x_{i}^{∗} = (\underline{M}^{−1}(|ř| + (\overline{r} − \underline{r})/2))_{i},
c_{i} = \frac{1}{2(\underline{M}^{−1})_{ii} − 1} ∈ (0, 1],

and ř, the midpoint vector, is that vector whose i-th component is (\underline{r}_{i} + \overline{r}_{i})/2.

In [16], Rohn also proves some interesting comparisons when the scheme is used to solve linear systems or invert interval matrices.

2. The theory. We extend the aforementioned techniques and results here. Theorem 2.1. Suppose that A ∈ 𝕀ℝ^{n×n} and b ∈ 𝕀ℝ^{n}. Assume D is a real

diagonal matrix with D_{ii} ≠ 0 for i = 1, 2, …, n, M = D^{−1}A, and r = D^{−1}b. Then the solution set of the system

Ax = b (15)

and the solution set of the preconditioned system

Mx = r (16)

are the same, i.e.

Σ(A, b) = Σ(M, r). (17)

Proof. The proof rests on the fact that subdistributivity does not occur, since D is diagonal. By the assumption on D, we have

D = ( \begin{matrix}D_{11} \\ & D_{22} \\ & & . & . & . \\ & & & & & D_{nn} \end{matrix} ) .

Thus,

D^{−1} = ( \begin{matrix}D_{11} & −1 \\ & & D_{22} & −1 \\ & & & & . & . & . \\ & & & & & & & D_{nn} & −1 \end{matrix} ) .

For any x ∈ Σ(A, b), there exists an A ∈ A and a b ∈ b such that Ax = b. (18) Multiplying both sides of the system (18) on the left by D^{−1}, we obtain the system

D^{−1}Ax = D^{−1}b.

Since D^{−1}A ∈ M and D^{−1}b ∈ r, we have x ∈ Σ(M, r). Hence Σ(A, b) ⊆ Σ(M, r). (19) Similarly for any x ∈ Σ(M, r), there exists an M ∈ M and a r ∈ r such that Mx = r. (20) Multiplying both sides of the system (20) on the left by D, we obtain the system

DMx = Dr.

Since D and D^{−1} are diagonal real matrices, we have

M_{ij} = \{ \begin{matrix}[D_{ii}^{−1}\overline{A}_{ij}, D_{ii}^{−1}\underline{A}_{ij}] & for D_{ii} < 0 \\ [D_{ii}^{−1}\underline{A}_{ij}, D_{ii}^{−1}\overline{A}_{ij}] & for D_{ii} > 0 \end{matrix} (21)

and

r_{j} = \{ \begin{matrix}[D_{ii}^{−1}\overline{b}_{j}, D_{ii}^{−1}\underline{b}_{j}] & for D_{ii} < 0 \\ [D_{ii}^{−1}\underline{b}_{j}, D_{ii}^{−1}\overline{b}_{j}] & for D_{ii} > 0 \end{matrix} (22)

for i, j = 1, 2, …, n. Multiplying the above two equations by D_{ii}, we obtain

D_{ii}M_{ij} = [\underline{A}_{ij}, \overline{A}_{ij}] = A_{ij}

and

D_{ii}r_{j} = [\underline{b}_{j}, \overline{b}_{j}] = b_{j}

6 S. Ning and R. B. Kearfott

for i, j = 1, 2, …, n. This means that

DM = A

and

Dr = b.

Hence DM ∈ A and Dr ∈ b. This implies that x ∈ Σ(A, b). Thus we have Σ(M, r) ⊆ Σ(A, b) (23) Combining (19) and (23), we have

Σ(A, b) = Σ(M, r).

Theorem 2.2. Let A_{n×n} be an H-matrix, b a right hand side,

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

and

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

for i = 1, 2, …, n. Then A^{H}b is contained in the vector x with components

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

for i = 1, 2, …, n. Moreover, if the midpoint matrix of A is diagonal, then A^{H}b = x. Proof. Since A is an H-matrix, 0 ∉ A_{ii}. First assume A_{ii} > 0. Let

D = ( Ǎ_{11} ; Ǎ_{22} ; ... ; Ǎ_{nn} ) , (25)

where Ǎ_{ii} = (\underline{A}_{ii} + \overline{A}_{ii})/2, and

(26) P = ( \begin{matrix}rad(A_{11}) & |A_{12}| & |A_{13}| & · & · & · & |A_{1n}| \\ |A_{21}| & rad(A_{22}) & |A_{23}| & · & · & · & |A_{2n}| \\ |A_{31}| & |A_{32}| & rad(A_{33}) & · & · & · & |A_{3n}| \\ · & · & · & · & & & · \\ · & · & · & & · & & · \\ · & · & · & & & · & · \\ |A_{n1}| & |A_{n2}| & |A_{n3}| & · & · & · & rad(A_{nn}) \end{matrix} ) .

Then A ⊆ D + [−1, 1]P. (27) Hence Σ(A, b) ⊆ Σ(D + [−1, 1]P , b). (28)

Solving Linear Interval Equations 7

Letting M = I + [−1, 1]D^{−1}P and r = D^{−1}b, it follows from Theorem 2.1 that Σ(D + [−1, 1]P , b) = Σ(M, r). (29) Therefore Σ(A, b) ⊆ Σ(M, r). (30) Since \underline{M} = I − D^{−1}P = D^{−1}(D − P) = D^{−1}⟨A⟩, (31) and ⟨A⟩ is an M-matrix, we have \underline{M}^{−1} = ⟨A⟩^{−1}D ≥ 0. (32) Theorem 1.3 then implies that M^{H}r = [\underline{x}, \overline{x}] (33) with

\underline{x}_{i} = \{ \begin{matrix}c_{i}e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} ≥ 0 \\ e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} < 0 \end{matrix} (34)

\overline{x}_{i} = \{ \begin{matrix}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} ≥ 0 \\ c_{i}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} < 0 \end{matrix} (35)

for i = 1, 2, …, n, where

(36) s^{(i)} = \{ \begin{matrix}\overline{r}_{i} & for j = i \\ max\{−\underline{r}_{j}, \overline{r}_{j}\} & for j ≠ i; j = 1, …, n \end{matrix}

(37) t^{(i)} = \{ \begin{matrix}\underline{r}_{i} & for j = i \\ min\{\underline{r}_{j}, −\overline{r}_{j}\} & for j ≠ i; j = 1, …, n \end{matrix}

and

c_{i} = \frac{1}{2(\underline{M}^{−1})_{ii} − 1}. (38)

Noting that s^{(i)} = |r| + (\overline{r}_{i} − |r_{i}|)e_{i} = D^{−1}(|b| + (\overline{b}_{i} − |b_{i}|)e_{i}), (39)

t^{(i)} = −|r| + (\underline{r}_{i} + |r_{i}|)e_{i} = D^{−1}(−|b| + (\underline{b}_{i} + |b_{i}|)e_{i}), (40) and

c_{i} = \frac{1}{2(⟨A⟩^{−1})_{ii}Ǎ_{ii} − 1} = \frac{\underline{A}_{ii} − α_{i}}{\overline{A}_{ii} + α_{i}} , (41)

8 S. Ning and R. B. Kearfott

we have

\begin{aligned}e_{i}^{T} \underline{M}^{−1}t^{(i)} & = e_{i}^{T} ⟨A⟩^{−1}(−|b| + (\underline{b}_{i} + |b_{i}|)e_{i}) \\ & = −u_{i} + (⟨A⟩^{−1})_{ii}(\underline{b}_{i} + |b_{i}|) \\ & = (−β_{i} + \underline{b}_{i})/(\underline{A}_{ii} − α_{i}) \end{aligned} (42)

and

\begin{aligned}e_{i}^{T} \underline{M}^{−1}s^{(i)} & = e_{i}^{T} ⟨A⟩^{−1}(|b| + (\overline{b}_{i} − |b_{i}|)e_{i}) \\ & = u_{i} + (⟨A⟩^{−1})_{ii}(\overline{b}_{i} − |b_{i}|) \\ & = (β_{i} + \overline{b}_{i})/(\underline{A}_{ii} − α_{i}). \end{aligned} (43)

By (34) and (35), we have

(44) x_{i} = \{ \begin{matrix}[(−β_{i} + \underline{b}_{i})/(\overline{A}_{ii} + α_{i}), (β_{i} + \overline{b}_{i})/(\underline{A}_{ii} − α_{i})] & for x_{i} > 0, \\ [(−β_{i} + \underline{b}_{i})/(\underline{A}_{ii} − α_{i}), (β_{i} + \overline{b}_{i})/(\underline{A}_{ii} − α_{i})] & for 0 ∈ x_{i}, \\ [(−β_{i} + \underline{b}_{i})/(\underline{A}_{ii} − α_{i}), (β_{i} + \overline{b}_{i})/(\overline{A}_{ii} + α_{i})] & for x_{i} < 0. \end{matrix}

That is,

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

It follows from (30) that A^{H}b ⊆ M^{H}r = x. (46)

If the midpoint of A is diagonal, then A = D + [−1, 1]P. (47) By (29), we have Σ(A, b) = Σ(M, r). (48) Hence A^{H}b = M^{H}r = x. (49)

For a general H-Matrix A, let E_{n×n} be a diagonal matrix with

E_{ii} = \{ \begin{matrix}−1 & for A_{ii} < 0, \\ 1 & for A_{ii} > 0. \end{matrix}

By Theorem 2.1, we have Σ(A, b) = Σ(EA, Eb). (50) Since ⟨EA⟩ = ⟨A⟩, |Eb| = |b|, (EA)_{ii} = E_{ii}A_{ii} > 0, and EA is also an H-matrix, using EA and Eb instead of A and b respectively, it follows from the above proof that

((EA)^{H}(Eb))_{i} ⊆ \frac{(Eb)_{i} + [−β_{i}, β_{i}]}{(EA)_{ii} + [−α_{i}, α_{i}]} (51)

Solving Linear Interval Equations 9

for i = 1, 2, …, n, and equality holds in Equation (51) if the midpoint of A is diagonal, since

\frac{(Eb)_{i} + [−β_{i}, β_{i}]}{(EA)_{ii} + [−α_{i}, α_{i}]} = \frac{(b)_{i} + [−β_{i}, β_{i}]}{(A)_{ii} + [−α_{i}, α_{i}]} = x_{i}. (52)

By (50) , (51) and (52), we have A^{H}b ⊆ x (53) and equality holds in Equation (53) when the midpoint of A is diagonal.

Theorem 2.3. Suppose that A ∈ 𝕀ℝ^{n×n} and b ∈ 𝕀ℝ^{n}. Assume A is inverse positive and \underline{b} = −\overline{b} ≠ 0. Then the hull of the system Ax = b (54) and the hull of the system Mx = r (55) are the same, i.e. A^{H}b = M^{H}r, (56) where M = Ǎ^{−1}A, r = Ǎ^{−1}b and Ǎ = (\underline{A} + \overline{A})/2. Thus, Hansen’s scheme of [6] (as in Theorem 1.3) gives exactly the hull in this case.

Proof. Theorem 1.2 implies A^{H}b = [\underline{A}^{−1}\underline{b}, \underline{A}^{−1}\overline{b}]. (57) Since Ǎ^{−1} ≥ 0, we have

M = Ǎ^{−1}A = I + [−1, 1]Ǎ^{−1}rad(A).

Hence,

M = [Ǎ^{−1}\underline{A}, Ǎ^{−1}\overline{A}].

From this, we have

\underline{M}^{−1} = \underline{A}^{−1}Ǎ = I + \underline{A}^{−1}(Ǎ − \underline{A}) ≥ 0.

It follows from Theorem 1.3 that M^{H}r = [\underline{x}, \overline{x}], (58) with

\underline{x}_{i} = \{ \begin{matrix}c_{i}e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} ≥ 0 \\ e_{i}^{T} \underline{M}^{−1}t^{(i)} & for \underline{x}_{i} < 0 \end{matrix} (59)

\overline{x}_{i} = \{ \begin{matrix}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} ≥ 0 \\ c_{i}e_{i}^{T} \underline{M}^{−1}s^{(i)} & for \overline{x}_{i} < 0 \end{matrix} (60)

for i = 1, 2, …, n, where s^{(i)}, t^{(i)} and c_{i} satisfy (9), (10) and (11). Since

r = Ǎ^{−1}b = Ǎ^{−1}b̌ + [−1, 1]Ǎ^{−1}rad(b) = [Ǎ^{−1}\underline{b}, Ǎ^{−1}\overline{b}],

we have \underline{r} = Ǎ^{−1}\underline{b}, \overline{r} = Ǎ^{−1}\overline{b}. Noting that \underline{b} = −\overline{b}, we have \underline{r} = −\overline{r}. The definition of s^{(i)} and t^{(i)} implies s^{(i)} = \overline{r} and t^{(i)} = \underline{r} for i = 1, 2, …, n. Since \underline{x} ≤ Ǎ^{−1}\underline{b} < 0 and \overline{x} ≥ Ǎ^{−1}\overline{b} > 0, we have

\underline{x} = \underline{M}^{−1}\underline{r} = (Ǎ^{−1}\underline{A})^{−1}(Ǎ^{−1}\underline{b}) = \underline{A}^{−1}\underline{b} (61)

and

\overline{x} = \underline{M}^{−1}\overline{r} = (Ǎ^{−1}\underline{A})^{−1}(Ǎ^{−1}\overline{b}) = \underline{A}^{−1}\overline{b}. (62)

Equations (57) and (58) then imply

A^{H}b = M^{H}r. (63)

We now turn our attention to the triangular factors of interval M-matrices, and to when the interval Gaussian elimination algorithm (IGA) gives the exact hull. The following lemma is from Fiedler and Pták [4].

Lemma 2.4. Suppose that A ∈ ℝ^{n} is an M-matrix and A = LU, where L is lower triangular, U is upper triangular, and L_{ii} = 1 for i = 1, 2, …, n. Then

U_{ij} = A_{ij} − ∑_{k=1}^{i−1} L_{ik}U_{kj}, for j ≥ i, (64)

L_{ij} = (A_{ij} − ∑_{k=1}^{j−1} L_{ik}U_{kj})/U_{jj}, for i > j, (65)

and U_{ii} > 0, U_{ji} ≤ 0, L_{ij} ≤ 0 for i > j. That is, the triangular factors of an M-matrix are themselves M-matrices.

[14, Theorem 4.5.8] implies the following theorem.

Theorem 2.5. (Barth and Nuding [2], Beeck [3]) Suppose that A = [\underline{A}, \overline{A}] ∈ 𝕀ℝ^{n×n} is an interval M-matrix. Let L = [\underline{L}, \overline{L}] be that lower triangular interval matrix with L_{ii} = [1, 1] for i = 1, 2, …, n, and let U = [\underline{U}, \overline{U}] be that upper triangular interval matrix, defined by

U_{ij} = A_{ij} − ∑_{k=1}^{i−1} L_{ik}U_{kj}, for j ≥ i (66)

(67) L_{ij} = (A_{ij} − ∑_{k=1}^{j−1} L_{ik}U_{kj})/U_{jj}, for i > j

Then A ⊆ LU, \underline{A} = \underline{L} \underline{U} and \overline{A} = \overline{L} \overline{U}. Hence L and U are interval M-matrices. Also A^{−1} = U^{−1}L^{−1}, that is, [\overline{A}^{−1}, \underline{A}^{−1}] = [\overline{U}^{−1}, \underline{U}^{−1}][\overline{L}^{−1}, \underline{L}^{−1}]. Moreover

A^{H}b ⊆ U^{−1}(L^{−1}b) = [\overline{A}^{−1}, \underline{A}^{−1}]b. (68)

If b ≤ 0, b ≥ 0 or 0 ∈ b, then IGA computes the exact hull A^{H}b, i.e.

A^{H}b = U^{−1}(L^{−1}b) = [\overline{A}^{−1}, \underline{A}^{−1}]b. (69)

The following theorem, a consequence of Theorem 1.2, gives additional conditions under which we may easily obtain the interval hull A^{H}b when A is inverse positive.

Theorem 2.6. Suppose A ∈ 𝕀ℝ^{n×n} is inverse positive, suppose b, x^{(0)} ∈ 𝕀ℝ^{n}, and suppose A^{H}b ⊆ x^{(0)}. For i, k = 1, 2, …, n, define A^{(1)}, A^{(2)} ∈ A by

A_{ik}^{(1)} = \overline{A}_{ik} if \underline{x}_{k}^{(0)} ≥ 0 and A_{ik}^{(1)} = \underline{A}_{ik} otherwise,

A_{ik}^{(2)} = \overline{A}_{ik} if \overline{x}_{k}^{(0)} ≤ 0 and A_{ik}^{(2)} = \underline{A}_{ik} otherwise.

Define

x = □ \{(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}\}
= ([\underline{x}_{1}, \overline{x}_{1}], …, [\underline{x}_{n}, \overline{x}_{n}])^{T} ,

where □{v, w} is the interval hull of the vectors v and w, i.e. the smallest interval vector that contains the set {v, w}. Then

x ⊆ A^{H}b ⊆ x^{(0)}.

In particular, if \underline{x}_{k}\underline{x}_{k}^{(0)} ≥ 0 and \overline{x}_{k}\overline{x}_{k}^{(0)} ≥ 0 for k = 1, 2, …, n, then

x = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = A^{H}b.

Proof. Since A^{(1)} ∈ A and A^{(2)} ∈ A, we have
(A^{(1)})^{−1}\underline{b} ∈ Σ(A, b) and (A^{(2)})^{−1}\overline{b} ∈ Σ(A, b).

Hence

x ⊆ A^{H}b ⊆ x^{(0)}. (70)

Let us define y by

y = [\underline{y}, \overline{y}] = A^{H}b.

Then Theorem 1.2 gives y = [(B^{(1)})^{−1}\underline{b}, (B^{(2)})^{−1}\overline{b}], where

B_{ik}^{(1)} = \overline{A}_{ik} if \underline{y}_{k} ≥ 0 and B_{ik}^{(1)} = \underline{A}_{ik} otherwise,

B_{ik}^{(2)} = \overline{A}_{ik} if \overline{y}_{k} ≤ 0 and B_{ik}^{(2)} = \underline{A}_{ik} otherwise,

for i, k = 1, 2, …, n. Also, (70) gives \underline{x}_{k} ≥ \underline{y}_{k} ≥ \underline{x}_{k}^{(0)} and \overline{x}_{k} ≤ \overline{y}_{k} ≤ \overline{x}_{k}^{(0)} for

k = 1, 2, …, n. If, in addition, \underline{x}_{k}\underline{x}_{k}^{(0)} ≥ 0 and \overline{x}_{k}\overline{x}_{k}^{(0)} ≥ 0 for k = 1, 2, …, n, then \underline{y}_{k}\underline{x}_{k}^{(0)} ≥ 0 and \overline{y}_{k}\overline{x}_{k}^{(0)} ≥ 0 for k = 1, 2, …, n. Thus, the definitions of A^{(1)} and B^{(1)} imply that

A_{ik}^{(1)} = B_{ik}^{(1)} if \underline{y}_{k} ≠ 0

for i, k = 1, 2, …, n, and

A^{(1)}\underline{y} = B^{(1)}\underline{y} = \underline{b}.

That is

(A^{(1)})^{−1}\underline{b} = (B^{(1)})^{−1}\underline{b} = \underline{y}.

Similarly, the definitions of A^{(2)} and B^{(2)} imply that

A_{ik}^{(2)} = B_{ik}^{(2)} if \overline{y}_{k} ≠ 0

for i, k = 1, 2, …, n, and

A^{(2)}\overline{y} = B^{(2)}\overline{y} = \overline{b}.

That is

(A^{(2)})^{−1}\overline{b} = (B^{(2)})^{−1}\overline{b} = \overline{y}.

Therefore,

A^{H}b = y = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = x,

thus completing the proof.

For example, if we know the box x^{(0)} obtained by IGA, [14, Proposition 3.6.6], Hansen’s technique of [6], or any other solution method for bounding the solution set, the box x of the above theorem is the hull (to within roundout error) whenever (A^{(1)})^{−1}\underline{b} ≤ (A^{(2)})^{−1}\overline{b}, \underline{x}_{k}\underline{x}_{k}^{(0)} ≥ 0, and \overline{x}_{k}\overline{x}_{k}^{(0)} ≥ 0 for every k = 1, 2, …, n.

3. Examples. The following examples illustrate the cases covered in the theory of the preceding section, and examine several cases outside this theory as well.

The computations in these examples were programmed in Fortran 90 using the system of [8]. The arithmetic was outwardly rounded interval arithmetic based on double precision IEEE arithmetic. This arithmetic was accessed with a modification of the Fortran 77 library INTLIB [7] that allows true directed roundings on Sun Sparc machines, as in [9]. The answers below are presented to three digits. The endpoints of the intervals in the answers are rounded using the default conversion routines in the printing functions of Fortran 90. Thus, the exhibited interval endpoints may differ from correctly outwardly rounded results by a unit in the last digit.

Example 3.1. Let the interval equation system Ax = b be given with

A = ( \begin{matrix}[4, 6] & [−1, 1] & [−1, 1] & [−1, 1] \\ [−1, 1] & [−6, −4] & [−1, 1] & [−1, 1] \\ [−1, 1] & [−1, 1] & [9, 11] & [−1, 1] \\ [−1, 1] & [−1, 1] & [−1, 1] & [−11, −9] \end{matrix} )

and

b = ( \begin{matrix}[−2, 4] \\ [1, 8] \\ [−4, 10] \\ [2, 12] \end{matrix} ) .

Then

Ǎ = ( \begin{matrix}5 \\ & −5 \\ & & 10 \\ & & & −10 \end{matrix} ) ,

Ǎ^{−1} = ( \begin{matrix}0.2 \\ & −0.2 \\ & & 0.1 \\ & & & −0.1 \end{matrix} ) ,

M = Ǎ^{−1}A = ( \begin{matrix}[0.8, 1.2] & [−0.2, 0.2] & [−0.2, 0.2] & [−0.2, 0.2] \\ [−0.2, 0.2] & [0.8, 1.2] & [−0.2, 0.2] & [−0.2, 0.2] \\ [−0.1, 0.1] & [−0.1, 0.1] & [0.9, 1.1] & [−0.1, 0.1] \\ [−0.1, 0.1] & [−0.1, 0.1] & [−0.1, 0.1] & [0.9, 1.1] \end{matrix} ) ,

and

r = Ǎ^{−1}b = ( \begin{matrix}[−0.4, 0.8] \\ [−1.6, −0.2] \\ [−0.4, 1] \\ [−1.2, −0.2] \end{matrix} ) .

Using interval Gaussian elimination, we obtain the box

( \begin{matrix}[−2.60, 3.10] \\ [−3.90, 1.50] \\ [−1.43, 2.15] \\ [−2.35, 0.60] \end{matrix} ) .

Using Theorem 2.1, we obtain A^{H}b = M^{H}r. Then, with Theorem 1.3, we compute the smallest box containing the solution set to be

A^{H}b = M^{H}r = ( \begin{matrix}[−2.50, 3.10] \\ [−3.90, 1.20] \\ [−1.40, 2.15] \\ [−2.35, 0.60] \end{matrix} ) .

Using Theorem 2.2, we obtain the same box as above.

Example 3.2. Let the interval equation system Ax = b be given with

A = ( \begin{matrix}[3.7, 4.3] & [−1.5, −0.5] & [0, 0] \\ [−1.5, −0.5] & [3.7, 4.3] & [−1.5, −0.5] \\ [0, 0] & [−1.5, −0.5] & [3.7, 4.3] \end{matrix} ) ,

b = ( \begin{matrix}[−14, 14] \\ [−9, 9] \\ [−3, 3] \end{matrix} ) .

Then

Ǎ = ( \begin{matrix}4 & −1 & 0 \\ −1 & 4 & −1 \\ 0 & −1 & 4 \end{matrix} ) ,

A is an M-matrix, and \underline{b} = −\overline{b}. Using interval Gaussian elimination, we obtain the box

( \begin{matrix}[−6.38, 6.38] \\ [−6.40, 6.40] \\ [−3.40, 3.40] \end{matrix} ) .

Using Theorem 1.3 and Theorem 2.2, we obtain the same box as above. Theorem 2.3 and Theorem 2.5 imply that in this case, we will obtain the exact hull of the original system, to within roundout error, by IGA, Hansen’s technique [6], or Rohn’s reformulation of Hansen’s technique.

Example 3.3. Let the interval equation system Ax = b be given, with A as in Example 3.2 and

b = ( \begin{matrix}[−14, 0] \\ [−9, 0] \\ [−3, 0] \end{matrix} ) .

Using interval Gaussian elimination, we obtain the box

( \begin{matrix}[−6.38, 0] \\ [−6.40, 0] \\ [−3.40, 0] \end{matrix} ) .

Theorem 2.5 implies that in this case IGA will give the smallest box containing the solution set. Using Hansen’s technique of [6] or Rohn’s reformulation of [16], we obtain the wider box

( \begin{matrix}[−6.38, 1.12] \\ [−6.40, 1.54] \\ [−3.40, 1.40] \end{matrix} ) .

Using Theorem 2.2, we obtain an even wider box:

( \begin{matrix}[−6.38, 1.67] \\ [−6.40, 2.77] \\ [−3.40, 2.40] \end{matrix} ) .

Note that Hansen’s technique and Rohn’s reformulation give the same results, since they both give the hull of the preconditioned systems.

Example 3.4. Let the interval equation system Ax = b be given with A as in Example 3.2 and

b = ( \begin{matrix}[0, 14] \\ [0, 9] \\ [0, 3] \end{matrix} ) .

Using IGA, we obtain the box

( \begin{matrix}[0, 6.38] \\ [0, 6.40] \\ [0, 3.40] \end{matrix} ) .

Theorem 2.5 implies that IGA computes the hull of the original system of equations in this case. Using Hansen’s technique of [6], we obtain the wider box

( \begin{matrix}[−1.12, 6.38] \\ [−1.54, 6.40] \\ [−1.40, 3.40] \end{matrix} ) .

Using Theorem 2.2, we obtain an even wider box:

( \begin{matrix}[−1.67, 6.38] \\ [−2.77, 6.40] \\ [−2.40, 3.40] \end{matrix} ) .

Example 3.5. Let the interval equation system Ax = b be given with A as in Example 3.2 and

b = ( \begin{matrix}[2, 14] \\ [−9, −3] \\ [−3, 1] \end{matrix} ) .

Using IGA, we obtain the box

( \begin{matrix}[−1.09, 4.29] \\ [−4.02, 1.24] \\ [−2.44, 0.773] \end{matrix} ) .

Using Hansen’s technique of [6], we obtain the wider box

( \begin{matrix}[−0.995, 5.01] \\ [−4.64, 1.52] \\ [−2.69, 1.38] \end{matrix} ) .

Using [14, Proposition 3.6.6], we obtain the box

A^{−1}b = [\overline{A}^{−1}, \underline{A}^{−1}]b = ( \begin{matrix}[−1.20, 4.69] \\ [−4.06, 1.73] \\ [−2.47, 1.18] \end{matrix} ) .

We may now use IGA, Hansen’s technique of [6], or [14, Proposition 3.6.6] to compute a box x^{(0)} for Theorem 2.6; the resulting box, regardless of which of these methods was used for x^{(0)}, is

x = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = A^{H}b = ( \begin{matrix}[−0.995, 4.29] \\ [−3.79, 1.24] \\ [−2.35, 0.773] \end{matrix} ) ,

since the signs of \underline{x} and \underline{x}^{(0)} are the same and the signs of \overline{x} and \overline{x}^{(0)} are the same.

Example 3.6. Let the interval equation system Ax = b be given with A as in Example 3.2 and

b = ( \begin{matrix}[2, 14] \\ [3, 9] \\ [−3, 1] \end{matrix} ) .

Using IGA, we obtain the box

( \begin{matrix}[0.517, 6.25] \\ [0.450, 6.07] \\ [−0.881, 2.73] \end{matrix} ) .

Using Hansen’s technique of [6], we obtain the wider box

( \begin{matrix}[−0.206, 6.25] \\ [−0.386, 6.07] \\ [−2.01, 2.73] \end{matrix} ) .

Using [14, Proposition 3.6.6], we obtain the box

A^{−1}b = [\overline{A}^{−1}, \underline{A}^{−1}]b = ( \begin{matrix}[0.356, 6.25] \\ [0.283, 6.07] \\ [−0.919, 2.73] \end{matrix} ) .

If x^{(0)} is obtained by IGA or [14, Proposition 3.6.6], Theorem 2.6, gives the box

x = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = A^{H}b = ( \begin{matrix}[0.523, 6.25] \\ [0.499, 6.07] \\ [−0.743, 2.73] \end{matrix} ) ,

since the signs of \underline{x} and \underline{x}^{(0)} are the same, and the signs of \overline{x} and \overline{x}^{(0)} are the same. However, if x^{(0)} is obtained by Hansen’s technique of [6], using Theorem 2.6, we obtain the box

x = [(A^{(1)})^{−1}\underline{b}, (A^{(2)})^{−1}\overline{b}] = ( \begin{matrix}[0.964, 6.25] \\ [1.04, 6.07] \\ [−0.387, 2.73] \end{matrix} ) ⊂ A^{H}b.

Example 3.7. Let the interval equation system Ax = b be given with

A = ( \begin{matrix}[15, 17] & [−3, 3.01] & [−3, 3.01] & [−3, 3.01] \\ [−3, 3.01] & [15, 17] & [−3, 2.99] & [−3, 2.99] \\ [−3, 2.99] & [−3, 2.99] & [15, 17] & [−3, 3.01] \\ [−3, 3.01] & [−3, 3.01] & [−3, 2.99] & [15, 17] \end{matrix} ) ,

b = ( \begin{matrix}[−6, −2] \\ [4, 5] \\ [−2, 4] \\ [8, 10] \end{matrix} ) .

By IGA, we obtain the box

( \begin{matrix}[−1.03, 0.495] \\ [−0.347, 0.974] \\ [−0.770, 0.917] \\ [0.150, 1.25] \end{matrix} ) .

Using Hansen’s technique of [6], we obtain the box

( \begin{matrix}[−1.03, 0.363] \\ [−0.223, 0.975] \\ [−0.752, 0.919] \\ [0.149, 1.25] \end{matrix} ) .

Examining the above examples, we see that sometimes IGA gives sharper bounds on the hull, while in other cases Hansen’s technique does. In some cases, the intersection of the results of IGA and Hansen’s technique is narrower than either result taken separately.

4. Acknowledgements. We wish to thank the referees for their deep and careful work. In particular, we wish to thank a referee for proposing (24), which resulted in a statement of Theorem 2.2 that is clearer and more general than our original one.

REFERENCES

[1] G. Alefeld and J. Herzberger, Introduction to Interval Computations, Academic Press, New York, 1983.

[2] W. Barth and E. Nuding, Optimale Lösung von Intervallgleichungssystemen, Computing, 12 (1974), pp. 117–125.

[3] H. Beeck, Zur Problematik der Hüllenbestimmung von Intervallgleichungssystemen, in Interval Mathematics, K. Nickel, ed., Springer-Verlag, New York, 1975, pp. 150–159.

[4] M. Fiedler and V. Pták, On Matrices with Non-Positive Off-Diagonal Elements and Positive Principal Minors, Czechoslovak Math. J., 12 (1962), pp. 382–400.

[5] E. R. Hansen, Global Optimization Using Interval Analysis, Marcel Dekker, Inc., New York, 1992.

[6] E. R. Hansen, Bounding the Solution of Interval Linear Equations, SIAM J. Numer. Anal., 29 (1992), pp. 1493–1503.

[7] R. B. Kearfott, M. Dawande, K.-S. Du and C. Hu, Algorithm 737: INTLIB: A Portable FORTRAN 77 Interval Standard Function Library, ACM TransṀathṠoftware, 20 (1994), pp. 447–459.

[8] R. B. Kearfott, A Fortran 90 Environment for Research and Prototyping of Enclosure Algorithms for Nonlinear Equations and Global Optimization, ACM TransṀathṠoftware, 21 (1995), pp. 63–78.

[9] O. Knüppel, PROFIL/BIAS — A Fast Interval Library, Computing, 53 (1994), pp. 277–287.

[10] J. Kuttler, A Fourth-Order Finite-Difference Approximation for the Fixed Membrane Eigenproblem, Math. Comp., 25 (1971), pp. 237–256.

[11] R. E. Moore, Methods and Applications of Interval Analysis, SIAM, Philadelphia, 1979.

[12] A. Neumaier, New Techniques for the Analysis of Linear Interval Equations, Linear Algebra Appl., 58 (1984), pp. 273–32.

[13] A. Neumaier, On the Comparison of H-Matrices with M-Matrices, Linear Algebra Appl. 83 (1986), pp. 135–141.

[14] A. Neumaier, Interval Methods for Systems of Equations, Cambridge University Press, Cambridge, England, 1990.

[15] H. Ratschek and J. Rokne, New Computer Methods for Global Optimization, Wiley, New York, 1988.

[16] J. Rohn, Cheap and Tight Bounds: The Recent Result by E. Hansen Can Be Made More Efficient, Interval Computations, (1993), pp. 13–21.

[17] J. Rohn, NP-Hardness Results for Linear Algebraic Problems with Interval Data, preprint,