Interval unions

This paper introduces interval union arithmetic, a new concept which extends the traditional interval arithmetic. Interval unions allow to manipulate sets of disjoint intervals and provide a natural way to represent the extended interval division. Considering interval unions lead to simplifications of the interval Newton method as well as of other algorithms for solving interval linear systems. This paper does not aim at describing the complete theory of interval union analysis, but rather at giving basic definitions and some fundamental properties, as well as showing theoretical and practical usefulness of interval unions in a few selected areas.


Introduction
Interval analysis is a branch of numerical analysis that was born in the 1960's. It consists of computing with intervals of reals instead of reals, providing a framework for handling uncertainties and verified computations (see e.g. [2,16,21,23] for a survey). Interval analysis is a key ingredient for numerical constraint satisfaction (see e.g. [14]) and global optimization (see e.g. [8]). Global optimization solvers like Gloptlab [5] and COCONUT [28,29] rely heavily on interval analysis to guarantee rigorous solutions, even non-rigorous solvers like BARON [27] and α-Branch and bound [1] use rigorous computations in some steps of the search.
In practice, interval arithmetic must be implemented using outward rounding in order to assure that the result of an interval calculation always contains the result of the corresponding real valued operation evaluated for each value(s) of the used interval(s). Interval arithmetic has been implemented in almost every programming language which is relevant for scientific computing, see for example Intlab [26] for Matlab, Filib++ [22] for C/C++, Interval [15] for Fortran and MathInterval [6] for Java.
This example shows the problem of interval arithmetic both from a theoretical and a computational point of view. For the theory of intervals it is an issue since the result of an elementary operation involving two intervals does not belong necessarily to the set of intervals 1 while for computations it is a problem since the interval division operator requires special treatment.
This paper extends the concept of interval arithmetic to interval unions. An interval union is a set of closed and disjoint intervals where the bounds of the extreme intervals can be ±∞. During the paper we demonstrate that interval unions generalize intervals and allow among others to represent the result of interval division in a natural way.
We show that the interval Newton method for univariate functions can be reformulated using interval union arithmetic. Our numerical experiments demonstrate that this new method, called interval union Newton method, produces much better results than the traditional approach.
Interval unions are also useful in algorithms for finding rigorous bounds on the solution set of interval linear systems. We give an example where the Gaussian elimination algorithm employed using interval union arithmetic gives a rigorous bound for an interval linear system that is not regular. Moreover, we show that the interval union Gauss-Seidel procedure can be used in the multivariate interval Newton method in order to produce multiple gaps instead of only one as suggested by Hansen [8].
Some of the theoretical results of interval analysis remain valid when we are dealing with interval unions. That is the case, e.g., for the fundamental theorem of interval arithmetic, and therefore the natural extension of real functions to interval unions is similar to the interval case. On the other hand, some inclusion results like the interval mean value theorem do not hold for interval unions, not even for the univariate case. During the paper it is shown that a large part of the interval union arithmetic can be easily implemented if we have an interval arithmetic library at our disposal.
A closely related concept to interval unions is that of multi-intervals, introduced by Yakovlev [32]. According to [31], they are defined as a union of closed intervals that are not necessarily disjoint, making them slightly more general from the interval unions of the present paper. Multi-interval arithmetic is (a not separately accessible) part of the publicly available software Unicalc [3,24] for solving constraint satisfaction problems and nonlinear systems of equations. Another variant are the discontinuous intervals by Hyvönen [10]. A discontinuous interval is the disjoint union of closed, half-open, or open intervals. For applications of discontinous intervals, see [11].
In Sect. 2 we present the basics of interval arithmetic. The section is mainly a revision of the traditional case in the extended context. Section 3 describes the generalization from intervals to interval unions, where the basic interval union operations are defined, isotonicity property shown, the fundamental theorem of interval union arithmetic is proven. In addition, in this section, hull and component-wise operations are also defined.
In Sect. 4 the interval union Newton method for univariate functions is presented. Similar as for the interval Newton method the aim is to enclose all roots of f (x) ∈ R subject to x ∈ X where both, R and X are interval unions. We show that the definition of Newton methods can be made through component-wise operations and compare our new approach with the traditional interval Newton algorithm in a set of 32 problems. Our experiment shows that interval union arithmetic can improve Newton methods significantly in the univariate case.
Finally in Sect. 5, interval union linear systems are studied and shown that the interval Gaussian elimination and Gauss-Seidel algorithms can be extended from intervals to interval unions. The advantages of replacing interval operations by interval unions in linear systems are demonstrated by performing tests on examples in low dimension.

