Real-root isolation

Short description: Methods for locating real roots of a polynomial

In mathematics, and, more specifically in numerical analysis and computer algebra, real-root isolation of a polynomial consist of producing disjoint intervals of the real line, which contain each one (and only one) real root of the polynomial, and, together, contain all the real roots of the polynomial.

Real-root isolation is useful because usual root-finding algorithms for computing the real roots of a polynomial may produce some real roots, but, cannot generally certify having found all real roots. In particular, if such an algorithm does not find any root, one does not know whether it is because there is no real root. Some algorithms compute all complex roots, but, as there are generally much fewer real roots than complex roots, most of their computation time is generally spent for computing non-real roots (in the average, a polynomial of degree n has n complex roots, and only log n real roots; see Geometrical properties of polynomial roots § Real roots). Moreover, it may be difficult to distinguish the real roots from the non-real roots with small imaginary part (see the example of Wilkinson's polynomial in next section).

The first complete real-root isolation algorithm results from Sturm's theorem (1829). However, when real-root-isolation algorithms began to be implemented on computers it appeared that algorithms derived from Sturm's theorem are less efficient than those derived from Descartes' rule of signs (1637).

Since the beginning of 20th century there is an active research activity for improving the algorithms derived from Descartes' rule of signs, getting very efficient implementations, and computing their computational complexity. The best implementations can routinely isolate real roots of polynomials of degree more than 1,000.

Specifications and general strategy

For finding real roots of a polynomial, the common strategy is to divide the real line (or an interval of it where root are searched) into disjoint intervals until having at most one root in each interval. Such a procedure is called root isolation, and a resulting interval that contains exactly one root is an isolating interval for this root.

Wilkinson's polynomial shows that a very small modification of one coefficient of a polynomial may change dramatically not only the value of the roots, but also their nature (real or complex). Also, even with a good approximation, when one evaluates a polynomial at an approximate root, one may get a result that is far to be close to zero. For example, if a polynomial of degree 20 (the degree of Wilkinson's polynomial) has a root close to 10, the derivative of the polynomial at the root may be of the order of [math]\displaystyle{ 10^{20}; }[/math] this implies that an error of [math]\displaystyle{ 10^{-10} }[/math] on the value of the root may produce a value of the polynomial at the approximate root that is of the order of [math]\displaystyle{ 10^{10}. }[/math] It follows that, except maybe for very low degrees, a root-isolation procedure cannot give reliable results without using exact arithmetic. Therefore, if one wants to isolate roots of a polynomial with floating-point coefficients, it is often better to convert them to rational numbers, and then take the primitive part of the resulting polynomial, for having a polynomial with integer coefficients.

For this reason, although the methods that are described below work theoretically with real numbers, they are generally used in practice with polynomials with integer coefficients, and intervals ending with rational numbers. Also, the polynomials are always supposed to be square free. There are two reasons for that. Firstly Yun's algorithm for computing the square-free factorization is less costly than twice the cost of the computation of the greatest common divisor of the polynomial and its derivative. As this may produce factors of lower degrees, it is generally advantageous to apply root-isolation algorithms only on polynomials without multiple roots, even when this is not required by the algorithm. The second reason for considering only square-free polynomials is that the fastest root-isolation algorithms do not work in the case of multiple roots.

For root isolation, one requires a procedure for counting the real roots of a polynomial in an interval without having to compute them, or, at least a procedure for deciding whether an interval contains zero, one or more roots. With such a decision procedure, one may work with a working list of intervals that may contain real roots. At the beginning, the list contains a single interval containing all roots of interest, generally the whole real line or its positive part. Then each interval of the list is divided into two smaller intervals. If one of the new intervals does not contain any root, it is removed from the list. If it contains one root, it is put in an output list of isolating intervals. Otherwise, it is kept in the working list for further divisions, and the process may continue until all roots are eventually isolated

Sturm's theorem

The first complete root-isolation procedure results of Sturm's theorem (1829), which expresses the number of real roots in an interval in terms of the number of sign variations of the values of a sequence of polynomials, called Sturm's sequence, at the ends of the interval. Sturm's sequence is the sequence of remainders that occur in a variant of Euclidean algorithm applied to the polynomial and its derivatives. When implemented on computers, it appeared that root isolation with Sturm's theorem is less efficient than the other methods that are described below. Consequently, Sturm's theorem is rarely used for effective computations, although it remains useful for theoretical purposes.

