Plantinga-Vegter algorithm takes average polynomial time

TL;DR

This paper provides a condition-based analysis showing that the Plantinga-Vegter algorithm for plane curves has an expected polynomial time complexity of O(d^7), contrasting with its worst-case exponential bound, supported by probabilistic and smoothed analysis.

Contribution

It introduces the first expected polynomial time complexity analysis for the PV Algorithm using probabilistic and smoothed analysis techniques.

Findings

Expected time complexity is bounded by O(d^7) for degree d curves.

Smoothed analysis yields similar polynomial bounds.

The approach combines probabilistic methods with condition numbers and continuous amortization.

Abstract

We exhibit a condition-based analysis of the adaptive subdivision algorithm due to Plantinga and Vegter. The first complexity analysis of the PV Algorithm is due to Burr, Gao and Tsigaridas who proved a $O\big(2^{\tau d^{4}\log d}\big)$ worst-case cost bound for degree $d$ plane curves with maximum coefficient bit-size $\tau$. This exponential bound, it was observed, is in stark contrast with the good performance of the algorithm in practice. More in line with this performance, we show that, with respect to a broad family of measures, the expected time complexity of the PV Algorithm is bounded by $O(d^7)$ for real, degree $d$, plane curves. We also exhibit a smoothed analysis of the PV Algorithm that yields similar complexity estimates. To obtain these results we combine robust probabilistic techniques coming from geometric functional analysis with condition numbers and the continuous amortization paradigm introduced by Burr, Krahmer and Yap. We hope this will motivate a fruitful exchange of ideas between the different approaches to numerical computation.

Equations

$$C_f(J)\text{: either }0\notin \square [hf](I)\text{ or }0\notin ⟨\square [h^{\prime }\partial f](J),\square [h^{\prime }\partial f](J)⟩$$

$$(\mathbb{E}|X{|}^p{)}^{\frac{1}{p}}\le K\sqrt{p}.$$

$$\max \left\{\mathbb{P}\left(|X-u|\le \epsilon \right)\mid u\in \mathbb{R}\right\}\le \rho \epsilon .$$

$$f:=\sum _{|\alpha |=d}{\left(\genfrac{}{}{0ex}{}{d}{\alpha }\right)}^{1/2}{c}_{\alpha }{X}^{\alpha }$$

$$d^{\frac{n^2+3n}{2}}\max \{1,a^n\}{2}^{\frac{n^2+16n\log (n)}{2}}(c_1{c}_{2}K\rho {)}^{n+1}$$

$$d^{\frac{n^2+5n}{2}}\max \{1,a^n\}{2}^{\frac{7n^2+9n\log (n)}{2}}(c_1{c}_{2}K\rho {)}^{n+1}$$

$$d^{\frac{n^2+3n}{2}}\max \{1,a^n\}{2}^{\frac{n^2+16n\log (n)}{2}}{\left(c_1{c}_{2}K\rho \right)}^{n+1}{\left(1+\frac{1}{\sigma }\right)}^{n+1}$$

$$d^{\frac{n^2+5n}{2}}\max \{1,a^n\}{2}^{\frac{7n^2+9n\log (n)}{2}}(c_1{c}_{2}K\rho {)}^{n+1}{\left(1+\frac{1}{\sigma }\right)}^{n+1}$$

$$2^{\mathcal{O}(n{d}^{n+1}(n\tau +nd\log (nd)+9n+d)\log a)}$$

$$x\mapsto {\left(\begin{array}{cc}1& x^T\end{array}\right)}^T/\sqrt{1+\Vert x{\Vert }^2}.$$

$$f^{\mathsf{h}}(\varphi (x))=f(x)/(1+\Vert x{\Vert }^2{)}^{d/2}.$$

$${\partial }_{\varphi (x)}{f}^h{\partial }_x\varphi =\frac{{\partial }_{x}f}{(1+\Vert x{\Vert }^2{)}^{d/2}}-\frac{d\cdot f(x)x^T}{(1+\Vert x{\Vert }^2{)}^{d/2+1}}$$

$$\Vert {\partial }_x\varphi \Vert =1/\sqrt{1+\Vert x{\Vert }^2}\text{ and }\Vert {\partial }_x{\varphi }^{-1}\Vert =1+\Vert x{\Vert }^2.$$

$$⟨f,g{⟩}_W:=\sum _{\alpha }{\left(\genfrac{}{}{0ex}{}{d}{\alpha }\right)}^{-1}{f}_{\alpha }{g}_{\alpha }$$

$$\hat{f}:x\mapsto f(x)/\left(\Vert f\Vert (1+\Vert x{\Vert }^2{)}^{(d-1)/2}\right)$$

$$\widehat{\partial f}:x\mapsto \partial f(x)/\left(d\Vert f\Vert (1+\Vert x{\Vert }^2{)}^{d/2-1}\right)$$

$$x\mapsto \tilde{\mathbf{f}}(x):=\mathbf{f}(x)/\left(\Vert \mathbf{f}\Vert (1+\Vert x{\Vert }^2{)}^{(d-1)/2}\right)$$

$$\mathbf{f}(x)/\left(\Vert \mathbf{f}\Vert (1+\Vert x{\Vert }^2{)}^{(d-1)/2}\right)=\sqrt{1+\Vert x{\Vert }^2}{\mathbf{f}}^{\mathsf{h}}(\varphi (x))/\Vert \mathbf{f}\Vert$$

$$\frac{{\mathbf{f}}^{\mathsf{h}}(\varphi (x))}{\Vert \mathbf{f}\Vert }\frac{x^T}{\sqrt{1+\Vert x{\Vert }^2}}+\sqrt{1+\Vert x{\Vert }^2}\frac{{\partial }_{\varphi (x)}\mathbf{f}}{\Vert \mathbf{f}\Vert }{\partial }_x\varphi .$$

