What is the best way to write a polynomial in the Bernstein basis?

Question

The Bernstein basis of polynomials of degree $n$ is the set of polynomials of the form

$$\binom{n}{k} t^{n-k}(1-t)^k$$

where $0 \leq k \leq n$.

What is the best way to transform a given polynomial into a list of coefficients for the Bernstein basis? i.e., I want to know if Mathematica has (or if not, the best way to write) a function which accepts as input a polynomial expression in $t$ of degree $n$ and produces as output a list of coefficients where the $(k+1)$-th element of the list corresponds to the basis element $\binom{n}{k}t^{n-k}(1-t)^k$ (apologies for the off-by-one indexing mess).

Answer 1

I'm not aware of any built in function, but here's one way to do it.

bb[n_, k_, t_] := Binomial[n, k] t^(n - k) (1 - t)^k;

A test cubic, and a sum of our basis functions:

p = -(1/12) + (5 x)/8 - (17 x^2)/12 + x^3;
s = Sum[a[k] bb[3, k, x], {k, 0, 3}]

Now solve for the coefficients:

Solve@Table[Coefficient[s, x, k] == Coefficient[p, x, k], {k, 0, 3}]

{{a[2] -> 1/8, a[3] -> -(1/12), a[0] -> 1/8, a[1] -> -(5/36)}}

Answer 2

In this paper of Farouki and Rajan, it is noted that there is a simple method to interconvert between the monomial and Bernstein forms: one merely needs to multiply an appropriate matrix with the coefficients.

Here, then, are the routines for generating the required matrices:

(* monomial to Bernstein conversion *)
pbmat[n_Integer?NonNegative] := 
 SparseArray[{{j_, k_} /; j >= k :> Binomial[n - k + 1, j - k]/Binomial[n, j - 1]},
             {n + 1, n + 1}]

(* Bernstein to monomial conversion *)
bpmat[n_Integer?NonNegative] := 
 SparseArray[{{j_, k_} /; j >= k :> (-1)^(j - k) Binomial[n, j - 1] Binomial[j - 1, k - 1]},
             {n + 1, n + 1}]

The matrices generated by pbmat[] and bpmat[] are necessarily inverses of each other; that is, pbmat[n].bpmat[n] == IdentityMatrix[n+1].

Now, an example:

With[{n = 4}, pbmat[n].CoefficientList[LaguerreL[n, x], x]]
   {1, 0, -1/2, -2/3, -5/8}

Check:

Assuming[0 < x < 1, 
         PiecewiseExpand[%.Table[BernsteinBasis[4, k, x], {k, 0, 4}]] == 
         LaguerreL[4, x]] // Simplify
   True

Another example:

cofs = {3/2, 1/3, 2/5, 1};
Assuming[0 < x < 1, 
         cofs.Table[BernsteinBasis[3, k, x], {k, 0, 3}] // PiecewiseExpand // Expand]
   3/2 - 7 x/2 + (37 x^2)/10 - (7 x^3)/10

Verify:

bpmat[3].cofs
   {3/2, -7/2, 37/10, -7/10}

---

One can exploit the recurrence relations satisfied by the binomial coefficients, so that the basis conversion matrices don't need to be generated in full. Here are conversion routines implementing this idea:

(* monomial to Bernstein conversion *)
pbconv[poly_, x_] /; PolynomialQ[poly, x] := Module[{bc, n},
                  n = Exponent[poly, x];
                  bc = CoefficientList[poly, x]/Binomial[n, Range[0, n]];
                  Do[bc[[j + 1]] += bc[[j]], {k, n, 1, -1}, {j, k, n}];
                  bc]

(* Bernstein to monomial conversion *)
bpconv[coeffs_?VectorQ] := Module[{n = Length[coeffs] - 1, p = coeffs},
                       Do[p[[j + 1]] -= p[[j]], {k, n, 1, -1}, {j, k, n}];
                       p = TakeWhile[p, # != 0 &];
                       p Binomial[n, Range[0, Length[p] - 1]]]

Using the same examples given above:

pbconv[LaguerreL[4, x], x]
   {1, 0, -1/2, -2/3, -5/8}

bpconv[{3/2, 1/3, 2/5, 1}]
   {3/2, -7/2, 37/10, -7/10}

---

Now, a warning. As noted by Farouki in this paper, the conversion between the Bernstein and monomial basis is ill-conditioned for sufficiently high degrees, since the condition number of the matrices involved increases sharply with the degree. Here is an illustration:

Table[Ceiling[Log10[LinearAlgebra`MatrixConditionNumber[pbmat[k]]]], {k, 15}]
   {1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8}

This, for instance, says that you stand to lose $8$ significant figures when converting a degree-$15$ polynomial in the monomial basis to Bernstein form.

Answer 3

Here is another solution, of course similar to that of wxffles, but using SolveAlways. The function takes any polynomial in one variable as argument and returns the list of coefficients.

bernsteincoefficients[pol_] := Module[{t, n, a},
  t = Variables[pol][[1]];
  n = Exponent[pol, t];
  Table[a[k], {k, 0, n}] /. 
    SolveAlways[Sum[a[k] Binomial[n, k] t^(n - k) (1 - t)^k, {k, 0, n}] == pol, t][[1]]
]