Certified counting of roots of random univariate polynomials

Certified counting of roots of random univariate polynomials

Joseph Cleveland, Jeffrey Dzugan, Jonathan D. Hauenstein,
Ian Haywood, Dhagash Mehta, Anthony Morse,
Leonardo Robol and Taylor Schlenk
Joseph Cleveland
Sam Houston State University
Jeffrey Dzugan
Department of Mathematics
Samford University
Jonathan D. Hauenstein
Department of Applied and Computational Mathematics and Statistics
University of Notre Dame
Notre Dame, IN 46556
USA
www.nd.edu/j̃hauenst Ian Haywood Department of Mathematics
North Carolina State University
Raleigh
North Carolina  27695
USA
Dhagash Mehta
Department of Applied and Computational Mathematics and Statistics
University of Notre Dame
Notre Dame, IN 46556
USA
www.nd.edu/d̃mehta Anthony Morse Department of Mathematics
State University of New York, Brockport
USA
Leonardo Robol Scuola Normale Superiore,
Piazza dei Cavalieri, 7 - 56127 Pisa
Italy.
Taylor Schlenk Department of Mathematics
University of North Dakota
Abstract.

A challenging problem in computational mathematics is to compute roots of a high-degree univariate random polynomial. We combine an efficient multiprecision implementation for solving high-degree random polynomials with two certification methods, namely Smale’s -theory and one based on Gerschgorin’s theorem, for showing that a given numerical approximation is in the quadratic convergence region of Newton’s method of some exact solution. With this combination, we can certifiably count the number of real roots of random polynomials. We quantify the difference between the two certification procedures and list the salient features of both of them. After benchmarking on random polynomials where the coefficients are drawn from the Gaussian distribution, we obtain novel experimental results for the Cauchy distribution case.

Research was supported by NSF grants DMS-1063010 and DMS-1262428, NSA project number H98230-12-1-0299, DARPA YFA, and Sloan Research Fellowship. The authors would like to thank Kambiz Farahmand, Antonio Lerario, Erik Lundberg and Yan Fyodorov for their feedback.

1. Introduction

A random univariate polynomial is a univariate polynomial whose coefficients are picked from a random distribution. Studying the roots of random polynomials is a classic problem in mathematics due to its numerous applications in many areas of mathematics and engineering, e.g., random matrix theory, hydrology, meteorology, aerodynamics, and structural engineering [13, 15]. Other situations also lead to solving random polynomials including solving polynomials where coefficients are subject to random noise; the solutions of differential equations with random coefficients; problems in spectral theory of random matrices; polynomial regression; and the development computer algebra methods [6, 44].

Recently, more applications have emerged from theoretical physics, especially in statistical physics, e.g., [16, 47], quantum mechanics, e.g., [10], and string theory, e.g., [3, 24]. The roots of certain random polynomials are also related to minimizing the energy of charged particles on a sphere [2], i.e., Smale’s 7th problem [49].

Formally, a random polynomial is given as

 (1) f(x)=N∑i=0aixi,

where is the degree of the polynomial and the coefficients are random numbers picked from a distribution. The Fundamental Theorem of Algebra yields that the number of complex roots of is , counting multiplicity. Since writing the exact roots of as a general algebraic expression of the coefficients for is not possible in general, due to the Abel-Ruffini theorem, we use numerical root-finding methods. There are numerous numerical methods available to find roots of univariate polynomials with floating point coefficients, e.g., companion matrices, the Aberth method, and the Jenkins-Traub method. All of these methods work well in certain cases (see Ref. [34] for an extensive review). There are also several symbolic methods to solve polynomials with rational coefficients, but the random polynomials under consideration will generally not have rational coefficients.

Approximating the roots of univariate polynomials with moderate degrees starts becoming an increasingly difficult task. However, in theory, numerical approximation methods typically have complexity around , but have the disadvantage of yielding only approximated results. This introduces the need to define what we actually consider a “good approximation” of a polynomial root. “Good” numerical methods often are backward-stable, that is, the results obtained at the end are exact results of a slightly perturbed input. Unfortunately, for ill-conditioned problems, even if they are posed in an exact setting, this might not be enough to provide bounds on the quality of the result.

Smale defined “good enough” as being in the quadratic convergence zone of Newton’s method for an exact root of the polynomial so that the exact root can be efficiently approximated to arbitrary accuracy. A point in such a quadratic convergence zone is called an approximate root and Smale showed how to certify using only data computed at the test point that it was an indeed an approximate root via -theory. An implementation of -theory is provided in alphaCertified [21, 22]. Other results of this kind also exist in the literature, e.g., [46, 52].

Other certification techniques are available for the univariate case, which, as we will see, can lead to better inclusion results. In this paper, we compare -theory with a global approach used in the MPSolve package [7, 8]. The MPSolve package implements an efficient multiprecision univariate polynomial solver which our tests show is efficient and effective at solving polynomials even as the degree increases.

