The multivariate bisection algorithm

The aim of this paper is the study of the bisection method in $\mathbb{R}^n$. In this work we propose a multivariate bisection method supported by the Poincar\'e-Miranda theorem in order to solve non-linear system of equations. Given an initial cube verifying the hypothesis of Poincar\'e-Miranda theorem the algorithm performs congruent refinements throughout its center by generating a root approximation. Throughout preconditioning we will prove the local convergence of this new root finder methodology and moreover we will perform a numerical implementation for the two dimensional case.


Introduction
The problem of finding numerical approximations to the roots of a non-linear system of equations was the subject of various studies, and different methodologies have been proposed between optimization and Newton's procedures. In [4] D. H. Lehmer proposed a method for solving polynomial equations in the complex plane testing increasingly smaller disks for the presence or absence of roots. In other work, Herbert S. Wilf developed a global root finder of polynomials of one complex variable inside any rectangular region using Sturm sequences [12].
The classical Bolzano's theorem or Intermediate Value theorem ensures that a continuous function that changes sign in an interval has a root, that is, if f : [a, b] → R is continuous and f (a)f (b) < 0 then there exists c ∈ (a, b) such that f (c) = 0. In the multidimensional case the generalization of this result is the known Poincaré-Miranda theorem that ensures that if we have f 1 , . . . , f n n continuous functions of the n variables x 1 , . . . , x n and the variables are subjected to vary between a i and −a i , then if f i (x 1 , . . . , a i , . . . , x n )f i (x 1 , . . . , −a i , . . . , x n ) < 0 for all x i then there exists c ∈ [−a i , a i ] n such that f (c) = 0. This result was announced for the first time by Poincaré in 1883 [8] and published in 1884 [9] with reference to a proof using homotopy invariance of the index. The result obtained by Poincaré has come to be known as the theorem of Miranda, who in 1940 showed that it is equivalent to the Brouwer fixed point [5]. For different proofs of the Poincaré-Miranda theorem in the n-dimensional case, see [3,10].  (Poincaré-Miranda theorem). Let K be the cube where ρ ≥ 0 and F = (f 1 , f 2 , . . . , f n ) : K → R n a continuous map on K. Also, let then the mapping F has at least one zero point r = (r 1 , r 2 , . . . , r n ) in K.
Throughout this paper we will recall the opposite signs condition (1.1) as the Poincaré-Miranda property P.M. The aim of this work is to develop a bisection method that allows us to solve the non-linear system of equations F(X) = 0, X = (x 1 , x 2 , . . . , x n ), using the above Poincaré-Miranda theorem. The idea of the algorithm will be similar to the classical one dimensional algorithm: we perform refinements of the cube domain in order to check the sign conditions on the parallel faces. In one dimension it is clear that an initial sign change in the border of an interval produces another sign change in a half partition of it, but in several dimensions we cannot guarantee that the Poincaré-Miranda conditions maintain after a refinement. Even if r is an exact solution, there may not be any such K (for which (1.1) holds). However, J. B. Kioustelidis [2] has pointed out that, forx close to a simple solution (where the Jacobian is nonsingular) of F(X), Miranda's theorem will be applicable to some equivalent system for suitable K. Therefore in case of a fail in the sign conditions with the original system, we should try to transform it. The idea will be to find an equivalent system through non-linear preconditioning where the equations are better balanced in the sense that the new system could be close to some hyperplane in order to improve the chances to check the sign conditions in some member of the refinement.
We will denote the infinity norm by x ∞ = max{|x 1 |, . . . , |x n |}, the Euclidean norm by x 2 = x 2 1 + · · · + x 2 n , and the 1-norm by x 1 = |x 1 | + · · · + |x n |. Given a vector norm on R n , the associated matrix norm for a matrix M ∈ R n×n is defined by It is known that in the case of the ∞-matrix norm it can be expressed as a maximum sum of its rows, that is if M = (m ij ) then M ∞ = max i=1,...,n n j=1 |m ij |, therefore it is easy to see that a sequence of matrices (M k ) k converge if and only if their coordinates converge. Since the domains involved are multidimensional cubes, the most proper norm to handle the distance will be the ∞-norm.
We will accept r as a root with a small tolerance level δ if F(r) p ≤ δ.

