Return to Question

The following Python code shows my attempt to implement the Holtz approximation scheme using the SymPy computer algebra library. It includes an implementation of the Lorentz operator as well as a method to find Bernstein polynomials for a function.

The following Python code (as well as lorentz.py) shows my attempt to implement the Holtz approximation scheme using the SymPy computer algebra library. It includes an implementation of the Lorentz operator as well as a method to find Bernstein polynomials for a function.

1. What are practical lower bounds for $s$, $\theta_{\alpha}$, and $D$ for the method above, given a function $f$?
2. What are other practical ways to apply the Holtz method to certain functions to find polynomials that converge to those functions?
3. To help me and others (especially programmers) understand and apply the Holtz method, can you provide examples of the method (in terms of Bernstein polynomial approximation schemes) for various functions or classes of functions (such as $C^2$ continuous, twice differentiable, or real analytic functions)?

TAURJ_CACHE = {}
TNJ_CACHE = {}
TAURJ_SYMBOL = symbols("T")

def taurj_x(r, j, n):
    if j == 0:
        return S(1)
    if j == 1:
        return S(0)
    if (r, j, n) in TAURJ_CACHE:
        return TAURJ_CACHE[r, j, n]
    ret = 0
    for l in range(2, j + 1):
        tnj = None
        if (n, l) in TNJ_CACHE:
            tnj = TNJ_CACHE[(n, l)]
        else:
            tnj = sum(
                (k - n * TAURJ_SYMBOL) ** l
                * binomial(n, k)
                * TAURJ_SYMBOL ** k
                * (1 - TAURJ_SYMBOL) ** (n - k)
                for k in range(0, n + 1)
            )
        ret -= tnj * taurj_x(r - l, j - l, n) / factorial(l)
    TAURJ_CACHE[(r, j, n)] = ret.simplify()
    return ret


def taurj(r, j, x, n):
    if j == 0:
        return 1
    if j == 1:
        return 0
    return taurj_x(r, j, n).subs(TAURJ_SYMBOL, x)
    
def lorentz_n_k(func, x, n, r, k, pt=None):
    # kth Bernstein coefficient of the Lorentz operator
    # of degree n and smoothness r at the point
    # pt (or at x if not given).
    if pt == None:
        pt = x
    ret = 0
    for j in range(0, r + 1):
        ret += diff(func, (x, j)).subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
    return ret  


def lorentz_n(func, x, n, r, pt=None):
    # Lorentz operator.
    # Create a polynomial that approximates func, which in turn uses
    # the symbol x.  The polynomial's degree is n, is assumed to
    # be at least r times differentiable, and is evaluated
    # at the point pt (or at x if not given).
    # NOTE: If r is 0, this is simply the Bernstein operator.
    if pt == None:
        pt = x
    ret = 0
    diffs = [diff(func, (x, j)) for j in range(0, r + 1)]
    for k in range(0, n + 1):
        rm = 0
        for j in range(0, r + 1):
            rm += diffs[j].subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
        ret += rm * bernstein_n_k(n, k, pt)
    return ret

 

def ncfo_lower(func, x, n, pt, alpha, theta, d):
    # NOTE: Assumes func is bounded away from 0 and 1,
    # but it's trivial to bound any function this way, such as
    # by doing a linear transformation such as
    # $f(x) * \frac{8}{10} + \frac{1}{10}$, then undoing
    # the transformation once the polynomials are found
    r = floor(alpha)
    return lorentz_n(func, x, n, r, pt) - d * bernstein_n(
        phi_n(n, theta, alpha, x), x, n
    )


def ncfo_upper(func, x, n, pt, alpha, theta, d):
    r = floor(alpha)
    return lorentz_n(func, x, n, r, pt) + d * bernstein_n(
        phi_n(n, theta, alpha, x), x, n
    )
      

1. What are practical lower bounds for $s$, $\theta_{\alpha}$, and $D$ for the method above, given a function $f$, with or without additional assumptions on $f$ (such as being $C^2$ continuous, concave, convex, piecewise polynomial, real analytic, or any combination of these)?
2. What are other practical ways to apply the Holtz method to certain functions to find polynomials that converge to those functions?
3. To help me and others (especially programmers) understand and apply the Holtz method, can you provide examples of the method (in terms of Bernstein polynomial approximation schemes) for various functions or classes of functions (such as $C^2$ continuous, concave, convex, piecewise polynomial, real analytic, or any combination of these)?