Notation
We mostly follow [18] for the notation of interval arithmetic. Throughout this paper R m×n denotes the vector space of all m × n matrices A with real entries A ik (i = 1, . . . , m, k = 1, . . . , n), and R n = R n×1 denotes the vector space of all column vectors v of length n and entries v i (i = 1, . . . , n). For vectors and matrices, the relations =, =, <, >, ≤, ≥ and the absolute value |A| of the matrix A are interpreted component-wise.
We write A T to represent the transpose of a matrix A and A −T is short for (A T ) −1 . The ith row vector of a matrix A is denoted by A i: and the jth column vector by A : j . For the n×n matrix A, diag(A) denotes the n-dimensional vector with diag(A) i = A ii .
The number of elements of the index set N is given by |N |. Let I ⊆ {1, . . . , m} and J ⊆ {1, . . . , n} be index sets and let n I := |I |, n J := |J |. For the n-dimensional vector x, x J denotes the n J -dimensional vector built from the components of x selected by the index set J . For the m × n matrix A, the expression A I : denotes the n I × n matrix built from the rows of A selected by the index sets I . Similarly, A :J denotes the m × n J matrix built from the columns of A selected by the index sets J .

Interval arithmetic
This section presents the basics of interval arithmetic. A comprehensive approach to this topic is given by [23]. We are mainly interested in extended interval arithmetic. i.e, when division by intervals containing 0 is allowed. Good references to extended interval arithmetic are [8] and [17].
Let a, a ∈ R with a ≤ a then a = [a, a] denotes a real interval with inf(a) = min(a) = a and sup(a) = max(a) = a. The set of nonempty compact real intervals is denoted by We extend the definition of real intervals by permitting the bounds of intervals to be one of the ideal points −∞ and ∞ and define IR as the set of closed real intervals. We write The width of the interval a ∈ IR \ {∅} is given by wid(a) := a − a, its mignitude by otherwise. and its magnitude by |a| := max(|a|, |a|). The midpoint of a ∈ IR isǎ := mid(a) := (a +a)/2 and the radius of a ∈ IR isâ := rad(a) := (a −a)/2. For a ∈ IR there is no natural definition of a midpoint. Moreover, ifǎ is well defined then a ∈ a ⇔ |a−ǎ| ≤â and we say that midrad(ǎ,â) is the midrad representation of interval a. For a set S the smallest box containing S is called the interval hull of S and denoted by S. An interval is called thin or degenerate if wid(a) = 0.
The inclusion relations are given as An interval vector x = [x, x] or box is the Cartesian product of the closed real intervals x i := [x i , x i ] ∈ IR. We write IR n to denote the set of all n-dimensional boxes. We also define the interval matrix A = [A, A] in a similar way and IR m×n denotes the set of all m × n interval matrices. Operations defined for intervals (like width, midpoint, radius, mignitude and magnitude) are defined component-wise when applied to boxes or matrices. Let a, b ∈ IR. The elementary real operations • ∈ {+, −, /, * , } are extended to the interval arguments a, b by defining the result of an elementary interval operation to be the set of real numbers which results from combining any two numbers contained in a and in b. Formally, This leads to operations on IR defined by a•b := (a• • b). The elementary operations are inclusion isotonic. That means: For a, b ∈ IR we get that As one can see in the cases marked with * , the result is not a single interval but the union of two disjoint ones. As shown in [25] the division defined by (1) is inclusion isotonic (also see [17]). Let x ∈ IR n and f : D ⊆ R n → R. We define rg•( f (x)) to be the set and call it the range of f over the box x. We extend the range to a function on IR by rg( f (x)) := rg•( f (x)), also called the range of f . We say that a function f : We already established that elementary interval operations are inclusion isotonic and it is also possible to construct interval functions with the isotonicity property for standard functions like exponential, logarithmic and trigonometric, see for example [26] or [6]. Moreover, it is easy to prove that the composition of inclusion isotonic functions is also inclusion isotonic. Formally we have The interval function f : IR n → IR is an interval extension of a function f : D ⊆ If f admits a closed form and can be expressed in terms of elementary operations and standard functions we call the interval function f given by replacing every real operation with its interval counterpart the natural extension. Using these definitions we can formulate the fundamental theorem of interval analysis and prove it as in [21]:

Proposition 2 (Fundamental theorem of interval analysis) If f is inclusion isotonic and is an interval extension of f
Interval arithmetic also allows to prove a general version of the mean value theorem for multivariate functions, see [23]: Proposition 3 (Interval mean value theorem) Let F : R n → R n be a differentiable function defined on a box x ⊂ R n . If F is an interval extension of F, J an interval extension of the Jacobian of F and z ∈ x then Proposition 3 leads to the following Taylor extension, see [23].
We define the set and call it the kth partial inverse image of f on y and for its interval hull we write

Interval unions 3.1 Motivation
The well known interval Newton iteration is the interval variant of Newton's method for finding the roots of a function f in a box x. If (2) is applied to an arbitrary univariate function f : R → R and the starting interval x 0 , the interval Newton method splits and contracts x 0 into several intervals enclosing the zeros of f over x 0 . By (1) the division operator applied to two intervals a, b ∈ IR in the cases marked by a * do not map into IR. To solve this issue one can either define / : IR×IR\{0} → IR or for the marked cases one could take the interval hull of the two resulting intervals. However, keeping the two disjoint intervals in the marked cases is the reason why (2) works properly if 0 ∈ f (x). Therefore, it is obvious to define a structure where the division operator and therefore the interval Newton method is defined in a consistent and natural way. It serves as a motivation to introduce interval unions and define operations similar to the interval versions.

Definition
Definition 1 Throughout this paper, interval unions are denoted by bold calligraphic letters. An interval union u of length l(u) := k is a finite set of k disjoint intervals. Since for all disjoint intervals the natural ordering exists we denote the elements of u by u i and write The set of all interval unions of length ≤ k is denoted by U k and U := k≥0 U k is the set of all interval unions. In addition to this U 0 = ∅ and we identify U 1 with IR.
Definition 2 Let u := (u 1 , . . . , u k ) ∈ U be an interval union. We will identify u with the subset k i=1 u i of R that u represents, so for a real number x we say Similarly, for the interval x

Definition 3
Let S be a finite set of intervals, the union creator U(S) is defined as the smallest interval union u that satisfies a ⊆ u for all a ∈ S.

Lemma 1 Let S be a set of intervals, the union creator is inclusion isotonic:
Proof Follows directly from the definition.