$${\partial }_{\varphi (x)}\mathbf{f}{\partial }_x\varphi ={\partial }_{\varphi (x)}\left({\mathbf{f}}_{|{\mathbb{S}}^n}\right){\partial }_x\varphi$$

$$J\mapsto \hat{f}(m(J))+(1+\sqrt{d})\sqrt{n}w(J)\left[-1/2,1/2\right]$$

$$J\mapsto\widehat{\partial f}(m(J))+\big{(}1+\sqrt{d-1}\big{)}\sqrt{n}\,w(J)\left[-1/2,1/2\right]^{n}$$

$$h(x)=\frac{1}{\Vert f\Vert (1+\Vert x\Vert {)}^{(d-1)/2}}\text{ and }{h}^{\prime }(x)=\frac{1}{d\Vert f\Vert (1+\Vert x\Vert {)}^{d/2-1}}$$

$$\mathrm{dist}\left(\square[hf](J),(hf)(m(J))\right)\leq\big{(}1+\sqrt{d}\big{)}\sqrt{n}\,w(J)/2$$

$$\mathrm{dist}\left(\square[h^{\prime}\partial f](J),(h^{\prime}\partial f)(m(J))\right)\leq\big{(}1+\sqrt{d-1}\big{)}nw(J)/2.$$

$$C_{f}^{\prime}(J)\,:\,\left\{\begin{array}[]{rl}\left|\widehat{f}(m(J))\right|>&(1+\sqrt{d})\sqrt{n}w(J)\\ &\text{or }\left\|\widehat{\partial f}(m(J))\right\|>\sqrt{2}(1+\sqrt{d-1})nw(J)\end{array}\right.$$

$$\kappa (F,y):=\Vert F\Vert /\sqrt{F(y{)}^2+\Vert {\partial }_y{F}_{|{\mathrm{T}}_y{\mathbb{S}}^n}{\Vert }^2/d}.$$

$$\Vert f\Vert /{\kappa }_{\mathsf{aff}}(f,x)=\mathrm{dist}(f,{\Sigma }_x)$$

$$|F(y)|/\Vert F\Vert \ge 1/\left(\sqrt{2}\kappa (F,y)\right)$$

$$\Vert {\partial }_y{F}_{|{\mathrm{T}}_y{\mathbb{S}}^n}\Vert /\left(\sqrt{d}\Vert F\Vert \right)\ge 1/\left(\sqrt{2}\kappa (F,y)\right).$$

Full text

computational algebraic geometry, numerical methods, adaptive subdivision methods, isotopy of curves, complexity

††conference: ; ; ††ccs: Mathematics of computing Computations on polynomials††ccs: Mathematics of computing Interval arithmetic††ccs: Theory of computation Design and analysis of algorithms††ccs: Theory of computation Computational geometry

1. Introduction

In 2004 Plantinga and Vegter proposed an algorithm for computing a regularly isotopic piecewise linear approximation of a curve or surface (Plantinga and Vegter, 2004). Their algorithm relied on a subdivision method enhanced with interval arithmetic to certificate the procedure (i.e., ensure its correctness) and in Section 7 of their paper they provided some examples with their approximations and the record of how many cubes (squares in the case of plane curves) were in the description of these approximations. This number of cubes appears to be proportional to, and dominate, the cost of the computation. The paper however, contained no complexity analysis and not even a formal setting fixing either the kind of functions implicitly defining the considered curves and surfaces or the arithmetic used.

An article doing so was published in 2017 by Burr, Gao and Tsigaridas (Burr et al., 2017). The functions this article deals with are polynomials with integer coefficients and smooth zero set. Consistently with this choice of data, the arithmetic is infinite precision. The main result in the paper is a worst-case complexity analysis for the number of cubes in the description of the approximation which, as we just mentioned, dominates the cost of the computation. The bounds proved for this quantity are shown to be optimal. Yet, these bounds are exponential (both in the degree of the input polynomial and in its logarithmic height), a fact that motivates the following comment at the end of the paper

Even though our bounds are optimal, in practice, these are quite pessimistic […]

The authors further observe that, following from their Proposition 5.2 (see Theorem 6.3 below) an instance-based analysis of the algorithm (i.e., one yielding a cost that depends on the input at hand) could be derived from the evaluation of a certain integral. And they conclude their paper by writing

Since the complexity of the algorithm can be exponential in the inputs [size], the integral must be described in terms of additional geometric and intrinsic parameters.

A number of features in this state of affairs suggest that a condition-based approach to the analysis of our quantity of interest could be useful. To begin with, the fact that a condition number is a perfect fit for the notion of an “additional geometric and intrinsic parameter.” To which we may add the fact that the obvious set of ill-posed inputs, the set of polynomials having a non-smooth zero set, is precisely the set of data which are not allowed as inputs in (Burr et al., 2017). Of course, such a condition-based analysis would drop the assumption of integer coefficients and replace it by that of real coefficients but this is a common practice for numerical algorithms and, as we will see, it pays off in our case as it yields small (i.e., polynomial) average complexity bounds for a large class of probability measures.

Although our approach follows the condition-based ideas of, e.g., (Demmel, 1988; Renegar, 1995; Bürgisser and Cucker, 2013; Cucker et al., 2008, 2009), the complexity analysis in this paper would have been impossible without the continuous amortization technique developed in the exact numerical context (Burr et al., 2009; Burr, 2016). We hope that this merging of techniques will start a fruitful exchange of ideas between different approaches to continuous computation.

1.1. Notation

Throughout the paper, we will assume some familiarity with the basics of differential geometry and with the sphere $\mathbb{S}^{n}$ as a Riemannian manifold. For scalar smooth maps $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$, we will write the tangent map at $x\in\mathbb{R}^{m}$ as $\partial_{x}f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ when we want to emphasize it as a linear map and as $\partial f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$, $x\mapsto\partial f(x)$, when we want to emphasize it as a smooth function. For general smooth maps $F:\mathcal{M}\rightarrow\mathcal{N}$, we will just write $\partial_{x}F:\mathrm{T}_{x}\mathcal{M}\rightarrow\mathrm{T}_{x}\mathcal{N}$ as the tangent map.

In what follows, $\mathcal{P}_{n,d}$ will denote the set of real polynomials in $n$ variables with degree at most $d$, $\mathcal{H}_{n,d}$ the set of homogeneous real polynomials in $n+1$ variables of degree $d$, and $\|~{}\|$ and $\langle\,~{},~{}\rangle$ will denote the usual norm and inner product in $\mathbb{R}^{m}$ as well as the Weyl norm and inner product in $\mathcal{P}_{n,d}^{m}$ and $\mathcal{H}_{n,d}^{m}$. Given a polynomial $f\in\mathcal{P}_{n,d}$, $f^{\mathsf{h}}\in\mathcal{H}_{n,d}$ will be its homogenization and $\partial f$ the polynomial map given by its partial derivatives. For details about the concrete definition of each of these notions, see Section 4. Additionally, $V_{\mathbb{R}}(f)$ and $V_{\mathbb{C}}(f)$ will be, respectively, the real and complex zero sets of $f$.

We will denote by $\mathcal{I}_{n}$ the set of $n$-cubes of $\mathbb{R}^{n}$ and, for a given $J\in\mathcal{I}_{n}$, $m(J)$ will be its middle point, $w(J)$ its width, and $\mathrm{vol}\,J=w(J)^{n}$ its volume.

Also, $\mathbb{P}(A)$ will denote the probability of the event $A$, $\mathbb{E}_{x\in K}g(x)$ the expectation of $g(x)$ when $x$ is sampled uniformly from $K$ and $\mathbb{E}_{y}g(x)$ the expectation of $g(y)$ with respect to a previously specified probability distribution of $y$.

Regarding complexity parameters, $n$ will be the number of variables, $d$ the degree bound, and $N=\binom{n+d}{n}$ the dimension of $\mathcal{P}_{n,d}$.

Finally, $\ln$ will denote the natural logarithm and $\log$ the logarithm in base $2$.

1.2. Outline

In Section 2, we discuss the $\mathsf{PV}$  Algorithm and the $n$-dimensional generalization of its subdivision method that we will analyze. In Section 3, we state the main complexity results of this paper. In Section 4, we present the geometric framework of polynomials we will work with. Following a common practice in condition-based analysis we use homogenization to get many of our results. In Section 5, we introduce the condition number along with some of its main properties. In Section 6, we present the existing results of complexity of the subdivision method of the $\mathsf{PV}$  Algorithm based on local size bound functions from (Burr et al., 2017) and we relate them to the local condition number. In Section 7, we rely on the bounds for the condition number obtained in Section 5 to derive average and smoothed complexity bounds under (quite) general randomness assumptions.

2. The $\mathsf{PV}$ Algorithm

Given a real smooth hypersurface in $\mathbb{R}^{n}$ described implicitly by a map $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and a region $[-a,a]^{n}$, the $\mathsf{PV}$  Algorithm constructs a piecewise-linear approximation of the intersection of its zero set $V_{\mathbb{R}}(f)$ with $[-a,a]^{n}$ isotopic to this intersection inside $[-a,a]^{n}$.

Let $\mathcal{I}_{m}$ be the set of $m$-cubes of $\mathbb{R}^{m}$. Recall that an interval approximation of a function $F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m^{\prime}}$ is a map $\square[F]:\mathcal{I}_{m}\rightarrow\mathcal{I}_{m^{\prime}}$ such that for all $J\in\mathcal{I}_{m}$, $F(J)\subseteq\square[F](J)$ (c.f. (Ratschek and Rokne, 1984)). We notice that if we see $J$ as error bounds for the midpoint $m(J)$, then $\square[F](J)$ is nothing more than error bounds for $F(m(J))$.