Descartes' rule of signs and its generalizations

Descartes' rule of signs asserts that the difference between the number of sign variations in the sequence of the coefficients of a polynomial and the number of its positive real roots is a nonnegative even integer. It results that if this number of sign variations is zero, then the polynomial does not have any positive real roots, and, if this number is one, then the polynomial has a unique positive real root, which is a single root. Unfortunately the converse is not true, that is, a polynomial which has either no positive real root or has a single positive simple root may have a number of sign variations greater than 1.

This has been generalized by Budan's theorem (1807), into a similar result for the real roots in a half-open interval (a, b]: If f(x) is a polynomial, and v is the difference between of the numbers of sign variations of the sequences of the coefficients of f(x + a) and f(x + b), then v minus the number of real roots in the interval, counted with their multiplicities, is a nonnegative even integer. This is a generalization of Descartes' rule of signs, because, for b sufficiently large, there is no sign variation in the coefficients of f(x + b), and all real roots are smaller than b.

Budan's may provide a real-root-isolation algorithm for a square-free polynomial (a polynomial without multiple root): from the coefficients of polynomial, one may compute an upper bound M of the absolute values of the roots and a lower bound m on the absolute values of the differences of two roots (see Properties of polynomial roots). Then, if one divides the interval [–M, M] into intervals of length less than m, then every real root is contained in some interval, and no interval contains two roots. The isolating intervals are thus the intervals for which Budan's theorem asserts an odd number of roots.

However, this algorithm is very inefficient, as one cannot use a coarser partition of the interval [–M, M], because, if Budan's theorem gives a result larger than 1 for an interval of larger size, there is no way for insuring that it does not contain several roots.

Vincent's and related theorems

Vincent's theorem (1834) provides a method for real-root isolation, which is at the basis of the most efficient real-root-isolation algorithms. It concerns the positive real roots of a square-free polynomial (that is a polynomial without multiple roots). If [math]\displaystyle{ a_1, a_2,\ldots, }[/math] is a sequence of positive real numbers, let

[math]\displaystyle{ c_k=a_1 + \cfrac{1}{a_2 + \cfrac{1}{a_3 + \cfrac{1}{ \ddots + \cfrac{1}{a_k} }}} }[/math]

be the kth convergent of the continued fraction

[math]\displaystyle{ a_1 + \cfrac{1}{a_2 + \cfrac{1}{a_3 + \cfrac{1}{ \ddots}}}. }[/math]

Vincent's theorem — Let [math]\displaystyle{ p_0(x) }[/math] be a square-free polynomial of degree n, and [math]\displaystyle{ a_1, a_2,\ldots, }[/math] be a sequence of real numbers. For i = 1, 2,..., consider the polynomial

[math]\displaystyle{ p_i(x)=x^np_{i-1}\left(a_i+1/x\right). }[/math]

Then, there is an integer k such that either [math]\displaystyle{ p_k(0)=0, }[/math] or the sequence of the coefficients of [math]\displaystyle{ p_k }[/math] has at most one sign variation.

In the first case, the convergent c_k is a positive root of [math]\displaystyle{ p_0. }[/math] Otherwise, this number of sign variations (either 0 or 1) is the number of real roots of [math]\displaystyle{ p_0 }[/math] in the interval defined by [math]\displaystyle{ c_{k-1} }[/math] and [math]\displaystyle{ c_k. }[/math]

For proving his theorem, Vincent proved a result that is useful on its own:

Vincent's auxiliary theorem — If p(x) is a square-free polynomial of degree n, and a, b, c, d are nonnegative real numbers such that [math]\displaystyle{ \left|\frac ac-\frac bd\right| }[/math] is small enough (but not 0), then there is at most one sign variation in the coefficients of the polynomial

[math]\displaystyle{ q(x)=(cx+d)^n p\left(\frac{ax+b}{cx+d}\right), }[/math]

and this number of sign variations is the number of real roots of p(x) in the open interval defined by [math]\displaystyle{ \frac ac }[/math] and [math]\displaystyle{ \frac bd. }[/math]

