Structured inversion of the Bernstein mass matrix

TL;DR

This paper develops fast, accurate algorithms for inverting the univariate Bernstein mass matrix using explicit formulas, matrix decompositions, and spectral methods, with analysis of conditioning and numerical stability.

Contribution

It introduces several novel approaches for inverting the univariate Bernstein mass matrix, including explicit inverse formulas and spectral decomposition, improving computational efficiency and accuracy.

Findings

Eigendecomposition can be constructed in O(n^2) operations.

The spectral method's accuracy is comparable to Cholesky decomposition.

Conditioning in the L^2 norm is less severe than in the standard 2-norm.

Abstract

Bernstein polynomials, long a staple of approximation theory and computational geometry, have also increasingly become of interest in finite element methods. Many fundamental problems in interpolation and approximation give rise to interesting linear algebra questions. Previously, we gave block-structured algorithms for inverting the Bernstein mass matrix on simplicial cells, but did not study fast alorithms for the univariate case. Here, we give several approaches to inverting the univariate mass matrix based on exact formulae for the inverse; decompositions of the inverse in terms of Hankel, Toeplitz, and diagonal matrices; and a spectral decomposition. In particular, the eigendecomposition can be explicitly constructed in $\mathcal{O}(n^2)$ operations, while its accuracy for solving linear systems is comparable to that of the Cholesky decomposition. Moreover, we study conditioning and accuracy of these methods from the standpoint of the effect of roundoff error in the $L^2$ norm on polynomials, showing that the conditioning in this case is far less extreme than in the standard 2-norm.

Equations

$$M\mathbf{c}=\mathbf{b},$$

$$B_i^n(x)=\left(\genfrac{}{}{0ex}{}{n}{i}\right)x^i(1-x{)}^{n-i}.$$

$$E_{ij}^{m,n}=\frac{\left(\genfrac{}{}{0ex}{}{m}{j}\right)\left(\genfrac{}{}{0ex}{}{n-m}{i-j}\right)}{\left(\genfrac{}{}{0ex}{}{n}{i}\right)},$$

$${\Vert L^n\Vert }_{L^2}^2=\frac{1}{2n+1}.$$

$$L^n(x)=\sum _{i=0}^n(-1{)}^{n+i}(\genfrac{}{}{0ex}{}{n}{i})B_i^n(x).$$

$$M_{ij}^n={\int }_0^1{B}_i^n(x)B_j^n(x)dx,$$

$$M_{ij}^n=\left(\genfrac{}{}{0ex}{}{n}{i}\right)\left(\genfrac{}{}{0ex}{}{n}{j}\right)\frac{(2n-i-j)!(i+j)!}{(2n+1)!}.$$

$$\pi (\mathbf{c})(x)=\sum _{i=0}^n{\mathbf{c}}_i{B}_i^n(x).$$

$${\int }_0^{1}p(x)q(x)dx=\sum _{i,j=0}^n{\mathbf{p}}_i{\mathbf{q}}_j{\int }_0^1{B}_i^n(x)B_j^n(x)dx={\mathbf{p}}^T{M}^n\mathbf{q}.$$

$$\Vert \mathbf{p}{\Vert }_{M^n}=\sqrt{{\mathbf{p}}^T{M}^n\mathbf{p}}$$

$$\Vert p{\Vert }_{L^2}=\Vert \mathbf{p}{\Vert }_{M^n}.$$

$$M^m={\left(E^{m,n}\right)}^T{M}^n{E}^{m,n}.$$