Assume that we have interval approximations of both $f$ and its tangent map $\partial f$ or, more generally, of $hf$ and $h^{\prime}\partial f$ for some positive maps $h,h^{\prime}:\mathbb{R}^{n}\rightarrow(0,\infty)$. The $\mathsf{PV}$ Algorithm on $[-a,a]^{n}$ will subdivide this region into smaller and smaller cubes until the condition

[TABLE]

is satisfied in each of the $n$-cubes $J$ of the obtained subdivision of $[-a,a]^{n}$. In Section 4, we will be more precise on the assumptions on our interval approximations and the functions $h$ and $h^{\prime}$ that we will use.

The procedure in Algorithm 2.1 is only the subdivision routine of the $\mathsf{PV}$  Algorithm but it dominates its complexity in the sense that the remaining computations do not add to the final cost estimates in Landau notation. Moreover, these additional computations have been implemented only for $n\leq 3$. So, proceeding as in (Burr et al., 2017), we will only analyze the complexity of the subdivision routine, keeping track of the dependency on $n$. Also as in (Burr et al., 2017), our complexity analysis will not deal with the precision needed for the algorithm.

3. Main result

In this section, we outline without proofs the main results of this paper. In the first part, we describe our randomness assumptions for polynomials. In the second one, we give precise statements for our bounds on the average and smoothed complexity of the $\mathsf{PV}$ Algorithm.

3.1. Randomness Model

Most of the literature on random multivariate polynomials considers polynomials with Gaussian independent coefficients and relies on techniques that are only useful for Gaussian measures. We will instead consider a general family of measures relying on robust techniques coming from geometric functional analysis. Let us recall some basic definitions.

1. P1

A random variable $X$ is called centered if $\mathbb{E}X=0$. 2. P2

A random variable $X$ is called subgaussian if there exist a $K$ such that for all $p\geq 1$,

[TABLE]

The smallest such $K$ is called the $\Psi_{2}$-norm of $X$. 3. P3

A random variable $X$ satisfies the anti-concentration property with constant $\rho$ if

[TABLE]

