Thurston, Selberg, and Random Polynomials, Part I.

Apart from everything else, you could always count on Bill Thurston to ask interesting questions. This is the first of a small number of posts which were motivated in part by figure two from this paper, and this accompanying MO question. I liked this problem enough to give it as a thesis problem to my student Zili Huang, and much of what I discuss below arose from this project.

Say that an algebraic integer \(\alpha\) is Perron if \(|\alpha| > |\sigma \alpha|\) for every conjugate \(\sigma \alpha\) of \(\alpha\). One immediately observes that \(\alpha\) must be real. Say that a monic polynomial is Perron if it is irreducible and has a Perron integer as a root. Thurston’s question is (roughly) to describe the distribution of Perron algebraic integers, especially those chosen in some (small) fixed interval in \(\mathbf{R}\). This question has several interpretations, but one experiment Thurston does is to take 20,000 monic polynomials of degree 21 with integer coefficients in \([-5,5]\), and plots the quantities \(\sigma \alpha /\alpha \in B(1)\) for all the conjugates of the 5,932 resulting Perron polynomials such that the corresponding Perron integer was in the interval \([1,2]\). The result is this:

The first observation is that this graph has (apart from some noise coming from real roots) rotational symmetry. The next observation is that the roots tend to be concentrated in a ring of some radius, which (from experiment) becomes more concentrated the more one restricts the range in \(\mathbf{R}\) of the Perron integers one is considering. The first question is: can one explain this graph, and does it reflect reality (that is, the actual distribution of Perron integers)?

The answers to these questions turn out to be: yes, and no. The first problem is that it is hard (a priori) to “randomly” generate Perron algebraic integers of large degree in \([1,2]\). Knowing a bound on the roots places a bound on the coefficients, but a randomly chosen polynomial with coefficients satisfying the required bounds will almost always have a root larger than 2. Thus Thurston “cheats” with his algorithm, making the coefficients of his polynomials very small in order to increase the probability that the largest root will also be small. (Full disclosure, Thurston makes no claims that his algorithm reflects reality, and explicitly asks whether it does so or not.) The issue is then whether this will skew the distribution of the roots. It turns out that it does! To explain why this might not be surprising, let’s talk about the size of the spaces over which Thurston is sampling. Let \(\Omega^P_{21}\) be the set of monic polynomials of degree 21 with real coefficients and with a unique largest real root \(\lambda \le 2\). Thurston is sampling over a space with \(11^{20}\) lattice points and volume \(10^{20}\). On the other hand, it turns out that the volume of \(\Omega^{P}_{21}\) is equal to

\(\displaystyle{\frac{2^{399}}{3^{24} 5^{12} 7^{10} 11^{11} 13^{9} 17^{5} 19^{3}}} \sim 2.249 \times 10^{60}.\)

So Thurston was only really sampling a \(10^{-40}\)th of the entire space! Thurston’s picture can be explained as follows: polynomials with (suitably) small coefficients (contingent on the initial and final coefficients not being too small) tend to have all their roots clustering uniformly around the disc of radius one. This follows in the radial direction by a famous theorem of Erdös and Turán, and for the absolute values it follows (in a related way using Jensen’s formula) from a paper of Hughes and Nikeghbali here. So the apparent “radius” in Thurston’s picture is just representing \(1/R\), where \(R\) is the approximate size of the Perron integers being considered. It turns out that, in reality, most of the conjugates of Perron integers have size comparible to the Perron integer itself. That is, the correct version of Thurston’s picture should show the roots clustering (roughly) uniformly around the boundary.

OK, now a pause when I look at Thurston’s graph and see that the radius is not something like a half as I claimed above, but something much smaller. So I just repeated Thurston’s experiment, and out of 20,000 monic polynomials with coefficients randomly chosen in [-5,5], only 1011 were Perron polynomials with largest root less than 2, and the resulting picture came out like this:

Here one really sees the (misleading) accumulation around the radius \(1/2\). I’m guessing that Thurston actually kept all polynomials whose largest root was in \([1,5]\), which would account for the larger success rate for choosing Perron polynomials as well as the smaller radius. This is also consistent with how Thurston describes the corresponding graphs in the MO question rather than in Figure 2 of his preprint.

So how does one study Perron integers? Let us re-wind slightly and discuss a more elementary problem. How does one count algebraic integers? The most natural way to count algebraic integers is to order them by height. However, Thurston’s problem clearly suggests a different measure, namely, to count by the size of the largest conjugate. This has a profound effect on some of the statistical properties under consideration. Roughly, algebraic integers ordered by height are much more likely to have a small number of “outliers” with large absolute value, whereas when one orders by the size of the largest conjugate, most of the other conjugates accumulate around the circle with radius the size of the largest root as the degree goes to infinity.

