Handy Polynomial Fitting with Bernstein Polynomials

[Thomas Klemm]

I assumed you used the Bernstein Expansion to approximate a continuous function.

Here's a Python program that does this expansion:

Code:

from math import factorial

def binominal(n, k):
    return factorial(n) // factorial(k) // factorial(n-k)

def bernstein(n, k, x):
    return binominal(n, k) * x**k * (1-x)**(n-k)

def expand(f, x, n):
    return sum([f(j/n) * bernstein(n, j, x) for j in range(n+1)])

We can now compare the different results:

Code:

n = 10
for j in range(n+1):
    x = j / n
    f = lambda x: x**2
    approx = lambda x: x**2 + x*(1-x)/n
    print("%.1f: %.4f %.4f %.4f" % (x, f(x), expand(f, x, n), approx(x)))

x:   f(x)   expand approx
0.0: 0.0000 0.0000 0.0000
0.1: 0.0100 0.0190 0.0190
0.2: 0.0400 0.0560 0.0560
0.3: 0.0900 0.1110 0.1110
0.4: 0.1600 0.1840 0.1840
0.5: 0.2500 0.2750 0.2750
0.6: 0.3600 0.3840 0.3840
0.7: 0.4900 0.5110 0.5110
0.8: 0.6400 0.6560 0.6560
0.9: 0.8100 0.8190 0.8190
1.0: 1.0000 1.0000 1.0000


We can see that the expansion and the approximation agree.
That's astonishing since the expansion is based on Bernstein polynomials of degree 10 whereas the approximation is only of degree 2.

We can also see that approximation isn't really good except at \(x=0\) and \(x=1\).

---

However you appear to calculate a best fit of a function using Bernstein polynomials of a certain degree.

Here's how we can do this with Python:

Code:

from numpy import array, dot
from numpy.linalg import solve

b = array([x**2 for j in range(m+1) for x in [j/m]])
A = array([[bernstein(n, k, x) for k in range(n+1)] for j in range(m+1) for x in [j/m]])

These are the intermediate results:

>>> b
array([0. , 0.01, 0.04, 0.09, 0.16, 0.25, 0.36, 0.49, 0.64, 0.81, 1. ])

>>> A
array([[1.  , 0.  , 0.  ],
       [0.81, 0.18, 0.01],
       [0.64, 0.32, 0.04],
       [0.49, 0.42, 0.09],
       [0.36, 0.48, 0.16],
       [0.25, 0.5 , 0.25],
       [0.16, 0.48, 0.36],
       [0.09, 0.42, 0.49],
       [0.04, 0.32, 0.64],
       [0.01, 0.18, 0.81],
       [0.  , 0.  , 1.  ]])


For a best fit we solve the system that we get when both A and b are multiplied by A.T, the transpose of A:

Code:

solve(A.T @ A, A.T @ b)

The result is:

array([0., 0., 1.])

That's not a surprise as \(B_{2,2}(x) = x^2\).
However we would get the same result with fewer or more points as long as they are exactly on the curve as in Valentin's example.

---

What's the benefit of using Bernstein polynomials of order 2 compared to the ordinary basis \(1, x, x^2\)?

And why would you need an additional constant?
Since \(\sum_{i=0}^n B_{i,n}=1\) the constant 1 is already present as a linear combination.

Kind regards
Thomas