The algorithm and its description
This section gives a step-by-step description of our algorithm, whose core lies in the classical bisection algorithm in one dimension. Definition 2.1. A 2 n -refinement of a cube K ⊂ R n is a refinement into 2 n congruent cubes Q = {K 1 , K 2 , . . . , K 2 n }.  We say that a 2 n -refinement of Q satisfies the Poincaré-Miranda condition if there exist K l ∈ Q such that F : K l → R n satisfies the condition of Theorem 1.1.
Given a system F(X) = 0, the preconditioned system is G(X) = M F(X) = 0 for some matrix M such that the jacobian at X 0 satisfies DG(X 0 ) = Id. Since DG(X 0 ) = M DF(X 0 ) it turns out that M = DF(X 0 ) −1 and it is clear that the preconditioned system is an equivalent system of F and both have the same roots. After preconditioning, the equations in G(X) = 0 are close to a hyperplane having equation x i = k i , where k i is some constant. This fact comes from the Taylor expansion of G around X 0 , indeed if X is close to X 0 then and therefore it is clear that the equations are close to some hyperplane. Moreover if X 0 is nearly a zero point of F then G(X) ≈ X − X 0 and therefore it will behave like the components of X −X 0 , and take nearly opposite values on the corresponding opposite faces of the cube.
] be the quarter of Q 1 where the conditions of Theorem 1.1 are satisfied, we chose the center of K 1 . If Q 1 does not satisfy the Poincaré-Miranda condition we preconditioning the system in c 1 setting and then we check again the sign conditions with the preconditioned system This recursion is repeated while the Poincaré-Miranda condition are satisfied, generating a sequence of equivalent systems and a decreasing cube sequence K k , such that where the vertices satisfy for each j = 1(1)n and where the length of the current interval [a j k , b j k ] is a half of the last iteration, The root's approximation after the k-th iteration will be and the method is stopped until the zero's estimates gives sufficiently accuracy or until the Poincaré-Miranda condition no longer holds.
Remark 2.2. The kth-preconditioning system G k (X) = (g 1 k (X), . . . , g n k (X)) can be expressed as DF(c k ) −1 F(X), indeed, by induction suppose that it is true for k − 1, then differencing and evaluating in c k we have Since we cannot always ensure that a refinement of a given cube will satisfy the Poincaré-Miranda condition, we cannot ensure the converge for any map that only has a sign change in a given initial cube. So, in case of a fail in the sign conditions in some step, we try to rebalance the system using preconditioning in the center of the current box recursion. The preconditioning allows us to increase the chances to be more often in the sign conditions and therefore keep going with the quadrisection procedure in order to get a better root's approximation. In [2], J. B. Kioustelidis found sufficient conditions for the validity of the Poincaré-Miranda Miranda condition for preconditioning system, there it was proved that the sign conditions are always valid if the center of the cube K is close enough to some root of F. So, if we start the multivariate bisection algorithm with an initial guess close to some root, Kioustelidis's theorem will guarantee the validity of Poincaré-Miranda in each step of our method allowing the local convergence of it.
In the next theorem we will prove the local convergence for the multivariate bisection algorithm when we preconditioning at each step.
] with ρ small enough satisfying the Poincaré-Miranda sign condition; assume that DF(X) is invertible for all X ∈ K 0 ; furthermore, suppose that we perform the preconditioning at each step. Then the multivariate bisection algorithm generates a sequence c k such that Proof.
(1) The Poincaré-Miranda sign conditions guarantee the existence of a root inside K 0 and given a refinement Q 1 of K 0 since ρ is small enough Item c of Theorem 2 in [2] guarantees the validity of Poincaré-Miranda sign conditions for a member of Q 1 . Performing successive refinements we will always find a member K k of the refinement Q k satisfying the sign conditions for the preconditioned system G k (X). For each j = 1(1)n the sequences (a j k ) k , (b j k ) k are monotone and bounded and therefore they converge. From equation (2.3) we have for each j = 1(1)n, and from the border conditions, Since the diameter of K k tends to zero by Cantor's intersection theorem the intersection of the K k contains exactly one point, and the equations (2.4) guarantee that p = (r 1 , . . . , r n ). Then, we can evaluate equations (2.5) in p = (r 1 , . . . , r n ), getting It is clear that then by the continuity of DF and the continuity of the inversion in the ∞-matrix norm we have for each X ∈ K 0 then we get the punctual convergence for each coordinate function From equations (2.4) we have . . , r n ) for each j = 1(1)n; therefore taking limit in equations (2.6) we get and finally it is clear that F(r 1 , . . . , r n ) = 0. (2) Let (c j k ) k , j = 1(1)n, be the coordinates of the sequence (c k ) k ; we have the following estimation: Indeed, since the sequences (a j k ) and (b j k ) are monotone and bounded by r j we get for each j = 1(1)n, On the other hand, Therefore, As in the classical one dimensional bisection algorithm, Item 2 of Theorem 2.3 gives a way to determine the number of iterations that the bisection method would need to converge to a root to within a certain tolerance. The number of iterations needed, k, to achieve the given tolerance δ is given by 3. Implementation, performance, and testing Throughout this section we will focus in the implementation and performance of the bisection algorithm in R 2 . The bisection algorithm was developed in Matlab in a set of functions running from a main function. In order to check the P.M. conditions for the function F = (f 1 , f 2 ) we need to compute the intervals , 2) and one way to achieve this is by using interval analysis (IA). IA was marked by the appearance of the book Interval Analysis by Ramon E. Moore in 1966 [6] and it gives a fast way to find an enclosure for the range of the functions. A disadvantage of IA is the well-known overestimation. If intervals are available then the P.M. follows from the condition sup{y :  Various interval-based software packages for Matlab are available; we have chosen the well-known INTLAB toolbox [11]. The toolbox has several interval class constructor for intervals, affine arithmetic, gradients, hessians, slopes and more. Ordinary interval arithmetic has sometimes problems with dependencies and wrapping effect given large enclosures of the range and therefore overestimating the sign behaviour. A way to fight with this is affine arithmetic. In affine arithmetic an interval is stored as a midpoint X 0 together with error terms E 1 , . . . , E k and it represents where U 1 , . . . , U k are parameters independently varying within [−1, 1]. In order to avoid an overestimation in the range enclosure of f i (F − i ) and f i (F + i ) we also compute the interval extension using the affine arithmetic.
Another way to improve the enclosure of the range and get sharper lower and upper bounds is through subdivision or refinements. In this methodology we perform subdivision of the domain and then we take the union of interval extensions over the elements of the subdivision; this procedure is called a refinement of [f ] over X. Let N be a positive integer; we define We have X i = ∪ N j=1 X i,j and w(X i,j ) = w(X i ) N and furthermore, The interval quantity is the refinement of [f ] over X.
The algorithms that we have performed to compute equations (3.3) and (3.4) combine all the above methodologies and were adapted from [7]. In the following steps we summarize the routines that we have performed. The mean value extension was implemented using an approximation of [ ∂fi ∂xi ] through the central finite difference of the natural interval extension of f , that is Algorithm 1 shows the routine for the mean value extension.