The problem of understanding algebraic integers of bounded size (where by bounded we mean a bound on the largest conjugate) amounts to understanding the lattice points in a certain region of \(\mathbf{R}^N\). Now as long as one fixes the degree and increases the bound, such counting problems (including this one) typically reduce to a volume problem. (One also uses the fact that almost all polynomials are irreducible, and that the regions are “nice” in some explicit way, i.e. not Cantor sets.) Moreover, the corresponding regions are essentially (up to a simple stretching) independent of the bound. Hence the key region to understand is the region \(\Omega_{N} \subset \mathbf{R}^N\) of monic degree \(N\) polynomials all of whose roots have absolute value at most one, and the region \(\Omega^{P}_N \subset \Omega_N\) consisting of such polynomials whose largest root is real. Of course, one is not only interested in the volumes of these regions, but also the integrals of various quantities. As an example, one can consider the integral

\(C_N(T,\alpha) = \displaystyle{\int_{\Omega_N} P(T) |a_N|^{\alpha – 1} dV}\)

where \(P \in \Omega_N\) represents the monic polynomial at any point, and \(a_N\) is the constant term. Evaluating this integral at \(\alpha = 1\) and taking the leading term (in \(T\)) recovers the volume. On the other hand, there are some other relations. A fairly simple computation shows that

\(\mathrm{Vol}(\Omega^P_N) = \displaystyle{\frac{4}{N(N+1)} C_{N-1}(1,1)},\)

which is how one can compute the left hand side exactly for any \(N\). In order to evaluate these integrals, it makes more sense to integrate not over the “coefficient space” of polynomials, but rather the “configuration space” of roots. The coefficient space is naturally stratified by the number of real and complex roots. For that reason, it makes sense to decompose \(\Omega_N\) as

\(\coprod_{R + 2 S = N} \Omega_{R,S}\)

where \(\Omega_{R,S}\) corresponds to polynomials whose roots all have absolute value at most one and have signature \((R,S)\) (since we are interested only in integrals, we elide issues concerning whether one wants these spaces to be open or closed or somewhere in between). As a special case, let’s think about the integral \(C_{N,0}(T,\alpha)\) where one restricts the integrand to \(\Omega_{N,0}\). The configuration space is simply \([-1,1]^N\). On the other hand, the map from configuration space to coefficient space is just given in terms of the symmetric polynomials, and the corresponding Jacobian matrix is the Vandermonde determinant. Hence, taking into account the action of \(S_N\) on the fibres, one finds that

\( C_{N,0}(T,\alpha) = \displaystyle{\frac{1}{N!} \int_{[-1,1]^N} \left| \prod x_i \right|^{\alpha – 1}
\prod (T – x_i) \prod |x_i – x_j| dx_1 \ldots dx_N}.\)

This is now very reminiscent of the classical Selberg Integral. There is some beautiful mathematics related to the Selberg integral; let me direct you here for a nice survey. The integrals arising here are, however, not quite Selberg integrals except for some very degenerate cases.

Once you start writing these integrals down, and computing some of them (by hook or crook), there are a number of problems which naturally come to mind. For example, what is the probability that a random polynomial all of whose roots have absolute value at most one is Perron? Well, by explicitly computing the ratio of the volume of \(\Omega^P_N\) to \(\Omega_N\), you find that the answer is \(1/N\) if \(N\) is odd and \(1/(N-1)\) if \(N\) is even (this checks out for \(N = 1,2\)). On the other hand, you might ask, given a polynomial all of whose roots have absolute value at most one, what the expected number of real roots?, or what the probability is (at least in even degree) that the polynomial has no real roots at all? Having asked these questions, it is then sensible to ask the same questions for other ways of choosing random polynomials. The classical way to choose a real random polynomial is to write

\(f(x) = a_N x^N + \ldots + a_0\)

where the \(a_i\) are independent normal variables with mean zero (this is the Kac ensemble). To what extent do the statistics of random polynomials with this measure mirror the constrained problem consisting of polynomials all of whose roots have absolute value at most one? Obviously, it depends on the type of problem one considers. The most classical problem for real polynomials concerns counting the expected number of real roots. A famous theorem of Kac says that, under the ensemble above, the expected number of real roots is approximately \(2/\pi \cdot \log(N)\). I recommend reading this paper for an introduction to the subject; I learnt these things from chatting with Peter Sarnak at the IAS.) The methods of Kac also show that the real roots concentrate for large N around \(– 1\) and \(+ 1\). In fact, the complex roots also concentrate along the unit circle as well. How does this compare to our constrained model? First of all, the real roots in Kac model either lie in \([-1,1]\) or in \([-\infty,1] \cup [1,\infty]\). Certainly our polynomials have no roots in the larger region. If one restricts the Kac polynomials to \([-1,1]\), then the expected number of real roots decreases to \(1/\pi \cdot \log(N)\). This is in some sense easy to see from the previous formula, because the map on coefficients \(a_k \rightarrow a_{N-k}\) is measure preserving and inverts the roots. In fact, a stronger result follows from Kac. If one takes an inteveral \([a,b]\) strictly contained inside \([-1,1]\), then the expected number of real roots in the polynomial for sufficiently large \(N\) converges to

