# Explicit and fast error bounds for polynomial approximation

## Main Question

This question is about finding explicit, calculable, and fast error bounds when approximating continuous functions with polynomials to a user-specified error tolerance.

In this question:

• A polynomial $$P(x)$$ is written in Bernstein form of degree $$n$$ if it is written as $$P(x)=\sum_{k=0}^n a_k {n \choose k} x^k (1-x)^{n-k},$$ where $$a_0, ..., a_n$$ are the polynomial's Bernstein coefficients.
• For the Bernstein polynomial of $$f(x)$$ of degree $$n$$, $$a_k = f(k/n)$$.

Let $$f:[0,1]\to [0,1]$$ be continuous and polynomially bounded (both $$f$$ and $$1-f$$ are bounded below by min($$x^n$$, $$n(1-x)^n$$) for some integer $$n$$), let $$r\ge 3$$, and denote the Bernstein polynomial of degree $$n$$ of a function $$g$$ as $$B_n(g)$$.

Given that $$f$$ has a continuous $$r$$-th derivative (or has a Lipschitz continuous $$(r-1)$$-th derivative), are there results that give a sequence of polynomials $$P_n$$ with the following error bound?

$$| f(x) - P_n(f)(x) | \le \epsilon(f, n, x) = O(1/n^{r/2}),$$ where:

1. $$\epsilon(f, n, x)$$ is a fully determined function, with all constants in the expression having a known exact value or upper bound.
2. $$P_n(f)(x)$$ is an approximating polynomial of degree $$n$$ that can be readily rewritten to a polynomial in Bernstein form with coefficients in $$[0, 1]$$. Preferably, $$P_n(f)$$ has the form $$B_n(W_n(f))$$ where $$W_n(f)$$ is easily computable from $$f$$ using rational arithmetic only (see "Remarks", later).

One way to answer this (see this question) is to find a sequence $$W_n(f)$$ and an explicit and tight upper bound on $$C_1>0$$ such that, for each integer $$n\ge 1$$ that's a power of 2— $$\max_{0\le k\le 2n}\left|\left(\sum_{i=0}^k \left(W_n\left(\frac{i}{n}\right)\right) {n\choose i}{n\choose {k-i}}/{2n \choose k}\right)-W_{2n}\left(\frac{k}{2n}\right)\right|\le \frac{C_1 M}{n^{r/2}},$$ where $$M$$ is the maximum absolute value of $$f$$ and its derivatives up to the $$r$$-th derivative (or, respectively, the maximum of $$|f|$$ and the Lipschitz constants of $$f$$ and its derivatives up to the $$(r-1)$$-th derivative).