The subgaussian property (P2) has other equivalent formulations. We refer the interested reader to (Vershynin, 2018).

Definition 3.1.

A dobro random polynomial $f\in\mathcal{H}_{n,d}$ with parameters $K$ and $\rho$ is a polynomial

[TABLE]

such that the $c_{\alpha}$ are independent centered subgaussian random variables with $\Psi_{2}$-norm $\leq K$ and anti-concentration property with constant $\rho$. A dobro random polynomial $f\in\mathcal{P}_{n,d}$ is a polynomial $f$ such that its homogenization $f^{\mathsf{h}}$ is so.

Some dobro random polynomials of interest are the following three.

1. N

A KSS random polynomial is a dobro random polynomial such that each $c_{\alpha}$ in (3.1) is Gaussian with unit variance. For this model we have $K\rho=1/\sqrt{2\pi}$. 2. U

A Weyl random polynomial is a dobro random polynomial such that each $c_{\alpha}$ in (3.1) have uniform distribution in $[-1,1]$. For this model we have $K\rho\leq 1$. 3. E

A $p$-random polynomial is a dobro random polynomial whose coefficients are independent identically distributed random variables with the density function $g(t)=c_{p}e^{-|t|^{p}}$ with $c_{p}$ being the appropriate constant and $p\geq 2$.

Remark 3.2*.*

When we are interested in integer polynomials, dobro random polynomials may seem inadequate. One may be inclined to consider random polynomials $f\in\mathcal{P}_{n,d}$ such that $c_{\alpha}$ is a random integer in the interval $[-2^{\tau},2^{\tau}]$, i.e., $c_{\alpha}$ is a random integer of bit-size at most $\tau$. As $\tau\to\infty$ and after we normalize the coefficients dividing by $2^{\tau}$, this random model converges to that of Weyl random polynomials.

To have a more satisfactory understanding of random integer polynomials, one has to consider random variables without a continuous density function. The techniques used in this note are already extended to include such random variables in the case of random matrices  (Rudelson and Vershynin, 2008; Vershynin, 2018). We hope to pursue this delicate case in a more general setting (including complete intersections) in future work.

3.2. Average and Smoothed Complexity

The following two theorems give bounds for, respectively, the average and smoothed complexity of Algorithm 2.1. In both of them $c_{1}$ and $c_{2}$ are, respectively, the universal constants in Theorems 7.2 and 7.4.

Theorem 3.1.

Let $f\in\mathcal{P}_{n,d}$, $\sigma>0$, and $g\in\mathcal{P}_{n,d}$ a dobro random polynomial with parameters $K$ and $\rho$. The expected number of $n$-cubes in the final subdivision of Algorithm 2.1 on input $(f,a)$ is at most

[TABLE]

if the interval approximations satisfy (4.4) and (4.5) and

[TABLE]

if they satisfy the hypothesis of (Burr et al., 2017).

Theorem 3.2.

Let $f\in\mathcal{P}_{n,d}$, $\sigma>0$, and $g\in\mathcal{P}_{n,d}$ a dobro random polynomial with parameters $K$ and $\rho$ . Then the expected number of $n$-cubes of the final subdivision of Algorithm 2.1 for input $(q_{\sigma},a)$ where $q_{\sigma}=f+\sigma\|f\|g$ is at most

[TABLE]

if the interval approximations satisfy (4.4) and (4.5) and

[TABLE]

if they satisfy the hypothesis of (Burr et al., 2017).

If we compare the results above with the worst-case bound of (Burr et al., 2017, Theorem 4.3), which is

[TABLE]

with $\tau$ being the largest bit-size of the coefficients of $f$, we can see that our probabilistic bounds are exponentially better: they may provide an explanation of the efficiency of the $\mathsf{PV}$  Algorithm in practice.

We note, however, that the bound in (Burr et al., 2017) and our bounds cannot be directly compared. Not only because the former is worst-case and the latter average-case (or smoothed) but because of the different underlying complexity settings: the bound in (Burr et al., 2017) applies to integer data, ours to real data. A first approach to bridge this difference relies on the approximation of distributions described in Remark 3.2. But, as mentioned there, this approach does not give completely satisfactory results. A more detailed study of how the $\mathsf{PV}$  Algorithm behaves on random integer polynomials is desirable.

4. Geometric framework

There is an extensive literature on norms of polynomials and their relation to norms of gradients in $\mathcal{H}_{n,d}$. We can use homogenization to carry these results from $\mathcal{H}_{n,d}$ to $\mathcal{P}_{n,d}$.

To be more precise, let $\phi:\mathbb{R}^{n}\rightarrow\mathbb{S}^{n}$, given by

[TABLE]

This gives a diffeomorphism between $\mathbb{R}^{n}$ and the upper half of $\mathbb{S}^{n}$, and we have

[TABLE]

Using the chain rule, we can see that

[TABLE]

where $\partial_{y}f^{\mathsf{h}}:\mathbb{R}^{n+1}\rightarrow\mathbb{R}$, $\partial_{x}\phi:\mathbb{R}^{n}\rightarrow\mathrm{T}_{x}\mathbb{S}^{n}=x^{\perp}$ and $\partial_{x}f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ are respectively the tangent maps of $f^{\mathsf{h}}$, $\phi$ and $f$.

As it will be useful later, we note that a direct computation shows

[TABLE]

An inner product on $\mathcal{H}_{n,d}$ with desirable geometric properties is known as the Weyl inner product and it is given by

[TABLE]

for $f=\sum_{\alpha}f_{\alpha}X^{\alpha},\,g=\sum_{\alpha}g_{\alpha}X^{\alpha}\in\mathcal{H}_{n,d}$. This product is extended to $\mathcal{P}_{n,d}$ using $\langle f,g\rangle:=\langle f^{\mathsf{h}},g^{\mathsf{h}}\rangle$ and to $\mathcal{P}_{n,d}^{k}$ using $\langle\mathbf{f},\mathbf{g}\rangle:=\sum_{i=1}^{k}\langle f_{i},g_{i}\rangle$.