TAURJ_CACHE = {}
TNJ_CACHE = {}
TAURJ_SYMBOL = symbols("T")

def taurj_x(r, j, n):
    if j == 0:
        return S(1)
    if j == 1:
        return S(0)
    if (r, j, n) in TAURJ_CACHE:
        return TAURJ_CACHE[r, j, n]
    ret = 0
    for l in range(2, j + 1):
        tnj = None
        if (n, l) in TNJ_CACHE:
            tnj = TNJ_CACHE[(n, l)]
        else:
            tnj = sum(
                (k - n * TAURJ_SYMBOL) ** l
                * binomial(n, k)
                * TAURJ_SYMBOL ** k
                * (1 - TAURJ_SYMBOL) ** (n - k)
                for k in range(0, n + 1)
            )
        ret -= tnj * taurj_x(r - l, j - l, n) / factorial(l)
    TAURJ_CACHE[(r, j, n)] = ret.simplify()
    return ret


def taurj(r, j, x, n):
    if j == 0:
        return 1
    if j == 1:
        return 0
    return taurj_x(r, j, n).subs(TAURJ_SYMBOL, x)
    
def lorentz_n_k(func, x, n, r, k, pt=None):
    # kth Bernstein coefficient of the Lorentz operator
    # of degree n and smoothness r at the point
    # pt (or at x if not given).
    if pt == None:
        pt = x
    ret = 0
    for j in range(0, r + 1):
        ret += diff(func, (x, j)).subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
    return ret  


def lorentz_n(func, x, n, r, pt=None):
    # Lorentz operator.
    # Create a polynomial that approximates func, which in turn uses
    # the symbol x.  The polynomial's degree is n, is assumed to
    # be at least r times differentiable, and is evaluated
    # at the point pt (or at x if not given).
    # NOTE: If r is 0, this is simply the Bernstein operator.
    if pt == None:
        pt = x
    ret = 0
    diffs = [diff(func, (x, j)) for j in range(0, r + 1)]
    for k in range(0, n + 1):
        rm = 0
        for j in range(0, r + 1):
            rm += diffs[j].subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
        ret += rm * bernstein_n_k(n, k, pt)
    return ret

def bernstein_n_k(n, k, x):
    return binomial(n, k) * x ** k * (1 - x) ** (n - k)

This question is a continuation of another question of mine, which seeks ways to compute polynomials that converge from above and below to $f$ in a manner that solves the Bernoulli factory problem for $f$. These polynomials form an approximation scheme for $f$. See that question for a formal statement of such approximation schemes.

This question is a continuation of another question of mine, which seeks ways to compute polynomials that converge from above and below to $f$ in a manner that solves the Bernoulli factory problem for $f$. These polynomials form an approximation scheme for $f$. See that question for more details, but here is a recap:

Formal statement: There exist two sequences of polynomials, namely—

• $g_{n}(\lambda): =\sum_{k=0}^{n}a(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$, and
• $h_{n}(\lambda): =\sum_{k=0}^{n}b(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$,

for every integer $n\ge1$, such that—

1. $0\le a(n, k)\le b(n, k)\le1$,
2. $\lim_{n}g_{n}(\lambda)=\lim_{n}h_{n}(\lambda)=f(\lambda)$ for every $\lambda\in[0,1]$, and
3. for every $m<n$, both $(g_{n} - g_{m})$ and $(h_{m} - h_{n})$ have non-negative coefficients once $g_{n}$, $g_{m}$, $h_{n}$, and $h_{m}$ are rewritten as degree-$n$ polynomials in Bernstein form,

where $f(\lambda)$ is continuous on $[0, 1]$ (Nacu and Peres 2005; Holtz et al. 2011), and the goal is to find the appropriate values for $a(n, k)$ and $b(n, k)$.

It is allowed for $a(n, k)\lt0$ for a given $n$ and some $k$, in which case all $a(n, k)$ for that $n$ are taken to be 0 instead. It is allowed for $b(n, k)\gt1$ for a given $n$ and some $k$, in which case all $b(n, k)$ for that $n$ are taken to be 1 instead.