A Faster Solution to Smale's 17th Problem I: Real Binomial Systems

TL;DR

This paper presents a faster deterministic algorithm for finding real roots of binomial systems, improving efficiency over previous methods for certain polynomial systems, and extends to more general probability measures.

Contribution

It introduces a new average-case polynomial-time algorithm for real roots of binomial systems, surpassing prior approaches that focused on complex roots and generic polynomial systems.

Findings

Deterministic algorithm with O(n^2(log n + log max d_i)) complexity for real roots

Applicable to Gaussian coefficients with arbitrary variance

Discusses limitations for systems with more terms

Abstract

Suppose $F:=(f_1,\ldots,f_n)$ is a system of random $n$-variate polynomials with $f_i$ having degree $\leq\!d_i$ and the coefficient of $x^{a_1}_1\cdots x^{a_n}_n$ in $f_i$ being an independent complex Gaussian of mean $0$ and variance $\frac{d_i!}{a_1!\cdots a_n!\left(d_i-\sum^n_{j=1}a_j \right)!}$. Recent progress on Smale's 17th Problem by Lairez --- building upon seminal work of Shub, Beltran, Pardo, B\"{u}rgisser, and Cucker --- has resulted in a deterministic algorithm that finds a single (complex) approximate root of $F$ using just $N^{O(1)}$ arithmetic operations on average, where $N\!:=\!\sum^n_{i=1}\frac{(n+d_i)!}{n!d_i!}$ ($=n(n+\max_i d_i)^{O(\min\{n,\max_i d_i)\}}$) is the maximum possible total number of monomial terms for such an $F$. However, can one go faster when the number of terms is smaller, and we restrict to real coefficient and real roots? And can one still maintain average-case polynomial-time with more general probability measures? We show the answer is yes when $F$ is instead a binomial system --- a case whose numerical solution is a key step in polyhedral homotopy algorithms for solving arbitrary polynomial systems. We give a deterministic algorithm that finds a real approximate root (or correctly decides there are none) using just $O(n^2(\log(n)+\log\max_i d_i))$ arithmetic operations on average. Furthermore, our approach allows Gaussians with arbitrary variance. We also discuss briefly the obstructions to maintaining average-case time polynomial in $n\log \max_i d_i$ when $F$ has more terms.

Equations

$$\Phi (t):=\mathbb{P}(Z\le t):=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^t{e}^{-\frac{s^2}{2}}ds=:{\int }_{-\infty }^t\varphi (s)ds.$$

$$F_Y(t):=\mathbb{P}(Y\le t)=\mathbb{P}(|Z|\le e^t)=\mathbb{P}(-e^{-t}\le Z\le e^t)=1-2\Phi (-e^t).$$

$${\rho }_Y(t):=\sqrt{\frac{2}{\pi }}{e}^t{e}^{-\frac{e^{2t}}{2}},-\infty <t<\infty .$$

$$a:=\mathbb{E}Y:=\sqrt{\frac{2}{\pi }}{\int }_{-\infty }^{\infty }t{e}^t{e}^{-\frac{e^{2t}}{2}}dt.$$

$$0<a<e\sqrt{\frac{2}{\pi }}+2<5.$$

$$\sqrt{\frac{\pi }{2}}a\le {\int }_0^{\infty }t{e}^{-t}dt+{\int }_0^{1}t{e}^{t}dt+{\int }_1^{\infty }t{e}^t{e}^{-\frac{e^{2t}}{2}}dt$$

$$\le 1+(e-1)+{\int }_e^{\infty }\log s{e}^{-\frac{s^2}{2}}ds\le e+\sqrt{2\pi }{\int }_e^{\infty }{s}^2\varphi (s)ds\le e+\sqrt{2\pi }.$$

$$W:=Y-a.$$

$$W_{\mathbf{a}}:=\sum _{i=1}^k{a}_i{Y}_i,V_{\mathbf{a}}:=e^{W_{\mathbf{a}}}\mathrm{and}{X}_{\mathbf{a}}:=\max \{V_{\mathbf{a}},V_{\mathbf{a}}^{-1}\}$$

$$\Vert W{\Vert }_p\simeq \Vert \Theta {\Vert }_p\simeq \Vert L{\Vert }_p.$$

$$\Vert W{\Vert }_p^p=\sqrt{\frac{2}{\pi }}{\int }_{-\infty }^{\infty }|t-a{|}^p{e}^t{e}^{-\frac{e^{2t}}{2}}dt=$$

$$\sqrt{\frac{2}{\pi }}{\int }_0^{\infty }|t+a{|}^p{e}^{-t}{e}^{-\frac{e^{-2t}}{2}}dt+\sqrt{\frac{2}{\pi }}{\int }_0^{\infty }|t-a{|}^p{e}^t{e}^{-\frac{e^{2t}}{2}}dt=:\sqrt{\frac{2}{\pi }}{I}_1+\sqrt{\frac{2}{\pi }}{I}_2.$$

$$I_1\ge {\int }_1^{\infty }|t+a{|}^p{e}^{-t}{e}^{-\frac{e^{-2t}}{2}}dt\ge {\int }_1^{\infty }{t}^p{e}^{-\frac{3t}{2}}dt={\left(\frac{2}{3}\right)}^{p+1}{\int }_{3/2}^{\infty }{s}^p{e}^{-s}ds=$$

$${\left(\frac{2}{3}\right)}^{p+1}\left((p+1)!-{\int }_0^{3/2}{s}^p{e}^{-s}ds\right)\ge \frac{\Vert \Theta {\Vert }_p^p}{4^p},$$

$$I_1\le {\int }_0^{\infty }|t+a{|}^p{e}^{-t}dt=\Vert \Theta +a{\Vert }_p^p\le 2\Vert \Theta {\Vert }_p^p,$$

$$\frac{\Vert \Theta {\Vert }_p^p}{4^p}\le I_1\le \Vert \Theta {\Vert }_p^p.$$

$$I_2\le {\int }_0^1|t-a{|}^p{e}^t{e}^{-\frac{e^{2t}}{2}}dt+{\int }_1^{\infty }|t-a{|}^p{e}^t{e}^{-\frac{e^{2t}}{2}}dt$$