We will use $\|~{}\|$ to denote both the Weyl norm for polynomials and the usual Euclidean norm in $\mathbb{R}^{n}$.

4.1. Lipschitz properties

Given $f\in\mathcal{P}_{n,d}$, let us consider the maps

[TABLE]

and

[TABLE]

which are just ”linearized” version of $f$ and its derivative. The intuition behind this fact is that for large values of $x$ a polynomial map of degree $d$ grows like $\|x\|^{d}$.

Proposition 4.1.

Let $\mathbf{f}\in\mathcal{P}_{n,d}^{k}$ be a polynomial map. Then the map

[TABLE]

is $(1+\sqrt{d})$-Lipschitz and, for all $x$, \big{\|}\tilde{\mathbf{f}}(x)\big{\|}\leq\sqrt{1+\|x\|^{2}}.

Proof.

For the Lipschitz property, it is enough to bound the norm of the derivative of the map by $1+\sqrt{d}$. Due to (4.1),

[TABLE]

and so the derivative equals

[TABLE]

Now, \big{\|}\mathbf{f}^{\mathsf{h}}(\phi(x))\big{\|}/\|\mathbf{f}\|\leq 1, by (Bürgisser and Cucker, 2013, Lemma 16.6), and

[TABLE]

with $\left\|\partial_{\phi(x)}\left(\mathbf{f}_{|\mathbb{S}^{n}}\right)\right\|\leq\sqrt{d}\|\mathbf{f}\|$, by the Exclusion Lemma (Bürgisser and Cucker, 2013, Lemma 19.22). Thus taking norms and using (4.3) finishes the proof.

The claim about \big{\|}\tilde{\mathbf{f}}(x)\big{\|} is just (Bürgisser and Cucker, 2013, Lemma 16.6). ∎

Corollary 4.2.

Let $f\in\mathcal{P}_{n,d}$. Then $\widehat{f}$ and $\widehat{\partial f}$ are Lipschitz with Lipschitz constants $(1+\sqrt{d})$ and $(1+\sqrt{d-1})$, respectively, and for all $x$, $|\widehat{f}(x)|,\|\widehat{\partial f}(x)\|\leq\sqrt{1+\|x\|^{2}}$.

Proof.

The claims about $\widehat{f}$ are immediate from Proposition 4.1. For the claims about $\widehat{\partial f}$, observe that $\partial f\in\mathcal{P}_{n,d-1}^{n}$ and that, by a direct computation, $\|\partial f\|\leq d\|f\|$. Thus Proposition 4.1 completes the proof. ∎

4.2. Interval approximations

The Lipschitz properties above (Corollary 4.2) ensure that the mapping

[TABLE]

is an interval approximation of $f\in\mathcal{P}_{n,d}$ and the mapping

[TABLE]

one of $\partial f$. Also, letting

[TABLE]

the mappings above become interval approximations $\square[hf]$ and $\square[h^{\prime}\partial f]$ of $hf$ and $h^{\prime}\partial f$, respectively, satisfying that, for each $J\in\mathcal{I}_{n}$,

[TABLE]

and

[TABLE]

With these conditions on our interval approximations, we can reformulate a weaker, but easier to check, condition $C_{f}^{\prime}(J)$.

Theorem 4.3.

Let $J\in\mathcal{I}_{n}$ and assume that our interval approximation satisfies (4.4) and (4.5). If the condition

[TABLE]

is satisfied, then $C_{f}(J)$ is true.

Lemma 4.4.

Let $x\in\mathbb{R}^{n}$. Then for all $v,w\in B_{\|x\|/\sqrt{2}}(x)$, $\langle v,w\rangle>0$.

Proof.

Let $W:=\{u\in\mathbb{R}^{n}\mid\langle x,u\rangle\geq\|x\|\|u\|/\sqrt{2}\}$ be the convex cone of those vectors $u$ whose angle $\widehat{x,u}$ with $x$ is at most $\pi/4$. For $v,w\in W$, $\widehat{v,w}\leq\widehat{v,x}+\widehat{x,w}\leq\pi/2$, by the triangle inequality. Thus $\cos\,\widehat{v,w}\geq 0$.

Now, $\mathrm{dist}(x,\partial W)=\min\{\|x-u\|\mid\langle x,u\rangle=\|x\|\|u\|/\sqrt{2}\}$ where the latter equals the distance of $x$ to a line having an angle $\pi/4$ with $x$, which is $\|x\|/\sqrt{2}$. Hence $B_{\|x\|/\sqrt{2}}(x)\subseteq W$ and we are done. ∎

Proof of Theorem 4.3.

When the condition on $\widehat{f}(m(J))$ is satisfied,  (4.4) guarantees that $0\notin\square[hf](J)$. Whenever the condition on $\widehat{\partial f}(m(J))$ is satisfied, (4.5) and Lemma 4.4 guarantee that $0\notin\langle\square[h^{\prime}\partial f](J),\square[h^{\prime}\partial f](J)\rangle$. Hence $C_{f}^{\prime}(J)$ implies $C_{f}(J)$ under the given assumptions. ∎

Remark 4.5*.*

The interval approximations in (Burr et al., 2017) are based on Taylor expansion at the midpoint, so they are different from ours. However, our complexity analysis also applies to the interval approximations considered in (Burr et al., 2017), see §6.2 below for the details.

5. Condition number

As other numerical algorithms in computational geometry, the $\mathsf{PV}$  Algorithm has a cost which significantly varies with inputs of the same size. One wants to explain this variation in terms of geometric properties of the input. Condition numbers allow for such an analysis.

Definition 5.1.

(Cucker et al., 2018; Bürgisser et al., 2018) Given $F\in\mathcal{H}_{n,d}$, the local condition number of $F$ at $y\in\mathbb{S}^{n}$ is

[TABLE]