Algorithm 1: Function meanValue(), computes the mean value extension
Data: f, X Result: returns F mv the value for the mean value extension form for f evaluated over the interval X.
The refinement procedure was implemented twice, one for the case of mean value extension and other for the affine arithmetic implementation. Algorithm 2 computes the mean value extension over an uniform refinement of the interval X with N subintervals and Algorithm 3 computes the natural extension using affine arithmetic.

MANUEL LÓPEZ GALVÁN
Algorithm 2: Function meanValueRefinement(), computes the refinement procedure using mean value extension Data: f, X, N Result: returns Y the value for the mean value extension form for f evaluated over a partition of X.
Xs(i) ← infsup(xi,xip1) // Interval class constructor for each subinterval.   In order to check the accuracy and performance of the algorithm, we test it through different systems of equations. We start testing the new algorithm in a transcendental system of equations used by Broyden in [1], The initial guess rectangle used for this problem was K 0 = [0.4, 0.55]×[3, 3.5] and the tolerance level was setted in δ = 10 −15 . The interval analysis refinement used to compute the sign along the edges is N = 3. Figure 3 illustrates the algorithm behaviour with the respective preconditioning procedure and Table 1 shows the numerical solution.  3.141592653589793 -0.643347227536409 1e-15 Table 1. Root's approximation, evaluation, iteration performed and initial guess for Broyden system.
A straightforward computation shows that the solution for the system is (x 1 , x 2 ) = (0.5, π) and therefore we can check the consistency of the error estimation given by Theorem 2.3 at each iteration,  In the following steps we test the algorithm in several other problems. We will see that in some systems the algorithm needs preconditioning in order to guarantee the P.M. conditions through the refinement. Let be the testing maps. In Table 2 we show the numerical performance for the testing maps. The method was implemented setting the tolerance level in δ = 10 −15 and the interval analysis refinement in N = 3.  Table 2. Root's approximation, evaluation, iteration performed and initial guess for testing maps. Figure 5 illustrates the algorithm behaviour for the testing maps with the refinement procedure. The systems of equations and their successive possible preconditionings are represented by a zero contour level on an mesh on the initial guess K 0 and the refinement procedure was illustrated using the rectangle Matlab's functions.