$$\le e{5}^p+\sqrt{2\pi }{\int }_e^{\infty }|\log s-a{|}^p\varphi (s)ds\le 5^p\left(e+\Vert Z{\Vert }_p^p\right)\le 5^p\Vert \Theta {\Vert }_p^p.$$

$$0\le I_2\le 5^p\Vert \Theta {\Vert }_p^p.$$

$${\Vert \sum _{i=1}^n{X}_i\Vert }_p\simeq \Vert |(X_1,\dots ,X_n)\Vert {|}_p,$$

$$a\Vert X_i{\Vert }_q\le \Vert {\tilde{X}}_i{\Vert }_q\le b\Vert X_i{\Vert }_q$$

$$a\Vert |(X_1,\dots ,X_n)\Vert {|}_p\le \Vert |({\tilde{X}}_1,\dots ,{\tilde{X}}_n)\Vert {|}_p\le b\Vert |(X_1,\dots ,X_n)\Vert {|}_p.$$

$$\mathbb{E}\eta \left(a{X}_i/t\right)\le \mathbb{E}\eta \left({\tilde{X}}_i/t\right)\le \mathbb{E}\eta \left(b{X}_i/t\right),1\le i\le n,$$

$$=\frac{1}{2}\sum _{k=0,k\mathrm{even}}^p\left(\genfrac{}{}{0ex}{}{p}{k}\right)\mathbb{E}\left(({\tilde{X}}_i/t{)}^k)+({\tilde{X}}_i/t{)}^k)\right)$$

$$\le \frac{1}{2}\sum _{k=0,k\mathrm{even}}^p\left(\genfrac{}{}{0ex}{}{p}{k}\right)\mathbb{E}\left((b{X}_i/t{)}^k)+(b{X}_i/t{)}^k)\right)$$

$$=\frac{1}{2}\sum _{k=0}^p\left(\genfrac{}{}{0ex}{}{p}{k}\right)\mathbb{E}\left((b{X}_i/t{)}^k)+(-b{X}_i/t{)}^k)\right)=\mathbb{E}\eta \left(b{X}_i/t\right).$$

$$c_{1}a\Vert ⟨\mathbf{X},\theta ⟩{\Vert }_r\le \Vert ⟨\tilde{\mathbf{X}},\theta ⟩{\Vert }_r\le c_{2}b\Vert ⟨\mathbf{X},\theta ⟩{\Vert }_r,$$

$$C_{1}p\Vert {\theta }^p{\Vert }_{\infty }+C_1\sqrt{p}\Vert {\theta }_p{\Vert }_2\le \Vert ⟨L,\theta ⟩{\Vert }_p\le C_{2}p\Vert {\theta }^p{\Vert }_{\infty }+C_2\sqrt{p}\Vert {\theta }_p{\Vert }_2.\text{ ■}$$

$$C_{1}p\Vert {\theta }^p{\Vert }_{\infty }+C_1\sqrt{p}\Vert {\theta }_p{\Vert }_2\le \Vert ⟨W,\theta ⟩{\Vert }_p\le C_{2}p\Vert {\theta }^p{\Vert }_{\infty }+C_2\sqrt{p}\Vert {\theta }_p{\Vert }_2.$$

$$\mathbb{P}\left(\left|\sum _{i=1}^n{\theta }_i\log |Z_i|-\sum _{i=1}^n{\theta }_{i}a\right|\ge t\right)$$

Full text

We show the answer is yes when $F$ is instead a binomial system — a case whose numerical solution is a key step in polyhedral homotopy algorithms for solving arbitrary polynomial systems. We give a deterministic algorithm that finds a real approximate root (or correctly decides there are none) using just $O(n^{2}(\log(n)+\log\max_{i}d_{i}))$ arithmetic operations on average. Furthermore, our approach allows Gaussians with arbitrary variance. We also discuss briefly the obstructions to maintaining average-case time polynomial in $n\log\max_{i}d_{i}$ when $F$ has more terms.

Smale’s 17th Problem, real roots, sparse polynomial, Newton iteration, approximate root

††journalyear: 2019††copyright: acmlicensed††conference: ISSAC ’19: ACM International Symposium on Symbolic and Algebraic Computation; July 15–18, 2019; Beijing, China††booktitle: ISSAC ’19: ACM International Symposium on Symbolic and Algebraic Computation, July 15–18, 2019, Beijing, China††price: 15.00††doi: 10.1145/1122445.1122456††isbn: 978-1-4503-9999-9/18/06††ccs: Theory of computation Algebraic complexity theory

1. Introduction

Polynomial system solving has occupied a good portion of research in algebraic geometry for centuries, and inspired numerous algorithms in engineering and optimization. In recent years, homotopy continuation (see, e.g., (MS, 87; LW, 91; Li, 97; SW, 05; BHSW, 13)) has emerged as one of the most practical and efficient approaches to leverage high performance computing for the approximation of roots of large polynomial systems. A refinement particularly useful for sparse systems is polyhedral homotopy (HS, 95; Ver, 10; LL, 11). To be brutally concise, polyhedral homotopy reduces the solution of an arbitrary polynomial system to (a) solving a finite collection of binomial systems to high precision and then (b) iterating a finite collection of rational functions.

It is thus important to have rigorous and, ideally, optimal complexity estimates for solving binomial systems. Since solving arbitrary polynomial systems is a numerical problem involving solutions of unknown minimal spacing, we will need to incorporate the cost of approximating well enough to distinguish distinct solutions. A recent and elegant way to handle this is via the notion of approximate root in the sense of Smale. In what follows, we use $|\cdot|$ for the standard $\ell_{2}$-norm on $\mathbb{C}^{n}$.

Definition 1.1.