The main inclusion result, originally proved in the paper by Tilli [52], roughly states that if the distance of an approximation from the root is less than , where is a constant depending on the separation of the roots, then Newton’s method converges quadratically from the start. The constant can be easily computed or bounded without computing the roots by means of Gerschgorin’s theorem applied to a particular companion matrix. One may view this as a global method since depends on all of the roots whereas -theory is a local method. Since we are focusing on random polynomials where, with probability , has distinct roots, the global root separation bound will be positive. Thus, this global method is applicable for the problems under consideration.

The output of the software MPSolve is well-suited for this random context: it can deliver good approximations using multiprecision arithmetic (given exact input, it obtains the roots to an arbitrary number of digits), can guarantee the inclusion radii, and can certify if a root is real. In this way, we were able to solve large samples of random polynomials and to certifiably count the number of real roots.

For concrete problems, we work with two different sets of random polynomials: coefficients are independently picked from the Gaussian distribution with mean zero and variance one, and from the standard Cauchy distribution. Since the Gaussian case is well studied with many analytical results available, we use this as our benchmark system. The Cauchy case is relatively under-studied and, except for the mean number of real roots, analytical results are usually not available. Hence, our numerical results for the Cauchy case can provide much needed input for potentially leading to new analytical results.

The rest of the paper is organized as follows. Section 2 describes numerical methods under consideration for approximating roots of univariate polynomials as well as discusses two certification methods. We also pictorially show the certifiable regions of a simple and exactly solvable polynomial equation. Section 3 describes our setup for solving and certifying the roots of random polynomials. Section 4 highlights some differences between the two certification methods, especially as the degree increases, and provides numerical results for various quantities associated to random polynomials.

2. Numerical Methods

One common approach for solving univariate polynomial equations is to use the companion matrix. This method is widely used in mathematical software such as Matlab. After describing the companion matrix method, we summarize Smale’s -theory followed by the certification approach based on Gerschgorin’s theorem used in MPSolve. The global nature of certification based on Gerschgorin’s theorem generally provides larger certification regions due to exploiting more information about the roots than the local -theory approach. This section ends with certification methods for detecting real roots.

2.1. The Companion Matrix Method

Suppose that for presented in (1), the leading coefficient is nonzero, i.e., is a univariate polynomial of degree . One can find the roots of by computing the eigenvalues of the companion matrix:

 (2) Cf:=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0000−a0/aN1000−a1/aN0100−a2/aN0010...........00…1−aN−1/aN⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,

since is the characteristic polynomial of . The eigenvalues of can be computed, for example, using the QR algorithm as implemented in LAPACK [23]. The roots() command in Matlab uses this companion matrix method.

2.2. Smale’s α-Theorem