Lemma 2 The interval hull of an union u ∈ U is given by
Proof Follows directly from Definition 1.
Definition 4 An interval union vector or box of dimension n is the cartesian product of n interval unions. We define U n k and U n as the set of all interval union vectors of dimension n and denote interval union boxes by lower case bold calligraphic letters like x or y. In a similar way we define interval union matrices as n × m arrays of interval unions. We introduce U n×m k and U n×m as the sets of interval union matrices of size n × m with the usual definition of the operations. Interval union matrices are given by capital bold calligraphic letters like A or B.
The interval union vector u ∈ U regarded as a subset of R n is always a finite set of boxes. More specifically, if u j has length k j we get the n j=1 k j disjoint boxes n j=1 u j, j , 1 ≤ j ≤ k j . We write for u ∈ IR n that u ∈ u iff u is one of these boxes.
Note that storing this set as an interval union vector requires just n j=1 k j intervals which is a clear advantage over storing all the individual boxes, especially in higher dimensions.
If u ∈ U k \{∅} we define the magnitude and mignitude of the interval union respectively by and u := min( u 1 , . . . , u k ).
We also define for u ∈ U k \{∅} the maximum, minimum and maximum width of Given the interval union u ∈ U k and a point x ∈ R we define the projection of x as follows Some functions defined for intervals do not extend naturally to interval unions. For such functions we present different definitions that can be useful in several contexts. Let u ∈ U k \{∅} be an interval union, we denote the component-wise midpoint and radius respectively byǔ c := (ǔ 1 , . . . ,ǔ k ) andû c := (û 1 , . . . ,û k ) whenever −∞ < u 1 ≤ u k < ∞. We denote the component-wise width and magnitude of u by wid(u) c := (wid(u 1 ), . . . , wid(u k )) and |u| c := (|u 1 |, . . . , |u k |) respectively. In some applications we also need to define operations over the hull of u. In such cases we add a subscript h to identify the hull operation. For example the hull mid-point operator and hull width of u are given byǔ h := mid( u) and wid(u) := wid( u).

Maximum length and filling gaps
The motivation from Sect. 3.1 hints at a problem which can arise when considering interval unions, since during iterative evaluations the number of intervals inside an union can grow uncontrollably. This can be easily anticipated if considering the task of finding zeros of a function having an infinite number of zeros in the starting box via the interval Newton method. Actually, this problem arises in several other interval methods where intervals unions could prove quite useful. We propose to solve the problem by restricting the maximum length of unions and by defining gap filling strategies.
Definition 5 Let u ∈ U be an interval union and let u i , u i+1 ∈ u. The open interval g i between the intervals u i and u i+1 is called the ith gap of u and is defined as (4) We denote by U k the set of all gap collections of size ≤ k and by U := i∈N U i the set of all gap collections. We will again identifyv ∈ U with the set v∈vv ⊆ R and write x ∈v, x ⊆v, and w ⊆v for x ∈ R n , x ∈ IR, and w ∈ U.

Lemma 3 Let u be an interval union of length k, and letû
Proof The result follows from Definition 5 and the strict inequality in (3).
Definition 7 Let u ∈ U k \U 1 and g ⊆û be a set of gaps of u. We define the gap filling F(u, g) ∈ U k−|g| as the unique interval union with F(u, g) =û\g and F(u, g) = u, i.e., we fill all the gaps from g in u.
If g i is the ith gap of u we get F(u, g i ) by setting u i := u i+1 and removing the interval u i+1 from u.
Lemma 5 For u ∈ U and g ⊆û we have Since F(u, g) = F(u, g\{g}) the general case follows by induction on the size of g.
Now we will introduce the concept of gap ordering to determine which gap to fill first. Usually, the width of the gap plays a part in that ordering (sometimes also a relative width with respect to the position of the interval along the real axis), and also the position of the gap might be interesting. Since we do not want to fix this ordering for developing the theory we will just assume that we are given a linear order on the set of all open intervals of R with the property that for arbitrary x ∈ IR every collection of disjoint open intervals contained in x has a maximal element w.r.t. . For example, two possible linear orders are given below.
Example 1 We say that g i g j if at least one of the following criteria is attained: where g i = (u i , u i+1 ) is given by Definition 4 and where C(u 1 , u 2 ) is the same as in the example above.
We must take care in both examples to avoid comparisons of type ∞ < ∞ when evaluating the width. In practice we use min(wid(g), M) for fixed M very large instead of wid(g).

Definition 8
The index set of the n smallest gaps of u (w.r.t. ) is defined by Similarly, the index set of the n largest gaps of u (w.r.t. ) is defined by For r ∈ {L , S} we denote by g r n (u) := {g i ∈û | i ∈ G r n (u)} the set of smallest respectively largest gaps of u. For convenience we define g r n (u) =û if n ≥ l(u) and g r n (u) = ∅ if n ≤ 0.