(Sma, 86; BCSS, 98)* Given any analytic function $F:\mathbb{C}^{n}\longrightarrow\mathbb{C}^{n}$, we define the Newton endomorphism of $F$ to be $N_{F}(z)\!:=\!z-F^{\prime}(z)^{-1}F(z)$, where we think of $F(z)$ as a column vector and we identify the derivative $F^{\prime}(z)$ with the matrix of partial derivatives $\left.\left[\frac{\partial f_{i}}{\partial x_{j}}\right]\right|_{x=z}$. We call $\zeta\!\in\!\mathbb{C}^{n}$ a non-degenerate root of $F$ if and only if $F^{\prime}(\zeta)$ is invertible. Given $z_{0}\!\in\!\mathbb{C}^{n}$, we then define its sequence of Newton iterates $(z_{n})_{n\in\mathbb{N}\cup\{0\}}$ via the recurrence $z_{n+1}\!:=\!N_{F}(z_{n})$ (for all $n\!\geq\!0$). We then call $z_{0}$ an approximate root of $F$ in the sense of Smale (with associated true root $\zeta$) if and only if $F$ has a non-degenerate root $\zeta\!\in\!\mathbb{C}^{n}$ satisfying $|z_{n}-\zeta|\!\leq\!\left(\frac{1}{2}\right)^{2^{n-1}}|z_{0}-\zeta|$ for all $n\!\geq\!1$. $\diamond$*

In essence, once one has an approximate root in the sense above, one can easily compute coordinates within any desired $\varepsilon\!>\!0$ of the coordinates of a true root, simply by computing $O\!\left(\log\log\frac{1}{\varepsilon}\right)$ Newton iterates. The special case $F(z_{1})\!:=\!z^{2}_{1}-2$ already shows that one needs $\Omega\!\left(\log\log\frac{1}{\varepsilon}\right)$ arithmetic operations to compute $\sqrt{2}$ within $\varepsilon$ (BMST, 97). So one can arguably consider an approximate root to be the gold standard for specifying a true root. In particular, one no longer has to worry about finding the minimal root spacing of $F$ (to get the right $\varepsilon$ for approximations within $\varepsilon$), since an approximate root in the sense of Smale is guaranteed to converge optimally fast to a unique true root.

Of course, this begs the question of how one can possibly find an approximate root. This is the crux of Smale’s 17${}^{\text{\text@underline{th}}}$ Problem (see (Sma, 98, 00) and Section 1.1 below), which was recently positively solved by Lairez (Lai, 17). (See also the seminal work of Beltran and Shub (Shu, 09; BS, 09), Beltran and Pardo (BP, 08, 09, 11) and Bürgisser and Cucker (BC, 12).) Roughly, Lairez’s discovery was an algorithm that, for a certain class of random polynomial systems, finds a single (complex) approximate root in polynomial-time on average. We now introduce some more terminology to be precise:

Definition 1.2.

Suppose $\mathcal{A}_{1},\ldots,\mathcal{A}_{n}\!\subset\!\mathbb{Z}^{n}$ are finite subsets and $\{c_{i,a}\;|\;i\!\in\!\{1,\ldots,n\}\text{ and }a\!\in\!\mathcal{A}_{i}\text{ for all }i\}$ is a collection of independent complex Gaussians with mean [math] and the variance of $c_{i,a}$ equal to $v_{i,a}$. Letting $a\!:=\!(a_{1},\ldots,a_{n})$, $x^{a}\!:=\!x^{a_{1}}_{1}\cdots x^{a_{n}}_{n}$, and $f_{i}(x)\!:=\!\sum_{a\in\mathcal{A}_{i}}c_{i,a}x^{a}$, we call $F\!:=\!(f_{1},\ldots,f_{n})$ an $n\times n$ random polynomial system with support $(\mathcal{A}_{1},\ldots,\mathcal{A}_{n})$. $\diamond$

Lairez’s Theorem .

(Lai, 17, Thm. 23)111We have paraphrased a bit: Lairez’s main theorem is stated in terms of homogeneous polynomials, and he counts square roots as arithmetic operations as well. Via the techniques of, say, (BP, 09), one can easily derive our affine statement. Following the notation above, let $d_{1},\ldots,d_{n}\!\in\!\mathbb{N}$, $\mathcal{A}_{i}\!:=\!\{(a_{1},\ldots,a_{n})\!\in\!(\mathbb{N}\cup\{0\})^{n}\;|\;\sum^{n}_{j=1}a_{j}\!\leq\!d_{i}\}$ for all $i$, and $v_{i,a}\!:=\!\frac{d_{i}!}{a_{1}!\cdots a_{n}!\left(d_{i}-\sum^{n}_{j=1}a_{j}\right)!}$. Then one can find a (complex) approximate root of $F$ using just $O(nd^{3/2}N(N+n^{3}))$ arithmetic operations on average, where $N\!:=\!\sum^{n}_{i=1}\frac{(d_{i}+n)!}{d_{i}!n!}$ and $d\!:=\!\max_{i}d_{i}$. $\blacksquare$

Note that restricting the support $(\mathcal{A}_{1},\ldots,\mathcal{A}_{n})$ is a way to consider sparsity for one’s polynomial system. In particular, one can think of Lairez’s Theorem as solving Smale’s 17${}^{\text{\text@underline{th}}}$ Problem in the “dense” case, since Lairez assumes that all monomial terms up to a given degree appear (with probability $1$) in each polynomial $f_{i}$. Indeed, one should note that Smale never specified what kind of probability measure one should use in his 17${}^{\text{\text@underline{th}}}$ Problem (Sma, 98, 00). So Smale’s 17${}^{\text{\text@underline{th}}}$ Problem actually includes sparse systems if some of the random coefficients have mean, and all higher moments, equal to [math]. Smale also observed that one can pose a more difficult analogue of his 17${}^{\text{\text@underline{th}}}$ problem over the real numbers.

Observe that $\sum^{n}_{i=1}\frac{(d_{i}+n)!}{d_{i}!n!}$ is exactly the maximal possible total number of monomial terms in an $n\times n$ polynomial system where $f_{i}$ has degree $d_{i}$. Note also that just evaluating a monomial of degree $d$ takes $\Omega(\log d)$ arithmetic operations: Simply consider the straight-line program complexity of the integer $2^{d}$ (see, e.g., (Bra, 39; dMS, 96; Mor, 97)). One should pay attention to the evaluation complexity of $F$ since Lairez’s algorithm uses Newton iteration, which in turn requires evaluating $F$ (and $F^{\prime}$) many times. So one can then naturally ask, in the spirit of real fewnomial theory (Kho, 91): Can one find a real approximate root of $F$ (or decide whether there is no real root) using, say, $(t\log d)^{O(1)}$ arithmetic operations on average, when $t$ is the total number of monomial terms of $F$ and $d\!:=\!\max_{i}d_{i}$? This would be a significant new speed-up. For instance, the special case $t\!=\!O(n)$ is already quite non-trivial since there are standard algebraic tricks (e.g., the bottom of the first page of (ES, 96)) to reduce arbitrary polynomial systems to trinomial systems.

