Effective Topological Degree Computation Based on Interval Arithmetic
We describe a new algorithm for calculating the topological degree where is a product of closed real intervals and is a real-valued continuous function given in the form of arithmetical expressions. The algorithm cleanly separates numerical from combinatorial computation. Based on this, the numerical part provably computes only the information that is strictly necessary for the following combinatorial part, and the combinatorial part may optimize its computation based on the numerical information computed before. We present computational experiments based on an implementation of the algorithm. In contrast to previous work, the algorithm does not assume knowledge of a Lipschitz constant of the function , and works for arbitrary continuous functions for which some notion of interval arithmetic can be defined.
The notion of topological degree was introduced by Jan Brouwer  and was motivated by questions in differential topology [26, 19]. The degree of a continuous function is an integer, describing some topological properties of it. Degree theory has many applications, including geometry , nonlinear differential equations [24, 11, 25, 6], dynamical systems , verification theory , fixed point theory  and others.
The presented algorithm is able to calculate the degree of any real-valued continuous function defined on a box such that and is given in the form of arithmetical expressions containing function symbols for which interval enclosures can be computed [28, 34]. Computational experiments show that for low-dimensional examples of simple functions (up to dimension 10) the algorithm terminates in reasonable time.111The program is accessible on topdeg.sourceforge.net. In addition to efficiency, the algorithm has several advantages over previous work that we now describe in more detail.
The idea of computing the degree algorithmically is not new. Since the seventies, many algorithms were proposed and implemented that calculate the degree of a function defined on a bounded set via a symbolical expression. However, all these methods have various restrictions. One of the first such methods was proposed by Erdelsky in 1973 . His assumption is that the function is Lipschitz, with a known Lipschitz constant. Thomas Neil published another method for automatic degree computation in 1975 . It is based on the approximation of a multidimensional integral of a function derived from and its partial derivatives. Here, the error analysis uses only probabilistic methods. Other authors constructed algorithms that cover the boundary with a large set of -simplices and use the information about the signs of on the vertices of these simplices to calculate the degree in a combinatorial way. However, the calculated result is proved to be correct only if a parameter is chosen to be sufficiently large [21, 36, 37]. Boult and Sikorski developed a different method for degree calculation in the eighties, but their algorithm also requires the knowledge of a Lipschitz constant for . Later, many algorithms arose where the degree was calculated recursively from partial information about on the boundary . For example, one has
where and is a -dimensional open neighborhood of in . Aberth described an algorithm using this formula, based on interval arithmetic . This method was not implemented and is rather a recipe than a precise algorithm. Later, Murashige published a method for calculating the degree that uses concepts from computational homology theory .
Although a broad range of ideas and methods for automatic degree computation has been implemented, the effectivity of these algorithms decreases fast with the dimension of . For example, in the Murashige homological method, computation of the degree of the identity function takes more than 100 seconds already in dimension 5 [29, Figure 3]. Other approaches were developed that calculate the degree of high-dimensional examples quickly, provided the functions are of some special type. For instance, there exist effective degree algorithms for complex functions [10, 22, 23].
Our approach is based on a formalization, extension, and implementation of the rough ideas of Oliver Aberth . In our setting, we assume that the function is real valued and continuous, and it is possible to implement an interval-valued function which computes box enclosures for the range of over a box. We don’t require the function to be differentiable and not even Lipschitz. This enables us to work with algebraic expressions containing functions such as , and , but also with any function that cannot be defined by algebraic expressions and only an algorithm is given that computes a superset of for any interval s.t. the measure of can be arbitrary small for small intervals . Throughout the paper, we assume that the domain of the function is a box (product of compact intervals), but the algorithm works without major changes for more general domains, such as finite unions of boxes with more complicated topology. This will be discussed at the end of Section 3.3.
From the algorithmic point of view, our algorithm consists of a numerical part, that provably computes only information that is strictly necessary for determining the degree, and a combinatorial part that computes the degree from this information. The separation of those two parts has the advantage that both can be used and improved independently. The first, numerical part covers the boundary of a -dimensional set with -dimensional regions where a particular component of has constant sign. The combinatorial part recursively gathers the information about the signs of the remaining components of on . All the sets are represented as lists of oriented boxes. They do not have to represent manifolds and we allow the boundary of these sets to be complicated (see Def. 2.4). In this setting, it is computationally nontrivial to identify the boundary of a -dimensional set embedded in and to decompose the boundary into a sum of “nice” sets. Instead of doing this, we calculate an “over-approximation” of that is algorithmically simpler and then prove that it has no impact on the correctness of the result. This involves some theoretical difficulties whose solution necessitates the development of several technical results.
Some interest in automatic degree computation is motivated by verification theory. Methods have been developed for automatic verification of the satisfiability of a system of nonlinear equations in variables, written concisely as , where is a continuous function. Most of these methods first find small boxes that potentially contain a root of and then try to formally prove the existence of a root in such a box [27, 33, 16, 15] using tests based on theorems such as the Kantorovich theorem, Miranda theorem, or Borsuk theorem. From those, the test based on Borsuk theorem is the most powerful [2, 15]. It can be easilly shown that the assumptions of Miranda theorem imply that and the assumption of Borsuk theorem imply that the degree is an odd number. It is well known that implies the existence of a root of in . An efficient test developed by Beelitz can verify that the degree is , if it is , and hence prove the existence of a solution . By not restricting oneself to degree but computing the degree in general, one can prove the existence of a root of in all cases that are robust in a certain sense [8, 14].
The second section contains the main definitions needed from topological degree theory—Theorem 2.9 is a fundamental ingredient of our algorithm. Section 3 describes the algorithm itself and its connection to Theorem 2.9. In Section 4, we present some experimental results. The last section contains the proof of two auxiliary lemmas that we need throughout the paper. These proofs do not involve deep ideas but are quite long and technical—hence the separate section at the end of the paper.
2 Mathematical Background
2.1 Definitions and Notation
In this section, we first summarize the definition and main characteristics of the topological degree on which there exists a wide range of literature, such as [13, 32]. Degree theory works with continuous maps between oriented manifolds, and in order to represent these topological objects on computers we will then introduce Definitions 2.1 to 2.6. Finally, the original Theorem 2.9 will be the main ingredient of our algorithm for computing the topological degree.
Let be open and bounded, continuous and smooth (i.e., infinitely often differentiable) in , . For regular values (i.e., values such that for all , ), the degree is defined to be
This definition can be extended for non-regular values in a unique way, such that for given and , —as a function in —is locally constant on the connected components of .
Here we give an alternative, axiomatic definition, that determines the degree uniquely. For any continuous function s.t. the degree is the unique integer satisfying the following properties [13, 32, 17]:
For the identity function , iff is in the interior of .
If then has a solution in .
If there is a continuous function (a “homotopy”) such that , then .
If and , then .
For given and , —as a function in —is constant on any connected component of .
This can be generalized to the case of a continuous function , where and are oriented manifolds of the same dimension and is compact. If is smooth, denotes the matrix of partial derivatives of some coordinate representation of and formula is still meaningful. For example, if is a scalar valued function from an oriented curve (i.e., an oriented set of dimension ) to and on the endpoints of , then is well defined. If is a function between two oriented manifolds without boundary, then the degree is defined to be for any .
A simple consequence of the degree axioms is that for a continuous , implies that . So we will be only interested in calculating .
We will represent geometric objects like manifolds, orientation, boundaries and functions in a combinatorial way, using the following definitions.
A -dimensional box (simply -box) in is the product of non-degenerate closed intervals and degenerate intervals (one-point sets). A sub-box of a -box is any -box s.t. . The diameter of a box is the width of its widest interval.
The orientation of a -box is a number from the set . An oriented box is a pair where is a box and its orientation. We say that is an oriented sub-box of an oriented box , if , the dimensions of and are equal and the orientations are equal.
Let be an oriented -box in with orientation . Let, for every , . Assume that the intervals are non-degenerate, , the other intervals are degenerate (one-point) intervals. For , the -dimensional boxes
are called faces of . Any sub-box of a face is called a sub-face of . If we choose the orientation of to be and the orientation of to be , then we call oriented faces of . An oriented sub-box of an oriented face is called oriented sub-face. The orientation of the oriented faces and sub-faces is called the induced orientation from the orientation of .
An oriented cubical set is a finite set of oriented boxes of the same dimension such that the following conditions are satisfied:
For each , the dimension of is at most .
Whenever is a -dimensional box, then the orientations of and are compatible. This means that has an opposite induced orientation as a sub-face of as the orientation induced from .
The dimension of an oriented cubical set is the dimension of any box it contains. If is an oriented cubical set, we denote by the set it represents (the union of all the oriented boxes contained in ).
An oriented cubical set is sketched in Figure 1. An immediate consequence of the definition is that each sub-face of a box in an oriented cubical set is a boundary sub-face of at most two boxes in .
Note that an oriented cubical set does not have to represent a manifold, because some boxes may have lower-dimensional intersection, like and in Figure 1.
An oriented boundary of an oriented -dimensional cubical set is any set of -dimensional oriented boxes , such that
Any two boxes in have intersection of dimension at most .
For each , and each -dimensional sub-box of , there exists exactly one box such that is an oriented sub-face of .
An oriented cubical set and its oriented boundary are sketched in Figure 2. Geometrically, this definition describes the topological boundary of an oriented cubical set and we denote the union of all oriented boxes in by . Clearly, , the meaning of the left hand side being the topological boundary of the set . Note that if is a -dimensional oriented cubical set and an oriented boundary of , then each sub-face of some box in s.t. is at most -dimensional, is a sub-face of exactly two boxes in with opposite induced orientation (see in Figure 1).
An oriented boundary of an oriented cubical set does not have to form an oriented cubical set, because the second condition of Definition 2.4 may be violated (for a counter-example, see Fig. 3 where the -boxes and have -dimensional intersection but not compatible orientations).
The notion of topological degree can be naturally generalized to oriented cubical sets. So, if is a continuous function from a -dimensional oriented cubical set to such that , then is well-defined, extending the definition of for oriented manifolds . 222For an oriented cubical set , one can define an oriented manifold for a small enough and define the degree to be
Finally, we will represent functions as algorithms that can calculate a superset of for any given box .
Let . We call a function interval-computable if there exists a corresponding algorithm that, for a given box with rational endpoints and positive diameter, computes a closed (possibly degenerate) interval such that
for every there is a such that for every box with , .
We call a function interval-computable iff each is interval-computable. In this case, the algorithm returns a tuple of intervals, one for each .
Usually such functions are written in terms of symbolic expressions containing symbols denoting certain basic functions such as rational constants, addition, multiplication, exponentiation, trigonometric function and square root. Then, can be computed from the expression by interval arithmetic [30, 28]. The interval literature usually calls an interval function fulfilling the first property of Definition 2.6 “enclosure”. Instead of the second property, it often uses a slightly stronger notion of an interval function being “Lipschitz continuous” [30, Section 2.1]. We will use interval computable functions and expressions denoting them interchangeably and assume that for an expression denoting a function , a corresponding algorithm is given.
2.2 Main Theorem
Now we define the combinatorial information we use to compute the degree, and prove that it is both necessary and sufficient for determining the degree.
A -dimensional sign vector is a vector from .
Let be a set of oriented -boxes. A sign covering of is an assignment of a -dimensional sign vector to each . For a sign covering and we will denote this sign vector by , and its -th component by .
A sign covering is sufficient if each sign vector contains at least one non-zero element.
A sign covering is a sign covering wrt. a function with components , if for every oriented box and for every , implies that has constant sign on .
In the following we will often recursively reduce proofs/algorithms for -dimension oriented cubical sets, to proofs/algorithms on their oriented boundary. Since—as we have already seen—an oriented boundary of an oriented cubical set does not necessarily have to form an oriented cubical set, we will need the following lemma that will allow us to decompose this oriented boundary again into oriented cubical sets:
Let be a -dimensional oriented cubical set, an oriented boundary of , a sufficient sign-covering of with respect to and assume that for each , has exactly one nonzero component. Let for each and . Then there exist oriented cubical sets and corresponding oriented boundaries s.t. the following conditions are satisfied:
For each , there exists such that ,
Each is a sub-face of some where .
The lemma is illustrated in Figure 3. The proof of this lemma is technical and we postpone it to the appendix in order to keep the text fluent.
Let be an oriented -dimensional cubical set, an oriented boundary of and a continuous function with components such that . Then a sign covering of wrt. determines the degree uniquely if and only if it is sufficient.
We first prove that sufficiency of the sign covering implies a unique degree. We proceed by induction on the dimension of . If is a -dimensional oriented cubical set , then is determined by the sufficient sign covering of wrt. . Let . For each box , choose an index such that . For all and , let . It follows from Lemma 2.8 that we may decompose into oriented cubical sets and oriented boundaries , such that for unique and each is a sub-face of some where . For each , define . Then and the degree is defined. Let and be arbitrary. It follows from [22, Theorem 2.2] and [36, Section 4.2] that
where is the -dimensional zero.
For each set from the sum on the right hand side, has sign on . Each is a sub-box of some where , so we may assign a new sign vector for by deleting the -th component from . In this way, we define a sufficient sign covering of wrt. and the degree can be then calculated recursively using .
Now assume that the sign covering of is not sufficient. We will prove that in this case, the degree is not uniquely determined.
Let be a -dimensional box such that . Choose to be arbitrary. We will construct a function such that the sign covering of is a sign covering with respect to and .
Denote the oriented manifold with boundary by . is a union of the oriented manifolds and , the boundaries and are equal with opposite orientations, homeomorphic to the sphere . The degree where is a map to the sphere. Let be such that , let and . We construct a map such that . The homotopy group for , so each map from a -sphere to the -sphere is homotopic to a constant map. Let us define on . Then is homotopic to a constant map. There exists a sub-box and a continuous extension of such that on and is constant on . Using the fact that , there exists a map of degree . It follows from the identity that we can extend to a map such that . Finally, extend to a map by on . Then
Let be the inclusion. Multiplying by some scalar valued function, we can obtain a function such that on . Extending to a continuous arbitrarily (this is possible due to Tietze’s Extension Theorem [7, Thm. 4.22],) we obtain a function such that the original sign covering is a sign covering of wrt. and . This completes the proof. ∎
3 Algorithm description
3.1 Informal Description of the Algorithm
We describe now our algorithm for degree computation of an interval computable function. If is an interval computable function nowhere zero on the boundary , then the corresponding interval computation algorithm from Definition 2.6 may be used to construct a sufficient sign covering of wrt. . This sign covering will be represented as a list of oriented boxes and sign vectors. The main ingredient of the algorithm is Equation from the proof of Theorem 2.9. For some index and sign , we select all the boxes with . From Lemma 2.8, we know that these boxes form some oriented cubical sets . Then a new list of -dimensional oriented boxes is constructed that covers the boundaries of . Possibly subdividing boxes in this new list, we assign -dimensional sign vectors to its elements in such a way that we obtain a sufficient sign covering of wrt. . Equation is used for a recursive dimension reduction.
We work with lists of oriented boxes and sign vectors rather than with sets, because it will be convenient for our implementation to allow an oriented box to be contained in a list multiple times. However, we will usually ignore the order of the list elements (i.e., the algorithm actually is based on multi-sets which we implement by lists). For two lists and , we denote by the concatenation of and and will also use the symbol for the concatenation of several lists. We use the notation if is contained in at least once. If is a sub-list of , we denote by the list with the sub-list omitted.
Now we define a version of the notion of sign covering based on lists:
A sign list (of dimension ) is a list of pairs consisting of
an oriented -box, and
a corresponding -dimensional sign vector.
A sign list is sufficient, iff each sign vector contains at least one non-zero element. A sign list of dimension is a sign list wrt. a function iff for each element and corresponding sign vector , for all , implies that has sign on .
By misuse of notation, we will sometimes refer to the elements of a sign list as pairs consisting of an oriented box and a sign vector, and sometimes just as an oriented box.
The basic ingredient of the algorithm is a recursive function with input a sufficient sign list and output an integer. This function involves no interval arithmetic and is purely combinatorial. For an input that is a sufficient sign list wrt. such that the boxes in form an oriented boundary of an oriented cubical set , this function returns . If the function input is a -dimensional sign list , then the output is returned. This is compatible with the the formula for the degree of a function on an oriented edge, .
If the input consists of oriented -boxes and sign vectors of dimension for , we choose and and compute a list of boxes (the selected boxes) having as the -th component of the sign vector. We split the boundary faces of all selected boxes until each face of a selected box is either contained in some non-selected box or the intersection of with each non-selected box is at most -dimensional. For each face of a selected box that is a sub-face of some non-selected box , we delete the -th entry from the sign vector of and assign this as a new sign vector to . The list of all such oriented -boxes and their sign vectors is denoted by . This is a sufficient sign list wrt. and is returned.
The choice of and has no impact on the correctness of the algorithm but can optimize its speed. We choose and in such a way that the number of selected boxes is minimal. See Section 4 for a more detailed discussion of this issue. The algorithm for calculating , and is displayed in Figure 4.
If the input of the Deg function is a sign list representing a boundary of an oriented cubical set, then the list of selected boxes is exactly the set from Lemma 2.8. We will prove in Section 3.3 that the list can be subdivided into where is the decomposition of into oriented cubical sets and contains each box the same number of times as , where represents the box with opposite orientation. We will prove that . Together with equation this implies
One example of a possible construction is displayed in Figure 5.
|:||algorithmic representation of an interval-computable -valued function|
|Output: the degree|
|boundary_info refineCov(, )|
For an interval-computable function (see Definition 2.6) and box , if , we can infer a non-zero sign vector entry (see Definition 3.1) for . Moreover, due to interval-computability, if is small enough and , then . Hence, a function with the following specification can be easily implemented by starting with the list of faces of , using to assign sign vectors to them so that the constructed sign list is wrt. , and recursively splitting the boxes in the list until the interval evaluation computes the necessary sign information for it to be sufficient.
|:||an -box in|
|:||algorithmic representation of an interval-computable -valued function|
|Sufficient sign list wrt. , covering the oriented boundary of .|
Now, finally, we can compute the degree from a sufficient sign covering.
|Input: : Sufficient sign list wrt. some function , covering the oriented boundary of .|
|else if then|
|for all do|
|bound list of the oriented faces of|
|split the boxes in bound until for all , either|
|– is a subset of some element of , or|
|– the intersection of with any element of has dimension smaller than|
|for every that is a subset of some box in do|
|is the sign vector of with omitted -th component|
Note that the input/output specification of the function Deg describes the behavior for calls from the outside. Recursive calls of the function Deg might violate the condition on the input—it might be a more complicated sign list. We will discuss details on the structure of that list and correctness of recursive calls in the following section.
3.3 Proof of Correctness
The algorithm first creates a sufficient sign list wrt. where is the input box. This sign list is then an input for the recursive function . We want to prove that if is a sufficient sign list wrt. covering the boundary of a box , then returns the degree .
To prove this, we will analyze the function body. When dealing with -dimensional sufficient sign lists, we always assume that some and has been chosen. Let be a sufficient sign list wrt. . We denote and the sub-list of selected and non-selected boxes. For each , the function refines the boundary until each is either a subset of some or has at most -dimensional intersection with each . For each that is a subset of a , it assigns to the sign vector with deleted -th coordinate. We denote the sub-list of all such constructed from by . The list constructed in the function body satisfies
and is returned.
In this section, we will suppose that some implementation of the algorithm is given. This includes a rule for the choices of , subdivision of the boundary faces of the selected boxes, order of the lists and and the choice of . We will show that if the sign list satisfies a certain regularity condition defined in Definition 3.3, then the function output is invariant with respect to some changes of the input list, including any change of order, merging and splitting some boxes or adding and deleting a pair of identical boxes with opposite orientation. This is shown in Lemma 3.5. Further, we show that the list can be decomposed into the sum of oriented cubical sets such that and such that the list constructed in the function body is a merging of and a set of pairs , so that . In Theorem 3.6 we prove that and combining this with equation in Theorem 2.9, we show that if is a sufficient sign list wrt. covering the boundary of a box , then returns the degree .
Let and be two sufficient sign lists wrt. . We say that is equivalent to and write , if can be created from by applying a finite number of the following operations:
Changing the order of the list,
Replacing some oriented box in one list by two boxes where is the splitting of into two oriented sub-boxes with equal sign vectors ,
Merging two oriented boxes , , that form a splitting of some box and have the same sign vector , to one list element ,
Adding or deleting a pair of oriented boxes and where is the box with opposite orientation (the sign vectors and do not have to be necessary equal in this case),
Changing the sign vectors of some oriented boxes so that the sign covering is still sufficient and wrt. .
Clearly, is an equivalence relation on sign lists.
Let be a -dimensional sufficient sign list wrt. . We say that is balanced, if each sub-face of some such that for each , either , or is at most -dimensional 333Here and represent just the box, without taking care of the orientation., satisfies
where is a sub-list of containing all s.t. is an oriented sub-face of .
In other words, is a sub-face of some oriented box in the same number of times as .
A sign list representing the oriented boundary of an -box is clearly balanced, because for each -dimensional sub-face of some that is small enough to have either lower-dimensional or full intersection with each , is an oriented sub-face of exactly one and is an oriented sub-face of exactly one . The following Lemma says that the property of being balanced is also preserved in the construction procedure. This implies that all input lists within the recursive function are balanced.
Let be a sufficient sign list wrt. that is balanced. Then the list created in the Deg function body is also balanced.
The proof of this is technical and we postpone it to the appendix.
Let be a balanced sufficient sign list wrt. and be equivalent to . Then .
We prove this by induction on the dimension of the sign list. If is a -dimensional sign list, then nontrivial merging and splitting of a box is impossible. Independence of order of the list follows from the formula and adding a pair to the list will add to the sum .
Assume that the lemma holds up to dimension . Let be a permutation (i.e. the same multiset with different order of elements) of a -dimensional sign list and be the list created for in the function body. Changing the order of the list possibly changes the order of and . However, if and only if and the same number of times. Further, and can be constructed from each other by a finite number of splitting, merging and sign vector changing operations, because both are sufficient sign list wrt. representing a sign covering of the set
So, and .
Further, let be created from by splitting or merging some oriented box and , resp. the list constructed in the Deg function body. If we split or merge a non-selected box, then will be equivalent to , because the equivalence class of depends only on the union of all non-selected boxes. Splitting a selected box into will result in splitting some elements of , possibly changing their sign-vectors (depending on the choice of in the algorithm) compatibly with and generate a finite number of new pairs and s.t. and . So, is again equivalent to and we can apply the induction.
Assume that we change the sign vector of an element in in such a way that we still have a sufficient sign list wrt. . If we change the sign vector of a box such that we don’t change a selected box to a non-selected or vice versa, then this change may only result in a possible change of sign vectors in wrt. (and possibly splitting and merging of the boxes in , if the sign vector change has an impact on the choice of in the algorithm). So, in this case, . Assume that we change the sign vector so that some will become selected. Denote to be the original sign list () and to be the new sign list in which and let , resp. be the corresponding sign lists created in the function body. First note that the sublists of containing all elements that are not sub-faces of and the sublist of containing all elements that are not sub-faces of , are equivalent, so we only have to analyze the changes caused by the changed sign-vector of . We claim that the sign list is equivalent to , where is a sign list covering a boundary of such that all are endowed with the old sign vectors with -th entry deleted. An implementation of the Deg function body will create, in the construction, a decomposition , where each oriented box in has at most -dimensional intersection with each and each oriented box in is a subset of some . It follows that each is contained in and the list contains one more time than . Further, due to the fact that is balanced, for each , there exist the same number of boxes in s.t. is an oriented sub-face of as boxes s.t. is an oriented sub-face of , being among the ’s. All such and ’s are in , being the only of these boxes contain in . This implies that the list is equivalent to a list containing each such one more time than . After deleting a finite number of pairs , is equivalent to a list containing one for each (it comes from for some containing a sub-face of . In , there is no such , because is not contained in any . Summarizing this, can be constructed from be adding a sign list covering and deleting a sign list covering . This is equivalent to adding a sign list covering all and we obtain that