Given $f\in\mathcal{P}_{n,d}$, the local affine condition number of $f$ at $x\in\mathbb{R}^{n}$ is $\kappa_{\sf aff}(f,x):=\kappa(f^{\mathsf{h}},\phi(x))$.

5.1. What does $\kappa_{\sf aff}$ measure?

The nearer the hypersurface $V_{\mathbb{R}}(f)$ is to having a singularity at $x\in\mathbb{R}^{n}$, the smaller are the boxes drawn by the $\mathsf{PV}$ Algorithm around $x$. A quantity controlling how close if $f$ to have a singularity at $x$ will therefore control the size of these boxes. This is precisely what $\kappa_{\sf aff}(f,x)$ does.

Theorem 5.2 (Condition Number Theorem).

Let $x\in\mathbb{R}^{n}$ and $\Sigma_{x}$ be the set of hypersurfaces having $x$ as a singular point. That is, $\Sigma_{x}:=\{g\in\mathcal{P}_{n,d}\mid g(x)=0,\,\partial_{x}g=0\}$. Then for every $f\in\mathcal{P}_{n,d}$,

[TABLE]

where the distance is induced by the Weyl norm of $\mathcal{P}_{n,d}$.

Proof.

This follows from (Bürgisser and Cucker, 2013, Proposition 19.6), (Bürgisser et al., 2018, Theorem 4.4) and the definition of $\kappa_{\sf aff}$. ∎

Theorem 5.2 provides a geometric interpretation of the local condition number, and a corresponding ”intrinsic” complexity parameter as desired by the authors of (Burr et al., 2017, 2018). The next result will be useful in the probabilistic analyses.

Corollary 5.3.

Let $x\in\mathbb{R}^{n}$ and let $P_{x}:\mathcal{P}_{n,d}\rightarrow\Sigma_{x}^{\perp}$ be the orthogonal projection onto the orthogonal complement of the linear subspace $\Sigma_{x}$. Then $\kappa_{\sf aff}(f,x)=\|f\|/\|P_{x}f\|$.

Proof.

We have that $\mathrm{dist}(f,\Sigma_{x})=\|P_{x}f\|$ since $\Sigma_{x}$ is a linear subspace. Hence Theorem 5.2 finishes the proof. ∎

We notice that the above expression should not come as a surprise, since $\kappa_{\sf aff}$ is define in a way that the denominator is the norm of a vector depending linearly of $f$.

5.2. A fundamental proposition

The following result plays a fundamental role in our development. Despite having been used in many occasions within various proofs, it wasn’t explicitly stated until recently.

Lemma 5.4.

Let $F\in\mathcal{H}_{n,d}$ and $y\in\mathbb{S}^{n}$. Then either

[TABLE]

or

[TABLE]

Proof.

This is (Bürgisser et al., 2018, Proposition 3.6). ∎

Proposition 5.5.

Let $f\in\mathcal{P}_{n,d}$ and $x\in\mathbb{R}^{n}$. Then either

[TABLE]

Proof.

Without loss of generality assume that $\|f\|=1$. Let $y:=\phi(x)$, $F:=f^{\mathsf{h}}$ and assume that the first inequality does not hold. Then, by (4.1), $|F(y)|\leq 1/\left(2\sqrt{2d}\,\kappa(F,y)\sqrt{1+\|x\|^{2}}\right)$.

By (4.2),  (4.3) and Lemma 5.4, we get

[TABLE]

We divide by $\sqrt{d}$ and use the triangle inequality to obtain

[TABLE]

Using (4.1) and our initial assumption on the second term in the sum, which we subtract, we get the desired inequality since $\|x\|<\sqrt{1+\|x\|^{2}}$. ∎

6. Adaptive Complexity Analysis

As stated in Section 2 (and as done in (Burr et al., 2017)), our complexity analysis will focus on the number of subdivisions steps of the subdivision routine of the $\mathsf{PV}$ Algorithm (Algorithm 2.1). That is, our cost measure will be the number of $n$-cubes in the final subdivision of $[-a,a]^{n}$.

We note that we do not deal with the precision needed to run the algorithm. This issue should be treated in the future.

6.1. Local size bound framework (Burr

et al., 2017)

The original analysis in (Burr et al., 2017) was based on the notion of local size bound.

Definition 6.1.

A local size bound for $f$ is a function $b_{f}:\mathbb{R}^{n}\rightarrow[0,\infty)$ such that for all $x\in\mathbb{R}^{n}$,

[TABLE]

Arguing as in (Burr et al., 2017, Proposition 4.1), one can easily get the following general bound.

Proposition 6.2.

(Burr et al., 2017)* The number of $n$-cubes of the final subdivision of Algorithm 2.1 on input $(f,a)$, regardless of how the subdivision step is done, is at most*

[TABLE]

The bound above is worst-case, it considers the worst $b_{f}(x)$ among the $x\in[-a,a]^{n}$. Continuous amortization developed by Burr, Krahmer and Yap (Burr et al., 2009; Burr, 2016), provides the following refined complexity estimate (Burr et al., 2017, Proposition 5.2) which is adaptative.

Theorem 6.3.

(Burr et al., 2009; Burr, 2016; Burr et al., 2017)* The number of $n$-cubes of the final subdivision of Algorithm 2.1 on input $(f,a)$ is at most*

[TABLE]

Moreover, the bound is finite if and only if the algorithm terminates. ∎

To effectively use Theorem 6.3 we need to explicit estimates for the local size bound.

6.2. Construction in (Burr

et al., 2017)

In (Burr et al., 2017), the authors use the following function

[TABLE]

where $g_{f}$ is the polynomial $\langle\partial f(X),\partial f(Y)\rangle$, to construct a local size bound.

Theorem 6.4.

(Burr et al., 2017)* Assume that the interval approximation satisfies the hypothesis of (Burr et al., 2017). Then*

[TABLE]

is a local size bound function for $f$. ∎