Due to conditioning of , the approximations of the roots of computed by, for example, LAPACK or Matlab, using double precision may not be in a quadratic convergence zone of Newton’s method. One way to certify they are indeed in a quadratic convergence zone is to use -theory. Although the approach works for more general situations, we focus on univariate polynomials of degree . Let be the first derivative of and be the derivative. A Newton iteration of starting at is defined by

 Nf(x∗):={x∗−f(x∗)/f′(x∗)if f′(x∗)≠0,x∗otherwise.

Clearly, if , then if and only if .

Newton’s method relies upon repeatedly applying Newton iterations. For simplicity, the Newton iteration will be denoted

 Nkf(x):=Nf∘⋯∘Nfk times(x).

After defining the notion of an approximate root, we state a theorem which can be used to certify that a point is indeed an approximate root without a priori knowledge about the exact root. A point is an approximate root of with associated root if, for each ,

 (3)

That is, and is in the quadratic convergence basin of for .

Theorem 1 ([9], pg. 160).

If and is a univariate polynomial of degree  such that and

 α(f,x∗)<13−3√174,

then is an approximate root of with associated root where

 α(f,x):=β(f,x)⋅γ(f,x),\omit\span\@@LTX@noalign\vskip3.0ptplus1.0ptminus1.0pt\omitβ(f,x):=|f(x)/f′(x)|, and\omit\span\@@LTX@noalign\vskip3.0ptplus1.0ptminus1.0pt\omitγ(f,x):=max2≤k≤N∣∣∣f(k)(x)k!⋅f′(x)∣∣∣1/(k−1).

Moreover, .

Given two approximate roots and , they must have distinct associated roots, via the triangle inequality, if

 |x1−x2|>2(β(f,x1)+β(f,x2)).

This and several other -theoretic algorithms are implemented in the software package alphaCertified [21, 22] which, among other things, has been recently used to certify solutions to multivariate problems arising in theoretical physics and chemistry [42, 43].

In the present article, we take advantage of the fact that random polynomials of degree have, with probability , distinct complex roots. This means that if we find approximate roots with distinct associated roots, we have accounted for all of the roots.

2.3. Approximate Roots using MPSolve

Another way to certify a point is an approximate root of a univariate polynomial is the following global approach.

Theorem 2 ([52]).

Let be a monic univariate polynomial of degree with distinct roots , i.e., with for . Suppose that and such that

 |x∗−ζi|≤|x∗−ζj|3(N−1)~{}~{}for all~{}j≠i.

Then, is an approximate root of with associated root .

This theorem involves computing the distance from a point to the a priori unknown roots of . The following lemma provides an effective bound.

Lemma 3.

If is a monic univariate polynomial of degree and such that , then contains at least one root of .

Proof.

See [20]. ∎

The following combines Theorem 2 and Lemma 3 to allow one to algorithmically verify that a given point is an approximate root given information about all of the roots. This type of global approach is used in the MPSolve package [7, 8].

Corollary 4.

Let be a monic univariate polynomial of degree with distinct roots and suppose that . If, for all ,

 (4) |xi−xj|>N(3N−1)(β(f,xi)+β(f,xj))

then each is an approximate root of with distinct associated roots.

Proof.

By Lemma 3, the ball centered at of radius contains at least one root. For , it follows immediately from (4) that . Thus, each contains a distinct root of , say . For ,

 |xi−ζj|≥|xi−xj|−|xj−ζj|>N(3N−2)(β(f,xi)+β(f,xj))−Nβ(f,xi)≥3N(N−1)(β(f,xi)+β(f,xj)).

Thus, for all ,

 |xi−ζj|3(N−1)≥N(β(f,xi)+β(f,xj))≥|xi−ζi|

and the statement follows from Theorem 2. ∎

The balls provides some information about the roots. Gerschgorin’s theorem is another way to obtain strict and useful inclusion results to analyze clusters of roots, such as the following statement.

Theorem 5.

Let be a matrix with . Let and be the ball centered at with radius . Then,

1. each eigenvalue of is contained in , and

2. every connected component of this union made up of balls contains exactly eigenvalues of (counted with multiplicity).

Suppose that is a monic univariate polynomial and are distinct. Define . Consider the following matrix:

 A:=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣w1+x1w2⋯wNw1w2+x2⋯wN⋮⋮⋱⋮w1w2⋯wN+xN⎤⎥ ⎥ ⎥ ⎥ ⎥⎦
Lemma 6.

The matrix is a companion matrix for , i.e., .

Proof.

Let , and . Recall that . We have then that so we can write:

 det(xI−A) =det((xI−D)(I−(xI−D)−1ewT)) =N∏i=1(x−xi)⋅(1−N∑i=1wix−xi) =N∏i=1(x−xi)+N∑i=1f(xi)∏j≠i(x−xjxi−xj)=:q(x).

The expression above proves that we have for every . Recalling that is a degree polynomial we can conclude that the relation must hold for every . In particular, for every , we know and so we can conclude that . ∎

Since we have that is a companion matrix for , we can state the following.

Corollary 7.

Let be distinct and let be the values defined above. The union of the balls contains all the roots of . Moreover, every connected component of this union made up of balls contains exactly roots.

One may wonder which one between Gerschgorin radii and Newton radii is more strict. It turns out that if each is a good approximation of a simple root (or more formally, when ), we have that they are exactly the same and are both equal to up to the first order expansion. This can be verified by the following simple computation under the assumption that with sufficiently small. In the cases of Gerschgorin radii, we have

 N⋅∣∣ ∣∣f(xi)∏j≠i(xi−xj)∣∣ ∣∣=N⋅|xi−ζi|⋅∣∣ ∣∣∏j≠ixi−xj−ϵjxi−xj∣∣ ∣∣≈N⋅|xi−ζi|

while, for Newton radii, we obtain

This implies that not only they are essentially equivalent, but also that they are optimal up to a factor of .

2.4. Detecting Real Roots

For a polynomial with real coefficients, it was observed in [21] that -theory can be used to certifiably decide the reality or nonreality of the associated root given an approximate root. This local certification test is based on the simple fact that a root of is real if and only if there is a real approximate root.

Corollary 7 can be used to develop a global certification test as developed in the following which uses and to denote the real and imaginary part, respectively, of a number .

Theorem 8.

Let be a monic univariate polynomial of degree with distinct roots and . Suppose that such that for all ,

1. ,

2. , and

3. .

Then, has a real root such that .

Proof.

To prove by contradiction, suppose that is the unique root of in the ball such that . Hence, the conjugate of , namely , is also a root of . So, there must exists such that . Observe that Condition implies that every circle centered in , cannot intersect the circle . In particular, this implies that the complex conjugate of cannot be a root of , unless is a real root. ∎

In the random setting where, with probability , the univariate polynomial of degree will have distinct roots, we are trying to yield a set of points consisting of approximate roots with distinct associated roots. When the roots have a large pairwise separation, numerical methods such as those implemented in MPSolve yield an effective tool for approximating such points. As the minimum of the pairwise distance decreases, one needs to use higher precision in order to yield approximations to separate the roots and achieve quadratic convergence of Newton’s method. In the limit when the minimum pairwise distance is zero, standard Newton’s method can only achieve linear convergence. In this case, for example, MPSolve can compute an approximation correct to digits, where the integer is selected by the user.

2.5. Certifiable Regions for Comparing the Methods

As a comparison of the two certification methods, we consider a “Wilkinson-type” polynomial which has exact integer roots in order to compare numerical results. Figure 1 summarizes results from -theory and Figure 2 summarizes results from MPSolve. In both of these plots, the white regions correspond to certifiable regions where the black points are approximations computed by the roots() command in Matlab and MPSolve, respectively. When using MPSolve, the roots are inside the certifiable regions while the roots() command, probably due to the limitation of using double precision computations, obtained roots outside of some -theoretic certifiable quadratic convergence basins.

In the rest of the figure the areas were filled with different colors based on the root to which the Newton method converges if using that point as the starting one.

We can note how the regions obtained through MPSolve’s certification procedure have a similar radius, while this is not true for alphaCertified. This is a direct consequence of the uniformity of the constants involved in the first strategy.

3. Numerical Setup

The following describes our numerical setup for generating random polynomials. In particular, we generate a random polynomial of degree by drawing the coefficients from a chosen random distribution such as the uniform distribution or the Gaussian distribution with specified mean and variance using Matlab which generates random numbers with 32-bit precision. Then, we use two numerical root-finding methods and two certification methods.

First, we approximate all of the roots using either the companion matrix method as implemented in Matlab’s roots() command111We have tried other root-finding methods such as Lindsey-Fox algorithm, as well the built-in root-finding methods in Maple and Mathematica. However, we find the same conclusions as with Matlab’s roots() command. or using MPSolve.

To use -theory to certify the roots, we compute a set of rational approximations of the coefficients of the polynomial such that, for each we have . With numerically estimated roots and rational coefficients in hand, we use alphaCertified using -bit floating point arithmetic to certify the approximated roots. If any of the roots is not certified, then we run several Newton iterations on the uncertified root until it finally reaches the quadratic convergence zone of some exact root. In our experiments, we obtained  certified distinct roots for the random sample and then counted the real roots.

As an alternative approach, we use the approach presented in Section 2.3 to obtain approximate solutions for the roots of a given polynomial. MPSolve provides a Matlab interface which allows the user to request an arbitrary number of digits and certify the result. This approach is valuable because it alleviates the numerical ill-conditioning issues that may be encountered when solving high degree polynomials with (possibly big) integer coefficients. For these cases, the Matlab’s roots() command, even if using a backward stable algorithm, can give very poor approximations of the roots. As a consequence, it is possible that alphaCertified may not be able to certify them due to large errors. This issue is solved by MPSolve using a global approach that is guaranteed to produce correct digits on the output (and certify they are approximate solutions). A suitable level of floating point multiprecision will be used by MPSolve internally to achieve this, but it may be not necessary to expose this to Matlab (unless more than digits are required).

This method is effective in our experiments. As an example, we solved random polynomials of degrees between and using MPSolve, which was able to certify all the roots and determine their reality. While in principle it may happen that a root cannot be classified via Theorem 8, this did not happen in our experiments.

4. Results

The following presents our results for the roots of random polynomials with coefficients taking indenpendent and identically distributed (i.i.d.) values from two distributions: the Gaussian distribution with mean and variance which is one of the most extensively studied cases and can serve as our benchmark as many analytical results are available here, and the Cauchy distribution for which only a handful of analytical results are available and hence our numerical results can provide valuable input to the theory.

4.1. Quantifying Certification Procedures

While trying to approximate the roots, we wished to ensure that our numerical approximations were in the quadratic convergence zone of Newton’s method. In Section 2, two different certification methods are discussed. The first, based on -theory, which requires that , is a local certification approach in that it requires only data from one point. The second, based on Gerschgorin’s theorem, which requires that the Newton correction is smaller than the relative distance between the roots divided by approximately , is a global certification approach in that it requires information about all roots simultaneously.

The second approach can be reformulated in the same framework as -theory if we replace with . One can now compare the mean value of at the roots of with with the smaller value leading to a larger certification region.

If we assume that the distance of the roots is proportional to , which is a realistic expectation for random polynomials, then . Experimentally, we have verified that in the -theory has an exponential growth, as can be clearly seen in Figure 3. This causes certification to be computationally prohibitive as increases. For example, using only double precision, Figure 3 suggests that degrees less than are feasible using double precision, but may be insufficient for degrees larger than due to large values of .

It is important to note that the comparison between these two strategies is not fair since is directly computable given a single approximation while needs to be bounded using some kind of inclusion theorem since it relies on the (unknown) distance between roots of . However, our numerical results should provide valuable feedback when comparing the two certification procedures. In particular, our experiments showed that the usage of the second strategy, using a global approach, is much more effective for the cases under consideration.

In Section 2.3, we showed that the Gerschgorin radii are suboptimal only up to a factor of . This additional factor of leads to a moderate overestimate of , namely .

4.2. Real Roots

For a first experiment, we analyzed the number of real roots of random polynomials. Many theoretical results are known, but further analysis could be carried out to obtain more information about the real roots. In particular, the following summarizes some experimental results that we have obtained by solving polynomials of degrees , , …, . We measured the mean and variance for our sample and checked our results against the theoretical predictions, when available.

4.2.1. Average Number of Real Roots

The average number of real roots is one of the most well-studied quantities in the context of random polynomials. We have estimates both for the Gaussian case and for the Cauchy case. We only summarize the main results here and leave more in-depth information to [15].

Theorem 9.

Let be a random polynomial of degree where are i.i.d. selected from a standard normal distribution. If is the number of real roots of , there exists a constant such that

 E(R(a0,…,an))=2πlogN+C0+2Nπ+O(1N2).

The constant is approximately equal to .

Theorem 10.

Let be a random polynomial of degree where are i.i.d. with a standard Cauchy distribution. If is the number of real roots of , then there exists constants , , and such that

 E(R(a1,…,aN))=Clog(N+1)+A0+A2(N+1)2+O(1N3).

The constants , , and are approximately equal to:

 C=0.7413,A0=0.559132,A2=0.230596.

The experimental results are summarized in Figure 4. This figure shows that the theoretical results agree with the experimental results, but the convergence can be quite slow. That is, even with sample polynomials, there are still some oscillations around the theoretical value that are visible from the plot.

4.2.2. Variance of the Number of Real Roots

From the experimental data, we can estimate the variance on the number of real roots. For the Gaussian distribution, we compare this with the following theoretical result obtained in [33].

Theorem 11.

Let be a random polynomial such that the random variables are i.i.d. and satisfy the following:

• , and

• the random variables do have a moment of order for some .

If , then

 limN→∞Var(#{real roots of f(x)=0})M(N)=1.

We can see from Figure 5 that the experimental estimate and theoretical value seems to be off by a constant. This is due to the fact that, unlike in the case of the expected value, we do not know all the terms of the expansion so we can only state that the limit of as .

Moreover, the plot seems to confirm that both the Gaussian and the Cauchy case have a logarithmic behavior as increases. Note that the hypotheses of Theorem 11 are not satisfied in case of the Cauchy random polynomial since the random variables do not have moments of any order.

4.2.3. Histograms of the Number of Real Roots

In Figure 5(b) and 5(a) we report the number of real roots for random polynomials of degree in form of an histogram. We analyze both the Gaussian and the Cauchy case. The plots were generated starting from a sample of random polynomials with the different distributions for the coefficients.

4.2.4. Average Number of Positive/Negative Real Roots

Motivated by the problem of finding the average number of positive and negative real eigenvalues of real matrices, which amounts to finding the average number of positive and negative real roots of the corresponding characteristic polynomials, we ask what is the average number of positive and negative real roots for our random polynomials. We expect that the mean of the number of positive roots is equal to the mean number of negative roots since both of our random distributions are symmetric with respect to the origin. That is, it is easy to check that the polynomial has roots the opposites of the ones of but the distribution of the coefficients is the same (since both in the Gaussian and in the Cauchy case multiplying by does not change the random variable).

Even though this is clear from an asymptotic point of view, it may be interesting to further examine the topic by looking at the variance of the number of positive and negative roots.

As a further experimental proof of the above statement, one can look at Figure 7 that shows the ratio between the number of positive and negative real roots for a sample of random polynomials. We have kept the number of samples low to show that there is still some floating around the value , but it is quite clear that the problem of counting the roots is symmetric in the change of sign.

4.2.5. Variance on the Number of Positive Roots

Figure 8 plots our experimentally computed variance on the number of positive real roots using random polynomials for each sample point. As far as we know, there are no theoretical estimates for the variance of the number of positive roots. Unfortunately, since the variance is increasing very slowly, is not easy to guess its asymptotic behavior.

4.3. Complex Roots

A natural extension of the questions that we have tried to address in the previous section is the analysis of the distribution of complex roots. We are considering real random distributions so we expect to have a certain number of real roots and some complex conjugate ones. The following present evidence about the distribution of the real and imaginary parts and on their modulus.

4.3.1. Distribution of the Complex Roots

Figures 912 plot the approximate density of the distribution of the real and imaginary parts of the roots of Cauchy and Gauss polynomials, respectively.

4.3.2. Magnitude of the Roots

We also analyze the distribution of the modulus of the complex roots. The results of our experiments are reported in Figure 12(a) and in Figure 12(b). Only the case with is reported since all our plots look very similar and the modulus of almost all the roots have modulus near .

4.4. Stationary Points of Univariate Random Potential Energy Landscape

In many areas of science and engineering, studying the surface defined by a multi-dimensional hypersurface described by a multivariate nonlinear function, called potential energy function, is a fundamental problem – called analyzing the potential energy landscape of the given system [51]. In particular, a great deal of research is invested in finding stationary points of potential energy landscapes of different multivariate systems. Finding stationary points of multivariate functions amounts to solving the system of equations obtained by all the first derivatives with respect to all the variables equated to zero. To obtain statistical information about generic potential energy landscapes, it is customary to explore the potential energy landscapes of the random potentials, i.e., multivariate functions whose coefficients are randomly picked from a random distribution [19].

Using our numerical setup, we mimic exploring the random potential energy landscapes of the univariate random potential , where

 (5) F(x)=D∑i=0aixi.

In other words, we solve . The number of maxima can be obtained by just dividing the number of real zeros of by , since the number of inflection points is of measure (see [15]).

4.4.1. Average number of real critical points

With our numerical setting we have performed a simulation on random polynomials to check the theoretical estimates in [15]. The results can be seen in Figure 13(a).

4.4.2. Variance of the number of maxima

With the same data obtained above we can also estimate the variance of the number of zeros, which is summarized in Figure 13(b). We note that the analysis for the minima is totally symmetric.

4.4.3. Average number of minima at which F(x)>0

We now measure the minima for which . In a string theory set up, such minima are called de Sitter minima. In this set up, evaluated at a minimum corresponds to the cosmological constant. Since the observations have revealed that the cosmological constant is always positive definite, the search is then restricted to minima with .

It is easy to see that in our case if and then by defining we have that and . In particular this implies that the mean number of minima points where are bound to be equal to the mean number of minima points where is negative due to symmetry reasons. So we can expect that the number of positive minima will be half of the negative ones.

The only question that we can ask is what can be said about the variance of this number.

The results of this analysis are reported in Figure 15. Similarly to the variance of positive roots, we see that the variance is slowly increasing with the degree and it is also slowly converging with the size of the sample.

4.4.4. The Histogram of F′′(x) at the Real Roots

For multivariate random potentials, an important physical quantity is the eigenvalues of the Hessian matrix evaulated at the real stationary points. The distribution of the eigenvalues of random multivariate potentials is well studied in theoretical and statistical physics. It is shown that in certain cases the distribution follows Wigner’s semicircle law. Here, for the univariate random potential case, we want to analyze the distribution of the evaluations of on the real roots. Due to the very high oscillations of the value of given the high degree of the polynomials, we were able to easily analyze only the case of degree . We decided to numerically estimate the distribution of instead of estimating to obtain a much more meaningful plot.

The expermiental results obtained for are reported in Figures 15(a) and 15(b).

5. Conclusion and Outlook

Finding roots of high-degree univariate polynomials with coefficients chosen from particular random distributions, though a very important task that is related to advances in many different areas of science, mathematics and engineering, is a highly nontrivial task. Though there are several numerical root-finding methods available, the numerical approximates found by these methods may turn out to be in the linear or even worse in the chaotic convergence region of the nearby exact root. In this work, we attempted to find certified roots of such high-degree univariate polynomials. In particular, we used two different certification methods to certify if the numerical approximates are within the quadratic convergence region, namely based on Smale’s -theory, a local method, and Gerschgorin’s theorem, a global method. We first find all the potential numerical approximates of a univariate polynomial using traditional methods such as the companion matrix method as implemented in Matlab’s roots() command, or the multiprecision method implemented in MPSolve. Then, we certify if the numerical approximates are certifiable in the above mentioned sense, using both the certification methods via alphaCertified and MPSolve. We compared these methods based on the local conditioning constant from -theory with arising from Gerschgorin’s theorem. In our experiment, the average of increases exponentially with the degree whereas increases polynomially in the degree.

One advantage of Smale’s -theory is that it is a local method. That is, it can be performed using data computed at one point. This allows for a quite independent and self-contained analysis of the method. On the other hand, relying on universal, worst-case methods often causes the bounds to grow very fast. Global methods, which are based on the minimum distance between the roots, can alleviate the fast growth at the expense of computing all of the roots.

We emphasize that, given a high-degree polynomial, multiprecision computations may be needed to perform the computations. This is also true in cases where only floating point output ( digits) is desired but, due to the ill-conditioning of the problem, an intermediate algorithm implemented in multiprecision might be able to obtain all the digits.

We also provided detailed experiments from the Cauchy case, which is less studied than the Gaussian case. Our experiments suggest that:

1. the available theoretical results are easily verifiable experimentally, as we have shown, for example, in Figure 4 for the expected value for the number of real roots;

2. there is some additional difficulties in the theoretical analysis of Cauchy random polynomials given by the fact that the Cauchy distribution does not have moments of any finite order, e.g., it is not possible to apply Maslova’s theorem, but experimental results show that the variance of the number of real roots has a behaviour compatible with the Maslova’s estimate.

In conclusion, these results suggest a generalization to the Cauchy case is possible.

References

• [1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972.
• [2] D. Armentano, B. Carlos, and M.  Shub. Minimizing the discrete logarithmic energy on the sphere: The role of random polynomials, Trans. of the American Mathematical Society, 363 (6) 2955 (2011).
• [3] A. Ashmore and Y. -H. He, Calabi-Yau Three-folds: Poincare Polynomials and Fractals, arXiv:1110.1612 [hep-th].
• [4] D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler, Numerically solving polynomial systems with Bertini, Software, environments, and tools 25, SIAM books, Philadelphia, PA, 2013.
• [5] D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler, Software for Numerical Algebraic Geometry, Chicago (2006), Software available at http://bertini.nd.edu.
• [6] A. T. Bharucha-Reid and M. Sambandham. Random polynomials. 1986.
• [7] D. A. Bini and G. Fiorentino. Design, analysis, and implementation of a multiprecision polynomial rootfinder. Numerical Algorithms 23.2-3 (2000): 127-173.
• [8] D. A. Bini and L. Robol. Solving secular and polynomial equations: A multiprecision algorithm. Journal of Computational and Applied Mathematics (2013).
• [9] L. Blum, F. Cucker, M. Shub, and S. Smale, Complexity and real computation, with a foreword by Richard M. Karp, Springer-Verlag, New York, 1998.
• [10] E. Bogomolny, O. Bohigas, and P. Leboeuf. Distribution of roots of random polynomials. Phys. Rev. Lett. 68 (18) 2726 (1992).
• [11] D. Bshouty and A. Lyzzaik, Problems and conjectures for planar harmonic mappings: in the Proceedings of the ICM2010 Satellite Conference: International Workshop on Harmonic and Quasiconformal Mappings (HQM2010), Special issue in: J. Analysis 18 (2010), 69-82.
• [12] D. Bshouty, W. Hengartner, and T. Suez, The exact bound on the number of zeros of harmonic polynomials, J. Anal. Math. 67 (1995), 207-218.
• [13] W. E. H.  Culling. Dimension and entropy in the soil covered landscape. Earth Surface Processes and Landforms, 13 (7), 619 (1988).
• [14] A. Edelman and E. Kostlan, How many zeros of a random polynomial are real?, Bull. Amer. Math. Soc. 32 (1995), 1-37.
• [15] K.  Farahmand. Topics in random polynomials (1998).
• [16] Y. V. Fyodorov, G. A. Hiary and J. P. Keating.Freezing Transition, Characteristic Polynomials of Random Matrices, and the Riemann Zeta Function. Phys. Rev. Lett. 108 (17) 170601 (2012).
• [17] L. Geyer, Sharp bounds for the valence of certain harmonic polynomials, Proc. AMS 136 (2008), 549-555.
• [18] A. Granville, I. Wigman, The distribution of the number of zeros of trigonometric polynomials, Amer. J. of Math. 133 (2011), 295-357.
• [19] B. Greene, D. Kagan, A. Masoumi, D. Mehta, E. J. Weinberg and X. Xiao, Tumbling through a landscape: Evidence of instabilities in high-dimensional moduli spaces, Phys. Rev. D 88 (2013), 026005.
• [20] P. Henrici, Applied and Computational Complex Analysis: Vol.: 1.: Power Series, Integration, Conformal Mapping, Location of Zeros. New York: Wiley, 1974.
• [21] J. D. Hauenstein and F.  Sottile, Algorithm 921: alphaCertified: Certifying solutions to polynomial systems, ACM TOMS 38 (2012), 28.
• [22] J.D. Hauenstein and F. Sottile, alphaCertified: Software for certifying numerical solutions to polynomial equations, Available at
www.math.tamu.edu/sottile/research/stories/alphaCertified.
• [23] E. Anderson et al. LAPACK Users’ Guide, SIAM (1999).
• [24] Y. -H. He, Polynomial Roots and Calabi-Yau Geometries, Adv. High Energy Phys. 2011, 719672 (2011)
• [25] Y.-H. He, D. Mehta, M. Niemerg, M. Rummel and A. Valeanu, Exploring the Potential Energy Landscape Over a Large Parameter-Space, J. High Energy Phys. 1307 (2013), 050.
• [26] C. Hughes, D. Mehta and J. -I. Skullerud, Enumerating Gribov copies on the lattice, Annals Phys.  331 (2013), 188-215.
• [27] M. Kastner and D. Mehta, Phase Transitions Detached from Stationary Points of the Energy Landscape, Phys. Rev. Lett.  107 (2011), 160602.
• [28] D. Khavinson, G. Swiatek, On a maximal number of zeros of certain harmonic polynomials Proc. AMS, 131 (2003), 409-414.
• [29] S-Y. Lee, A. Lerario, E. Lundberg, Remarks on Wilmshurst’s theorem, to appear in Indiana U. Math. J.
• [30] W. V. Li, A. Wei (2009), On the expected number of zeros of random harmonic polynomials, Proc. AMS 137 (2009), 195-204.
• [31] M. Maniatis and D. Mehta, Minimizing Higgs potentials via numerical polynomial homotopy continuation, Eur. Phys. J. Plus 127 (2012), 91.
• [32] D. Martinez-Pedrera, D. Mehta, M. Rummel and A. Westphal, Finding all flux vacua in an explicit example, J. High Energy Phys. 2013 (2013), 110.
• [33] N.B. Maslova. On the variance of the number of real roots of random polynomial, Theory Prob. Appl., 19:35-52, 1974.
• [34] J. M. McNamee and V. Pan. Numerical methods for roots of polynomials. Vol. 16. Newnes, (2013).
• [35] D. Mehta, Ph.D. Thesis, The Uni. of Adelaide, Australasian Digital Theses Program (2009).
• [36] D. Mehta, A. Sternbeck, L. von Smekal and A. G. Williams, Lattice Landau Gauge and Algebraic Geometry, Proceedings of Science QCD-TNT09 (2009), 025.
• [37] D. Mehta, Finding all the stationary points of a potential energy landscape via numerical polynomial homotopy continuation method, Phys. Rev. E 84 (2011), 025702.
• [38] D. Mehta, Numerical Polynomial homotopy continuation method and string vacua, Adv. High Energy Phys.  2011 (2011), 263937.
• [39] D. Mehta, Y.-H. He and J. D. Hauenstein, Numerical algebraic geometry: A new perspective on string and gauge theories, J. High Energy Phys. 1207 (2012), 018.
• [40] D. Mehta, J. D. Hauenstein and M. Kastner, Energy landscape analysis of the two-dimensional nearest-neighbor model, Phys. Rev. E 85 (2012), 061103.
• [41] D. Mehta, D. A. Stariolo and M. Kastner, Energy landscape of the finite-size spherical three-spin glass model, Phys. Rev. E 87 (2013), 052143.
• [42] D. Mehta, J. D. Hauenstein and D. J. Wales, Certifying the Potential Energy Landscape, J. Chem. Phys., 138 (2013), 171101.
• [43] D. Mehta, J. D. Hauenstein and D. J. Wales, Certification and the Potential Energy Landscape. J. Chem. Phys. 140, 224114 (2014).
• [44] V. Y. Pan. Solving a Polynomial Equation: Some History and Recent Progress. SIAM Review 39 (2) 187 (1997).
• [45] R. Peretz, J. Schmid, On the zero sets of certain complex polynomials, Proceedings of the Ashkelon Workshop on Complex Function Theory (1996), 203-208, Israel Math. Conf. Proc. 11, Bar-Ilan Univ. Ramat Gan, 1997.
• [46] J. Renegar. On the worst-case arithmetic complexity of approximating zeros of polynomials. Journal of Complexity 3.2 (1987): 90-113.
• [47] G. Schehr and S. N.  Majumdar. Statistics of the number of zero crossings: from random polynomials to the diffusion equation, Phys. Rev. Lett. 99 (6) 060603 (2007).
• [48] T. Sheil-Small, Complex Polynomials, Cambridge University Press, 2002.
• [49] S. Smale, Mathematical problems for the next century, The Math. Intelli. 20 (2) 7 (1998).
• [50] A. J. Sommese and C. W. Wampler, The numerical solution of systems of polynomials arising in Engineering and Science, World Scientific Publishing Company, 2005.
• [51] D.J. Wales, Energy Landscapes, Cambridge University Press, 2004.
• [52] P. Tilli, Convergence conditions of some methods for the simultaneous computation of polynomial zeros. Calcolo 35.1 (1998): 3-15.
• [53] A. S. Wilmshurst, The valence of harmonic polynomials, Proc. AMS 126 (1998), 2077-2081.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters