Real Stability Testing
aa r X i v : . [ c s . D S ] O c t Real Stability Testing ∗ Prasad RaghavendraEECSUC Berkeley Nick RyderMathematicsUC Berkeley Nikhil SrivastavaMathematicsUC BerkeleyOctober 4, 2016
A univariate polynomial with real coe ffi cients is called real-rooted if all of its roots are real. Multivari-ate generalizations of this concept, known as hyperbolic and real stable polynomials, were defined inthe 50’s and in the 80’s in the context of Partial Di ff erential Equations and Control Theory, respec-tively , and have since made contact with several areas of mathematics. In particular, a polynomial p ∈ R [ x , . . . , x n ] is called real stable if it has no zeros with all coordinates in the open upper half ofthe complex plane. These polynomials have played a central role in several recent advances in theo-retical computer science and combinatorics — for instance, [AG14, MSS15a, MSS13, GSS11, BBL09].Each of these works relies in a critical way on (1) understanding which polynomials are real stable(2) understanding which linear operators preserve real-rootedness and real stability. Motivated by(1) and (2), this paper studies the following two fundamental algorithmic problems: Problem 1.
Given a bivariate polynomial p ∈ R n [ x , y ], is p real stable? Problem 2.
Given a linear operator T : R n [ x ] → R m [ x ], does T preserve real-rootedness?Problem 1 was solved in the univariate case by C. Sturm in 1835 [Stu35], who described a nowwell-known method that can be turned into a strongly polynomial quadratic time algorithm giventhe coe ffi cients of p [BPR05]. However, we are unaware of any algorithm (polynomial time or not)for the bivariate case, or for Problem 2.The main result of this paper is a strongly polynomial time algorithm that solves Problem 1. Theorem 1.1 (Main) . Given the coe ffi cients of a bivariate polynomial p ∈ R n [ x , y ] , there is a deterministicalgorithm which decides whether or not p is real stable in at most O ( n ) arithmetic operations, assumingexact arithmetic. Part of the motivation for solving Problem 1 is the following theorem of Borcea and Branden,which shows that Problem 2 can be reduced to Problem 1.
Theorem 1.2 (Borcea-Branden [BB09]) . For every linear transformation T : R n [ x ] → R m [ x ] , there is abivariate polynomial p ∈ R max( n , m ) [ x , y ] such that T preserves real-rootedness if and only if p is real stable.Moreover, the coe ffi cients of p can be computed from the matrix of T in linear time. ∗ This research was supported by NSF Grants CCF-1553751, NSF CCF-1343104 and NSF CCF-1407779. See [Pem12] for a more detailed history. We use R n [ x , . . . , x k ] to denote the vector space of real polynomials in x , . . . , x k of degree at most n in each variable. ff erentiation. Recent applications such as [MSS15b] rely on finding more elaboratedi ff erential operators with this property.We now describe the main ideas in our algorithm. It turns out that testing bivariate real stabilityis equivalent to testing whether a certain two parameter family of polynomials is real rooted. Itis not clear how to check this continuum of real-rootedness statements in polynomial, or evenin exponential time. To circumvent this, we use a deep convexity result from the theory ofhyperbolic polynomials to reduce the two parameter family to a one parameter family of degree n polynomials, whose coe ffi cients are themselves polynomials of degree n in the parameter. We thenuse a characterization of real-rootedness as postive semidefiniteness of certain moment matricesto further reduce this to checking that a finite number of univariate polynomials are nonnegative an interval. Finally, we solve each instance of the nonnegativity problem using Sturm sequencesand a bit of algebra.The set of polynomials nonnegative on an interval forms a closed convex cone, so the last step ofour algorithm may be viewed as a strongly polynomial time membership oracle for this cone. Wewould not be surprised if such a result is already known (at least as folklore) but we were unable tofind a concrete reference in the literature, so this component of our method may be of independentinterest.We see this result as being both mathematically fundamental, as well as useful for researcherswho work with stable polyomials, particularly since many of their known applications so far(e.g.[MSS15a]) put special emphasis on properties of bivariate restrictions. More speculatively, itis possible that being able to test membership in the set of real stable polynomials is a step towardsbeing able to optimize over them. The paper [Hen10] studied the problem of testing whether a bivariate polynomial is real zero (aspecial case of real stability). It reduced that problem to testing PSDness of a one-parameter familyof matrices which it then suggested could be solved using semidefinite programming, but withoutquite proving a theorem to that e ff ect. This work is partly inspired by ideas in [Hen10].The paper [KPV15] gives semidefinite programming based algorithms that can test whether certainrestricted classes of multia ffi ne polynomials are real stable (in more than 2 variables).The problem of certifying that a univariate polynomial is nonnegative is typically stated (forinstance, in lecture notes) as being the solution to a semidefinite program. If one were able to workout the appropriate error to which the SDP has to be solved, this could give a weakly polynomialtime algorithm for nonnegativity, which we suspect must be known as folklore. The paper [PP08]analyzes a semidefinite programming based algorithm in the special case when the polynomial isnondegenerate in an appropriate sense. 2 .2 Acknowledgments We thank Eric Hallman, Jonathan Leake, and Bernd Sturmfels for valuable discussions. We alsothank Didier Henrion for valuable correspondence regarding other algorthimic approaches tothese problems.
We recall below the definition of a real stable polynomial in an arbitrary number of variables.
Definition 2.1.
A polynomial p ∈ R [ x , . . . , x n ] is called real stable if it is identically zero or if p ( z , . . . , z n ) , z i ) > i = , . . . , n . Equivalently, p is real stable if and only ifthe univariate restrictions t p ( te + x , te + x , . . . , te n + x n )are real rooted whenever e , . . . , e n > x , . . . , x n ∈ R .The equivalence between the two formulations above is an easy exercise. Note that a univariatepolynomial is real stable if and only if it is real rooted. Note that we consider the zero polynomialto be real-rooted.We will frequently use the elementary fact that a limit of real-rooted polynomials is real-rooted,which follows from Hurwitz’s theorem (see, e.g. [Wag11, Sec. 2]), or from the argument principle.Real Stable polynomials are closely related to the following more general class of polynomials. Definition 2.2.
A homogeneous polynomial p ∈ R [ x , . . . , x n ] is called hyperbolic with respect to apoint e = ( e , . . . , e n ) ∈ R n if p ( e ) > t p ( te + x )are real rooted for all x ∈ R n . The connected component of { x ∈ R n : p ( x ) , } containing e is calledthe hyperbolicity cone of p with respect to e , and will be denoted K ( p , e ).Perhaps the most familiar example of a hyperbolic polynomial is the determinant of a symmetricmatrix: X det( X )for real symmetric X , which is hyperbolic with respect to the identity matrix since the characteristicpolynomial of a symmetric matrix is always real rooted. The corresponding hyperbolicity cone isthe cone of positive semidefinite matrices.The most important theorem regarding hyperbolic polynomials says that hyperbolicity cones are always convex, and that hyperbolicity at one point in the cone implies hyperbolicity at every otherpoint. Thus, hyperbolic polynomials and hyperbolicity cones may be viewed as generalizingdeterminants and PSD cones. Theorem 2.3 (Garding [Gar59]) . If p ∈ R [ x , . . . , x n ] is hyperbolic with respect to e ∈ R n then:1. K ( p , e ) is an open convex cone.2. p is hyperbolic with respect to every point y ∈ K ( p , e ) . Some works (e.g. [BB09]) consider only nonzero polynomials to be stable, while others [Wag11] include the zeropolynomial. We find the latter convention more convenient.
Theorem 2.4 (Borcea-Branden [BB09]) . A nonzero bivariate polynomial p ( x , y ) of total degree at most mis real stable if and only if its homogenizationp H ( x , y , z ) : = z m p ( x / z , y / z ) is hyperbolic with respect to every point in R > × { } = { ( e , e ,
0) : e , e > } . Thus, real stable polynomials enjoy the strong structural properties guaranteed by Theorem 2.3 aswell, and we exploit these in our algorithm.
In this section we use the properties of hyperbolic polynomials to reduce real stability of a bivariatepolynomial to testing real rootedness of a one parameter family of polynomials.
Theorem 3.1 (Reduction to One-Parameter Family) . A nonzero bivarite polynomial p ∈ R n [ x , y ] ofisreal stable if and only if following two conditions hold:1. The one-parameter family of univariate polynomials q γ ∈ R [ t ] given by,q γ ( t ) = p ( γ + t , t ) ∈ R [ t ] are real rooted for all γ ∈ R .2. The univariate polynomial t p H ( t , − t , is strictly positive on the interval (0 , ,Proof. ( real-stability of p = ⇒ (1) & (2))By Theorem 2.4, p H is hyperbolic with respect to the positive orthant R > × { } . Since (1 , , ∈ R > × { } , this implies that for all ( x , y , z ) ∈ R , q ( t ) = p H ( x + t , y + t , z )is real-rooted. Setting x = γ , y = z = q γ ( t ) = p H ( γ + t , t , = p ( γ + t , t ) isreal-rooted for all γ ∈ R which is condition (1). Finally, since { ( t , − t , | t ∈ (0 , } ⊂ R > × { } , and p H is hyperbolic with respect to R > × { } , it follows that p H ( t , − t , > t ∈ (0 , = ⇒ real-stability of p )First, we claim that the polynomial p H is hyperbolic with respect to (1 , , p H (1 / , / , > p H (1 , , >
0. It remains to show that q x , y , z ( t ) = p H ( x + t , y + t , z ) is real-rooted for all ( x , y , z ) ∈ R . First, consider the case of ( x , y , z ) ∈ R with z , ∀ ( x , y , z ) ∈ R with z , , p H ( x + t , y + t , z ) is real-rooted4 ⇒ ∀ ( x , y , z ) ∈ R with z , , p H ( xz + tz , yz + tz ,
1) is real-rooted ⇐⇒ ∀ ( x , y , z ) ∈ R with z , , p H ( xz + t , yz + t ,
1) is real-rooted (replacing t / z with t ) ⇐⇒ ∀ ( x , y ) ∈ R , p H ( x + t , y + t ,
1) is real-rooted ⇐⇒ ∀ ( x , y ) ∈ R , p H ( x + t , t ,
1) is real-rooted (replacing t with t − y ) ⇐⇒ ∀ γ ∈ R , p ( γ + t , t ) is real-rootedBy Hurwitz’s theorem, the limit of any sequence of real-rooted polynomials is real-rooted. There-fore, if q x , y , z ( t ) is real-rooted for all ( x , y , z ) ∈ R with z , q x , y , z ( t ) is real-rooted for all( x , y , z ) ∈ R .Given that p H is hyperbolic with respect to e = ( , , K ( p H , e ) is a convexcone containing (1 , , { x | p ( x ) , } containing (1 , ,
0) contains the open line segment from (1 , ,
0) to (0 , , R × { } ⊆ K ( p H , e ). By Theorem 2.4, this implies that p is real-stable. (cid:3) Thus, our algorithmic goal is reduced to testing whether a one-parameter family is real-rooted,and whether a given univariate polynomial is positive on an interval. We solve these problems inthe sequel.
In this section we present two algorithms for testing real-rootedness of a one-parameter family ofpolynomials. Both algorithms reduce this problem to verifying nonnegativity of a finite number ofpolynomials on the real line. The first algorithm produces n polynomials of degree roughly O ( n ),and has the advantage of being very simple, relying only on elementary techniques and standardalgorithms such as fast matrix multiplication and the discrete Fourier transform. The secondalgorithm produces n polynomials of degree roughly O ( n ) and runs significantly faster, but usessomewhat more specialized (but nonetheless classical) machinery from the theory of resultants. O ( n + ω ) Algorithm
The first algorithm is based on the observation that real-rootedness of a single polynomial isequivalent to testing positive semidefiniteness of its moment matrix, which in turn is equivalentto testing nonnegativity of the elementary symmetric polynomials of that matrix. In the moregeneral case of a one-parameter family, the latter polynomials turn out to be polynomials ofbounded degree in the parameter, and it therefore su ffi ces to verify that these are nonnegativeeverywhere.We begin by recalling the Newton Identities, which express the moments of a polynomial in termsof its coe ffi cients. Lemma 4.1 (Newton Identities) . Ifp ( x ) = n X k = ( − k x n − k c k = c n Y i = ( x − x i ) ∈ R [ x ]5 ith c , is a univariate polynomial with roots x , . . . , x n , then the momentsm k : = n X i = x ki satisfy the recurrence: m k = ( − k − c k c + k − X i = ( − k − + i c k − i c m i k n , m k = k − X k − n ( − k − + i c k − i c m i k > n , m = n . The following consequences of Lemma 4.1 will be relevant to analyzing our algorithm.
Corollary 4.2.
1. The moments m , . . . , m n − of a degree n polynomial can be computed from itscoe ffi cients in O ( n ) arithmetic operations.2. Suppose p ( x ) = P nk = ( − k x n − k c k ( γ ) is a polynomial whose coe ffi cients are polynomialsc ( γ ) , . . . , c n ( γ ) ∈ R d [ γ ] in a parameter γ . Then the moments of p are given bym k ( γ ) = r k ( γ ) / c ( γ ) k , for some polynomials r k ∈ R dk [ γ ] .Proof. The first claim follows because each application of the recurrence requires at most n arith-metic operations. For the second claim, observe that each ratio c k − i ( γ ) / c ( γ ) is a rational functionwith a numerator of degree at most d and denominator c ( γ ). Thus, each application of the recur-rence increases the degree of the numerator by at most d and introduces an additional c in thedenominator. (cid:3) As a subroutine, we will also need the following standard result in linear algebra.
Theorem 4.3 (Keller-Gehrig [KG85]) . Given an n × n complex matrix A, there is an algorithm whichcomputes the characteristic polynomial of A in time O ( n ω log n ) . We now specify the algorithm and prove its correctness.
Theorem 4.4.
A polynomial p γ ( x ) = P nk = ( − k x n − k c k ( γ ) is real-rooted for all γ ∈ R if and only if thepolynomials q , . . . , q n output by SimpleRR are nonnegative on R . Moreover, SimpleRR runs in time ˜ O ( dn + ω + d n ) .Proof. We first show correctness. Let m k ( p ) denote the k th moment of the roots of a polynomial. BySylvester’s theorem [BPR05, Theorem 4.58], a real polynomial p γ ( x ) = n X k = ( − k x n − k c k ( γ )is real-rooted if and only if the corresponding moment matrix M ( γ ) k , l : = m k + l − ( p γ )6s positive semidefinite. Since ν is even and c has real coe ffi cients, we have for every γ ∈ R that isnot a root of c : M ( γ ) (cid:23) ⇐⇒ c ( γ ) ν M ( γ ) = H ( γ ) (cid:23) . Since c has only finitely many roots and a limit of PSD matrices is PSD, we conclude that M ( γ ) (cid:23) ∀ γ ∈ R ⇐⇒ H ( γ ) (cid:23) ∀ γ ∈ R . Note that by Corollary 4.2 the entries of H ( γ ) are polynomials of degree at most d ( ν + n −
2) in γ .We now recall a well-known (e.g., [HJ12]) characterization of positive semidefiniteness as asemialgebraic condition: an n × n real symmetric matrix A is PSD if and only if e k ( A ) > k = , . . . , n , where e k ( A ) = X | S | = k det( A S , S )is the sum of all k × k principal minors of A . Thus, p γ is real-rooted for all γ ∈ R if and only if thepolynomials q k ( γ ) : = e k ( H ( γ ))for k = , . . . , n are nonnegative on R .Since each q k is a sum of determinants of order at most n in H ( γ ) it has degree at most n in the entriesof H ( γ ), and we conclude that q , . . . , q n ∈ R N [ γ ]. Thus, the q k can be recovered by interpolatingthem at the N th roots of unity. Since the k th elementary symmetric function of a matrix is thecoe ffi cient of z n − k in its characteristic polynomial, this is precisely what is achieved in Step 2.For the complexity analysis, it is clear that Step 1 takes O ( dn ) time. Constructing each Hankelmatrix H ( s i ) takes time O ( dn + n ) by Corollary 4.2, and computing its elementary symmetricfunctions via the characteristic polynomial takes time O ( n ω log n ), according to Theorem 4.3. Thus,the total time for each iteration is O ( n ω log n + dn ), so the time for all iterations is O ( dn + ω log n + d n ) . The final step requires O ( N log N ) time for each e k using fast polynomial interpolation via thediscrete Fourier transform, for a total of O ( dn log n ). Thus, the total running time is ˜ O ( dn + ω + d n ) , suppressing logarithmic factors. (cid:3) Here is a short proof: A is PSD i ff det( zI − A ) has only nonnegative roots. Since A is symmetric we know the rootsare real. We now observe that a real-rooted polynomial has nonnegative roots if and only if its coe ffi cients alternate insign. lgorithm SimpleRR Input: ( n +
1) univariate polynomials c , . . . , c n ∈ R d [ γ ] with c . Output: n univariate polynomials q , . . . , q n ∈ R n d [ γ ]1. Let ν be the first even integer greater than or equal to n and let N = nd (2 n − + ν ) = O ( dn ).Let s , . . . , s N ∈ C be the N th roots of unity.2. For each i = , . . . N : • Compute the n × n Hankel matrix H ( s i ) with entries H ( s i ) k , l : = c ( s i ) ν m k + l − ( p s i ) , by applying the Newton identities (Lemma 4.1). • Compute the characteristic polynomialdet( zI − H ( s i )) = n X k = ( − k z n − k e k ( H ( s i ))using the Keller-Gehrig algorithm (Theorem 4.3).3. For each k = , . . . , n : Use the points e k ( H ( s )) , . . . , e k ( H ( s N )) to interpolate the coe ffi cients ofthe polynomial q k ( γ ) : = e k ( H ( γ )) . Output q , . . . , q n . O ( n ) Algorithm Using Subresultants
The algorithm of the previous section is based on the generic fact that a matrix is PSD if and only ifits elementary symmetric polynomials are nonnegative. In this section we exploit the fact that ourmatrices have a special structure – namely, they are moment matrices – to find a di ff erent finite setof polynomials whose nonnegativity su ffi ces to certify their PSDness. These polynomials are called subdiscriminants , and turn out to be related to another class of polynomials called subresultants , forwhich there are known fast symbolic algorithms.Let M p denote the n × n moment matrix corresponding to a polynomial p of degree n . Recall that M p = VV T where V is the Vandermonde matrix formed by the roots of p . Let ( M p ) i denote theleading principal i × i minor of M p . We define subdiscriminants of a polynomial, and then showtheir relation to the leading principal minors of the moment matrix. For the remainder of thissection it will be more convenient to use the notation p ( x ) = n X k = a k x k for the coe ffi cients of a polynomial, with roots x , . . . , x n and a n , Definition 4.5.
The k th subdiscriminant of a polynomial p is defined as sDisc k ( p ) = a k − n X S ⊂{ ,..., n } , | S | = k Y { i , j }⊂ S ( x i − x j ) emma 4.6. The leading principal minors of the moment matrix are multiples of the subdiscriminants, ( M p ) i = a − kn sDisc k ( p ) = X S ⊂{ ,..., n } , | S | = k Y { i , j }⊂ S ( x i − x j ) Proof.
Let V i = . . . x . . . x n ... ... ... x i − . . . x i − n . Then ( M p ) i = det( V i V Ti ). By Cauchy-Binet, this determinant is the sum over the determinants ofall submatrices of size i × i . These submatrices are exactly the Vandermonde matrices formed bysubsets of the roots of size i . Then the identity follows from the formula for the determinant of aVandermonde matrix. (cid:3) Equipped with this we can provide an alternative characterization of real rootedness. Define thesign of a number, denoted sgn to be + − Lemma 4.7. p is real-rooted if and only if the sequence sgn ( sDisc ( p )) , . . . , sgn ( sDisc n ( p )) is first ’s andthen ’s.Proof. Note that since a n , sgn ( Disc k ) = sgn ( a − k ) n Disc k ) = sgn (( M p ) i ). It is clear fromthe definition of the subdiscriminants that if p is real-rooted with k distinct roots then sDisc i ispositive if i k and sDisc i = i > k .Conversely, given a polynomial p with k distinct roots, then if i > k we have all the minors of size i in V Ti contain two identical rows, and hence V Ti does not have full rank, so V i V Ti is singular. Let x , x , . . . , x j be the real distinct roots of p and y , ¯ y , . . . , y l , ¯ y l be the distinct complex conjugatepairs of p where j + l = k . Suppose the multiplicities of x i are n i and y i are m i . Then the top left k × k submatrix of M p is = X i n i x i ... x k − i h x i · · · x k − i i + X i m i y i ... y k − i h y i · · · y k − i i + y i ... ¯ y ik − h y i · · · ¯ y ik − i = X i n i x i ... x k − i h x i · · · x k − i i + X i m i y i ) Im( y i ) ... ... Re( y k − i ) Im( y k − i ) " − y i ) Im( y i ) ... ... Re( y k − i ) Im( y k − i ) T This shows that this submatrix is positive definite if and only if the distinct roots are all real. Notethat by Sylvester’s criterion this submatrix is positive definite if and only if all the leading principalminors of size k are positive. (cid:3) We now obtain a formula for the subdiscriminants of a polynomial in terms of its coe ffi cients. Theconnection is provided by another family of polynomials called the subresultants .9 efinition 4.8. Let p = P nk = a k x k where a n ,
0. The kth subresultant of p , denoted sRes k ( p , p ′ ) isthe determinant of the submatrix obtained from the first 2 n − − k columns of the following(2 n − − k ) × (2 n − − k ) matrix: a n · · · · · · · · · · · · a . . . . . . ... . . . a n · · · · · · · · · · · · a ... na n · · · · · · · · · a ... . . . . . . . . . . . . . . . . . . ... na n · · · · · · · · · a · · · We will use two properties of subresultants. The first is a good bound on their degree as aconsequence of the determinantal formula above. The second is quick algorithm to compute them.We refer the reader to [BPR05] for a more detailed discussion of subresultants.In this paper we will only be interested in subresultants of a polynomial with its derivative We areinterested in this because of its relation to our leading principal minors:
Lemma 4.9 ([BPR05] Proposition 4.27) . Let p ( x ) = P nk = a k x k where a n , sRes k ( p , p ′ ) = a n sDisc n − k ( p ) Corollary 4.10.
Since the first column of the determinant used to define the subresultant is divisible by a n ,we get sDisc k ( p ) is a polynomial in our coe ffi cients a n , . . . , a of degree at most n. The benefit of studying the principal minors instead of the coe ffi cients of the characteristic polyno-mial for our moment matrix is that we can use an algorithm from subresultant theory to quicklycalculate all the minors at once. Theorem 4.11 ([BPR05] Algorithm 8.21) . There exists an algorithm which, given a polynomial p of degreen returns a list of all of its subresultants sRes k ( p , p ′ ) for k = , . . . , n in O ( n ) time.Remark . Many computer algebra systems (e.g., Mathematica, Macaulay2) have built-in e ffi cientalgorithms to compute subresultants.We now combine the above facts to obtain a crisp condition for real-rootedness of a one-parameterfamily. Recall that by Theorem 3.1, we are interested in testing when a family of polynomials p γ ( x )are real-rooted for all γ ∈ R , where p γ ( x ) = n X k = a k ( γ ) x k with c k ∈ R n [ γ ]. Let c m ( γ ) be the highest coe ffi cient that is not identically zero. We are onlyinterested in the case when m > Proposition 4.13.
If p γ ( x ) = P nk = x k c k ( γ ) with c k ∈ R d [ γ ] , then sDisc k ( p γ ) is a polynomial in γ of degreeat most dn.Proof. From our previous lemma, we know that sDisc k is a polynomial in the coe ffi cients of p ofdegree at most 2 n . Since each of these coe ffi cients c k ( γ ) is a polynomial in γ of degree at most d ,our result follows. (cid:3)
10e now extend our characterization of real-rootedness in terms of the signs of the principal minorsof a fixed polynomial to a characterization for coe ffi cients which are polynomials in γ . Theorem 4.14. p γ ( x ) is real-rooted for all γ ∈ R if and only if there exists a k such that sDisc i ( p γ ) is anonnegative polynomial which is not identically zero for all i k and sDisc i ( p γ ) is identically zero fori > k.Proof. First suppose that p γ ( x ) is real rooted for all γ ∈ R . Observe that c m ( γ ) vanishes for at mostfinitely many points Z . Moreover, the degree m discriminant of p γ is a polynomial in γ , and iszero for at most finitely many points — call them Z . Thus, for γ < Z ∪ Z , we know that p γ hasexactly m distinct real roots, so by Lemma 4.7 sDisc i ( p γ ) is strictly positive for i m and zero for i > m on this set. By continuity this implies that sDisc i ( p γ ) is nonnegative and not identically zeroon R for i m , and sDisc i ( p γ ) is identically zero for i > m , as desired.To prove the converse, note that for i k , sDisc i ( p γ ( t )) is not identically zero, and hence thereare finitely many γ away from which sDisc i ( p γ ) is positive for all i k , and then all zero. ByLemma 4.7 we get that p γ ( x ) is real rooted for all these γ . Since real-rootedness is preserved bytaking limits (by Hurwitz’s theorem), we conclude that p γ ( x ) is real rooted for all γ ∈ R . (cid:3) Combining these observations, and using the O ( n ) time algorithm to compute the subdiscrimi-nants, we arrive at the following O ( n ) time algorithm for computing all the subdiscriminants. Algorithm FastRR
Input: ( n +
1) univariate polynomials c , . . . , c n ∈ R d [ γ ] with c . Output: n univariate polynomials q , . . . , q n ∈ R dn [ γ ]1. Find distinct points γ , . . . , γ dn ∈ such that c m ( γ i ) , γ i use the subresultant algorithm (Theorem 4.11) to compute all of the sRes k ( p γ i ),with k = , . . . , n .3. Use the above values to compute 2 dn di ff erent values q k ( γ ) , . . . , q k ( γ dn ) for each of thepolynomials q k ( γ ) : = sDisc k ( p γ ) = c m ( γ ) − sRes m − k ( p γ )) , k = , . . . , n .4. Use fast interpolation to compute the coe ffi cients of q , . . . , q n .Output q , . . . , q n . Theorem 4.15. FastRR runs in O ( n ) time.Proof. Since c n ( γ ) is of degree at most d we can test 2 dn + d points to find 2 dn points on which c n ( γ ) doesnt vanish. Each evaluation takes O ( d ) times, so total this takes O ( d n ) time. To compute sRes k ( p γ i ) for each 0 k n − i dn takes O ( dn ) time by Theorem 4.11. Then to scaleall the subresultants, since we have O ( dn ) data points and have already computed c n ( γ i ) takes O ( dn ) time. Finally, since the degrees of the q k are at most 2 dn , the total time to interpolate all ofthem is O ( dn log n ). (cid:3) Univariate Nonnegativity Testing
In this section, we describe an algorithm to test non-negativity of a univariate polynomial over thereal line.Let p ∈ R [ x ] denote a univariate polynomial of degree d . The goal of the algorithm is to test if p ( x ) > x ∈ R . A canonical approach for the problem would be to use a Sum-of-Squaressemidefinite program to express p as a sum of squares of low-degree polynomials. Unfortunately,the resulting algorithm is not a symbolic algorithm, i.e., its runtime is not strongly polynomial inthe degree d , since semidefinite programming is not known to be strongly polynomial.We will now describe a strongly polynomial time algorithm to test non-negativity of the polynomial p . Our starting point is an algorithm to count the number of real roots of a polynomial using Sturmsequences. We refer the reader to Basu et al. [BPR05] for a detailed presentation of Sturm sequencesand algorithms to compute them. For our purposes, we will need the following lemma. Lemma 5.1.
Given a univariate polynomial p ∈ R [ x ] , the algorithm based on computing Sturm sequencesuses O (deg( p ) ) arithmetic operations to determine the number of real roots of p. The polynomial p is positive, i.e., p ( x ) > x ∈ R , if and only if it has no real roots. Therefore,Lemma 5.1 yields an algorithm to test positivity using in O ( d ) arithmetic operations. To test non-negativity, the only additional complication stems from the roots of the polynomial p . We beginwith a simple observation. Fact 5.2.
If p ∈ R [ x ] is monic then p ( x ) > for all x ∈ R if and only if p has no real roots of odd multiplicity. Definition 5.3.
A square-free decomposition of a polynomial p ∈ R [ x ] of degree d , is a set ofpolynomials { a , . . . , a d } ∈ R [ x ] such that p ( x ) = d Y i = a i ( x ) i , and each a i has no roots with multiplicity greater than one. Alternately, for each i ∈ [ d ], a i ( x ) i consists of all roots of p with multiplicity exactly i .Square-free decompositions can be computed e ffi ciently using gcd computations. Yun [Yun76]carries out a detailed analysis of square-free decomposition algorithms. In particular, he showsthat an algorithm due to Musser can be used to compute square-free decompositions at the cost ofconstantly many gcd computations.Now, we are ready to describe an algorithm to test non-negativity.12 lgorithm Nonnegative Input
A monic polynomial p ∈ R [ x ], deg( p ) = dGoal Test if p ( x ) > x ∈ R .1. Using Musser’s algorithm, compute the square-free decomposition of p given by, p = Y i ∈ [ d ] a ii where a i ∈ R [ x ] has no roots with multiplicity greater than 1.2. For each i ∈ [ ⌈ d ⌉ ] • Using Sturm sequences, test if a i − has real roots. If a i − has real roots p is NOTnon-negative. Runtime.
Let T cd ( d ) denote the time-complexity of computing the gcd of two univariate polyno-mials of degree d . The runtime of Musser’s square-free decomposition algorithm is within constantfactors of T cd ( d ). Let S real ( ℓ ) denote the time-complexity of determining if a degree ℓ polynomialhas no real roots. Observe that X i deg( a i ) deg( p ) = d Since S real ( ℓ ) is super-linear in ℓ , we have P i ∈ [ d ] S real ( a i ) S real ( d ) . The run-time of the algorithmis given by O ( T cd ( d ) + S real ( d )). Using Sturm sequences, S real ( d ) = O ( d ) elementary operations onreal numbers (see [BPR05]). Using Euclid’s algorithm, T cd ( d ) = O ( d ) elemenetary operations onreal numbers. This yields an algorithm for non-negativity that incurs at most O ( d ) elementaryoperations. Finally, we combine the ingredients from sections 3, 4, and 5 to obtain the proof of our maintheorem.
Proof of Theorem 1.1.
Given the coe ffi cients of p , we can compute the coe ffi cients of the one-parameter family in (1) of Theorem 3.1 in time at most O ( n ). By Theorem 4.15, FastRR producesthe polynomials q , . . . , q n in time O ( n ). We check that some final segment of these polynomials areidentically zero by evaluating each one at O ( n ) points. These polynomials have degree O ( n ), so Nonnegative requires time O ( n ) to check nonnegativity of each remaining one, for a total runningtime of O ( n ).For part (2) of Theorem 3.1, we simply use a Sturm sequence to ensure that there are no roots in(0 , (cid:3) The algorithm in this paper o ff ers a starting point in the area of polynomial time algorithms forreal stability. In addition to the obvious possibility of improving the running time to say O ( n ) orbelow, several natural open questions remain: • Can the algorithm be generalized to 3 or more variables? The bottleneck to doing this isthat we do not know how to check real rootedness of 2-parameter families, or equivalently,nonnegativity of bivariate polynomials. 13
Is there an algorithm for testing whether a given polynomial is hyperbolic with respect to some direction, without giving the direction as part of the input? • Is there an algorithm for testing stability of bivariate polynomials with complex coe ffi cients?Perhaps leaving the realm of strongly polynomial time algorithms, the major open question inthis area is the following: a famous theorem of Helton and Vinnikov [HV07] asserts that everybivariate real stable polynomial can be written as p ( x , y ) = det( xA + yB + C )for some positive semidefinite matrices A , B and real symmetric C . Unfortunately, their proof doesnot give an e ffi cient algorithm for finding these matrices. Can the ideas in this paper, perhaps viausing SDPs to find sum-of-squares representations of certain nonnegative polynomials derivedfrom p , be used to obtain such an algorithm? References [AG14] Nima Anari and Shayan Oveis Gharan,
The kadison-singer problem for strongly rayleighmeasures and applications to asymmetric tsp , arXiv preprint arXiv:1412.1143 (2014).[BB09] Julius Borcea and Petter Brändén,
The lee-yang and pólya-schur programs. i. linear operatorspreserving stability , Inventiones mathematicae (2009), no. 3, 541–569.[BBL09] Julius Borcea, Petter Brändén, and Thomas Liggett,
Negative dependence and the geometryof polynomials , Journal of the American Mathematical Society (2009), no. 2, 521–567.[BPR05] Saugata Basu, Richard Pollack, and Marie-Francoise Roy, Algorithms in real algebraicgeometry , vol. 20033, Springer, 2005.[Gar59] Lars Garding,
An inequality for hyperbolic polynomials , Journal of Mathematics and Me-chanics (1959), no. 6, 957–965.[GSS11] Shayan Oveis Gharan, Amin Saberi, and Mohit Singh, A randomized rounding approachto the traveling salesman problem , Foundations of Computer Science (FOCS), 2011 IEEE52nd Annual Symposium on, IEEE, 2011, pp. 550–559.[Hen10] Didier Henrion,
Detecting rigid convexity of bivariate polynomials , Linear Algebra and itsApplications (2010), no. 5, 1218–1233.[HJ12] Roger A Horn and Charles R Johnson,
Matrix analysis , Cambridge university press,2012.[HV07] J William Helton and Victor Vinnikov,
Linear matrix inequality representation of sets , Com-munications on pure and applied mathematics (2007), no. 5, 654–674.[KG85] Walter Keller-Gehrig, Fast algorithms for the characteristics polynomial , Theoretical com-puter science (1985), 309–317.[KPV15] Mario Kummer, Daniel Plaumann, and Cynthia Vinzant, Hyperbolic polynomials, inter-lacers, and sums of squares , Mathematical Programming (2015), no. 1, 223–245.14MSS13] Adam Marcus, Daniel A Spielman, and Nikhil Srivastava,
Interlacing families i: Bipartiteramanujan graphs of all degrees , Foundations of Computer Science (FOCS), 2013 IEEE 54thAnnual Symposium on, IEEE, 2013, pp. 529–537.[MSS15a] Adam W Marcus, Daniel A Spielman, and Nikhil Srivastava,
Interlacing families ii: Mixedcharacteristic polynomials and the kadison–singer problem , Annals of Mathematics (2015),no. 1, 327–350.[MSS15b] ,
Interlacing families iv: Bipartite ramanujan graphs of all sizes , Foundations ofComputer Science (FOCS), 2015 IEEE 56th Annual Symposium on, IEEE, 2015, pp. 1358–1377.[Nui69] Wim Nuij,
A note on hyperbolic polynomials , Mathematica Scandinavica (1969), no. 1,69–72.[Pem12] Robin Pemantle, Hyperbolicity and stable polynomials in combinatorics and probability , arXivpreprint arXiv:1210.3231 (2012).[PP08] Helfried Peyrl and Pablo A Parrilo,
Computing sum of squares decompositions with rationalcoe ffi cients , Theoretical Computer Science (2008), no. 2, 269–281.[Stu35] Charles Sturm, Mémoire sur la résolution des équations numériques , 1835.[Wag11] David Wagner,