Looking at the definition of $\mathcal{C}(f,x)$ in (Burr et al., 2017) one can see that $1/\mathcal{C}$ measures how near is $x$ of being a singular zero of $f$. This is similar to $1/\kappa_{\sf aff}$ which, by Theorem 5.2, measures how near is $f$ of having a singular zero at $x$. The following result relates these two quantities.

Theorem 6.5.

Let $d>1$ and $f\in\mathcal{P}_{n,d}$. Then, for all $x\in\mathbb{R}^{n}$,

[TABLE]

Proof.

Note that Corollary 4.2 holds over the complex numbers as well. Due to this and the fact that $V_{\mathbb{C}}(f)=V_{\mathbb{C}}(\widehat{f})$, we have that

[TABLE]

Now, if $\sqrt{2}(1+\sqrt{d-1})\,\mathrm{dist}((y_{1},y_{2}),(x,x))<\|\widehat{\partial f}(x)\|$, then $\sqrt{2}(1+\sqrt{d-1})\|y_{i}-x\|<\|\widehat{\partial f}(x)\|$. Thus, by Corollary 4.2, $\sqrt{2}\|\widehat{\partial f}(y_{i})-\widehat{\partial f}(x)\|<\|\widehat{\partial f}(x)\|$ and so, by Lemma 4.4, $0\neq\langle\widehat{\partial f}(y_{1}),\widehat{\partial f}(y_{2})\rangle$. Hence

[TABLE]

The bound now follows from Proposition 5.5, together with $2^{3(n-1)}d+\sqrt{n}\leq 2^{3n-2}d$ and

[TABLE]

for which we use that

[TABLE]

∎

The main difference between $C(f,x)$ and $\kappa(f,x)$ is that $C(f,x)$ is a non-linear quantity and is hard to compute, while the local condition number $\kappa(f,x)$—as indicated in Corollary 5.3—is a linear quantity and is rather easy to compute. Theorem 6.6 below and the complexity analysis in Section 7 show that the local condition number $\kappa_{\sf aff}(f,x)$ is easily amenable to the adaptive complexity analysis techniques developed by Burr, Krahmer and Yap (Burr et al., 2009; Burr, 2016).

6.3. Condition-based complexity

The following result expresses a local size bound in terms of the local condition number $\kappa_{\sf aff}(f,x)$ directly, without using the construction in (Burr et al., 2017).

Theorem 6.6.

Assume that the interval approximation satisfies (4.4) and (4.5). Then

[TABLE]

is a local size bound for $f$.

Proof.

Let $x\in\mathbb{R}^{n}$. As, by Theorem 4.3, $C_{f}^{\prime}(J)$ implies $C_{f}(J)$, it is enough to compute the minimum volume of $J\in\mathcal{I}_{n}$ containing $x$ such that $C_{f}^{\prime}(J)$ is false. This will still give a local size function for $f$.

Since $x\in J$, $\|x-m(J)\|\leq\sqrt{n}w(J)/2$. Hence, by Corollary 4.2 and Proposition 5.5, either

[TABLE]

or

[TABLE]

This means that $C^{\prime}_{f}(J)$ is true if either

[TABLE]

or

[TABLE]

Hence we get that $C^{\prime}_{f}(J)$ is true when both conditions are satisfied and the inequality $1+\sqrt{d}\leq 2\sqrt{d}$ finishes the proof. ∎

Using the results above, we get the following theorem exhibiting a condition-based complexity analysis of Algorithm 2.1.

Theorem 6.7.

The number of $n$-cubes in the final subdivision of Algorithm 2.1 on input $(f,a)$ is at most

[TABLE]

if the interval approximation satisfies (4.4) and (4.4), and at most

[TABLE]

if the interval approximation satisfies the hypothesis of (Burr et al., 2017).

Proof.

This is just Theorems 6.3, 6.6 and 6.5 combined with the fact that the integral $\int_{[-a,a]^{n}}\,\kappa_{\sf aff}(f,x)^{n}\,\mathrm{d}x$ is nothing more than $(2a)^{n}\,\mathbb{E}_{x\in[-a,a]^{n}}\left(\kappa_{\sf aff}(f,x)^{n}\right)$. ∎

We observe that in contrast with the complexity analyses (of condition numbers closely related to $\kappa_{\sf aff}$) in the literature (see, e.g., (Cucker et al., 2008, 2009, 2012; Cucker et al., 2018; Bürgisser et al., 2018, 2018)), the bounds in Theorem 6.7 depend on $\mathbb{E}_{x\in[-a,a]^{n}}\left(\kappa_{\sf aff}(f,x)^{n}\right)$ and not on $\max_{x\in[-a,a]}\kappa_{\sf aff}(f,x)^{n}$. Whereas the former has finite expectation (over $f$), the latter has not. This shows that condition-based analysis combined with adaptive complexity techniques such as continuous amortization may lead to substantial improvements.

7. Probabilistic analyses

In this section, we prove Theorems 3.1 and 3.2 stated in Section 3.

7.1. Average Complexity Analysis

The following theorem is the main technical result from which the average complexity bound will follow.

Theorem 7.1.

Let $f\in\mathcal{P}_{n,d}$ be a dobro random polynomial with parameters $K$ and $\rho$. For all $x\in\mathbb{R}^{n}$ and $t\geq e^{n}$,

[TABLE]

where $c_{1}$ and $c_{2}$ are, respectively, the universal constants of Theorems 7.2 and 7.4.

The proof of Theorem 7.1 relies on two basic results from geometric functional analysis.

Theorem 7.2.

(Vershynin, 2018, Theorems 2.6.3 and 3.1.1)* There is a universal constant $c_{1}\geq 1$ with the following property. For all random vectors $X\>=(X_{1},\ldots,X_{N})^{T}$ with each $X_{i}$ centered and sub-Gaussian with $\psi_{2}$-norm $K$, and for all $t\geq c_{1}K\sqrt{N}$ the following inequality is satisfied*