Our first main theorem thus solves a special case of a refined version of Smale’s 17${}^{\text{\text@underline{th}}}$ Problem, and serves as a starting point for a deeper study of the randomized complexity of solving arbitrary real sparse polynomial systems.

Theorem 1.3.

Suppose $A\!=\![a_{i,j}]\!\in\!\mathbb{Z}^{n\times n}$ has nonzero determinant, and all the entries of $A$ have absolute value at most $d$. Suppose also that $c_{i,j}$ is an independent real Gaussian with mean [math] and fixed (but otherwise arbitrary) variance, for each $(i,j)\!\in\!\{1,\ldots,n\}\times\{0,1\}$. Let $F\!:=\!(f_{1},\ldots,f_{n})$ with $f_{i}(x)\!:=\!c_{i,0}+c_{i,1}\cdot x^{a_{1,1}}_{1}\cdots x^{a_{1,n}}_{n}$. Then, on average, one can find a real approximate root of $F$ (or correctly determine there are no real roots) using just $O(n^{2}\log(nd))$ arithmetic operations and $O(n^{\omega+1}\log^{2}(dn))$ bit operations, where $\omega$ is any upper bound on the matrix multiplication exponent.

We prove Theorem 1.3 in Section 4. The best current upper bound on $\omega$, as of January 2019, is $2.372873$ (Vas, 14). A fundamental ingredient behind our proof of Theorem 1.3 is a hybrid algorithm of Ye enabling the quick approximation of rational powers of a real number (Ye, 94), combined with some classic results on fast linear algebra over $\mathbb{Z}$ (Smi, 61; Sto, 00). A final key ingredient is estimating the expected value of linear combinations of logarithms of absolute values of standard real Gaussians (Proposition 3.9 in Section 3.3 below). We were unable to find any explicit asymptotics for such expectations, so we derive these from scratch in the latter half of Section 2 and Section 3.

We will explain some of the subtleties behind extending Theorem 1.3 to systems with arbitrary supports in Section 1.2 below. First, however, let us briefly review the original statement of Smale’s 17${}^{\text{\text@underline{th}}}$ Problem.

1.1. Quick Review of Smale’s 17${}^{\text{\text@underline{th}}}$ Problem

Smale’s 17${}^{\text{\text@underline{th}}}$ Problem (Sma, 98, 00) elegantly summarizes the subtleties behind polynomial system solving:

Can a zero of $n$ complex polynomial equations in $n$ unknowns be found approximately, on the average, in polynomial-time with a uniform algorithm? [Emphases added.]

We clarify the notion of “polynomial-time” below. As motivation, let us first see how the emphasized terms highlight fundamental difficulties in polynomial system solving:

• “a zero”: We can not expect a fast algorithm approximating all the roots since, for $n\!\geq\!2$, there may be infinitely many. In which case, for $d_{1}\!\geq\!3$ (e.g., the case of elliptic curves (ST, 94)), the roots will likely not admit a rational parametrization. When there are only finitely many roots, systems like$(x^{2}_{1}-1,\ldots,x^{2}_{n}-1)$ show that the number of roots can be exponential in $n$.
• “found approximately”: Even restricting to integer coefficients, the number of digits of accuracy needed to separate distinct roots can be exponential in $n$, e.g.,

$((2x_{1}-1)(3x_{1}-1),x_{2}-x^{2}_{1},\ldots,x_{n}-x^{2}_{n-1})$

has roots with $n^{\text{\text@underline{th}}}$ coordinates $\frac{1}{2^{2^{n-1}}}$ and $\frac{1}{3^{2^{n-1}}}$. So, especially for irrational coefficients, we need a more robust notion of approximation than digits of accuracy. (Hence’s Smale’s definition of approximate root from (Sma, 86).)

• “on the average”: Restricting to integer coefficients, distinguishing between a system having finitely many or infinitely many roots is ${\mathbf{NP}}$-hard (see, e.g., (Pla, 84; Koi, 97)). Furthermore, as already long known in the numerical linear algebra community (e.g., results on the distribution of eigenvalues of random matrices (Ede, 88; TV, 09)), even if the number of roots is finite, the accuracy needed to separate distinct roots can vary wildly as a function of the coefficients. So averaging over all inputs allows us to amortize the complexity of potentially intractable instances.

The original statement of Smale’s 17${}^{\text{\text@underline{th}}}$ Problem measures time (or complexity) as the total number of (a) (exact) field operations over $\mathbb{R}$, (b) comparisons over $\mathbb{R}$, and (c) bit operations (Sma, 98). (The underlying computational model is a BSS machine over $\mathbb{R}$ (BCSS, 98), which is essentially a classical Turing machine (Pap, 95; AB, 09; Sip, 12), augmented so that it can perform any field operation or comparison over $\mathbb{R}$ in one time step.) Polynomial-time was then meant as polynomial in the number of (nonzero) coefficients of $F$. Smale interpreted the number of coefficients (which can be as high as $\sum^{n}_{i=1}\binom{d_{i}+n}{n}$ for $F$ as specified above) as the input size.

Remark 1.4.

The precise probability distribution over which one averages was never specified in Smale’s original statement (Sma, 98, 00). In all the literature so far on the problem (see, e.g., (Sma, 98, 00; BP, 08, 09; Shu, 09; BS, 09; BP, 11; BC, 12; Lai, 17)), the Bombieri-Weyl measure was used: This is the choice of variances involving multinomial coefficients written earlier. $\diamond$

While the Bombieri-Weyl measure satisfies some very nice group invariance properties (see, e.g., (Kos, 93; SS92b, ; BSZ, 00; FLL, 14)), there is currently no widely-accepted notion of a “natural” probability distribution for a random polynomial. For instance, there are several different distributions of interest already for the matrix eigenvalue problem (see, e.g., (Ede, 88; Roj, 96; ABF, 11)). More to the point, much work has gone into finding useful properties of the roots of random polynomials that are distribution independent (see, e.g., (B-RS, 86; TV, 13)).