For working with real numbers, one may always choose c = d = 1, but, as effective computations are done with rational numbers, it is generally convenient to suppose that a, b, c, d are integers.

The "small enough" condition has been quantified independently by Nikola Obreshkov, and Alexander Ostrowski:

Obreschkoff–Ostrowski theorem: in blue and yellow, the regions of the complex plane where there should be no root for having 0 or 1 sign variation; on the left the regions excluded for the roots of p, on the right, the regions excluded for the roots of the transformed polynomial q; in blue, the regions that are excluded for having one sign variation, but allowed for zero sign variations.

Theorem (Obreschkoff–Ostrowski) — The conclusion of Vincent's auxiliary result holds if the polynomial p(x) has at most one root α + iβ such that

[math]\displaystyle{ \left(\alpha-\frac ac\right)\left(\alpha -\frac bd\right) +\beta^2\le \frac 1\sqrt 3\left|\beta\left(\frac ac-\frac bd\right)\right|. }[/math]

In particular the conclusion holds if

[math]\displaystyle{ \left|\frac ac -\frac bd\right| \lt \frac {\operatorname{sep}(p)}{2\sqrt 3}, }[/math]

where sep(p) is the minimal distance between two roots of p.

For polynomials with integer coefficients, the minimum distance sep(p) may be lower bounded in terms of the degree of the polynomial and the maximal absolute value of its coefficients; see Properties of polynomial roots § Root separation. This allows the analysis of worst-case complexity of algorithms based on Vincent's theorems. However, Obreschkoff–Ostrowski theorem shows that the number of iterations of these algorithms depend on the distances between roots in the neighborhood of the working interval; therefore, the number of iterations may vary dramatically for different roots of the same polynomial.