Definition 9
The length restriction mapping k : U → U k is given by k (u) := F(u, G S l(u)−k (u)), i.e., we fill the l(u) − k smallest gaps of u, and we do not change u if l(u) ≤ k.

Arithmetic for interval unions
In this section, similarly to interval arithmetic, basic set and elementary operations as well as properties like inclusion isotonicity are defined and explained for interval unions. Most of the theory translates nicely from intervals to interval unions, but some properties do not: e.g., due to the lack of convexity it is not possible to prove a mean value theorem for interval unions. (i) The union operation for u and x is defined as u ∪ x := U(u ∪ {x}). Obviously, we have (i ) The union operation for u and s is defined by (ii) The intersection operation for u and x is defined as u ∩x := U({u 1 ∩x, . . . , u k ∩ x}). We have (ii ) The intersection operation for u and S is defined by Note that there is a slight ambiguity in the notation, as u ∪ s can also denote the union of the two sets of intervals u and s. However, there will be no confusion between these two concepts, as the same real set is represented. (i) For the union operation defined by (6) we have x ∈ u ∪ x iff x ∈ u or x ∈ x.
(i ) For the union operation defined by (7) we have x ∈ u ∪ s iff x ∈ u or x ∈ s.
(ii) For the intersection operation defined by (8) we have x ∈ u ∩ x iff x ∈ u and x ∈ x. (ii ) For the intersection operation defined by (9) we have x ∈ u ∩ s iff x ∈ u and x ∈ s. Note that for the interval division operator (1) the above definition gives a natural embedding of the problematic cases into the set of interval unions: for arbitrary a, b ∈ IR we have Proof The union creator U is inclusion isotonic by Lemma 1. Interval operations are inclusion isotonic by Sect. 2, therefore the composition of them is also inclusion isotonic.
In addition to the usual definition of elementary operations we also introduce component-wise operations that will be useful in the context of interval union linear systems.
In the following we will fix a "cutoff" c ∈ R for filling the gaps as described before in Definition 10.

Definition 14
Let u ∈ U n be an interval union vector and s ∈ U an interval union, and let f : D ⊆ R n → R. For fixed l > 1 we define the range of length l of f over u (w.r.t. x) as and the kth partial inverse image of length l of f on u and s as As in the interval case, we call a function f : U n → U inclusion isotone if u ⊆ u ⇒ f(u ) ⊆ f(u). Moreover, we say f : U n → U is the interval union extension of a function f : We also refer to interval union extensions only as extensions when there is no possibility of misunderstandings. As in the interval case we can define a natural interval union extension for functions composed of elementary operations and standard function only by replacing real operations by their interval union counterparts. The following proposition states that the fundamental theorem of interval analysis can be naturally extended to interval unions.

Proposition 4 If f is inclusion isotonic and the interval union extension of f : R n → R then f rg (u) ⊆ f(u).
Proof immediately from the application of the fundamental theorem of interval analysis to every component u i of u = (u 1 , . . . , u k ).
On the other hand, due to the lack of convexity when working with interval unions we are not able to prove the interval union mean value theorem. For example, consider f (x) = x 2 and the interval union u = ([−3, −1], [1,3]). If we take x = −2 ∈ [−3, −1] and y = 2 ∈ [1, 3] then there is no ξ ∈ u such that 4 = 4 − 8ξ , and hence the statement fails even for univariate functions.

Interval union Newton method
In this section we consider the problem of rigorously enclosing all solutions of where f : R → R is a differentiable function. In Sect. 4.1 we review the interval Newton method for the case where r is set to be zero and x is a closed and bounded interval. In Sect. 4.2 we formulate the interval union Newton operator. Numerical experiments comparing both approaches are presented in Sect. 4.3.