$${\Delta }_{ij}^n=\{\begin{array}{ll}\left(\genfrac{}{}{0ex}{}{n}{i}\right),& i=j;\\ 0,& \text{otherwise}.\end{array}$$

$${\tilde{M}}_{ij}^n=\frac{(2n-i-j)!(i+j)!}{(2n+1)!},$$

$$M^n={\Delta }^n{\tilde{M}}^n{\Delta }^n.$$

$${\lambda }_i^n=\frac{(n!{)}^2}{(n+i+1)!(n-i)!},$$

$${\left(M^n\right)}_{ij}^{-1}=\frac{(-1{)}^{i+j}}{\left(\genfrac{}{}{0ex}{}{n}{i}\right)\left(\genfrac{}{}{0ex}{}{n}{j}\right)}\sum _{k=0}^n(2k+1-i+j){\left(\genfrac{}{}{0ex}{}{n+1}{i-k}\right)}^2{\left(\genfrac{}{}{0ex}{}{n+1}{j+k+1}\right)}^2.$$

$${\int }_0^1{B}_i^n(x)d_j^n(x)dx=\{\begin{array}{ll}1,& i=j;\\ 0,& \text{otherwise}.\end{array}$$

$${\int }_0^1{B}_i^n(x)d_j^n(x)dx=\sum _{k=0}^n{\mathbf{d}}_k^{n,j}{\int }_0^1{B}_i^n(x)B_k^n(x)dx={\left(M^n{\mathbf{d}}^{n,j}\right)}_i,$$

$${\left(M^n\right)}_{ij}^{-1}=\frac{(-1{)}^{i+j}}{\left(\genfrac{}{}{0ex}{}{n}{i}\right)\left(\genfrac{}{}{0ex}{}{n}{j}\right)}\sum _{k=0}^n(2k+1)\left(\genfrac{}{}{0ex}{}{n+k+1}{n-j}\right)\left(\genfrac{}{}{0ex}{}{n-k}{n-j}\right)\left(\genfrac{}{}{0ex}{}{n+k+1}{n-i}\right)\left(\genfrac{}{}{0ex}{}{n-k}{n-i}\right).$$

$$\frac{u(s)v(t)-u(t)v(s)}{s-t}=\sum _{i,j=0}^n{b}_{ij}{s}^i{t}^j.$$

$$b_{ij}=\sum _{k=0}^{\min \left\{i,n-j\right\}}\left(u_{j+k+1}{v}_{i-k}-u_{i-k}{v}_{j+k+1}\right).$$

$$H^{-1}=\frac{1}{v_{n+1}}\mathrm{Bez}(u,v),$$

$$0={\int }_0^{1}p(x)L^n(x)dx.$$

$$0=(E^{n-1,n}\mathbf{p}{)}^T{M}^n\mathbf{\Pi }(L^n)={\lambda }_n^n{\mathbf{q}}^T\mathbf{\Pi }(L^n)=\frac{(n!{)}^2}{(2n+1)!}{\mathbf{q}}^T\mathbf{\Pi }(L^n).$$

$${\mathbf{y}}_i^n=(-1{)}^{n+i}(n+1)(\genfrac{}{}{0ex}{}{n+1}{i}).$$

$${\mathbf{y}}^n=(2n+1)\mathbf{\Pi }(L^n)+E^{n-1,n}{\mathbf{y}}^{n-1}.$$

$${\left[(2n+1)\mathbf{\Pi }(L^n)+E^{n-1,n}{\mathbf{y}}^{n-1}\right]}_i$$

$$=(2n+1)(-1{)}^{n+i}(\genfrac{}{}{0ex}{}{n}{i})+i(-1{)}^{n+i}\left(\genfrac{}{}{0ex}{}{n}{i-1}\right)-(n-i)(-1{)}^{n+i}(\genfrac{}{}{0ex}{}{n}{i})$$

$$=(-1{)}^{n+i}\left[(n+i+1)\left(\genfrac{}{}{0ex}{}{n}{i}\right)+i\left(\genfrac{}{}{0ex}{}{n}{i-1}\right)\right]$$

$$=(-1{)}^{n+i}(n+1)(\genfrac{}{}{0ex}{}{n+1}{i})$$

$$={\mathbf{y}}_i^n.$$

Full text

1 Introduction

Bernstein polynomials, long used in splines, computer-aided geometric design, and computer graphics were first introduced more than a century ago [5]. More recently, they have also been considered as a tool for high-order approximation of partial differential equations via the finite element method. Tensorial decomposition of the simplicial Bernstein polynomials under the Duffy transform [6] has led to fast algorithms with comparable complexity to spectral element methods in rectangular geometry [1]. Moreover, Bernstein properties also lead to special recursive blockwise-structure for finite element matrices [12, 13]. Bernstein polynomials are not only a tool for standard $C^{0}$ finite element spaces, but also can be used to represent bases for polynomial differential forms discretizing the de Rham complex [2, 11].

In this paper, we focus on the inversion of the one-dimensional mass matrix. In [12], we gave recursive, block-structured fast algorithms for mass matrices on the $d$-simplex. These rely on inversion of the one-dimensional mass matrix as an essential building block, but we did not consider fast algorithms or the underlying stability issues for that case. We turn to these issues here. In addition to studying the inverse matrix and fast algorithms for applying it, we also give an argument that the effective conditioning for Bernstein polynomials is less extreme than it might appear. This amounts to introducing a nonstandard matrix norm in which to consider the problem. In fact, our main result on conditioning mass matrices – that in the “right” norm the condition number is the square root of the standard matrix 2-norm – is generic and not limited to the particularities of either univariate polynomials or the Bernstein basis. We believe this sheds important light on the observation that, despite massive condition numbers, Bernstein polynomials of degrees greater than ten can still give highly accurate finite element methods [3, 11]. It does also, however, suggest an eventual limit to the degree of polynomial absent algorithmic advances.

The mass (or Gram) matrix arises abstractly in the case of finding the best approximation onto a subspace of a generic Hilbert space. Given a Hilbert space $H$ with inner product $(\cdot,\cdot)$ and a finite-dimensional subspace $V^{N}$ with basis $\{\phi_{i}\}_{i=0}^{N}$, the problem of optimally approximating $u\in H$ by $\sum_{i=0}^{N}c_{i}\phi_{i}$ with respect to the norm of $H$ requires solving the linear system

[TABLE]

where $M_{ij}=(\phi_{i},\phi_{j})$ and $\mathbf{b}_{i}=(u,\phi_{i})$.

While choosing an orthogonal basis trivializes the inversion of $M$, matrices of this type also arise in other settings that constrain the choice of basis. For example, finite element methods typically require “geometrically decomposed” bases [4] that enable enforcement of continuity across mesh elements. Even when discontinuous bases are admissible for finite element methods, geometric decomposition can be used to simplify boundary operators [12]. Moreover, splines and many geometric applications make use of the Bernstein-Bezier basis as local representations.

2 Notation and Mass matrix

For integers $n\geq 0$ and $0\leq i\leq n$, we define the Bernstein polynomial $B^{n}_{i}(x)$ as

[TABLE]

Each $B^{n}_{i}(x)$ is a polynomial of degree $n$, and the set $\left\{B^{n}_{i}(x)\right\}_{i=0}^{n}$ forms a basis for the vector space of polynomials of degree $n$.

If $m\leq n$, then any polynomial expressed in the basis $\left\{B^{m}_{i}(x)\right\}_{i=0}^{m}$ can also be expressed in the basis $\left\{B^{n}_{i}(x)\right\}_{i=0}^{n}$. We denote by $E^{m,n}$ the $(n+1)\times(m+1)$ matrix that maps the coefficients of the degree $m$ representation to the coefficients of the degree $n$ representation. It is remarked in [7] that the entries of $E^{m,n}$ are given by

[TABLE]

with the standard convention that $\binom{n}{i}=0$ for $i<0$ or $i>n$.

We will also need reference to the Legendre polynomials, mapped from their typical home on $[-1,1]$ to $[0,1]$. We let $L^{n}(x)$ denote the Legendre polynomial of degree $n$ over $[0,1]$, scaled so that $L^{n}(1)=1$ and

[TABLE]

It was observed in [7] that $L^{n}(x)$ has the following representation:

[TABLE]

The Bernstein mass matrix $M^{n}$ is given by

[TABLE]

which can be exactly computed [10] as

[TABLE]

The mass matrix, whether with Bernstein or any other polynomials, plays an important role in connecting the $L^{2}$ topology on the finite-dimensional space to linear algebra. To see this, we first define mappings connecting polynomials of degree $n$ to $\mathbb{R}^{n+1}$. Given any $\mathbf{c}\in\mathbb{R}^{n+1}$, we let $\pi(\mathbf{c})$ be the polynomial expressed in the Bernstein basis with coefficients contained in $\mathbf{c}$:

[TABLE]

We let $\mathbf{\Pi}$ be the inverse of this mapping, sending any polynomial of degree $n$ to the vector of $n+1$ coefficients with respect to the Bernstein basis.

Now, let $p(x)$ and $q(x)$ be polynomials of degree $n$ with expansion coefficients $\mathbf{\Pi}(p)=\mathbf{p}$ and $\mathbf{\Pi}(q)=\mathbf{q}$. Then the $L^{2}$ inner product of $p$ and $q$ is given by the $M^{n}$-weighted inner product of $\mathbf{p}$ and $\mathbf{q}$, for

[TABLE]

Similarly, if

[TABLE]

is the $M^{n}$-weighted vector norm, then we have for $p=\pi(\mathbf{p})$,

[TABLE]

It was shown in [10] that if $m\leq n$, then

[TABLE]

Further, in [10] it was observed that the mass matrix is, up to a row and column scaling, a Hankel matrix (constant along antidiagonals). Let $\Delta^{n}$ be the $(n+1)\times(n+1)$ diagonal matrix with binomial coefficients on the diagonal:

[TABLE]

Then, let $\widetilde{M}^{n}$ be defined by

[TABLE]

so that

[TABLE]

Since $\widetilde{M}^{n}$ depends only on the sum $i+j$ and not $i$ and $j$ separately, it is a Hankel matrix. Therefore, $\widetilde{M}^{n}$ and hence $M^{n}$ can be applied to a vector in $\mathcal{O}(n\log n)$ operations via a circulant embedding and fast Fourier transforms [16], although the actual point where this algorithm wins over the standard one may be at a rather high degree.

In [12], we gave a spectral characterization of the mass matrix on simplices of any dimension. For the one-dimensional case, the eigenvalues and eigenvectors are:

Theorem 2.1**.**

The eigenvalues of $M^{n}$ are $\{\lambda_{i}^{n}\}_{i=0}^{n}$, where

[TABLE]

and the eigenvector of $M^{n}$ corresponding to $\lambda_{i}^{n}$ is $E^{i,n}\mathbf{\Pi}(L^{i})$.

3 Characterizing the inverse matrix

In this section, we present several approaches to constructing and applying the inverse of the mass matrix.

3.1 A formula for the inverse

The following result (actually, a generalization) is derived in [14]:

Theorem 3.1**.**

[TABLE]

This result is obtained by characterizing (possibly constrained) Bernstein dual polynomials via Hahn orthogonal polynomials. Much earlier, Farouki [7] gave a characterization of the standard dual polynomials by more direct means. Lu [15] builds on the work in [14] to give a recursive formula for applying this inverse.

In this section, we summarize a proof of Theorem 3.1 based on Farouki’s representation of the dual basis in Subsection 3.1.1, and then give an alternative derivation of the result based on Bézoutians in Subsection 3.1.2. Since the mass matrix is diagonal in the Legendre basis, we can also characterize the inverse operator by converting from Bernstein to Legendre representation, multiplying by a diagonal, and then converting back. This approach, equivalent to the spectral decomposition, is given in Subsection 3.1.3.

3.1.1 Inversion via the dual basis

The dual basis to $\left\{B^{n}_{i}(x)\right\}_{i=0}^{n}$ is the set of polynomials $\left\{d^{n}_{i}(x)\right\}_{i=0}^{n}$ satisfying

[TABLE]

If we consider $\mathbf{d}^{n,j}=\mathbf{\Pi}(d^{n}_{j})$, then

[TABLE]

and hence (18) implies that $\left(M^{n}\right)^{-1}_{ij}=\mathbf{d}^{n,j}_{i}$. Following Farouki’s representation of the dual basis coefficients [7], this gives us that

[TABLE]

Applying a few binomial identities gives the representation in Theorem 3.1.

3.1.2 Inversion via Bezoutians

In this section, we detail an alternate derivation of Theorem 3.1 via Bézout matrices, which we will use in Section 3.2 to decompose the inverse in terms of diagonal, Hankel, and Toeplitz matrices.

If $u(t)=\sum_{i=0}^{n+1}u_{i}t^{i}$ and $v(t)=\sum_{i=0}^{n+1}v_{i}t^{i}$ are polynomials of degree less than or equal to $(n+1)$ expressed in the monomial basis, then the Bézout matrix generated by $u$ and $v$, denoted $\operatorname{Bez}(u,v)$, is the $(n+1)\times(n+1)$ matrix with entries $b_{ij}$ satisfying

[TABLE]

By comparing coefficients, we have that a closed form expression for the entries is given by

[TABLE]

Heinig and Rost [9] gave a formula for the inverse of a Hankel matrix in terms of a Bézout matrix.

Theorem 3.2**.**

If $H$ is an $(n+1)\times(n+1)$ nonsingular Hankel matrix, and $\widehat{H}$ is any $(n+2)\times(n+2)$ nonsingular Hankel extension of $H$ obtained by appending a row and column to $H$, then

[TABLE]

where $u$ is the polynomial whose coefficients are given by the last column of $H^{-1}$, and $v$ is the polynomial whose coefficients are given by the last column of $\widehat{H}^{-1}$.

Since the matrices $H$ and $\widehat{H}$ are of different sizes, the corresponding polynomials have different degrees. We reconcile this difference by elevating by one degree in the monomial basis (that is, appending a zero to the vector of coefficients).

The next few results are dedicated to finding the $u^{n}$ and $v^{n+1}$ that correspond to $\widetilde{M}^{n}$. We begin by showing that the null space of $\left(E^{n-1,n}\right)^{T}$ is spanned by the eigenvector corresponding to $\lambda^{n}_{n}$. We then use this to project our candidate $\mathbf{y}^{n}$ for the last column of $\left(M^{n}\right)^{-1}$ onto the null space and its orthogonal complement. This decomposition gives a recurrence relation for the $\mathbf{y}^{n}$, which will be the foundation for an inductive proof that $\mathbf{y}^{n}$ is the last column of $\left(M^{n}\right)^{-1}$. The coefficients of $u^{n}$ are then obtained via (15).

Lemma 3.2.1**.**

The null space of $\left(E^{n-1,n}\right)^{T}$ is spanned by $\mathbf{\Pi}(L^{n})$.

Proof.

By (3), we have that $\left(E^{n-1,n}\right)^{T}$ is upper bidiagonal with nonzero entries on the main diagonal and the superdiagonal, and so the null space of $\left(E^{n-1,n}\right)^{T}$ is one-dimensional. Therefore, it is enough to show that $\mathbf{\Pi}(L^{n})$ belongs to the null space of $\left(E^{n-1,n}\right)^{T}$.

Let $\mathbf{q}=E^{n-1,n}\mathbf{p}$ be in the range of $E^{n-1,n}$. This means that $p=\pi(\mathbf{p})$ is a polynomial of degree at most $n-1$, and hence

[TABLE]

Therefore, Theorem 2.1 and (9) imply that

[TABLE]

This implies that $\mathbf{\Pi}(L^{n})$ is orthogonal in the Euclidean inner product to the range of $E^{n-1,n}$, and so $\mathbf{\Pi}(L^{n})$ belongs to the null space of $\left(E^{n-1,n}\right)^{T}$. ∎

Lemma 3.2.2**.**

For $n\geq 0$, let $\mathbf{y}^{n}$ be given by

[TABLE]

Then, for $n\geq 1$,

[TABLE]

Proof.

By (3) and (5), we have that for each $0\leq i\leq n$,

[TABLE]

∎

Proposition 3.2.1**.**

Let $\mathbf{y}^{n}$ be as in Lemma 3.2.2. Then

[TABLE]

where

[TABLE]

Proof.

We show the result by induction on $n$. Clearly, the result is true for $n=0$. So suppose $M^{n-1}\mathbf{y}^{n-1}=\mathbf{e}^{n-1}$ for some $n\geq 1$. We show that $M^{n}\mathbf{y}^{n}=\mathbf{e}^{n}$ by showing that $M^{n}\mathbf{y}^{n}-\mathbf{e}^{n}$ belongs to both the null space of $\left(E^{n-1,n}\right)^{T}$ and its orthogonal complement with respect to the Euclidean inner product.

Multiplying (27) on the left by $M^{n}$ gives us that

[TABLE]

Since the null space of $\left(E^{n-1,n}\right)^{T}$ is spanned by $\mathbf{\Pi}(L^{n})$, multiplying the previous equation on the left by $\left(E^{n-1,n}\right)^{T}$ and using (12) and the induction hypothesis gives us that

[TABLE]

Since $\mathbf{e}^{n-1}=\left(E^{n-1,n}\right)^{T}\mathbf{e}^{n}$, the previous equation implies that $M^{n}\mathbf{y}^{n}-\mathbf{e}^{n}$ belongs to the null space of $\left(E^{n-1,n}\right)^{T}$.

On the other hand, (30) implies that

[TABLE]

By (5), we have that

[TABLE]

and $\mathbf{\Pi}(L^{n})_{n}=1$. Therefore,

[TABLE]

Since $\mathbf{\Pi}(L^{n})$ is orthogonal in the Euclidean inner product to the range of $E^{n-1,n}$,

[TABLE]

and hence $\mathbf{\Pi}(L^{n})^{T}(M^{n}\mathbf{y}^{n}-\mathbf{e}^{n})=0$. Therefore, $M^{n}\mathbf{y}^{n}-\mathbf{e}^{n}$ also belongs to the orthogonal complement of the null space of $\left(E^{n-1,n}\right)^{T}$.

∎

Since $\Delta^{n}_{nn}=1$, $M^{n}\mathbf{y}^{n}=\mathbf{e}^{n}$ if and only if $\widetilde{M}^{n}\Delta^{n}\mathbf{y}^{n}=\mathbf{e}^{n}$. Therefore, the coefficients of $u^{n}$ are given by

[TABLE]

We now find the coefficients of $v^{n+1}$. By Theorem 3.2, the coefficients are given by the last column of the inverse of any nonsingular Hankel extension of $\widetilde{M}^{n}$. We use this freedom to choose an extension such that the coefficients of $v^{n+1}$ are given by elevating $u^{n}$ by one degree in the monomial basis and then reversing the order of the entries. To describe this process, we introduce the following notation: given a vector $\mathbf{x}$ of length $(n+1)$, denote by $\mathbf{x}^{P}$ the vector of length $(n+2)$ given by

[TABLE]

where

[TABLE]

is the $(n+2)\times(n+2)$ exchange matrix.

Lemma 3.2.3**.**

If $H$ is an $(n+1)\times(n+1)$ nonsingular Hankel matrix of the form

[TABLE]

and $\mathbf{x}$ satisfies

[TABLE]

with $\mathbf{x}_{0}\neq 0$, then there exists an $(n+2)\times(n+2)$ Hankel extension $\widehat{H}$ of $H$ such that

[TABLE]

Proof.

Since $\mathbf{x}^{P}$ is the action of $\begin{pmatrix}\mathbf{x}&0\end{pmatrix}^{T}$ under the exchange matrix $P$, we can instead apply the exchange matrix to $\widehat{H}$ and show that there exist $\alpha$ and $\beta$ such that

[TABLE]

The first $n$ equations are satisfied by assumption. Therefore, choosing

[TABLE]

and

[TABLE]

gives the desired result.

∎

Therefore, the coefficients of $v^{n+1}$ are given by

[TABLE]

and hence Theorem 3.2 and (22) imply that

Theorem 3.3**.**

[TABLE]

Applying the identity $\left(M^{n}\right)^{-1}=\left(\Delta^{n}\right)^{-1}\left(\widetilde{M}^{n}\right)^{-1}\left(\Delta^{n}\right)^{-1}$ gives the result from Theorem 3.1.

3.1.3 Bernstein-Legendre conversion

We can use Theorem 2.1 to construct the spectral (or eigenvector) decomposition

[TABLE]

where $\Lambda^{n}$ contains the eigenvalues and $Q^{n}$ is an orthogonal matrix of eigenvectors. Since $M^{n}$ is symmetric and all of its eigenvalues are distinct, the eigenvectors are orthogonal. Therefore, all that remains is to normalize the eigenvectors in the 2-norm.

Proposition 3.3.1**.**

[TABLE]

Proof.

The result follows from (4) and (9):

[TABLE]

∎

Proposition 3.3.1 implies that

[TABLE]

and

[TABLE]

Equation (48) characterizes $Q^{n}$ and suggests an algorithm for its construction – begin with $\mathbf{\Pi}(L^{i})$ and elevate it to degree $n$. This requires $\mathcal{O}(n^{3})$ operations – there are $\mathcal{O}(n)$ columns and each one needs to be elevated $\mathcal{O}(n)$ times at $\mathcal{O}(n)$ operations per elevation. However, we can also adapt the classical 3-term recurrence for the Legendre polynomials to give a $\mathcal{O}(n^{2})$ process for constructing $Q^{n}$.

We begin the process by constructing the degree $n$ representation of $L^{0}$ – which is simply the vector of ones. We can similarly construct the coefficients for $L^{1}$. Then, we proceed inductively. Assuming that $L^{i}$ has been constructed (for $1\leq i<n$), the critical step in the recurrence is to multiply $L^{i}(x)$ by $x$. Given any polynomial of degree $n$

[TABLE]

we have

[TABLE]

where $\tilde{p}_{0}=0$ and $\tilde{p}_{i}=\tfrac{ip_{i-1}}{n+1}$ for $1\leq i\leq n+1$. This $\mathcal{O}(n)$ computes $xp(x)$ in the degree $n+1$ basis, but we need $xL^{i}(x)$ in the Bernstein basis of degree $n$. That is, if $\tilde{\mathbf{p}}$ holds the degree $n+1$ Bernstein coefficients, we need to find $\mathbf{q}$ such that

[TABLE]

which is a consistent system of $n+2$ equations in $n+1$ unknowns. Several $\mathcal{O}(n)$ algorithms for this are possible, but Gaussian elimination on the (tridiagonal) normal equations seems stable in practice. Algorithm 1 summarizes building the entire matrix $Q^{n}$ using the three-term recurrence and the tridiagonal degree reduction, where $E$ denotes the matrix $E^{n,n+1}$.

Consequently, $(M^{n})^{-1}=Q^{n}\left(\Lambda^{n}\right)^{-1}\left(Q^{n}\right)^{T}$ can be applied to a vector by multiplying by two dense orthogonal matrices and a diagonal scaling.

3.2 Decomposing the inverse

The proof of Theorem 3.3 suggests the following decomposition of $(M^{n})^{-1}$ into diagonal, Toeplitz, and Hankel matrices:

Theorem 3.4**.**

Let $T^{n}$ and $\widetilde{T}^{n}$ be the Toeplitz matrices given by

[TABLE]

and let $H^{n}$ and $\widetilde{H}^{n}$ be the Hankel matrices given by

[TABLE]

Then

[TABLE]

Proof.

By Theorem 3.3, we have that

[TABLE]

Therefore, $\left(\widetilde{M}^{n}\right)^{-1}=\widetilde{T}^{n}H^{n}-T^{n}\widetilde{H}^{n}$, and so the result follows from (15). ∎

This result implies a superfast $\mathcal{O}(n\log n)$ algorithm, but our numerical experiments suggest that it is highly unstable.

4 Applying the inverse

Now, we describe several approaches to applying $M^{-1}$ to a vector.

Cholesky factorization

$M^{n}$ is symmetric and positive definite and therefore admits a Cholesky factorization

[TABLE]

where $L$ is lower triangular with positive diagonal entries [17]. Widely available in libraries, computing $L$ requires $\mathcal{O}(n^{3})$ operations, and each of the subsequent triangular solves require $\mathcal{O}(n^{2})$ operations to perform. Our numerical results below suggest that it is one of the more stable and accurate methods under consideration, although our technique based on the eigendecomposition has similar accuracy and $\mathcal{O}(n^{2})$ complexity for the startup phase.

Exact inverse

In light of Theorem 3.1, we can directly form $M^{-1}$. Since the formula for each entry requires a sum, forming the inverse requires $\mathcal{O}(n^{3})$ operations. The inverse matrix can then be applied to any vector using the standard algorithm and in $\mathcal{O}(n^{2})$ operations. This has the same startup and per-solve complexity as the Cholesky factorization, although the constants are different.

Spectral decomposition

In Section 3.1.3, we showed how to compute the eigendecomposition of $M^{n}$ and hence express its inverse as

[TABLE]

The inverse can be applied by two (dense) matrix multiplications and a diagonal scaling, requiring $\mathcal{O}(n^{2})$ operations. Thanks to Algorithm 1, we have only an $\mathcal{O}(n^{2})$ startup phase.

DFT-based application

Theorem 3.4 implies that we can multiply by the inverse of $M^{n}$ in $\mathcal{O}(n\log n)$ operations. We invert $\Delta^{n}$ onto a given vector, and then all of the Toeplitz and Hankel matrices can be applied via circulant embedding before inverting $\Delta^{n}$ onto the result. Moreover, the operations can be fused together so that some FFTs are re-used and FFT/inverse FFT pairs cancel in (53). However, our numerical results reveal this approach to be quite unstable in practice, becoming wildly inaccurate long before one could hope to win from the super-fast algorithm. Therefore, we do not go into much detail on this approach.

5 Conditioning and accuracy

As a corollary of Theorem 2.1, $M^{n}$ is terribly ill-conditioned as $n$ increases. For each $n$, the minimal eigenvalue occurs when $i=n$, which is

[TABLE]

Similarly, the maximal eigenvalue occurs when $i=0$:

[TABLE]

Given these extremal eigenvalues, the 2-norm condition number is

[TABLE]

While this conditioning seems spectacularly bad, in practice, the Bernstein basis seems to give entirely satisfactory results at moderately high orders of approximation [11]. Here, we hope to give at least a partial explanation of this phenomenon. Recall that solving $M^{n}\mathbf{c}=\mathbf{b}$ exactly yields the Bernstein expansion coefficients of the best $L^{2}$ polynomial approximation to a function $f$ whose moments against the Bernstein basis are contained in $\mathbf{b}$. Consequently, if our solution process yields some $\mathbf{c}+\delta\mathbf{c}$ instead of $\mathbf{c}$, the $L^{2}$ norm of the polynomial encoded by $\delta\mathbf{c}$ has more direct relevance than the size of $\mathbf{c}$ in the Euclidean norm. On the other hand, the perturbation $\delta\mathbf{b}$ exists as an array of numbers (say, the roundoff error in computing the moments of $f$ by numerical integration), and we continue to use the Euclidean norm here.

In standard floating point analysis, suppose we have computed the solution to a perturbed solution

[TABLE]

and while a backward stable algorithm yields a small $\delta\mathbf{b}$, the perturbation in the solution $\delta\mathbf{c}$ might itself be significant. Equivalently, the perturbation $\delta\mathbf{c}$ satisfies

[TABLE]

While classical error perturbation and conditioning analysis would estimate $\|\delta\mathbf{c}\|$ in the Euclidean norm, we wish to consider $\|\delta\mathbf{c}\|_{M^{n}}$, which is the $L^{2}$ norm of the perturbation in the computed polynomial.

Taking the inner product of (60) with $\delta\mathbf{c}$, we have:

[TABLE]

so that

[TABLE]

In other words, for a perturbation $\delta\mathbf{b}$ of fixed 2-norm, the perturbation $\delta\mathbf{c}$ is much smaller measured in the $M^{n}$ norm than in the 2-norm – the amplification factor is only the square root as large.

This discussion suggests the following matrix norm. We think of $M^{n}$ as an operator mapping $\mathbb{R}^{n+1}$ onto itself. However, we equip the domain with the $M^{n}$-weighted norm and the range with the Euclidean norm. Then, we define the operator norm

[TABLE]

and going the opposite direction,

[TABLE]

We will also need the matrix norm

[TABLE]

using the $M^{n}$-weight in both the domain and range.

In light of (64), we can interpret (62) as

Proposition 5.0.1**.**

The norm of $(M^{n})^{-1}$ satisfies

[TABLE]

We can, by the spectral decomposition, give a similar bound for $M^{n}$ in the $M\rightarrow 2$ norm:

Proposition 5.0.2**.**

The norm of $M^{n}$ satisfies

[TABLE]

Proof.

Since $M^{n}$ is symmetric and positive-definite, it has a well-defined positive square root via the spectral decomposition. For $\mathbf{c}\in\mathbb{R}^{n+1}$, we write

[TABLE]

Now, we can characterize the weighted norm of the square root matrix via the spectral decomposition $M^{n}=Q^{n}\Lambda^{n}(Q^{n})^{T}$.

Note that for any $\mathbf{c}\in\mathbb{R}^{n+1}$, we have its $M^{n}$ norm as

[TABLE]

Letting $\mathbf{d}=(Q^{n})^{T}\mathbf{c}$,

[TABLE]

Now, we consider the $M^{n}$ norm of $\sqrt{M^{n}}$. Again, for any $\mathbf{c}\in\mathbb{R}^{n+1}$, we have

[TABLE]

so that with $\mathbf{d}=(Q^{n})^{T}\mathbf{c}$, we have

[TABLE]

Consequently,

[TABLE]

∎

Theorem 5.1**.**

The condition number of solving $M^{n}\mathbf{c}=\mathbf{b}$ measuring $\mathbf{c}$ in the $M^{n}$-norm and $\mathbf{b}$ in the 2-norm satisfies

[TABLE]

A plot of the condition numbers up to degree 20 in both norms is shown in Figure 1. This discussion indicates that, when measuring the relevant $L^{2}$ rather than the Euclidean norm of the solution process, we can expect much better results than the alarming condition number in (58) suggests. It is important to note that nothing in this discussion, other than the eigenvalues of $M^{n}$, is particular to the univariate mass matrix. Consequently, this discussion can inform the accuracy of the multivariate mass inversion process in [12] as well as the preconditioners for global mass matrices given in [3].

6 Numerical results

Now, we consider the accuracy of the methods described above on a range of problems. We consider the best $L^{2}$ approximation of two smooth functions, $f(x)=1/(1+396(x-0.5)^{2})$ and $f(x)=0.01+x/(x^{2}+1)$ – see Figure 2 for plots of these. Both functions are smooth. However, the first function has large derivatives and so requires high order of approximation to obtain small error. The second is not symmetric about $x=0.5$ but is otherwise relatively easy to approximate with a polynomial.

Solving (1) accurately, the best polynomial approximation of degree $n$ gives an $L^{2}$ error decreasing exponentially as a function of $n$. To demonstrate this, Figures 3(a) and 4(a) show the $L^{2}$ error between the two functions and the $L^{2}$ projection computed with each of the four methods described above. However, at various degrees, each of the methods deviate from yielding exponential convergence as roundoff error accumulates. We note that the DFT-based method seems to be the worst, followed by multiplication by the exactly computed inverse.

We also compute the best approximation using the Legendre polynomials (which only involves numerical integration) and compute the $L^{2}$ difference between that and the polynomial produced by each solution technique. Equivalently, this is the relative $M^{n}$-norm error between the exact and computed solution to (1). These plots appear in Figures 3(b) and 4(b), futher demonstrating the instability of the DFT-based approach.

Differences between the three more stable methods begin to appear when we consider the Euclidean norm of the error (Figures 3(c) and 4(c)) and residual (Figures 3(d) and 4(d)) for each of our solution methods. The DFT-based approach again performs much worse than the other methods, although the exact inverse gives 2-norm error comparable to the spectral and Cholesky approaches. The residual plots show that the latter two methods give very small residual errors, not much above machine precision. This strong backward stability of the Cholesky decomposition is a known property [8], and we suspect (but have not proven) that it holds for our spectral decomposition as well.

To illustrate that these properties of the solution algorithm do not seem to depend on approximating a smooth solution, we also chose random solution vectors, computed the right-hand side by matrix multiplication, and then applied our four solution algorithms to attempt to recover the solution. In Figure 5, we see exactly the same behavior – the DFT-based method behaving very badly, the Cholesky and spectral methods seemingly backward stable, and the exact matrix inverse somewhere in between.

These numerical results highlight several points. First, although finite-dimensional norms are equivalent, the bounding constants can dramatically vary as a function of the matrix size. Consequently, we can see comparable Euclidean norm errors but very different $M^{n}$-norm errors for, say the matrix inverse and Cholesky methods. Second, the ill conditioning of our method is real, although we in fact see that the $M^{n}$ norm errors (Figures 3(b) and 4(b)) in our solution process for the (at least empirically) backward-stable methods are in fact quite a bit lower than the 2-norm errors as predicted by the conditioning analysis in Section 5. It seems that the empirical performance of the spectral method and Cholesky factorization are the best, although they have slightly different associated costs. The spectral decomposition can be computed more quickly ($\mathcal{O}(n^{2})$ versus $\mathcal{O}(n^{3})$ for Cholesky), but the per-solve costs is higher – two dense matrix multiplications rather than two triangular solves.

7 Conclusions

We have studied several algorithms for the inversion of the univariate Bernstein mass matrix. Our fast algorithm for inversion based on the spectral decomposition seems very stable in practice and has accuracy comparable to the Cholesky decomposition, while a superfast algorithm based on the discrete Fourier transform is unfortunately unstable. Moreover, we have given a new perspective on the conditioning of matrices for polynomial projection indicating that the problems are better-conditioned with respect to the $L^{2}$ norm of the output than the Euclidean norm. In the future, we hope to expand this perspective to other polynomial problems and continue the development of fast and accurate methods for problems involving Bernstein polynomials.