Newton's comparison
In this section we are going to compare the performance of our method against the classical multivariate Newton algorithm. It is well-known that the classical Newton's method has several numerical problems when the systems are illconditioned, i.e., systems having a "nearly singular" Jacobian at some iterate, getting slower rate of convergence and large numerical errors. The main advantage of our methodology is that it is not always necessary to perform the preconditioning at each step and therefore it can skip the ill-conditioned problem.
Let us consider the following system, F(x, y) = (x 2 + y 2 − 1, xy − x 2 ). This function has a zero in x * = ( ) and a straightforward computation shows that the system is ill-conditioned for any value sufficiently close to the origin. Newton's method cannot be initialized in (0, 0), however our algorithm can be evaluated in (0, 0) avoiding the singularity of the Jacobian.  Figure 5. The first row illustrates the algorithm procedure for F 1 , F 2 , the second for F 3 , F 4 and the third for F 5 , F 6 . The solid red line represents the first coordinate, while the dashed blue line represents the second coordinate for the equivalent system G k (X) = 0. Figure 6 and Table 3 show the behaviour of the reduction of error norms c k −x * and F(c k ) through iteration steps for the Bisection and Newton methods.

Newton Bisection
Step  Table 3. Error norms F(c k ) , c k − x * for the Newton and the Bisection methods after the iteration step for both cases.
The initial guess used to perform the Newton's method was (10 −8 , 10 −9 ) which is inside and very close to the vertex of the initial square K 0 = [0, 1] × [0, 1] used to perform the Bisection method. As depicted in Figure 6 the Bisection method is highly superior to the Newton method in respect to achieve convergence. Bisection's method achieves a norm error of 1,31e-09 after 28 iteration, however Newton's method achieves a norm error of 1,22e-01 demanding much higher numerical effort.
As a second example, we try the system F(x, y) = (x 2 − 4y + y 2 − 1, 2x − y 2 ). In Figure 7 and Table 4 we compare the Newton method and the Bisection; the starting square was K 0 = [1.0000001, 1.98] × [1,2] and the initial guess for Newton was the lower left vertex (1.0000001, 1) of K 0 .  2,33e+00 4,06e-08 Table 4. Error norms F(c k ) for the Newton and the Bisection method after the iteration step for both cases.
The condition number of the Jacobian matrix is close to 4e7 and therefore we are present to a ill-conditioned problem induced by a bad initial guess given by (1.0000001, 1). However, the bisection method can start in this same vertex and can skip the Jacobian illness without preconditioning in the first steps. Bisection's method achieves a norm error of 4,06e-08 after 24 iteration, however Newton's method achieves a norm error of 2,33e+00.

Conclusion
In this work we have clarified how a multidimensional bisection algorithm should be performed extending the idea of the classic one dimensional bisection algorithm. Due to the preconditioning at each step we could prove a local convergence theorem and we also found an error estimation. Interval analysis allowed a fast and reliable way of computing the Poincaré-Miranda conditions and the numerical implementation showed that the method has a very good accuracy similar with the classic methods like Newton or continuous optimization. We also have compared the performance of the Bisection method against the classical Newton method and we found that our methodology improves the speed of convergence in ill-conditioned problems.