Interval Newton method
Let x be a bounded interval and f : R → R a differentiable function. We are interested in enclosing all solutions of Interval Newton methods to solve this problem are based on the interval mean value theorem applied to (13). Formally, if y ∈ x such that f (y) = 0 then for any fixed x ∈ x. Therefore, the solution set of the problem can be given as (14) regardless of the choice of x. The usual interval Newton method fixes x as the midpoint of x and generates a sequence of nested intervals such that

The operator N (x) is called interval Newton function and is given by
Algorithms based on the interval Newton method can be divided into two groups depending on whether or not they rely on extended division, i.e. splitting intervals after the division into the two unconnected result intervals. Some authors like Moore [21] and Alefeld [2] only apply the interval Newton operator to boxes where 0 / ∈ f (x). More sophisticated algorithms like those proposed by Kearfott [17] and Hansen [8] allow division by intervals containing zero and process each box resulting from the division separately.
The simplest interval Newton method with extended division for enclosing all solutions of (13) is given in Algorithm 1. The algorithm takes the interval x and applies the interval Newton operator to it. If the resulting intervals are not empty or too thin then they are split, an interval to be processed is chosen and the iteration continues. The proof of finiteness and rigorousness of the interval Newton algorithm is given in [17]. For multivariate versions of this algorithm see [8,9].

Input:
The interval x 0 , the interval extensions f and f of f and f respectively and the narrow component tolerance > 0. Output: A list of intervals C with x ∈ C ⇒ wid(x) < and the guarantee that for all y ∈ x 0 with f (y) = 0 there exists at least one interval x ∈ C such that y ∈ x. 1: W ← x 0 ; 2: while W = ∅ do 3: x ← get_first(W); 4: x ←x; 5: The list of intervals C returned by the algorithm need not be disjoint. Moreover, it is possible that the algorithm saves an interval x in C even when it contains no root of f . The only guarantee we have is that when y ∈ x satisfies f (y) = 0 then y ∈ x i ⊆ C for some i.

Interval union Newton method
Let x and r be interval unions with p and q elements respectively. Applying the interval mean value theorem to each pair of intervals in x and r gives the solution set of (12) where S x (x, r) := {y ∈ x | ∃r ∈ r, f * ∈ f(x) and g * ∈ f (x) such that f * + g * (y − x) = r } for any fixed x ∈ x. Therefore, we can solve (12) by applying Algorithm 1 p ×q times. However, the interval union arithmetic provides a more natural approach, without the need of running multiple instances of the same algorithm. Let u k = (u k 1 , . . . , u k n ) be an interval union, f a differentiable function and f and f interval union extensions of f and of its derivative f . The interval union Newton iteration is given by where N (x) is the interval Newton function. Note that the interval union Newton iteration is rigorous since it is a component-wise application of the interval mean value theorem. Algorithm 2 uses (16) to enclose all solutions of (12). It also needs the auxiliary function checkAndRemove which is given in Algorithm 3. In the next section we perform numerical experiments to compare the performance of Algorithm 1 with Algorithm 2.

Algorithm 2 Interval union Newton algorithm
Input: The interval union u 0 , the interval union extensions f and f of f and f and the narrow component tolerance > 0. Output: The interval union s = (x i ) with wid(x i ) < and the guarantee that for all y ∈ u 0 with f (y) = 0 there exist an x i such that y ∈ x i . 1: u ← u 0 ; 2: while u = ∅ do 3: Newton operator 4: x ← ∅; 5: if wid(x i ) < then Solution test 8: S ← x i ; 9: else 10: x ← checkAndRemove(x i , , f)