The meaning of uniform algorithm is more technical and is formalized in (BCSS, 98) (see also (Pap, 95; AB, 09; Sip, 12) for the classical Turing case). Roughy, uniformity refers to having an implementation that can handle all input sizes, as opposed to having a different implementation for each input size.

1.2. Current Obstructions to Fully Incorporating Sparsity

As we’ll see from the proof of our main theorem, solving an $n\times n$ system of Gaussian random binomials of degree $d$ can be reduced to solving $n$ univariate binomials of degree $(nd)^{O(n)}$, where the underlying coefficients are no longer Gaussian but have reasonably estimable means. Algebraically, this will imply that the underlying field extension (where one adjoins the coordinates of the solutions to the field generated by the coefficients) is always a radical extension.

A natural next step then is to consider $n\times n$* unmixed $(n+1)$-nomial systems*:

$(c_{1,0}+c_{1,1}x^{a_{1}}+\cdots+c_{1,n}x^{a_{n}},\ldots,c_{n,0}+c_{n,1}x^{a_{1}}+\cdots+c_{n,n}x^{a_{n}})$, where $a_{i}\!:=\!(a_{1,i},\ldots,a_{n,i})$ for all $i$. Via Gauss-Jordan Elimination, one can reduce such a system to a binomial system without affecting the roots. Unfortunately, if one starts with a system of the form above, with Gaussian $c_{i,j}$, the resulting binomial system no longer has Gaussian coefficients. So one needs to consider binomial systems with coefficient distributions more general than Gaussian, and we do this in a sequel to this paper.

Going a bit farther, $n\times n$ unmixed $(n+2)$-nomial systems yield an interesting complication: The underlying field extensions need no longer be radical, even if $n\!=\!1$. A simple example is $x^{5}_{1}-2x_{1}+10$, which has Galois group $S_{5}$ over $\mathbb{Q}$. However, earlier results from (RY, 05) indicate that it should be possible to find real approximate roots quickly on average, at least for univariate trinomials. (One should also observe Sagraloff’s recent dramatic speed-ups for the worst-case arithmetic complexity of approximating real roots of univariate sparse polynomials (Sag, 14).) We conjecture that finding a real approximate root (or determining that there are no real roots) for a real Gaussian $n\times n$ unmixed $(n+2)$-nomial system is still possible in time $(n\log d)^{O(1)}$ on average, and hope to address this problem in the future.

2. Background

In what follows, for any $n\times n$ matrix $A\!\in\!\mathbb{Z}^{n\times n}$, we define $x^{A}$ to be the vector of monomials $\left(x^{a_{1,1}}_{1}\cdots x^{a_{n,1}}_{n},\ldots,x^{a_{1,n}}_{1}\cdots x^{a_{n,n}}_{n}\right)$. We call the substitution $x\!=\!z^{A}$ a monomial change of variables. The following proposition is elementary.

Proposition 2.1.

We have that $x^{AB}\!=\!(x^{A})^{B}$ for any $A,B\!\in\!\mathbb{Z}^{n\times n}$. Also, for any field $K$, the map defined by $m(x)\!=\!x^{U}$, for any unimodular matrix $U\!\in\!\mathbb{Z}^{n\times n}$, is an automorphism of $(K^{*})^{n}$. $\blacksquare$

Our main approach to solving binomial systems is to reduce them to systems of the form $(x^{d_{1}}_{1}-c_{1},\ldots,x^{d_{n}}_{n}-c_{n})$ via a monomial change of variables, and then prove that the distortion of the $c_{i}$ resulting from perturbing the original coefficients is controllable. Later on, we will also detail how a Gaussian distribution on the original coefficients implies that the $c_{i}$ still have well-behaved distributions. But now we will focus on quantifying our monomial changes of variables.

2.1. Linear Algebra Over $\mathbb{Z}$

Definition 2.2.

Let $\mathbb{G}\mathbb{L}_{n}(\mathbb{Z})$ denote the set of all matrices in $\mathbb{Z}^{n\times n}$ with determinant $\pm 1$ (the set of unimodular matrices). Given any $M\!\in\!\mathbb{Z}^{n\times n}$, we call any identity of the form $UMV\!=\!S$ with $U,V\!\in\!\mathbb{G}\mathbb{L}_{n}(\mathbb{Z})$ and $S$ diagonal a Smith factorization. In particular, if $S\!=\![s_{i,j}]$ and we require additionally that $s_{i,i}\!\geq\!0$ and $s_{i,i}|s_{i+1,i+1}$ for all $i\!\in\!\{1,\ldots,n\}$ (setting $s_{n+1,n+1}\!:=\!0$), then $S$ is uniquely determined and is called the Smith normal form of $M$. $\diamond$

Theorem 2.3.

(Sto, 00, Ch. 6 & 8, pg. 128)* For any $A\!=\![a_{i,j}]\!\in\!\mathbb{Z}^{n\times n}$, a Smith factorization of $A$ yielding the Smith normal form of $A$ can be computed within

$O\!\left(n^{\omega+1}\log^{2}(n\max_{i,j}|a_{i,j}|)\right)$

bit operations. Furthermore, the entries of all matrices in the underlying factorization have bit size $O(n\log(n\max_{i,j}|a_{i,j}|))$. $\blacksquare$*

2.2. From Approximate Roots of Univariate Binomials to Systems

We begin with an important observation from the middle author’s doctoral dissertation, building upon earlier work of Smale (Sma, 86) and Ye (Ye, 94).

Lemma 2.4.

(Phi, 16, Thm. 4.10)*

Let $d\!\in\!\mathbb{N}$ satisfy $d\!\geq\!2$, $c\!>\!0$, and $f(x_{1})\!:=\!x^{d}_{1}-c$. Then we can find an approximate root of $f$ using $O\!\left(\log(d)+\log\log\max\{c,c^{-1}\}\right)$ field operations over $\mathbb{R}$. $\blacksquare$*

Since a monomial change of variables enables us to replace an arbitrary binomial system by a simpler, diagonal system of univariate binomials, it’s enough to bound how the coefficients are distorted under such a change of variables. The following lemma gives us the bounds we need.

Lemma 2.5.

