Computing the Bézier control points of the lagrangian interpolant in arbitrary dimension

The Bernstein-Bézier form of a polynomial is widely used in the fields of computer aided geometric design, spline approximation theory, and, more recently, in high order finite element methods for the solution of partial differential equations. However, if one wishes to compute the classical Lagrange interpolant relative to the Bernstein basis, then the resulting Bernstein-Vandermonde matrix is found to be highly ill-conditioned (though not as severe as in the Vandermonde case). In the univariate case of degree n, Marco and Martínez [Linear Algebra Appl., 422(2007), pp. 616-628] showed, using an approach based on Neville elimination, that one can obtain an O(n²) algorithm for solving the system, which also exploits the total positivity of the Bernstein basis. Remarkable as it may be, the Marco-Martínez algorithm has some drawbacks: The derivation of the algorithm is quite technical; the interplay between the ideas of total positivity and Neville elimination are not part of the standard armory of many nonspecialists; and the algorithm does not seem to extend to higher dimensional simplices. The present work addresses these issues. An alternative algorithm for solving the univariate Bernstein-Vandermonde linear system is presented that has the same complexity as the Marco-Martínez algorithm and whose stability does not seem to be in any way inferior; a simple derivation using only the basic theory of Lagrange interpolation (at least in the univariate case); and a natural generalization to the multivariate case.