\( \displaystyle{\frac{1}{\pi} \int^{b}_{a} \frac{1}{1 -T^2}}.\)

This gives another strong indication of how the roots are concentrating at the points +1 and -1. OK, so now let us return to our constrained model consisting of monic polynomials all of whose roots have absolute value at most one. How many real roots does one expect such a polynomial to have? There’s a natural map

\(\Omega_{N-1} \times [-1,1] \rightarrow \Omega_{N}\)

which sends \(P(x)\) to \(P(x)(x-T)\). The Jacobian of this matrix turns out to be equal to \(|P(T)|\). On the other hand, the map is not one to one, rather, the image of \(\Omega_{R,S}\) has multiplicity \(R\). Hence, if \(Z(P)\) denotes the number of real roots of the polynomial \(P\), then

\(\displaystyle{\int_{\Omega_N} Z(P)} = \int_{0}^{1} \int_{\Omega_{N-1}} |P(T)| dV\)

The left hand side (after dividing by the volume) gives the expected number of real roots. So one is again reduced to a Selberg type integral. In this case, one apparently has (based on some Zagier-like integral guessing mojo, but unfortunately not yet Zagier-like integral proving mojo) for \(N = 2m\),

\( \displaystyle{\frac{1}{D_N}
\int_{\Omega_N} |P(T)| =
\frac{1}{2^{2m}
{2m \choose m}}
\left( \sum_{k=0}^{m} \frac{2m-2k+1}{2m+1}
{2m-2k \choose m-k}
{2k \choose k}
T^{2k} \right)
\left( \sum_{k=0}^{m}
{2m-2k \choose m-k}
{2k \choose k}
T^{2k} \right)},\)

and there is a similar formula for \(N = 2m+1\). After some analysis to estimate the resulting integral of the RHS from \(T = -1\) to \(1\), it turns out that, for large \(N\), the expected number of real roots is approximately

\( \displaystyle{\frac{1}{\pi} \log N},\)

whic is exactly in accordance with the Kac model! Indeed, if one restricts to real roots in an interval \([a,b]\) strictly in \([-1,1]\), then one also obtains the same integral formula as in the Kac ensemble. So, somewhat surprisingly to me, the number of real roots in \([-1,1]\) behaves in a very similar way whether one considers Kac polynomials or monic polynomials all of whose roots have absolute value at most one.

What then of the other problems? Given a polynomial in the Kac model of even degree \(2N\), what is the probability that is has no roots in the interval \([-1,1]\)? This problem was explicitly addressed by Dembo, Poonen, Shao, and Zeitouni here, where they show (under less restrictive hypotheses) that this occurs with probability \(O(N^{-b/2 + o(1)})\) for some universal constant \(b/2\) which they do not determine, although they estimate based on numerical evidence that \(b/2 = 0.38 \pm 0.015\). What happens in our constrained model? Once more it comes down to a Selberg-like integral, this time computing the ratio of volumes:

\(\displaystyle{\frac{\displaystyle{\int_{\Omega_{0,N}} dV}}{\displaystyle{\int_{\Omega_{2N}} dV}}}\)

It turns out that one can compute this explicitly as a product of factorials. Moreover, one can compute the exact asymptotic in this case as \(N \rightarrow \infty\), and the resulting probability is

\(\displaystyle{ \frac{2C}{\sqrt{2 \pi} (2N)^{3/8}}, \ \text{where} \ C = 2^{-1/24} e^{-3/2 \cdot \zeta'(-1)} = 1.24512 \ldots }\)

(It may be hard to read in the exponent, but that is the derivative of the Riemann zeta function \(\zeta'(-1)\) at \(-1\). That may seem strange, but in fact this is a fairly typical constant that comes up in asymptotics of the Barnes-G function, which is exactly the type of expression (a product of factorials) which turns up in the evaluation of the relevant integrals.) Now the result of DPSZ does not apply in our case (where the coefficients are a long way from being independent), but given the similarity in the distribution of real roots between our polynomials and the Kac model, we naturally make the following conjecture:

Conjecture: The constant \(b/2\) is equal to \(3/8\).

Optimistically, one might even try to prove this conjecture by showing that the statistics of our collection of polynomials mirror those of the Kac polynomials for sufficiently large \(N\).

Next time: we discuss a more concrete relationship between random polynomials and our models in terms of limits of gap probabilities. But let me also leave you with the following teaser question: What is the probability that the largest root of a polynomial of degree \(N\) is real?