*Suppose $c_{1},\ldots,c_{n}\!\in\!\mathbb{C}^{*}$ and $A\!\in\!\mathbb{Z}^{n\times n}$ has columns $a_{1},\ldots,a_{n}$ and entries of absolute value at most $d$. Also let $\sigma\!:=\!\max_{i}\{|\log|c_{i}||\}$, let $UAV\!=\!S$ be the Smith Factorization of $A$, and let $(\gamma_{1},\ldots,\gamma_{n})\!:=\!(c_{1},\ldots,c_{n})^{V}$. Then the following bounds hold:

     1. $\max_{i}|\log|\gamma_{i}||\!\leq\!n^{4+3n/2}d^{3n}\sigma$.

     2. If $\zeta\!=\!(\zeta_{1},\ldots,\zeta_{n})\!\in\!{(\mathbb{C}^{*})}^{n}$ is a true root of $F$ then

       $\max_{i}|\log|\zeta_{i}||\!\leq\!n^{O(n)}d^{O(n)}\sigma$. $\blacksquare$*

2.3. Logs of Absolute Values of Gaussians

We now finally address the change in probability distribution resulting from replacing a Gaussian coefficient by a monomial in several other Gaussians. Our derivation is, necessarily, a bit long. So the hurried reader can jump to Propositions 3.3 and 3.9, respectively in Sections 3 and 3.3 below.

Let $\Theta,\Theta_{1},\Theta_{2},\ldots$ be independent exponential random variables, i.e., $F_{\Theta}(t):=\mathbb{P}(\Theta\leq t)=1-e^{-t}$. Let $L,L_{1},L_{2},\ldots$ be independent symmetric exponential random variables, i.e., the density of $L$ is given by $\rho_{L}(t)=\frac{1}{2}e^{-|t|},\ -\infty<t<\infty$, and similarly for the $L_{i}$.

Let $Z,Z_{1},Z_{2},\ldots$ be independent standard real Gaussian random variables, i.e.,

[TABLE]

Let $Y,Y_{1},Y_{2},\ldots$ be independent random variables such that $Y_{i}:=\log|Z_{i}|$. We have that

[TABLE]

Taking derivatives we get $F_{Y}^{\prime}(t):=2e^{t}\phi(e^{t})$, which implies that the density of $Y$, $\rho_{Y}$, is

[TABLE]

We use $\mathbb{E}$ to denote expectation, $\mathbb{P}$ to denote probability, and define

[TABLE]

Note that

[TABLE]

Indeed, $a=\sqrt{\frac{2}{\pi}}\left(\int_{0}^{\infty}te^{-t}e^{-\frac{e^{-2t}}{2}}dt+\int_{0}^{\infty}te^{t}e^{-\frac{e^{2t}}{2}}dt\right)$. So clearly $a>0$ and also

[TABLE]

[TABLE]

We define a new, centered (i.e., mean [math]) random variable via

[TABLE]

We write $A\simeq B$ to indicate that there exist positive constants $c_{1},c_{2}$ with $c_{1}A\leq B\leq c_{2}A$. Let ${\bf a}:=(a_{1},\ldots,a_{k})\!\in\!\mathbb{R}^{k}$ and define

[TABLE]

Using the notation $\|R\|_{p}:=E(R^{p})^{1/p}$, we will prove the following fact:

Lemma 2.6.

Let $W$ be the centered random variable defined in (6) and let $p\geq 2$. Then

[TABLE]

Proof.

We have that

[TABLE]

[TABLE]

Note that for $t\geq 1$, $e^{-\frac{e^{-2t}}{2}}\geq e^{-\frac{t}{2}}$. So using the above we have that

[TABLE]

[TABLE]

since $\|\Theta\|_{p}^{p}=(p+1)!$. Moreover

[TABLE]

by Minkowski inequality the fact that $p\geq 2$ and (4). So we have shown that

[TABLE]

Moreover, using again (4),

[TABLE]

[TABLE]

So we have shown that

[TABLE]

Combining (8) and (9) we get that $\frac{\|\Theta\|_{p}}{4}\leq\|W\|_{p}\leq 6\|\Theta\|_{p}$. Finally, it is straightforward to check that $\|\Theta\|_{p}\simeq\|L\|_{p}$ for all $p>0$. $\blacksquare$

2.4. A Tool for Linear Combinations of Logs of Absolute Values of

Gaussians

We are going to use the following fundamental result of Latala:

Theorem 2.7.

(Lat, 97, Thm. 2 & Rem. 2)* Let $X_{1},\ldots,X_{n}$ be centered independent random variables and $p\geq 2$. Then*

[TABLE]

*where $\||(X_{1},\ldots,X_{n})\||_{p}$ is defined to be

$\inf\left\{t>0:\sum_{i=1}^{n}\log\left(\mathbb{E}\frac{\left|\frac{X_{i}}{t}+1\right|^{p}+\left|\frac{-X_{i}}{t}+1\right|^{p}}{2}\right)\leq p\right\}$. $\blacksquare$*

We will also need the following fact:

Lemma 2.8.

Let $X_{1},\ldots,X_{n}$ be independent random variables and let $\tilde{X}_{1},\ldots,\tilde{X}_{n}$ be another sequence of independent random variables. Fix $p\geq 2$ be an even integer and assume that there are $a,b>0$ such that

[TABLE]

for all $1\leq q\leq p$ and for all $1\leq i\leq n$. Then we have that

[TABLE]

Proof.

We will first prove the following

Claim: Under the assumptions of the Lemma we have that for every $t>0$

[TABLE]

where $\eta(x):=\frac{1}{2}\left(|x+1|^{p}+|1-x|^{p}\right)$.

Indeed, $\mathbb{E}\eta\!\left(\tilde{X}_{i}/t\right)=\frac{1}{2}\sum_{k=0}^{p}{p\choose k}\mathbb{E}\left((\tilde{X}_{i}/t)^{k})+(-\tilde{X}_{i}/t)^{k})\right)$

[TABLE]

[TABLE]

[TABLE]

The proof of the other side inequality in (13) is identical. Equation (12) then follows immediately from the claim and the definition of $\||(X_{1},\ldots,X_{n})\||_{p}$. $\blacksquare$

Our preceding lemma leads to the following:

Corollary 2.9.

