Today I demonstrated the Bernstein Operator

\(B_nf(x) = \sum\limits_{k=0}^n \dbinom{n}{k} x^k (1-x)^{n-k}. \)

This is known to converge for all continuous functions on [0,1], but its approximating properties are not very good. E.g., we have

\(B_nf(x) – f(x) = \dfrac{x(1-x)}{n} \qquad\text{for } f(x) = x^2.\)

To define such an operator in Euler Math Toolbox, I programmed

>function bernstein (x,f$,n) ...
$## Compute B_nf(x) for scalar x or row vectors x.
$  i=(0:n)'; B=bin(n,i)*x^i*(1-x)^(n-i);
$  return f$(i/n)'.B;
$endfunction

Then I proceed to use this operator in the following way.

>function f(x) := sqrt(x+1)
>n=10;
>x=0:0.01:1; y=bernstein(x,"f",n);
>shrinkwindow(>smaller); plot2d(x,y-f(x),title="n = "+n):

This plots the following error function.

The topic of this posting is to discuss if this is really a good way of doing things.

For the non-initiated, the shortness of the code for the function bernstein() must be frightening. So let me introduce the code step by step. The variable i is a column vector with elements 0 to n. Let us assume that x is a row vector. Then e.g. x^i is a matrix with rows 1, x, x^2 to x^n, where x^2 is the matrix row containing x(1)^2 to x(m)^2. Likewise for the other factors. This is combined by multiplication element for element. We get the following matrix.

\(B = \left(\dbinom{n}{i} x_j^i (1-x_j)^{n-i}\right)_{0 \le i \le n,\, 1 \le j \le m}\)

If we multiply this using matrix multiplication with the row vector

\(\left(f\left(\dfrac{0}{n}\right),f\left(\dfrac{1}{n}\right),\ldots,f\left(\dfrac{n}{n}\right)\right)\)

we get the row vector of values of the Bernstein approximation.

Even now, that I have explained it, you might have a hard time to get the details. My feeling and my experience tell me that the matrix languages, if used on this level, are for experts only. You can get rid of a lot of loops, but at the expense of having to think in row and column matrices and element-wise operators. For the casual user like a student, this requires some extra learning.

How could we get the very same result with less metal effort? One small step would be to define the operator in a less ambitious way.

>function map B(x,n) ...
$  k=0:n;
$  return sum(f(k/n)*bin(n,k)*x^k*(1-x)^(n-k));
$endfunction
>plot2d("B(x,10)-f(x)",0,1):

This is the same result, but we need only understand a row vector. The „map“ keyword allows to map B(x,n) to a vector of x values.

>B(0:0.1:1,10)
  [1, 1.04786, 1.09395, 1.13841, 1.1814, 1.22303, 1.26341,
  1.30263, 1.34079, 1.37796, 1.41421]

The effort to compute this is higher, however. We recompute f(k/n) and the binomial constants for all values of x. The matrix version calls the function f only n+1 times independent on the length of the row vector x.

The next step would be to avoid the matrix language completely. Then we have to write loops for the summations. Moreover, we would like to compute the vector of function values and the binomials beforehand. For a final optimization, we would have to compute x^i and (1-x)^i with one loop and combine both later. The binomials should be computed by looping too.

It becomes quite obvious now that an optimized version in C would take quite a lot of brain activity too, and we would not be sure if we get it right in the first shot. If you do not believe it try it. You can use Python or C via Euler Math Toolbox if you like.

Summarizing, the easiest way lies in the middle. It will not be completely optimized, but it will also not be obscure.