[TABLE]

Definition 7.3.

The concentration function of a random vector $X\in\mathbb{R}^{k}$ is the function $\mathcal{L}_{X}(\varepsilon):=\max_{u\in\mathbb{R}^{k}}\mathbb{P}\left(\lVert X-u\rVert\leq\varepsilon\right).$

Theorem 7.4.

(Rudelson and Vershynin, 2015, Corollary 1.4)* There is a universal constant $c_{2}\geq 1$ with the following property. For every random vector $X=(X_{1},\ldots,X_{N})^{T}$ with independent random variables $X_{i}$, and every $k$-dimensional linear subspace $S$ of $\mathbb{R}^{N}$ we have*

[TABLE]

where $P_{k}$ is the orthogonal projection onto $S$. ∎

Remark 7.5*.*

In (Livshyts et al., 2016) and references therein one can find information about the optimal value of the absolute constant $c_{2}$ in Theorem 7.4.

Remark 7.6*.*

We notice that for a dobro random polynomial $f$ with parameters, $K$ and $\rho$ we have $K\rho\geq\frac{1}{4}$ (Ergür et al., 2018a, (1)). Actually, the product $K\rho$ is invariant under scaling in the following sense; for $t>0$, $tf$ is again a dobro polynomial with parameters $tK$ and $\rho/t$. Hence, for the sake of simplicity and without loss of generality, we will assume $c_{1}c_{2}K\rho\geq 1$. Moreover, since $\kappa_{\sf aff}$ is scale invariant, we can assume, again without loss of generality, that $c_{1}K\geq 1$.

Proof of Theorem 7.1.

$\kappa_{\sf aff}(f,x)=\|f\|/\|P_{x}f\|$ by Corollary 5.3. So, by an union bound, for all $u,s>0$,

[TABLE]

By Theorems 7.2 and 7.4, we have

[TABLE]

We set $u=c_{1}K\sqrt{N\ln(s)}$ and use $s^{-N}\leq s^{-(n+1)}$, so we get

[TABLE]

By substituting $s=t^{\frac{1}{n}}$ we are done. ∎

Combining Theorem 6.7 with the next theorem, we get the proof of Theorem 3.1.

Theorem 7.7.

Let $f\in\mathcal{P}_{n,d}$ be a dobro random polynomial with parameters $K$ and $\rho$. Then

[TABLE]

where $c_{1}$ and $c_{2}$ are the universal constants of Theorems 7.2 and 7.4.

Proof.

By the Fubini-Tonelli theorem,

[TABLE]

so it is enough to have a uniform bound for

[TABLE]

Now, by Theorem 7.1, this is bounded by

[TABLE]

After the change of variables $t=e^{ns}$ the integral becomes

[TABLE]

where $\Gamma$ is Euler’s Gamma function. Using the Stirling estimates for it, we obtain

[TABLE]

and $N\leq(2d)^{n}$. Combining all these inequalities, we obtain the desired upper bound. ∎

7.2. Smoothed Complexity Analysis

The tools used for our average complexity analysis yield also a smoothed complexity analysis (see (Spielman and Teng, 2002) or (Bürgisser and Cucker, 2013, §2.2.7)). We provide this analysis following the lines of (Ergür et al., 2018b),

The main idea of smoothed complexity is to have a complexity measure interpolating between worst-case complexity and average-case complexity. More precisely, we are interested in the maximum—over $f\in\mathcal{P}_{n,d}$—of the average cost of Algorithm 2.1 with input

[TABLE]

where $g\in\mathcal{P}_{n,d}$ is a dobro random polynomials with parameters $K$ and $\rho$ and $\sigma\in(0,\infty)$. Notice that the perturbation $\sigma\|f\|g$ of $f$ is proportional to both $\sigma$ and $\|f\|$.

Lemma 7.8.

Let $q_{\sigma}$ be as in (7.2). Then for $t>1+\sigma\sqrt{N}$

[TABLE]

and for every $x\in\mathbb{R}^{n}$,

[TABLE]

where $P_{x}$ is as in Corollary 5.3. ∎

Proof.

By the triangle inequality we have $\mathbb{P}(\lVert q_{\sigma}\rVert\geq t\lVert f\rVert)\leq\mathbb{P}(\lVert g\rVert\geq(t-1)/\sigma)$. Then we apply Theorem 7.2 which finishes the proof of the first claim. The second claim is a direct consequence of Theorem 7.4. ∎

As in the average case, this leads to a tail bound.

Theorem 7.9.

Let $q_{\sigma}$ be as in (7.2). Then for $\sigma>0$ and $t\geq e^{n}$, $\mathbb{P}\left(\kappa_{\sf aff}(q_{\sigma},x)^{n}\geq t\right)$ is bounded by

[TABLE]

where $c_{1}$ and $c_{2}$ are, respectively, the universal constants of Theorems 7.2 and 7.4.

Proof.

We proceed as in the proof of Theorem 7.1, but with Lemma 7.8 using $u=\|f\|(\sigma c_{1}K\sqrt{N\ln(t)}+1)$. This gives the desired bound arguing as in that proof after noticing that

[TABLE]

which holds since $c_{1}K\sqrt{N\ln(t)}\geq 1$. ∎

Finally, the following theorem, together with Theorem 6.7, gives the proof of Theorem 3.2.

Theorem 7.10.

Let $q_{\sigma}$ be as in (7.2). Then for all $\sigma>0$, $\mathbb{E}_{q_{\sigma}}\mathbb{E}_{x\in[-a,a]^{n}}$ is bounded by

[TABLE]

where $c_{1}$ and $c_{2}$ are the universal constants of Theorems 7.2 and 7.4.

Proof.

The proof is as that of Theorem 7.7, but using Theorem 7.9 instead of Theorem 7.1. Actually, the integrand one ends up with is the same, up to a constant, so the calculation of the integral is the same up to that constant. ∎