James V. Uspensky gave a bound on the length of the continued fraction (the integer k in Vincent's theorem), for getting zero or one sign variations:

Theorem (Uspensky) —  Let p(x) be a polynomial of degree n, and sep(p) be the minimal distance between two roots of p. Let

[math]\displaystyle{ \varepsilon = \left(n+\frac 1n\right)^\frac 1{n-1}-1. }[/math]

Then the integer k, whose existence is asserted in Vincent's theorem, is not greater than the smallest integer h such that

[math]\displaystyle{ F_{h-1}\operatorname{sep}(p)\gt 2\quad\text{and}\quad F_{h-1}F_k\operatorname{sep}(p)\gt \frac 1\varepsilon, }[/math]

where [math]\displaystyle{ F_h }[/math] is the hth Fibonacci number.

Continued fraction method

The use of continued fractions for real-root isolation has been introduced by Vincent, although he credited Joseph-Louis Lagrange for this idea, without providing a reference. For making an algorithm of Vincent's theorem, one must provide a criterion for choosing the [math]\displaystyle{ a_i }[/math] that occur in his theorem. Vincent himself provided some choice (see below). Some other choices are possible, and the efficiency of the algorithm may depend dramatically on these choices. Below is presented an algorithm, in which these choices result from an auxiliary function that will be discussed later.

For running this algorithm one must work with a list of intervals represented by a specific data structure. The algorithm works by choosing an interval, removing it from the list, adding zero, one or two smaller intervals to the list, and possibly outputs an isolation interval.

For isolating the real roots of a polynomial p(x) of degree n, each interval is represented by a pair [math]\displaystyle{ (A(x), M(x)), }[/math] where A(x) is a polynomial of degree n and [math]\displaystyle{ M(x)=\frac{px+r}{qx+s} }[/math] is a Möbius transformation with integer coefficients. One has

[math]\displaystyle{ A(x)=p(M(x)), }[/math]

and the interval represented by this data structure is the interval that has [math]\displaystyle{ M(\infty)=\frac pq }[/math] and [math]\displaystyle{ M(0)=\frac rs }[/math] as end points. The Möbius transformation maps the roots of p in this interval to the roots of A in (0, +∞).

The algorithm works with a list of intervals that, at the beginning, contains the two intervals [math]\displaystyle{ (A(x)=p(x), M(x)=x) }[/math] and [math]\displaystyle{ (A(x)=p(-x), M(x)=-x), }[/math] corresponding to the partition of the reals into the positive and the negative ones (one may suppose that zero is not a root, as, if it were, it suffices to apply the algorithm to p(x)/x). Then for each interval (A(x), M(x)) in the list, the algorithm remove it from the list; if the number of sign variations of the coefficients of A is zero, there is no root in the interval, and one passes to the next interval. If the number of sign variations is one, the interval defined by [math]\displaystyle{ M(0) }[/math] and [math]\displaystyle{ M(\infty) }[/math] is an isolating interval. Otherwise, one chooses a positive real number b for dividing the interval (0, +∞) into (0, b) and (b, +∞), and, for each subinterval, one composes M with a Möbius transformation that maps the interval onto (0, +∞), for getting two new intervals to be added to the list. In pseudocode, this gives the following, where var(A) denotes the number of sign variations of the coefficients of the polynomial A.

function continued fraction is
    input: P(x), a square-free polynomial,
    output: a list of pairs of rational numbers defining isolating intervals
    /* Initialization */
        L := [(P(x), x), (P(–x), –x)]                /* two starting intervals */
        Isol := [ ]
    /* Computation */
    while L ≠ [ ] do
        Choose (A(x), M(x)) in L, and remove it from L
        v := var(A)
        if v = 0 then exit                /* no root in the interval */
        if v = 1 then                     /* isolating interval found */
            add (M(0), M(∞)) to Isol
            exit
        b := some positive integer
        B(x) = A(x + b)
        w := v – var(B)
        if B(0) = 0 then                         /* rational root found */
            add (M(b), M(b)) to Isol
            B(x) := B(x)/x
        add (B(x),  M(b + x)) to L           /* roots in (M(b), M(+∞)) */
        if w = 0 then exit                  /* Budan's theorem */ 
        if w = 1 then                       /* Budan's theorem again */ 
            add (M(0), M(b)) to Isol
        if w > 1 then
            add ( A(b/(1 + x)),  M(b/(1 + x)) )to L      /* roots in (M(0), M(b)) */

The different variants of the algorithm depend essentially on the choice of b. In Vincent's papers, and in Uspensky's book, one has always b = 1, with the difference that Uspensky did not use Budan's theorem for avoiding further bisections of the interval associated to (0, b)

The drawback of always choosing b = 1 is that one has to do many successive changes of variable of the form x → 1 + x. These may be replaced by a single change of variable x → n + x, but, nevertheless, one has to do the intermediate changes of variables for applying Budan's theorem.

A way for improving the efficiency of the algorithm is to take for b a lower bound of the positive real roots, computed from the coefficients of the polynomial (see Properties of polynomial roots for such bounds).

Bisection method

The bisection method consists roughly of starting from an interval containing all real roots of a polynomial, and divides it recursively into two parts until getting eventually intervals that contain either zero or one root. The starting interval may be of the form (-B, B), where B is an upper bound on the absolute values of the roots, such as those that are given in Properties of polynomial roots § Bounds on (complex) polynomial roots. For technical reasons (simpler changes of variable, simpler complexity analysis, possibility of taking advantage of the binary analysis of computers), the algorithms are generally presented as starting with the interval [0, 1]. There is no loss of generality, as the changes of variables x = By and x = –By move respectively the positive and the negative roots in the interval [0, 1]. (The single changes variable x = (2By – B) may also be used.)

The method requires an algorithm for testing whether an interval has zero, one, or possibly several roots, and for warranting termination, this testing algorithm must exclude the possibility of getting infinitely many times the output "possibility of several roots". Sturm's theorem and Vincent's auxiliary theorem provide such convenient tests. As the use Descartes' rule of signs and Vincent's auxiliary theorem is much more computationally efficient than the use of Sturm's theorem, only the former is described in this section.

The bisection method based on Descartes' rules of signs and Vincent's auxiliary theorem has been introduced in 1976 by Akritas and Collins under the name of Modified Uspensky algorithm, and has been referred to as the Uspensky algorithm, the Vincent–Akritas–Collins algorithm, or Descartes method, although Descartes, Vincent and Uspensky never described it.

The method works as follows. For searching the roots in some interval, one changes first the variable for mapping the interval onto [0, 1] giving a new polynomial q(x). For searching the roots of q in [0, 1], one maps the interval [0, 1] onto [0, +∞]) by the change of variable [math]\displaystyle{ x\to\frac{1}{x+1}, }[/math] giving a polynomial r(x). Descartes' rule of signs applied to the polynomial r gives indications on the number of real roots of q in the interval [0, 1], and thus on the number of roots of the initial polynomial in the interval that has been mapped on [0, 1]. If there is no sign variation in the sequence of the coefficients of r, then there is no real root in the considered intervals. If there is one sign variation, then one has an isolation interval. Otherwise, one splits the interval [0, 1] into [0, 1/2] and [1/2, 1], one maps them onto [0, 1] by the changes of variable x = y/2 and x = (y + 1)/2. Vincent's auxiliary theorem insures the termination of this procedure.