Numerical experiments
We compare interval and interval union Newton methods for univariate functions using a Java implementation that is part of JGloptlab [6]. We used 32 test functions listed in Table 1 most of them taken from [4]. For each function we consider the natural extensions for both f and f . In our implementation, we have followed the pseudo-codes of Algorithms 1 and 2 precisely, without any additional acceleration or optimization.
For each function f i we seek the enclosure of all solutions f i (x) = 0 where x ∈ x and x is a bounded interval. The narrow component tolerance is set to 10 −7 and the maximum number of function evaluations is set to 100,000. If we are unable to reduce the width of every component of the solution set below before the maximum number of function evaluations is reached we relax the tolerance parameter by a factor of 10 and restart the process. Table 1 shows the test functions while Tables 2 and 3 present the results of the experiment. A supplementary table comparing other aspects of both algorithms can be found in http://www.mat.univie.ac.at/~dferi/research/UnionsTests.
The test results are given for both the interval Newton (Algorithm 1 in column INewton) and for the interval union Newton (Algorithm 2 in column IUNewton).
In particular, Table 2 shows for each test function (func) the number of boxes possibly containing solutions found (Sol), the number of function evaluations needed to enclose all solutions (FunEv) and the narrow component tolerance (Wid) used.
It is clear from that interval union arithmetic significantly increases the efficiency of the Newton method. Table 2 shows that both the number of function evaluations and the number of boxes possibly containing solutions is smaller when using the interval union Newton method. Moreover, the tolerance achieved with the interval union method is, in every case, at least as small as the tolerance achieved with interval Newton method. Table 3 shows that the interval union Newton method is faster than the interval approach in 37 % of instances and requires less storage memory in 53 % of cases.

Systems of interval union equations
This section extends the concept of interval linear systems to interval unions. The algorithms used to solve interval linear systems can be naturally adapted to the interval union case with a few modifications. The basic definitions of interval union linear systems are given in Sect. 5.1, the Gaussian elimination and the Gauss-Seidel algorithm Table 1 The test functions f 1f 32

Basics
Let A ∈ U n×n be an interval union matrix and b ∈ U n an interval union vector. An interval union linear system of equations is the family of linear systems given bỹ The solution set of interval union linear systems is the union of solution sets from every combination of interval matrices and vectors contained in A and b, formally we have

Definition 15
The set S := {x ∈ R n |Ãx =b for someÃ ∈ A andb ∈ b} is the solution set of (17).
If A ∈ U n×n 1 and b ∈ U n 1 then problem (17) reduces to a typical interval linear system. Finding the interval hull of the solution set is N P−Hard for general interval linear systems and therefore it is also N P−Hard to find the interval hull of S.
We say that a square interval matrix A is regular if every matrix A ∈ A is nonsingular. In the same way, the interval union matrix A is regular if every real matrix A ∈ A with A ∈ A is non-singular. The interval matrix A ∈ U n×n 1 is diagonally Table 2 Comparison between the interval and the interval union Newton methods The interval union matrix A is diagonally dominant if relation (18) remains valid when we replace interval operations with interval union operations. In general, algorithms for solving interval linear systems of equations benefit greatly from preconditioning. We say that the interval linear system A x = b is preconditioned if where M is a real matrix. Typically M =Ǎ −1 is chosen, but some authors suggests better strategies for choosing M, see for example [17]. Similarly, algorithms for solving interval union linear systems may also take advantage of preconditioning, however, the choice of the preconditioning matrix is harder than in the interval case. The study of this topic will be addressed in a future work.