Let ${\bf X}:=(X_{1},\ldots,X_{n})$, ${\bf\tilde{X}}:=(\tilde{X}_{1},\ldots,\tilde{X}_{2})$ be two centered random vectors with independent coordinates and let $\theta=(\theta_{1},\ldots,\theta_{n})\in\mathbb{R}^{n}$. We assume that (11) holds true. Then for every $1\leq r\leq p$,

[TABLE]

where $c_{1},c_{2}>0$ are universal constants.

Proof.

The result follows from Theorem 2.7 and Lemma 2.8 applied to the random variables $\theta_{i}X_{i}$ and $\theta_{i}\tilde{X_{i}}$. $\blacksquare$

3. Additional Probabilistic Estimates

Let ${\bf W}:=(W_{1},\ldots,W_{n})$ be the centered random vector with independent entries that are logs of absolute values of real standard Gaussians. Let ${\bf L}:=(L_{1},\ldots,L_{n})$. Let $\theta\in S^{n-1}$ (Here $S^{n-1}$ is the unit sphere in dimension $n$.) The next theorem below is a special case of a more general result of Gluskin and Kwapien (GK, 95). Let us introduce some notation. Let $x\in\mathbb{R}^{n}$. We write $x^{\ast}$ for the non-increasing rearrangement of the vector $(|x_{1}|,\ldots,|x_{n}|)$. Given any $1\leq s\leq n$ and a vector $x$ we denote $x^{s}$ the vector with entries $x_{i}^{\ast}$ for $i\leq s$ and [math] otherwise and by $x_{s}$ the vector with entries [math] for $i\leq s$ and entries $x^{\ast}$ for $i>s$.

Theorem 3.1.

*(Special case of (GK, 95))

There are constants $C_{1},C_{2}>0$ such that for every $n\geq 1$, $p\geq 1$, and every $\theta\in S^{n-1}$, one has that*

[TABLE]

Lemma 2.6, Theorem 3.1, and Corollary 2.9 together imply the following

Proposition 3.2.

There exists two constants $C_{1},C_{2}>0$ such that for every $n\geq 1$, $p\geq 1$ and every $\theta\in S^{n-1}$, one has that

[TABLE]

The above result gives very precise estimates about the concentration of the function

[TABLE]

for all $t$. A less precise but simpler to use statement than Theorem 3.1 is the following estimate: For every $\theta\in S^{n-1}$

[TABLE]

Using the above we arrive at the following

Proposition 3.3.

Let $Z_{1},\ldots,Z_{n}$ be independent standard real Gaussian random variables, $\theta\in S^{n-1}$, and $a$ as defined in (3). Then the following holds:

[TABLE]

for $t\!>\!0$, where $C,C^{\prime}>0$ are absolute constants.

Proof.

By (16) we have that $\|\langle W,\theta\rangle\|_{p}\leq C_{2}\sqrt{p}$ if $p\leq\|\theta\|_{\infty}^{-2}$ and $\|\langle W,\theta\rangle\|_{p}\leq C_{2}p\|\theta\|_{\infty}$ otherwise. Using Markov’s Inequality we get that $\mathbb{P}\!\left(|\langle W,\theta\rangle|\geq eC_{2}\sqrt{p}\right)\leq e^{-p},\ {\rm if}\ p\leq\|\theta\|_{\infty}^{-2}$, or (if we will set $eC_{2}\sqrt{p}=t$), for $t\geq C_{3}$,

[TABLE]

and $\mathbb{P}\!\left(|\langle W,\theta\rangle|\geq eC_{2}p\|\theta\|_{\infty}\right)\leq e^{-p},\ {\rm if}\ p\geq\|\theta\|_{\infty}^{-2}$ or (if we will set $eC_{2}p\|\theta\|_{\infty}=t$), for $t\geq\frac{C_{4}}{\|\theta\|_{\infty}}$,

[TABLE]

Combining (19) and (20) and adjusting the constants properly we get (17). $\blacksquare$

3.1. On the Expectation of $\log\log$

Let $Z,Z_{i}$ be independent real standard Gaussian random variables and let $d$ be a positive integer.

Let $X:=\max\{|Z|,|Z|^{-1}\}$. We have that

[TABLE]

Indeed, for $t\geq 0$, we have

$\mathbb{P}\!\left(\log\log\{eX\}\geq t\right)=\mathbb{P}\!\left(X\geq e^{e^{t}-1}\right)=\mathbb{P}\!\left(|Z|\geq e^{e^{t}-1}\ {\rm or}\ |Z|\leq e^{-(e^{t}-1)}\right)$

[TABLE]

$\leq 4\frac{1}{\sqrt{2\pi}}e^{-(e^{t}-1)}$. So, we get that

[TABLE]

Indeed, since $\log\log\{eX\}\geq 0$, using (21),

[TABLE]

[TABLE]

In what follows we assume that

[TABLE]

We will use the following elementary inequality:

[TABLE]

Since $eX\geq e$ and $d/e\geq e$, using (24) and (22) we get

[TABLE]

[TABLE]

[TABLE]

Moreover, since $X\geq 1$, $\log\log\{dX\}\geq\log\log{d}$ and we conclude that

[TABLE]

3.2. Log-Concavity

A Borel measure $\mu$ in $\mathbb{R}^{n}$ is called log-concave if for every compact sets $A,B$ and $\lambda\in(0,1)$ one has

[TABLE]

Theorem 3.4 (Borell (Bor, 75)).

Let $\mu$ be a Borel measure in $\mathbb{R}^{n}$ that gives positive mass to some open ball. Then $\mu$ is $log$-concave if and only if has a density $\rho_{\mu}$ that is a $log$-concave function i.e. $\rho_{\mu}$ is non-negative, supported on a convex set and

[TABLE]

Theorem 3.5 (Prékopa).

Sum of independent $log$-concave random variables is $log$-concave.

Proposition 3.6.

Let $\mu$ be a $log$-concave probability measure and let $K$ be a symmetric closed convex set in $\mathbb{R}^{n}$. Then if $\delta:=\mu(K)\geq\frac{1}{2}$ for every $t>1$ we have that

[TABLE]

Corollary 3.7.

Let $X$ be a $log$-concave random variable with mean [math] and variance $\gamma^{2}$. Then

[TABLE]

Proof.