Except for the initialization, all these changes of variable consists of the composition of at most two very simple changes of variable which are the scalings by two x → x/2 , the translation x → x + 1, and the inversion x → 1/x , the latter consisting simply of reverting the order of the coefficients of the polynomial. As most of the computing time is devoted to changes of variable, the method consisting of mapping every interval to [0, 1] is fundamental for insuring a good efficiency.

Pseudocode

The following notation is used in the pseudocode that follows.

• p(x) is the polynomial for which the real roots in the interval [0, 1] have to be isolated
• var(q(x)) denotes the number of sign variations in the sequence of the coefficients of the polynomial q
• The elements of working list have the form (c, k, q(x)), where
  • c and k are two nonnegative integers such that c < 2^k, which represent the interval [math]\displaystyle{ \left[\frac c{2^k}, \frac{c+1}{2^k}\right], }[/math]
  • [math]\displaystyle{ q(x)=2^{kn} p\left(\frac{x+c}{2^k}\right), }[/math] where n is the degree of p (the polynomial q may be computed directly from p, c and k, but it is less costly to compute it incrementally, as it will be done in the algorithm; if p has integer coefficients, the same is true for q)

function bisection is
    input: p(x), a square-free polynomial, such that p(0) p(1) ≠ 0, 
                      for which the roots in the interval [0, 1] are searched
    output: a list of triples (c, k, h), 
                      representing isolating intervals of the form [math]\displaystyle{ \left[\frac c{2^k}, \frac{c+h}{2^k}\right] }[/math]
    /* Initialization */
    L := [(0, 0, p(x))] /* a single element in the working list L */
    Isol := [ ]
    n := degree(p)
    /* Computation */
    while L ≠ [ ] do
        Choose (c, k, q(x)) in L, and remove it from L
        if q(0) = 0 then
            q(x) := q(x)/x
            n := n – 1                /* A rational root found */
            add (c, k, 0) to Isol
        v := [math]\displaystyle{ \operatorname{var}((x+1)^n q(1/(x+1))) }[/math]
        if v = 1 then                /* An isolating interval found */
            add (c, k, 1) to Isol
        if v > 1 then                /* Bisecting */
            add (2c, k + 1, [math]\displaystyle{ 2^n q(x/2) }[/math]  to L
            add (2c + 1, k + 1, [math]\displaystyle{ 2^n q((x + 1)/2) }[/math]  to L
    end

This procedure is essentially the one that has been described by Collins and Akritas. The running time depends mainly on the number of intervals that have to be considered, and on the changes of variables. There are ways for improving the efficiency, which have been an active subject of research since the publication of the algorithm, and mainly since the beginning of the 21st century.

Recent improvements

Various ways for improving Akritas–Collins bisection algorithm have been proposed. They include a method for avoiding storing a long list of polynomials without losing the simplicity of the changes of variables, the use of approximate arithmetic (floating point and interval arithmetic) when it allows getting the right value for the number of sign variations, the use of Newton's method when possible, the use of fast polynomial arithmetic, shortcuts for long chains of bisections in case of clusters of close roots, bisections in unequal parts for limiting instability problems in polynomial evaluation.

All these improvement lead to an algorithm for isolating all real roots of a polynomial with integer coefficients, which has the complexity (using soft O notation, Õ, for omitting logarithmic factors)

[math]\displaystyle{ \tilde O (n^2(k+t)), }[/math]

where n is the degree of the polynomial, k is the number of nonzero terms, t is the maximum of digits of the coefficients.

The implementation of this algorithm appears to be more efficient than any other implemented method for computing the real roots of a polynomial, even in the case of polynomials having very close roots (the case which was previously the most difficult for the bisection method).

Sources

• Alesina, Alberto; Massimo Galuzzi (1998). "A new proof of Vincent's theorem". L'Enseignement Mathématique 44 (3–4): 219–256. http://retro.seals.ch/cntmng?type=pdf&rid=ensmat-001:1998:44::149&subp=hires. Retrieved 2018-12-16. 
• Akritas, Alkiviadis G. (1986). "There's no "Uspensky's Method"". Proceedings of the fifth ACM Symposium on Symbolic and Algebraic Computation (SYMSAC '86). Waterloo, Ontario, Canada. pp. 88–90. http://dl.acm.org/citation.cfm?id=32457. 
• Akritas, Alkiviadis G.; Strzeboński, A. W.; Vigklas, P. S. (2008). "Improving the performance of the continued fractions method using new bounds of positive roots". Nonlinear Analysis: Modelling and Control 13 (3): 265–279. doi:10.15388/NA.2008.13.3.14557. http://www.lana.lt/journal/30/Akritas.pdf. 
• Akritas, Alkiviadis G.; Strzeboński, Adam W. (2005). "A Comparative Study of Two Real Root Isolation Methods". Nonlinear Analysis: Modelling and Control 10 (4): 297–304. doi:10.15388/NA.2005.10.4.15110. http://www.lana.lt/journal/19/Akritas.pdf. 
• Collins, George E.; Akritas, Alkiviadis G. (1976). "Polynomial Real Root Isolation Using Descartes' Rule of Signs". SYMSAC '76, Proceedings of the third ACM symposium on Symbolic and algebraic computation. Yorktown Heights, NY, USA: ACM. pp. 272–275. doi:10.1145/800205.806346. 
• Kobel, Alexander; Rouillier, Fabrice; Sagraloff, Michael (2016). "Computing real roots of real polynomials ... and now for real!". Waterloo, Canada. doi:10.1145/2930889.2930937. 
• Obreschkoff, Nikola (1963) (in de). Verteilung und Berechnung der Nullstellen reeller Polynome. Berlin: VEB Deutscher Verlag der Wissenschaften. p. 81. 
• Ostrowski, A. M. (1950). "Note on Vincent's theorem". Annals of Mathematics. Second Series 52 (3): 702–707. doi:10.2307/1969443. 
• Rouillier, F.; Zimmerman, P. (2004). "Efficient isolation of polynomial's real roots". Journal of Computational and Applied Mathematics 162 (1): 33–50. doi:10.1016/j.cam.2003.08.015. Bibcode: 2004JCoAM.162...33R. 
• Sagraloff, M.; Mehlhorn, K. (2016). "Computing real roots of real polynomials". Journal of Symbolic Computation 73: 46–86. doi:10.1016/j.jsc.2015.03.004. 
• Tsigaridas, Elias P.; Emiris, Ioannis Z. (2006). "Algorithms - ESA 2006, 14th Annual European Symposium, Zurich, Switzerland, September 11-13, 2006, Proceedings". in Azar, Yossi; Erlebach, Thomas. 4168. Springer. pp. 817–828. doi:10.1007/11841036_72. 
• Uspensky, James Victor (1948). Theory of Equations. New York: McGraw–Hill Book Company. https://www.google.com/search?q=uspensky+theory+of+equations&btnG=Search+Books. 
• Vincent, Alexandre Joseph Hidulphe (1834). "Mémoire sur la résolution des équations numériques" (in fr). Mémoires de la Société Royale des Sciences, de L' Agriculture et des Arts, de Lille: 1–34. http://gallica.bnf.fr/ark:/12148/bpt6k57787134/f4.image.r=Agence%20Rol.langEN. 
• Vincent, Alexandre Joseph Hidulphe (1836). "Note sur la résolution des équations numériques". Journal de Mathématiques Pures et Appliquées 1: 341–372. http://sites.mathdoc.fr/JMPA/PDF/JMPA_1836_1_1_A28_0.pdf. 
• Vincent, Alexandre Joseph Hidulphe (1838). "Addition à une précédente note relative à la résolution des équations numériques". Journal de Mathématiques Pures et Appliquées 3: 235–243. http://math-doc.ujf-grenoble.fr/JMPA/PDF/JMPA_1838_1_3_A19_0.pdf. Retrieved 2018-12-16.