Then $$| f(x) - B_n(W_n(f))(x) | \le \frac{C_1}{1-\sqrt{2/2^{r+1}}}\frac{M}{n^{r/2}}=O(1/n^{r/2})$$ (see Lemma 3 in "Proofs for Polynomial-Building Schemes), although this is only guaranteed to work for power-of-2 values of $$n$$. For example, $$W_n$$ can be $$2f-B_n(f)$$ and $$r$$ can be 3 or 4, or $$W_n$$ can be $$B_n(B_n(f))+3(f-B_n(f))$$ and $$r$$ can be 5 or 6.

### Motivation

My motivation for this question is to implement "approximate Bernoulli factories", or algorithms that toss heads with a probability equal to a polynomial in Bernstein form that comes within $$\epsilon$$ of a continuous function $$f(x)$$. This involves finding a reasonably small degree $$n$$ of that polynomial, then the algorithm works as follows:

1. Flip the coin $$n$$ times, count the number of heads as $$h$$.
2. With probability equal to the $$h$$-th Bernstein coefficient, return heads; otherwise tails.

Note that the algorithm requires finding only one Bernstein coefficient per run. And for ordinary Bernstein polynomials, finding it is trivial — $$f(h/n)$$ — but the degree $$n$$ can be inordinate due to Bernstein polynomials' slow convergence; for example, if $$\epsilon=0.01$$ and $$f$$ is Lipschitz with constant 1, the required polynomial degree is 11879.

• Approximating $$f$$ with a rational function is also interesting, but is outside the scope of this question.
• Exact Bernoulli factories require a slightly different approach to finding the polynomials; see another question of mine.

### Polynomials with faster convergence than Bernstein polynomials

As is known since Voronovskaya (1932), the Bernstein polynomials converge uniformly to $$f$$, in general, at a rate no faster than $$O(1/n)$$, regardless of $$f$$'s smoothness, which means that it won't converge in a finite expected running time. (See also a related question by Luis Mendo on ordinary Bernstein polynomials.)

But Lorentz (1966, "The degree of approximation by polynomials with positive coefficients") has shown that if $$f(x)$$ is positive (the case that interests me) and has $$k$$ continuous derivatives, there are polynomials with non-negative Bernstein coefficients that converge to $$f$$ at the rate $$O(1/n^{k/2})$$ (and thus can be faster than the $$O(1/n^{2+\epsilon})$$ needed for a finite expected running time, depending on $$f$$).*

Thus, people have developed alternatives, including iterated Bernstein polynomials, to improve the convergence rate. These include:

• Micchelli, C. (1973). The saturation class and iterates of the Bernstein polynomials. Journal of Approximation Theory, 8(1), 1-18.
• Guan, Zhong. "Iterated Bernstein polynomial approximations." arXiv preprint arXiv:0909.0684 (2009).
• Güntürk, C. Sinan, and Weilin Li. "Approximation with one-bit polynomials in Bernstein form" arXiv preprint arXiv:2112.09183 (2021).
• The "Lorentz operator": Holtz, Olga, Fedor Nazarov, and Yuval Peres. "New coins from old, smoothly" Constructive Approximation 33, no. 3 (2011): 331-363.
• Draganov, Borislav R. "On simultaneous approximation by iterated Boolean sums of Bernstein operators." Results in Mathematics 66, no. 1 (2014): 21-41.

Usually, papers like those express a bound on the error when approximating a function with polynomials as follows: $$| f(x) - P_n(f)(x) | \le c_n \epsilon(f, n, x),$$ where $$\epsilon(f, n, x)$$ is a fully determined function, $$c_n>0$$ is a constant that may depend on $$n$$, and $$P_n(f)(x)$$ is an approximating polynomial of degree $$n$$.

There are results where the error bound $$\epsilon(.)$$ is in $$O(1/n^{k/2})$$, but in all those results I've seen so far (e.g., Theorem 4.4 in Micchelli; Theorem 5 in Güntürk and Li), $$c_n$$ is unknown, and no upper bound for $$c_n$$ is given by the results in the papers above, so that the error bound is unimplementable and there is no way of knowing beforehand whether $$P_n$$ will come close to $$f$$ within a user-specified error tolerance. (There is also a separate matter of rewriting the polynomial to its Bernstein form, but this is much more manageable, especially with the Lorentz operator.)

### Remarks

• I prefer approaches that involve only rational arithmetic and don't require transcendental or trigonometric functions to build the Bernstein-form polynomials.

• Unlike with rational arithmetic (where arbitrary precision is trivial thanks to Python's fractions module), transcendental and trig. functions require special measures to support arbitrary accuracy, such as constructive/recursive reals — floating point won't do for my purposes.
• In general, "rounding" a polynomial's coefficients or "nodes" to rational numbers will add a non-trivial error that, for my purposes, has to be accounted for in any error bound.

* If the polynomials are not restricted in their coefficients, then the rate $$O(1/n^k)$$ is possible (e.g., DeVore and Lorentz 1993). But this result is not useful to me since my use case (approximate Bernoulli factories) requires the polynomials to have Bernstein coefficients in $$[0, 1]$$.

If $$f, f', \dotsc, f^{(\nu-1)}$$ are all absolutely continuous on $$[-1, 1]$$, and $$f^{(\nu)}$$ has bounded variation $$V$$ on $$[-1, 1]$$, where $$\nu \ge 1$$, then the polynomial interpolant $$p_n$$ of degree $$n > \nu$$ through the Chebyshev points of the second kind, $$p_n(x_j) = f(x_j), \quad x_j = \cos \tfrac{j\pi}{n}, \quad j=0,\dotsc,n,$$ satisfies the bound $$\sup_{-1 \le x \le 1} | f(x) - p_n(x)| \le \frac{4V}{\pi \nu (n - \nu)^\nu}.$$ See Theorem 7.2, equation (7.5) of [1].

It should be easy to map the result to [0, 1]. If all you know is that $$f$$ has $$k$$ continuous derivatives, then you can take $$\nu = k-1$$, provided $$k \ge 2$$. The convergence rate $$O(n^{-(k-1)})$$ is faster than you asked for when $$k > 2$$ (and equal at $$k=2$$).

I am unclear exactly what you mean by being "readily rewritten" to the Bernstein polynomial basis. One explicit possibility: you can use a discrete Fourier transform (specifically, the FFT) to recover the coefficients of the Chebyshev series of $$p_n$$ from the values at the Chebyshev points above (again, see [1]), then use an explicit expansion of the Chebyshev polynomials in the Bernstein basis (e.g., the one on the Wikipedia page for the Bernstien polynomials).

[1] Trefethen, Lloyd N., Approximation theory and approximation practice, Other Titles in Applied Mathematics 128. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM) (ISBN 978-1-611972-39-9/pbk). 305 p. (2013). ZBL1264.41001.

• I thought of Chebyshev approximation then conversion to the Bernstein basis (e.g., using the paper of Rababah 2003). But what I had in mind when I spoke of polynomials that are "readily convertible" to Bernstein form includes polynomials formed with the Lorentz operator (Holtz et al. 2011), which doesn't rely on trigonometric functions, but only on rational arithmetic. Other examples include the iterated Bernstein operator in the Guan paper, or the one-bit polynomials in Güntürk and Li (which are already in Bernstein form). Jun 13, 2022 at 23:47
• For the Lorentz operator, see also my Python implementation. Jun 13, 2022 at 23:53
• You could always choose rational nodes close to the Chebyshev nodes, though if you want an explicit bound you'd need then to bound the resulting perturbation. Other nodes, like the Gauss-Legendre quadrature nodes, would also give polynomial interpolants with good approximation properties. The key property of the nodes is the asymptotic density is proportional to $1/\sqrt{1-x^2}$. I don't know if there any explicitly given rational interpolation nodes with this asymptotic density. Jun 14, 2022 at 21:07
• Maybe you should clarify your requirement (3) that you want the coefficients to lie in [0, 1]? Jun 14, 2022 at 21:07
• On your first comment: At this point I would rather see if there are upper bounds for the constant factors given in the papers I listed in my question (e.g., Theorem 4.4 in Micchelli; Theorem 5 in Güntürk and Li; numerous lemmas in Holtz et al.). My implementation attempts have shown that Chebyshev polynomials are far from being "readily convertible" to Bernstein form (first vs. second kind; $[0, 1]$ vs. $[-1, 1]$; Chebyshev points vs. Chebyshev coefficients; etc.) Thus, I will consider Chebyshev polynomials as a last resort. Jun 14, 2022 at 21:21

After analyzing the proof of Güntürk and Li (2021), Theorem 3.3, I believe I found explicit error bounds for the Micchelli–Felbecker polynomials (iterated Bernstein polynomials) when $$f(x)$$ has a given number of continuous derivatives. In the table below, let—

• $$||f(x)||_{C^k} = \max(\max_{0\le x \le 1} |f(x)|, \max_{0\le x \le 1} |f^{(k)}(x)|),$$
• $$I$$ be the identity operator, and
• $$B_n$$ be the Bernstein operator of degree $$n$$.
No. of continuous derivatives Polynomial Error bound
3 $$I-(I-B_n)^2$$ 0.3489 $$||f(x)||_{C^3}/n^{3/2}$$
4 $$I-(I-B_n)^2$$ 0.275 $$||f(x)||_{C^4}/n^2$$
5 $$I-(I-B_n)^3$$ 0.7284 $$||f(x)||_{C^5}/n^{5/2}$$
6 $$I-(I-B_n)^3$$ 0.9961 $$||f(x)||_{C^6}/n^3$$

Providing the full proof for these error bounds is a bit tedious, so here is a sketch. The proof involves finding upper bounds for binomial moments (discussed later), then plugging them in to estimates for the Bernstein polynomial approximation error (denoted as $$(B_n-I)(f)$$, $$G_{n,r+1}$$, and $$(B_n-I)^{\lceil (r+1)/2 \rceil}(f)$$ in the proof of Theorem 3.3) as well as derivatives for a function denoted as $$F_{n,\alpha}$$, along with the bound, mentioned in the proof of Theorem 3.3, that $$||(B_n-I)^k|| \le 2^k$$ for every $$k\ge 1$$. See later for Python code that calculates these error bounds.

I would appreciate any corrections.

EDIT (Aug. 6): It appears that the exact definition of the norm $$||f(x)||_{C^k}$$ matters a great deal. The paper Güntürk and Li (2021) defined this norm (in the univariate case) as— $$\max(||f(x)||_\infty, ||f^{(k)}(x)||_\infty),$$ where $$||f(x)||_\infty$$ is the essential supremum. Thus, this definition looks only at the function and its $$k$$-th derivative, and not any derivatives in between. However, with this definition, the results above fail in the case of the polynomial $$2x(1-x)$$, which is a quadratic polynomial whose third and higher derivatives vanish, as well as $$(\sin(x)+2x(1-x))/2$$, which is a convex combination of a concave function and a quadratic polynomial.

There are other definitions for the norm $$||f(x)||_{C^k}$$ that may work better, including:

• $$\sum_{0\le i\le k} ||f^{(i)}(x)||_\infty,$$ (see the Encyclopedia of Math and these lecture notes).
• $$\max_{0\le i\le k} ||f^{(i)}(x)||_\infty,$$ (Chichilnisky, G.,1986, Topological complexity of manifolds of preferences; Petersen, Philipp Christian. "Neural network theory." University of Vienna (2020)).

I don't know whether these results will hold, even for quadratic polynomials, with either definition of the $$C^k$$ norm, though.

END EDIT

Remark 3.4 in Güntürk and Li mentions that Theorem 3.3 works for functions with Lipschitz continuous 2nd, 3rd, 4th, or 5th derivative rather than continuous 3rd, 4th, 5th, or 6th derivative, respectively, after replacing the $$C^k$$ norm with the Lipschitz $$C^{k-1}$$ norm and making other "natural modifications". Assuming that the bounds above are true, I don't know whether they remain true under these weaker assumptions, but I conjecture that they do.

Finding these bounds relies, in part, on finding upper bounds for the $$d$$-th central moments of the binomial($$n$$, $$p$$) distribution. This is discussed in Molteni (2022); before Molteni there were almost no works on upper bounds for those moments.

When $$d$$ is even, the moment is no greater than $$A_d n^{d/2}$$, where $$A_d$$ is a constant that depends on $$d$$ (DeVore and Lorentz 1993). The goal here is to find this constant.

The $$d$$-th moment is upper bounded by $$\frac{d!}{(d/2)!8^{d/2}} n^{d/2}$$ for even $$r\le 44$$ (Molteni 2022).

• The 3rd moment is bounded by $$\frac{\sqrt{3}}{18\sqrt{n}} n^{3/2}$$, so is bounded by $$\frac{\sqrt{3}}{18} n^{3/2} < (963/10000) n^{3/2}$$ for every $$n\ge 1$$.
• Proof: The critical points of the moment are at $$p=0$$, $$p=1$$, $$p=1/2-\sqrt{3}/6$$, and $$p=1/2+\sqrt{3}/6$$. The moment equals 0 at the points 0 and 1, so that leaves the last two. Since the odd moments are antisymmetric (Skorski 2020), it's enough to take the third critical point, where the moment is positive. By inspection, the third moment at that critical point is decreasing for every $$n\ge 1$$.
• The 5th moment is bounded by $$0.083 n^{5/2}$$.
• Proof: Follows from evaluating the moment for each $$1\le n \le 303$$ at their critical points (resulting in $$< 0.083 n^{5/2}$$ for every such $$n$$), and for each $$n\ge 304$$, applying the bound $$n/4+2 {n \choose 2}$$ based on the upper bound given in sec. 3.1 of Skorski (2020). This latter bound, when divided by $$n^{5/2}$$, decreases with $$n$$, so it's enough to take its value at 304, namely $$92188 < 0.05722/n^{5/2}$$ for every $$n\ge 304$$.

The following is Python code I used to calculate the error bounds for the iterated Bernstein polynomials. It uses the SymPy computer algebra library.

def tnr(n,r):
if r%2==0 and r<=44:
return (factorial(r)/(factorial(r//2)*8**(r//2)))*n**(r//2)
if r==1: return 0
if r==3: return (S(963)/10000)*sqrt(n**3)
if r==5: return (S(83)/1000)*sqrt(n**5)
return 2*factorial(S(r)/2)*sqrt(n**r)
raise ValueError

def gnr1(n,r,derivs):
return (Max(derivs[0],derivs[r+1])/n**(r+1))*sqrt(tnr(n,2*(r+1)))/factorial(r+1)

def bnerror(n,r,derivs):
return sum(derivs[i]*tnr(n,i)/(n**i*factorial(i)) for i in range(2,r+1))+gnr1(n,r,derivs)

def fnrderivs(n,r,alpha,derivs):
d=[]
for beta in range(0,(r+1-alpha)+1):
d.append((Max(derivs[0],derivs[r+1])/(n**floor(alpha/2)*factorial(alpha)))*\
sum(binomial(beta,g)*tnr(n,alpha)*factorial(alpha)/factorial(alpha-g) \
for g in range(0,min(alpha,beta)+1)))
return d

def bnr(n,rr,derivs,r=None,gnd=True): # r=s-1, s is no. of cont. deriv.
s=len(derivs)-1 # No. of continuous derivatives
if r==None: r=s-1
if rr==1: return bnerror(n,r,derivs)
gn=gnr1(n,r,derivs)
return sum(bnr(n,rr-1,fnrderivs(n,r,i,derivs))/n**((i+1)//2) \
for i in range(2,r+1))+gn*2**(rr-1)

n=symbols('n',nonnegative=True,integer=True)
d0,d1,d2,d3,d4,d5,d6=symbols('d0 d1 d2 d3 d4 d5 d6',real=True,positive=True)
print(bnr(n,2,[d0,d1,d2,d3]).simplify().n())
print(bnr(n,2,[d0,d1,d2,d3,d4]).simplify().n())
print(bnr(n,3,[d0,d1,d2,d3,d4,d5]).simplify().n())
print(bnr(n,3,[d0,d1,d2,d3,d4,d5,d6]).simplify().n())


## References

• C.S. Güntürk, W. Li, "Approximation of functions with one-bit neural networks", arXiv:2112.09181 [cs.LG], 2021.
• Skorski, Maciej. "Handy formulas for binomial moments." arXiv preprint arXiv:2012.06270 (2020).
• DeVore, R.A., Lorentz, G.G., Constructive approximation, 1993.
• Molteni, Giuseppe. "Explicit bounds for even moments of Bernstein’s polynomials." Journal of Approximation Theory 273 (2022): 105658.
• Can you state any of these bounds as self-contained propositions, that can be understood without having to read the whole question or the whole answer? Jul 25, 2022 at 9:46