Algorithms
Let A be an interval union matrix and b an interval union vector. We present two methods to enclose the solution set S given by Definition 15. The algorithms discussed here can be easily generalized to the case where A is not square.
Interval Gaussian elimination, as described in [17], is obtained by just replacing real operations with interval ones in the Gaussian elimination algorithm. The interval version of the algorithm also allows to perform partial or full pivoting using the mignitude for element comparison. As proved in [17], the fundamental theorem of interval arithmetic guarantees that if x is the interval vector obtained with interval Gaussian elimination then S ⊆ x. Since the fundamental theorem of interval union arithmetic is already proved, the same conclusion holds if we replace all real operations with interval union counterparts in the Gaussian elimination. Moreover, the definition of the mignitude for interval unions allows the same pivoting strategies as in the interval case.
Consider A and b of the form A = a 11 a 12 a 21 a 22 The interval union Gaussian elimination with backward substitution and without pivoting gives q = −a 21 a 11 , It is trivial to generalize the Gaussian elimination to higher dimensions, but the two dimensional case is good enough to show some interesting properties of the Gaussian elimination applied to interval union systems. Let us first assume that every entry of (19) is an interval instead of an interval union. In this case, if 0 ∈ a 11 then the interval Gaussian elimination will fail even with extended division. However, as demonstrated on Example 3 below, using interval union arithmetic we may obtain useful bounds for x 1 and x 2 even if 0 ∈ a 11 .
Even for systems with 0 / ∈ a 11 the union Gaussian elimination may give us sharper bounds for x 1 and x 2 than the interval Gauss-Seidel algorithm. This is demonstrated on Example 4 below, where by using interval union Gaussian elimination we obtain bounds almost as sharp as solving several interval linear sub-systems separately.
During the interval Newton method, in each iteration, we have to solve an interval linear system of form A(x − x) = b where x is the box currently processed, A is the interval matrix given by evaluating the Jacobian of the function f over x and b is usually set to −f(x). The usual approach this system is the interval Gauss-Seidel algorithm that is based on the so called Gauss-Seidel operator where The interval Gauss-Seidel algorithm applies Eq. (20) as long the bounds of the processed box are improved. In practice, we iterate as long as the difference between the largest widths of x k+1 and x k is bigger than a given tolerance , see Algorithm 4. Note that Algorithm 4 does not update the variables x i when 0 ∈ a ii . When this happens several authors (see [8,17]) suggest a second step of the Gauss-Seidel algorithm which is based on the extended interval division (1). The second step consists of applying Eq. (20) to all indices i for which 0 ∈ a ii and then save the largest gap produced by the interval division. Then two boxes that are identical in every entry except for the one with the largest gap are returned.
Based on Algorithm 4 the interval union version of the Gauss-Seidel elimination can be formulated, where the interval union version of the Gauss-Seidel operator (20) is applied to every equation. The interval union Gauss-Seidel procedure differs from Algorithm 4 in steps 1, 5 and 16. Steps 1 and 16 must be modified to use the component-wise interval union midpoint instead of the interval midpoint, since this is necessary in order to guarantee that the interval union fundamental theorem holds for r i .
As a natural consequence, Algorithm 4 with interval unions returns an interval union vector which stores not only the boxes with the largest gap but all gaps. This simple modifications can lead to significant improvements over the interval Newton procedures for multivariate functions.

Algorithm 4 Interval Gauss-Seidel
Input: The interval matrix A, the interval vectors b and x and the tolerance ≥ 0 Output: The interval vector y such that S ⊆ y ⊆ x or a proof that x ∩ S = ∅. 1: y ← x and x ←x; 2: while true do 3: for i = 1, . . . , n do 4: if 0 / ∈ a ii then 5: j =i a i j (y j − x j ));

Examples
We conclude the section by showing some advantages of using interval union arithmetic to solve interval or interval union linear systems.
Now as (21) suggests (and shown in Fig. 1, left) S may be split into two disjoint sets, and we see that the Gaussian elimination with interval unions provided useful information about S even though A is not regular. The solution set of Ax = b is the union of each interval linear system A i x = b for i = 1, . . . , 4. Figure 1, right shows the result of applying the interval union Gaussian elimination and the interval union Gauss-Seidel algorithm to Ax = b as well as the interval hull of each interval linear system. Note again that the Gauss-Seidel procedure overestimates the bounds of the interval hull while Gaussian elimination give us a sharp enclosure of the four sets. The reason for this is that every interval matrix A ∈ A is regular and diagonally dominant. Our final example shows how the multivariate interval Newton method can benefit from interval union analysis.

Example 4 Now let
Example 5 Assume that we want to enclose the solution set of We use the Gauss-Seidel algorithm applied to the interval Newton operator and precondition the Jacobian matrix by the inverse of its midpoint as described by Hansen [8]. It gives Despite the significant improvement in the resulting box, the result is still not optimal. Applying the interval union Gauss-Seidel algorithm we have