Let $A:=\{|x|\leq 2\gamma\}$. Then, by Chebychev’s inequality we have that $\mathbb{P}(A)=\delta\geq\frac{3}{4}$. By Proposition 3.6 we get that

[TABLE]

$\blacksquare$

3.3. Final Estimates

Recall that if $Z$ is a standard Gaussian then $Y:=\log|Z|$ has density

[TABLE]

Let $a:=\mathbb{E}[Y]$ and $\tau^{2}$ be the variance of $Y$.

We have the following

Proposition 3.8.

Let ${\bf a}\in\mathbb{R}^{k}$ and assume that $\sum_{i=1}^{k}a_{i}=0$. Then $W_{{\bf a}}$ is a $log$-concave random variable with expectation [math] and variance $\gamma^{2}:=\|{\bf a}\|_{2}^{2}\tau^{2}$. Then we have

[TABLE]

Moreover,

[TABLE]

Proof.

Note that $v(t):e^{\frac{2t}{2}}-t$ is a convex function so by Borell’s theorem $Y$ is a log-concave random variable. We have that

[TABLE]

and since $Y_{i}$ are independent

[TABLE]

So, we can estimate as follows:

[TABLE]

[TABLE]

[TABLE]

as long $e^{t}-1\geq\gamma$, where we have also used Corollary 3.7. Finally, since $eX_{\bf a}\geq e$, $\log\log\{X_{\bf a}\}\geq 0$, we have that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\blacksquare$

Working as in Subsection 3.1 we arrive at the following key result:

Proposition 3.9.

Let ${\bf a}\in\mathbb{R}^{k}$ and assume that $\sum_{i=1}^{k}a_{i}=0$. We have that

[TABLE]

Proof.

Since $eX\geq e$ and $d/e\geq e$, using (24), (30) we get

[TABLE]

[TABLE]

[TABLE]

Moreover, since $X_{\bf a}\geq 1$, $\log\log\{dX_{\bf a}\}\geq\log\log{d}$ we get (31). $\blacksquare$

4. The Proof of Theorem 1.3

First note that the $c_{i,j}$ are all nonzero with probability $1$, so we may assume (since we are considering average-case complexity) that all the $c_{i,j}$ are nonzero. In which case, we can focus solely on roots in ${(\mathbb{R}^{*})}^{n}$.

Now note that by Proposition 2.1, we can easily decide whether our input binomial system $F$ has a real root: If $F$ is diagonal, i.e., if $F\!=\!(c_{1,0}+c_{1,1}x^{d_{1}}_{1},\ldots,c_{n,0}+c_{n,1}x^{d_{n}})$, then $F$ has a real root if and only if $c_{i,0}c_{i,1}\!<\!0$ for all $i$ with $d_{i}$ even. (In which case, each orthant of $\mathbb{R}^{n}$ contains at most $1$ root of $F$.) If $F$ is not diagonal, then after computing a Smith factorization $UAV\!=\!S$ (which accounts for our stated bit complexity bound, thanks to Theorem 2.3), we can reduce to the diagonal case and simply check $n$ inequalities. If there are no real roots, no further work needs to be done.

So let us now assume that there are real roots. Without loss of generality, we may assume there is a root in the positive orthant $\mathbb{R}^{n}_{+}$. This will be the root we will try to approximate. So we may now assume that we are trying to approximate the roots of $G\!:=\!(z^{s_{1,1}}_{1}-\gamma_{1},\ldots,z^{s_{n,n}}_{n}-\gamma_{n})$ where

$(\gamma_{1},\ldots,\gamma_{n})\!:=\!(-c_{1,0}/c_{1,1},\ldots,-c_{n,0}/c_{n,1})^{V}$

lies in $\mathbb{R}^{n}_{+}$, and the $s_{i,i}$ are the diagonal entries of the Smith normal form $S$ of $A$. In particular, we need to approximate the unique root $\mu$ of $G$ in $\mathbb{R}^{n}_{+}$ well enough so that $\zeta\!:=\!\mu^{U}$ is an approximate root of $F$.

Thanks to Lemmata 2.4 and 2.5, a quick derivative calculation tells us that it suffices to find an approximate root of $G$. (One needs some extra precision to ensure that $\zeta$ is an approximate root of $F$ but the bounds from Lemma 2.5 easily imply that the necessary extra work is negligible compared to our stated arithmetic complexity bound.) So it suffices to compute an upper bound on the expectation of $\sum^{n}_{i=1}\left[\log(|s_{i,i}|)+\log\log(e\max\{|\gamma_{i}|,|\gamma^{-1}_{i}|\})\right]$. We are almost done, save for the fact that the $\gamma_{i}$ are monomials in real Gaussians that need not have variance $1$.

However, we can precede our construction of $G$ with another renormalization to reduce to the variance $1$ case: Observe that $x^{A}\!=\!c$ if and only if $(rx)^{A}\!=\!r^{A}c$, for any $r\!\in\!{(\mathbb{R}^{*})}^{n}$. So if we take $r$ to be a suitable matrix power of a vector of ratios of variances, we can replace our original binomial system $F$ by a new binomial system $\tilde{F}$ with all coefficients being standard real Gaussians (and new root a rescaling of our old root). In particular, we merely take $r\!:=\!(v_{1,1}/v_{1,0},\ldots,v_{n,1}/v_{n,0})^{A^{-1}}$. Lemmata 2.4 and 2.5 once again imply that the cost of the necessary increase in precision to convert an approximate root of $\tilde{F}$ to an approximate root of $F$ is negligible.

We now conclude via Proposition 3.9 and Theorem 2.3: Our desired expectation is at most

  $\sum^{n}_{i=1}\left[n\log(nd)+0+2+\log(2)+\log\left(1+\tau\sqrt{n}e^{O(n\log(nd))}\right)\right]$.

The last quantity is clearly $O(n^{2}\log(nd)+n\log(\sqrt{n}n\log(nd)))$ or, more simply, $O(n^{2}\log(nd))$. $\blacksquare$

Acknowledgements

We humbly thank REU students Caleb Bugg and Paula Burkhardt for important discussions on preliminary versions of this work. We also thank the NSF for their support through grants CCF-1409020, DMS-1460766, DMS-1757872, and DMS-1812240.