# Recursive construction of biorthogonal polynomials for handling polynomial
regression

Laura Rebollo-Neira<EMAIL_ADDRESS>Department of Applied
Mathematics and Data Science, College of Engineering and Physical Sciences,
Aston University, B4 7ET, Birmingham, UK  Jason Laurie<EMAIL_ADDRESS>Department of Applied Mathematics and Data Science, College of Engineering and
Physical Sciences, Aston University, B4 7ET, Birmingham, UK

###### Abstract

An adaptive procedure for constructing a series of biorthogonal polynomials to
a basis of monomials spanning the same finite-dimensional inner product space
is proposed. By taking advantage of the orthogonality of the original basis,
our procedure circumvents the well-known instability problem arising from the
matrix inversion involved in classical polynomial regression. Moreover, the
recurrent generation of biorthogonal polynomials in our framework facilitates
the upgrading of all polynomials to include one additional element in the set
whilst also allowing for a natural downgrading of the polynomial regression
approximation. This is achieved by the posterior removal of any basis element
leading to a straightforward approach for reducing the approximation order. We
illustrate the usefulness of this approach through a series of examples where
we derive the resulting biorthogonal polynomials from Legendre, Laguerre, and
Chebyshev orthogonal bases.

Biorthogonal polynomials, polynomial regression, biorthogonal representation
of orthogonal Projections

##  I Introduction

Polynomial regression is becoming increasingly important for predictive
modelling [1, 2, 3], data analysis [4], signal and images processing [5],
machine learning and neuron networks [7, 6, 8]. Thus, a renewed interest in
efficient and effective techniques suitable for tackling common polynomial
regression problems has arisen. Fundamentally, polynomial regression problems
concern the modelling of a process represented by a function $f$ by a
$k$-order polynomial $f_{k}$, i.e.

$\displaystyle f_{k}(x)=c_{0}+c_{1}x+\cdots+c_{k}x^{k},\quad\text{for
$x\in{\cal{I}}\subseteq\mathbb{R}$},$

and finding the coefficients $c_{n}$ for $n=0,\dots,k$ by minimising the
$L^{2}$-norm of the approximation error $\|f-f_{k}\|_{2}$. Most applications
using polynomial regression, e.g. such as the method of moments used in Ref.
[4], approach the problem through the traditional method of solving the least
squares equation via the normal equations. However, such approaches require
the inversion of a Gram matrix of Hankel-type which are notorious for being
very badly conditioned [9], and can only provide a reliable solution for low-
order polynomial approximations. This is unfortunate, because the Weierstrass
approximation theorem states that every continuous function defined on a
compact interval can be uniformly approximated by a polynomial up to any
desired precision, even if the order of the polynomial might be high.
Moreover, because monomials on a compact interval are among the simplest
functions to evaluate, integrate, and differentiate, polynomial approximations
have enormous practical relevance in terms of both analytical and numerical
computations of integrals and derivatives. This article introduces an adaptive
biorthogonalisation technique which, by leveraging the orthogonality of
classic orthogonal polynomials (OP), avoids the need for a matrix inversion
and provides the coefficients of the polynomial approximation by simple
computation of inner products. Additional advantages stem from the recursive
nature of the approach that provides an efficient strategy for adapting the
order (or accuracy) of the polynomial approximation i.e., it is convenient to
iteratively increase or decrease the polynomial order to easily adjust the
residual error of the approximation.

The general theory of biorthogonal polynomials sequences (BPS) is not as
extensively developed as the theory of orthogonal polynomials sequences (OPS).
Most of the initial research on the topic was conducted in the 1970s and 1980s
with the contributions of Szegiö [10], and Iserle and Nørsett [11, 12, 13]. In
particular, the main focus being on the existence conditions, zeros, and the
Rodriguez-type representation of BPS. However, since then, there have been few
applications or further studies on BPS. We believe that their limited current
applicability might be partly due to the absence of convenient methodologies
for their practical construction. Therefore, our objective is to innovate new
methodologies for building BPS, enhancing comprehension, and expanding the
range of potential applications across diverse domains.

Since several definitions of biorthogonal polynomials (BP) have been used, let
us specify the definition we refer to here: Two polynomial bases
$\\{p_{n}\\}_{n=0}^{k}$ and $\\{q_{n}\\}_{n=0}^{k}$ are biorthogonal on an
interval ${\cal{I}}\subseteq\mathbb{R}$ if the inner product

$\displaystyle\left\langle
p_{n},q_{m}\right\rangle=\int_{\mathcal{I}}\,p_{n}(x)q_{m}(x)\,d\mu=\delta_{m,n},$

where $d\mu=w(x)dx$ is a measure with $w(x)$ an acceptable weight function,
and $\delta_{m,n}$ is the Kronecker delta-function. More precisely, we refer
to the classic-like biorthogonality by the set of polynomials
$\left\\{\beta_{n}^{k}\right\\}_{n=0}^{k}$ being biorthogonal to the set of
monomials $\\{x^{n}\\}_{n=0}^{k}$ with the
${\mbox{\rm{Span}}}\left(\\{x^{n}\\}_{n=0}^{k}\right)={\mbox{\rm{Span}}}\left(\\{\beta_{n}^{k}\\}_{n=0}^{k}\right)$
for each $k\in\mathbb{N}$. Here, the superscript of $\beta_{n}^{k}$ indicates
that all the elements $n=0,\dots,k$ must be re-calculated if the number $k$ of
elements in the basis changes. This is a major difference with OP for which if
the order $k$ is increased, say to $k+1$, only a single new independent
element needs to be added to the original orthogonal set, without modification
of the previously calculated elements.

While recursive equations for generating biorthogonal functions were proposed
in Refs. [14, 15, 16], the general nature of that framework misses the main
advantage of dealing with the particular case of monomials, namely, the wealth
of available results concerning the construction of OP. The central
contribution of this work is to produce an adaptive methodology for
constructing classic-like BP giving rise to representations of orthogonal
projections for the subspace spanned by the monomials. In particular, we
illustrate the method by deriving BP based on orthogonal bases built from
Legendre, Laguerre, and Chebyshev polynomials.

The paper is organised as follows: Sec. II provides the necessary definitions
and background material for the construction of classic-like BP. Sec. III
proves the biorthogonal property of the proposed construction as well as their
convenient adaption for upgrading and downgrading orthogonal projections.
Examples involving biorthogonality with respect to particular weight functions
and domains are also produced in this section. The conclusions and final
remarks are presented in Sec. IV.

##  II Preliminary Considerations

In this section we introduce the necessary background for a self-contained
derivation of the proposed construction of BP. Since our purpose is to develop
a methodology suitable for practical applications, we consider only
biorthogonal representations of orthogonal projections, and therefore,
throughout the paper we restrict ourselves to consider only finite-dimensional
inner product spaces, which we denote as $V({\cal{I}},w)$ to indicate that the
inner product $\left\langle f,g\right\rangle$ between two real functions in
$V({\cal{I}},w)$ is defined as

$\displaystyle\left\langle f,g\right\rangle=\int_{{\cal{I}}}f(x)g(x)w(x)\,dx,$

where $w(x)$ is a weight function which satisfies $w(x)\geq 0$ for all
$x\in{\cal{I}}$. For $k\in\mathbb{N}$, the subspace
$V_{k}({\cal{I}},w)\subseteq V({\cal{I}},w)$ is specifically defined as the
vector space spanned by the set of monomials $\\{x^{n}\\}_{n=0}^{k}$, i.e.

$\displaystyle
V_{k}={\mbox{\rm{Span}}}\left(\\{x^{0},x,\ldots,x^{k}\\}\right),\quad\text{for
$x\in{\cal{I}}$}.$

#### Type A and B orthonormal polynomial sequences

We introduce two definitions for the construction of orthonormal polynomial
sequences (OPS). These two types, type A and type B, are general forms based
on formulas used to construct sequences of orthogonal polynomials. Before, we
define the two types, we highlight a clarification of our notation to the
reader. Unless otherwise specified, throughout the paper, superscripts are
used as labels. The only exception is for the variable $x$ where the
superscript indicates the variable $x$ raised to a power
$i\in\mathbb{N}\cup\\{0\\}$, i.e., $x^{i}$.

######  Definition 1 (Type A Orthonormal Polynomial Sequence).

The sequence of functions generated by

$\displaystyle p_{n}(x)=\sum_{i=0}^{n}a_{i}^{n}x^{i},\quad\text{with
$a_{n}^{n}\neq 0$}$ (1)

such that $\left\langle p_{n},p_{m}\right\rangle=\delta_{n,m}$ is called a
Type A orthonormal polynomial sequence.

For the convenience of some particular cases we differentiate the above
definition with a second type, type B of OPS.

######  Definition 2 (Type B Orthonormal Polynomial Sequence).

The sequence of functions generated by

$\displaystyle p_{n}(x)=\sum_{i=0}^{\lfloor
n/2\rfloor}a_{i}^{n}x^{n-2i},\quad\text{with $a_{0}^{n}\neq 0$},$ (2)

with $\lfloor n/2\rfloor$ indicating the integer part of $n/2$, i.e. the
greatest integer less than or equal to $n/2$, that satisfy $\left\langle
p_{n},p_{m}\right\rangle=\delta_{n,m}$, is called a Type B orthonormal
polynomial sequence.

The polynomials $\left\\{p_{n}(x)\right\\}_{n=0}^{k}$ of Type A or B, form an
orthonormal set that the span the same subspace as the set of monomials
$\\{x^{n}\\}_{n=0}^{k}$ with respect to the inner product in $V({\cal{I}},w)$
so that

$\displaystyle
V_{k}={\mbox{\rm{Span}}}\left(\\{x^{n}\\}_{n=0}^{k}\right)={\mbox{\rm{Span}}}\left(\\{p_{n}\\}_{n=0}^{k}\right)\subseteq
V({\cal{I}},w).$

######  Property 1.

If $n\leq m$, then an OPS satisfies

$\displaystyle\left\langle
x^{n},p_{m}\right\rangle=\begin{cases}\displaystyle\frac{\delta_{n,m}}{a_{m}^{m}}&\text{if
$p_{m}$ is Type A,}\\\ \displaystyle\frac{\delta_{n,m}}{a_{0}^{m}}&\text{if
$p_{m}$ is Type B}.\end{cases}$ (3)

###### Proof.

Since $x^{n}\in{\mbox{\rm{Span}}}\left(\\{p_{i}\\}_{i=0}^{n}\right)$ it can be
expressed in terms of the first $n$ orthonormal polynomials as

$\displaystyle x^{n}=\sum_{i=0}^{n}d_{i}^{n}p_{i}(x),$

where $d_{i}^{n}\in\mathbb{R}$ for $i=0,\dots,n$. Hence, for any $m$
satisfying $n<m$ we have that $\left\langle x^{n},p_{m}\right\rangle=0$ for
both Type A and Type B OPS. For $n=m$, if $p_{m}$ is of Type A, using
definition (1), we can express

$\displaystyle a_{m}x^{m}=p_{m}-\sum_{i=0}^{m-1}a_{i}^{m}x^{i}$ (4)

so that, from the previous case, it follows that $\left\langle
x^{m},p_{m}\right\rangle=1/{a_{m}^{m}}$.

Similarly, if $p_{m}$ is of Type B, with $n=m$, then using definition (2) we
have

$\displaystyle a_{0}^{m}x^{m}=p_{m}-\sum_{i=0}^{\lfloor
m/2\rfloor}a_{i}^{m}x^{m-2i},$ (5)

and subsequently, $\left\langle x^{m},p_{m}\right\rangle=1/{a_{0}^{m}}$. This
proves Property 1. ∎

#### Orthogonal Projections

In what follows, it is important to define some properties of subspaces that
allows us to use results regarding orthogonal projections.

######  Definition 3 (Direct sum of subspaces).

A vector space $V$ is said to be the direct sum of two subspaces, say $V_{k}$
and $U$, written as

$\displaystyle V=V_{k}\oplus U$

if each $f\in V$ has a unique representation

$\displaystyle f=f_{k}+g,\quad\text{where $f_{k}\in V_{k}$ and $g\in U$}.$

If $V$ is an inner product space $V({\cal{I}},w)$ a case of much interest
concerns complementary subspaces which are orthogonal.

######  Definition 4 (Orthogonal complement subspace).

A vector space $U$ is said to be the orthogonal complement of $V_{k}$ in
$V({\cal{I}},w)$ if

$\displaystyle U=V_{k}^{\perp}=\left\\{g\in V\ :\ \left\langle
g,f_{k}\right\rangle=0,\forall f_{k}\in V_{k}\right\\}.$ (6)

######  Definition 5.

(Orthogonal projector) Let two vectors spaces $V_{k}$ and $V_{k}^{\perp}$ be
orthogonal complements of each other. Then, given any $f\in V$, the map
$\hat{\operatorname{P}}_{V_{k}}:V\to V_{k}$ defined as

$\displaystyle\hat{\operatorname{P}}_{V_{k}}f=f_{k}\in V_{k},$

is called the orthogonal projection (or orthogonal projector) of $V$ onto
$V_{k}$.

######  Property 2.

The orthogonal projector $\hat{\operatorname{P}}_{V_{k}}:V\to V_{k}$ has the
following properties:

  * (i)

$\hat{\operatorname{P}}_{V_{k}}$ is idempotent, i.e.
$\hat{\operatorname{P}}^{2}_{V_{k}}=\hat{\operatorname{P}}_{V_{k}}$.

  * (ii)

If $V=V_{k}\oplus V_{k}^{\perp}$ then
$\hat{\operatorname{P}}_{V_{k}}f_{k}=f_{k}$ if $f_{k}\in V_{k}$ and
$\hat{\operatorname{P}}_{V_{k}}g=0$ if $g\in V_{k}^{\perp}.$

As a consequence of Property 2(ii) we can express
$V_{k}^{\perp}=\hat{\operatorname{I}}_{V}V-\hat{\operatorname{P}}_{V_{k}}V=(\hat{\operatorname{I}}_{V}-\hat{\operatorname{P}}_{V_{k}})V=\hat{\operatorname{P}}_{V_{k}^{\perp}}V$,
where $\hat{\operatorname{I}}_{V}$ is the identity operator in $V$.
Subsequently, it follows that every $f\in V$ can be decomposed in terms of two
contributions
$f=\hat{\operatorname{P}}_{V_{k}}f+\hat{\operatorname{P}}_{V_{k}^{\perp}}f$.

######  Property 3 (Self-adjoint operator).

It readily follows that $\hat{\operatorname{P}}_{V_{k}}$ satisfies

$\displaystyle\left\langle
f,\hat{\operatorname{P}}_{V_{k}}g\right\rangle=\left\langle\hat{\operatorname{P}}_{V_{k}}f,g\right\rangle.$

An operator fulfilling the above property is said to be self-adjoint.

The next two theorems are well-known fundamental results concerning the
optimal approximation in a finite-dimensional inner product subspace, also
known as the least squares approximation. For pedagogical reasons, we present
a proof tailored to our framework.

######  Theorem 1.

Let $V({\cal{I}},w)$ be an inner product space and $V_{k}$ a finite-
dimensional subspace of $V({\cal{I}},w)$. Then, the least squares
approximation in $V_{k}$ of any $f\in V$ can be obtained by the projection of
$f$ onto the subspace $V_{k}$, i.e. $\hat{\operatorname{P}}_{V_{k}}f$.

###### Proof.

Recall that the optimal approximation that minimises the $L^{2}$-norm in
$V_{k}$ for any $f\in V$ satisfies

$\displaystyle\|f-f^{\ast}\|_{2}=\min_{f_{k}\in V_{k}}\|f-f_{k}\|_{2}.$

For any $f_{k}\in V_{k}$ we can write that
$f=f+\hat{\operatorname{P}}_{V_{k}}f-\hat{\operatorname{P}}_{V_{k}}f$ to
express $\|f-f_{k}\|_{2}^{2}$ in a convenient form

$\displaystyle\|f-f_{k}\|_{2}^{2}$ $\displaystyle=\left\langle
f-f_{k},f-f_{k}\right\rangle=\left\langle
f-f_{k}+\hat{\operatorname{P}}_{V_{k}}f-\hat{\operatorname{P}}_{V_{k}}f,f-f_{k}+\hat{\operatorname{P}}_{V_{k}}f-\hat{\operatorname{P}}_{V_{k}}f\right\rangle.$

Since $\left\langle
f_{k},f-\hat{\operatorname{P}}_{V_{k}}f\right\rangle=\left\langle
f_{k},\hat{\operatorname{P}}_{V_{k}^{\perp}}f\right\rangle=0$ we have

$\displaystyle\|f-f_{k}\|_{2}^{2}=\left\langle
f-\hat{\operatorname{P}}_{V_{k}f},f-\hat{\operatorname{P}}_{V_{k}}f\right\rangle+\left\langle\hat{\operatorname{P}}_{V_{k}}f-f_{k},\hat{\operatorname{P}}_{V_{k}}f-f_{k}\right\rangle.$

i.e.

$\displaystyle\|f-f_{k}\|_{2}^{2}=\|f-\hat{\operatorname{P}}_{V_{k}}f\|_{2}^{2}+\|\hat{\operatorname{P}}_{V_{k}}f-f_{k}\|_{2}^{2}.$

Hence, $\|f-f_{k}\|_{2}$ is minimised when
$f^{\ast}=f_{k}=\hat{\operatorname{P}}_{V_{k}}f$. ∎

The next theorem prescribes a method of constructing an orthogonal projector
onto the finite-dimensional subspace
$V_{k}={\mbox{\rm{Span}}}\left(\\{x^{n}\\}_{n=0}^{k}\right)$, with respect to
the inner product defined on $V({\cal{I}},w)$.

######  Theorem 2.

Let $\\{p_{n}\\}_{n=0}^{k}$ be an orthonormal basis for the inner product
subspace $V_{k}({\cal{I}},w)$. Then the orthogonal projector operator
$\hat{\operatorname{P}}_{V_{k}}:V\to V_{k}$ can be expressed as

$\displaystyle\hat{\operatorname{P}}_{V_{k}}=\sum_{i=0}^{k}p_{i}\left\langle\cdot,p_{i}\right\rangle,$
(7)

where $\left\langle\cdot,p_{i}\right\rangle$ performs an inner product with
the basis polynomial $p_{i}$, i.e. for any $f\in V$

$\displaystyle\hat{\operatorname{P}}_{V_{k}}f=\sum_{i=0}^{k}p_{i}\left\langle
f,p_{i}\right\rangle.$

###### Proof.

We need to show that $\hat{\operatorname{P}}_{V_{k}}$ as defined in Eq. (7)
fulfils

  * (i)

$\hat{\operatorname{P}}_{V_{k}}^{2}=\hat{\operatorname{P}}_{V_{k}}$.

  * (ii)

$\hat{\operatorname{P}}_{V_{k}}f_{k}=f_{k}$ for all $f_{k}\in V_{k}$ and
$\hat{\operatorname{P}}_{V_{k}}f_{k}^{\perp}=0$ for all $f_{k}^{\perp}\in
V_{k}^{\perp}$.

Since $\left\langle p_{j},p_{i}\right\rangle=\delta_{j,i}$ it follows from Eq.
(7), that

$\displaystyle\hat{\operatorname{P}}_{V_{k}}^{2}=\hat{\operatorname{P}}_{V_{k}}\hat{\operatorname{P}}_{V_{k}}=\sum_{j=0}^{k}\sum_{i=0}^{k}p_{i}\left\langle
p_{j},p_{i}\right\rangle\left\langle\cdot,p_{j}\right\rangle=\sum_{j=0}^{k}p_{j}\left\langle\cdot,p_{j}\right\rangle=\hat{\operatorname{P}}_{V_{k}},$

which proves (i).

For (ii) we use $\\{p_{n}\\}_{n=0}^{k}$ to write an arbitrary vector in
$V_{k}$ as the linear combination
$\displaystyle{f_{k}=\sum_{n=0}^{k}s_{n}p_{n}}$, for some real coefficients
$s_{n}$ for $n=0,\dots,k$. Then

$\displaystyle\hat{\operatorname{P}}_{V_{k}}f_{k}=\sum_{n=0}^{k}\sum_{i=0}^{n}p_{i}\left\langle
s_{n}p_{n},p_{i}\right\rangle=\sum_{n=0}^{k}\sum_{i=0}^{n}s_{n}p_{i}\delta_{n,i}=\sum_{n=0}^{k}s_{n}p_{n}=f_{k}.$

Since, $f_{k}^{\perp}\in V_{k}^{\perp}$ we have that $\left\langle
f_{k}^{\perp},p_{n}\right\rangle=0$ for $n=0,\dots,k$, and it follows that
$\sum_{n=0}^{k}p_{n}\left\langle f_{k}^{\perp},p_{n}\right\rangle=0$.
Subsequently, we conclude that for $f\in V({\cal{I}},w)$ the operation
$\hat{\operatorname{P}}_{V_{k}}f=\sum_{n=0}^{k}p_{n}\left\langle
f,p_{n}\right\rangle$ is the orthogonal projection of $f$ onto
$V_{k}({\cal{I}},w)$. ∎

We are now ready to construct BPS from given OPS of either Type A or B.

##  III Biorthogonal representation of orthogonal projectors

The next theorem establishes two distinctive features of the proposed BPS. (i)
Their construction readily follows from known OPS of either Type A or B. (ii)
The BPS provide us with the biorthogonal representation for the orthogonal
projection operator $\hat{\operatorname{P}}_{V_{k}}$.

######  Theorem 3.

Let $\\{p_{n}\\}_{n=0}^{k}$ be a finite OPS. Then the orthogonal projector
operator $\hat{\operatorname{P}}_{V_{k}}$ can be represented as

$\displaystyle\hat{\operatorname{P}}_{V_{k}}=\sum_{i=0}^{k}x^{i}\left\langle\cdot,\beta_{i}^{k}\right\rangle,$
(8)

where the functions $\beta_{n}^{k}(x)$ are defined as

  * (i)

$\displaystyle\beta_{n}^{k}(x)=\sum_{j=n}^{k}a_{n}^{j}\
p_{j}(x),\,x\in{\cal{I}},$ (9)

if $\\{p_{n}\\}_{n=0}^{k}$ is an OPS of Type A as defined in Eq. (1), or

  * (ii)

$\displaystyle\beta_{n}^{k}(x)=\sum_{i=0}^{\lfloor(k-n)/2\rfloor}a_{i}^{2i+n}\
p_{2i+n}(x),$ (10)

if $\\{p_{n}\\}_{n=0}^{k}$ is an OPS of Type B as defined in Eq. (2).

###### Proof.

For $n=k=0$ we have $\beta_{0}^{0}=p_{0}$ and $x^{0}=p_{0}$. Hence, operator
$\hat{\operatorname{P}}_{V_{0}}$ as given by Eq. (8) for $k=0$ is the
orthogonal projection operator onto $V_{0}$. Assuming that
$\hat{\operatorname{P}}_{V_{k}}$ in Eq. (8) is the orthogonal projection
operator onto $V_{k}$ we will show that then
$\hat{\operatorname{P}}_{V_{k+1}}$ is the orthogonal projection operator onto
$V_{k+1}$.

Case (i). Since
${\displaystyle{\beta_{n}^{k+1}=\beta_{n}^{k}+a_{n}^{k+1}p_{k+1}}}$, we have

$\displaystyle\sum_{n=0}^{k+1}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle$
$\displaystyle=\sum_{n=0}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle+x^{k+1}\left\langle\cdot,\beta_{k+1}^{k+1}\right\rangle$
$\displaystyle=\sum_{n=0}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k}\right\rangle+\sum_{n=0}^{k}x^{n}a_{n}^{k+1}\left\langle\cdot,p_{k+1}\right\rangle+x^{k+1}a_{k+1}^{k+1}\left\langle\cdot,p_{k+1}\right\rangle$
$\displaystyle=\hat{\operatorname{P}}_{V_{k}}+\sum_{n=0}^{k+1}x^{n}a_{n}^{k+1}\left\langle\cdot,p_{k+1}\right\rangle$
$\displaystyle=\hat{\operatorname{P}}_{V_{k}}+p_{k+1}\left\langle\cdot,p_{k+1}\right\rangle\quad{\text{using
\eqref{eq:OPS}}}$ $\displaystyle=\hat{\operatorname{P}}_{V_{k+1}}.$

Case (ii). From the definition of $\beta_{n}^{k}$ given by Eq. (12) we have
that

$\displaystyle\beta_{n}^{k+1}=\begin{cases}\displaystyle\beta_{n}^{k}+a^{k+1}_{(k-n+1)/2}\
p_{k+1}&\text{for $k-n+1$ even},\\\ \beta_{n}^{k}&\text{for $k-n+1$ odd},\\\
\beta_{k+1}^{k+1}=a^{k+1}_{0}\,p_{k+1}&\text{for $n=k+1$}.\end{cases}$

For $n\in\mathbb{N}$ let us define two sets of numbers

$\displaystyle{\cal{S}}_{e}$ $\displaystyle=\\{n\leq k\ :\ k-n+1\
\text{even}\\},\quad\text{and}\quad{\cal{S}}_{o}=\\{n\leq k\ :\ k-n+1\
\text{odd}\\}.$

Then,

$\displaystyle\sum_{n=0}^{k+1}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle$
$\displaystyle=\sum_{n=0}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle+x^{k+1}\left\langle\cdot,\beta_{k+1}^{k+1}\right\rangle=\sum_{n\in{\cal{S}}_{o}}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle+\sum_{n\in{\cal{S}}_{e}}x^{n}\left\langle\cdot,\beta_{n}^{k+1}\right\rangle+x^{k+1}\left\langle\cdot,\beta_{k+1}^{k+1}\right\rangle$
$\displaystyle=\sum_{n=0}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k}\right\rangle+\sum_{n\in{\cal{S}}_{e}}x^{n}\left\langle\cdot,a_{(k-n+1)/2}^{k+1}p_{k+1}\right\rangle+x^{k+1}\left\langle\cdot,a_{0}^{k+1}p_{k+1}\right\rangle$
$\displaystyle=\hat{\operatorname{P}}_{V_{k}}+\sum_{i=0}^{\lfloor
k/2\rfloor}x^{k+1-2i}a_{i}^{k+1}\left\langle\cdot,p_{k+1}\right\rangle=\hat{\operatorname{P}}_{V_{k}}+p_{k+1}\left\langle\cdot,p_{k+1}\right\rangle=\hat{\operatorname{P}}_{V_{k+1}}.$

∎

######  Remark III.1.

It is worth noticing that Property 3 implies that

$\displaystyle\hat{\operatorname{P}}_{V_{k}}=\sum_{i=0}^{k}x^{i}\left\langle\cdot,\beta_{i}^{k}\right\rangle=\sum_{i=0}^{k}\beta_{i}^{k}\left\langle\cdot,x^{i}\right\rangle.$

######  Property 4.

The set $\\{\beta^{k}_{n}\\}_{n=0}^{k}$ defined in (9) or in (12) is
biorthogonal to the set of monomials $\\{x^{n}\\}_{n=0}^{k}$ on the interval
$x\in{\cal{I}}\subseteq\mathbb{R}$, i.e.

$\displaystyle\left\langle\beta_{n}^{k},x^{m}\right\rangle=\delta_{n,m}.$

###### Proof.

We consider only case (i), in which we will prove the result by induction. The
proof of case (ii) follows a similar approach. For $n=k=0$ the property holds,
because $\beta^{0}_{0}=a_{0}^{0}p_{0}$, $x^{0}=p_{0}/a_{0}^{0}$, and
$\left\langle p_{0},p_{0}\right\rangle=1$. For the induction step within the
proof, we assume that

$\displaystyle\left\langle\beta_{n}^{k},x^{m}\right\rangle=\delta_{n,m},\quad\text{for
$n=0,\dots,k$ and $m=0,\dots,k$}.$ (11)

We must consider four specific cases separately, to show that

$\displaystyle\left\langle\beta_{n}^{k+1},x^{m}\right\rangle=\delta_{n,m},\quad\text{for
$n=0,\dots,k+1$ and $m=0,\dots,k+1$}.$

The four cases to consider are

  1. (i)

$n=0,\dots,k$ and $m=0,\dots,k$.

  2. (ii)

$n=k+1$ and $m=0,\dots,k$.

  3. (iii)

$n=0,\dots,k$ and $m=k+1$.

  4. (iv)

$n=k+1$ and $m=k+1$.

Case (i): Using (9) we write

$\displaystyle\left\langle\beta_{n}^{k+1},x^{m}\right\rangle=\left\langle\beta_{n}^{k},x^{m}\right\rangle+a_{n}\left\langle
p_{k+1},x^{m}\right\rangle,$

so that from the assumption
$\left\langle\beta_{n}^{k},x^{m}\right\rangle=\delta_{n,m}$, with
$n=0,\dots,k$ and $m=0,\dots,k$, and Eq. (3), it follows that
$\left\langle\beta_{n}^{k+1},x^{m}\right\rangle=\delta_{n,m}$ for
$n=0,\dots,k$ and $m=0,\dots,k$.

Case (ii):
$\left\langle\beta_{k+1}^{k+1},x^{m}\right\rangle=a_{k+1}^{k+1}\left\langle
p_{k+1},x^{m}\right\rangle=0$ from (3), since $m<k+1$.

Case (iii): Using Eqs. (9), (4), and (11) we have

$\displaystyle\left\langle\beta_{n}^{k+1},x^{k+1}\right\rangle$
$\displaystyle=\left\langle\beta_{n}^{k},x^{k+1}\right\rangle+a_{n}^{k+1}\left\langle
p_{k+1},x^{k+1}\right\rangle=\left\langle\beta_{n}^{k},\frac{p_{k+1}}{a_{k+1}^{k+1}}\right\rangle-\left\langle\beta_{n}^{k},\sum_{i=0}^{k}\frac{a_{i}^{k+1}}{a_{k+1}^{k+1}}x^{i}\right\rangle+\frac{a_{n}^{k+1}}{a_{k+1}^{k+1}}$
$\displaystyle=0-\frac{a_{n}^{k+1}}{a_{k+1}^{k+1}}+\frac{a_{n}^{k+1}}{a_{k+1}^{k+1}}=0.$

Case (iv) The result follows by setting $m=k+1$ in Case (ii):

$\displaystyle\left\langle\beta_{k+1}^{k+1},x^{k+1}\right\rangle=a_{k+1}^{k+1}\left\langle
p_{k+1},x^{k+1}\right\rangle=\frac{a_{k+1}}{a_{k+1}}=1.$

From the results of all four cases (i), (ii), (iii), and (iv) we can conclude
by mathematical induction that Eq. (11) holds. ∎

####  III.0.1 Application to downgrading orthogonal projectors

Let us suppose that from knowing the orthogonal projection of $f$,
$\hat{\operatorname{P}}_{V_{k}}f=f_{k}$, one is interested in producing the
orthogonal projection onto the subspace left by removing one element,
$x^{\ell_{j}}$, from the set $\\{x^{n}\\}_{n=0}^{k}$, i.e. the orthogonal
projection on the new subspace spanned by the set

$\displaystyle{\cal{D}}_{k\backslash{\ell_{j}}}=\\{x^{0},x^{1},\ldots,x^{\ell_{j}-1},x^{\ell_{j}+1},\ldots,x^{k}\\}.$

We use the notation
$V_{k\backslash{\ell_{j}}}={\mbox{\rm{Span}}}\left({\cal{D}}_{k\backslash{\ell_{j}}}\right)$
to indicate the subspace spanned by this set. The downgrading of the projector
$\hat{\operatorname{P}}_{V_{k}}\to\hat{\operatorname{P}}_{V_{k\backslash{\ell_{j}}}}$,
is simpler using BP than having to downgrade an OPS. It can be achieved
through the backward adaptive biorthogonalisation formula derived in Ref. [15]
and extended in Ref. [16]:

$\displaystyle\beta_{n}^{k\backslash\ell_{j}}(x)=\beta_{n}^{k}(x)-\beta_{\ell_{j}}^{k}(x)\frac{\left\langle\beta_{\ell_{j}}^{k},\beta_{n}^{k}\right\rangle}{\left\langle\beta_{\ell_{j}}^{k},\beta_{\ell_{j}}^{k}\right\rangle},\quad\text{for
$n=0,\dots,\ell_{j-1},\ell_{j+1},\dots,k$ and $x\in{\cal{I}}$}.$ (12)

######  Property 5.

Let the set $\\{\beta_{n}^{k}\\}_{n=0}^{k}$ satisfy the following two
properties:

  * (i)

is biorthogonal to the set $\\{x^{n}\\}_{n=0}^{k}$.

  * (ii)

provides a representation of $\hat{\operatorname{P}}_{V_{k}}$ of the form give
by Eq. (8),

then the set of functions $\\{\beta_{n}^{k\backslash\ell_{j}}\\}_{n=0}^{k}$
defined as in Eq. (12),

  * (iii)

is biorthogonal to the set ${\cal{D}}_{k\backslash{\ell_{j}}}$.

  * (iv)

provides a representation of
$\hat{\operatorname{P}}_{V_{k\backslash{\ell_{j}}}}$, the orthogonal
projection onto
$V_{k\backslash{\ell_{j}}}={\mbox{\rm{Span}}}\left(\\{{\cal{D}}_{k\backslash{\ell_{j}}}\\}\right)$.

###### Proof.

The proof is presented in Appendix A. ∎

######  Remark III.2.

Since the approximation error arising from removing the element $x^{\ell_{j}}$
from the projection of $f$ onto the subspace $V_{k}$ given by
$\hat{\operatorname{P}}_{V_{k}}f$ is (c.f. (IV))

$\displaystyle\hat{\operatorname{P}}_{V_{k}}f-\hat{\operatorname{P}}_{V_{k\backslash{\ell_{j}}}}f=\frac{\beta_{\ell_{j}}^{k}\left\langle
f,\beta_{\ell_{j}}^{k}\right\rangle}{\|\beta_{\ell_{j}}^{k}\|^{2}}.$

It follows that by removing the element $x^{\ell_{j}}$ such that

$\displaystyle\frac{\left|\left\langle
f,\beta_{\ell_{j}}^{k}\right\rangle\right|^{2}}{\|\beta_{\ell_{j}}^{k}\|^{2}}\,\quad\text{is
minimum}$

the remaining downgraded approximation minimises the new approximation error
norm.

###  III.1 Biorthogonal polynomials from Type A classic OPS

We illustrate our methodology by considering the particular families of
biorthorthogonal polynomials which are generated by classical OPS. In this
subsection we consider only those OPS of Type A (c.f. Eq. (1)). These include
OPS composed of

  * (i)

Legendre polynomials $\\{\tilde{P}_{n}(x)\\}_{n=0}^{k}$ defined on the
interval ${\cal{I}}=[0,b]$, with inner product weight $w(x)=1$.

  * (ii)

Laguerre polynomials $\\{L_{n}(x)\\}_{n=0}^{k}$ defined on the interval
${\cal{I}}=[0,\infty)$, with inner product weight $w(x)=\exp(-x)$.

####  III.1.1 Biorthogonal polynomials derived from Legendre polynomials

In this case the orthogonal polynomials $\\{p_{n}(x)\\}_{n=0}^{k}$ that
generate the biorthogonal polynomials Eq. (9) are the Legendre polynomials
$\\{\tilde{P}_{n}(x)\\}_{n=0}^{k}$ defined on the interval ${\cal{I}}=[0,b]$.
These polynomials on $[0,b]$ can be obtained by the formula:

$\displaystyle\tilde{P}_{j}(x)=(-1)^{j}\frac{\sqrt{2j+1}}{\sqrt{b}}\sum_{i=0}^{j}(-1)^{i}\binom{j}{i}\binom{j+i}{i}\frac{x^{i}}{b^{i}},$

which is in the exact form defined by being Type A. Here, we note that $b^{i}$
indicates the $i$-th power of $b$. For this particular case, the coefficients
$a_{i}^{j}$ that are required to obtain the functions $\beta_{n}^{k}$ as
defined in Eq. (9) are

$\displaystyle
a_{i}^{j}=(-1)^{i+j}\frac{\sqrt{2j+1}}{b^{i}\sqrt{b}}\binom{j}{i}\binom{j+i}{i}.$

Accordingly, the BP generated by Legendre polynomials on $[0,b]$ adopt the
form

$\displaystyle\beta_{n}^{k}$
$\displaystyle=\sum_{j=n}^{k}\frac{{2j+1}}{{b}}(-1)^{j+n}\binom{j}{n}\binom{j+n}{n}\sum_{i=0}^{j}(-1)^{i+j}\binom{j}{i}\binom{j+i}{i}\frac{x^{i}}{b^{n+i}}$
$\displaystyle=\frac{(-1)^{n}}{b^{n+1}}\sum_{j=n}^{k}\sum_{i=0}^{j}(-1)^{i}\frac{{2j+1}}{b^{i}}\binom{j}{n}\binom{j+n}{n}\binom{j}{i}\binom{j+i}{i}x^{i}.$

Subsequently, the coefficients $c_{n}^{k}$ calculated by

$\displaystyle c_{n}^{k}$
$\displaystyle=\int_{0}^{b}f(x)\beta_{n}^{k}(x)\,dx=\frac{(-1)^{n}}{b^{n+1}}\sum_{j=n}^{k}\sum_{i=0}^{j}\frac{{2j+1}}{b^{i}}(-1)^{i}\binom{j}{n}\binom{j+n}{n}\binom{j}{i}\binom{j+i}{i}\int_{0}^{b}f(x)x^{i}\,dx,$
(13)

provide the orthogonal projection of $f$ onto $V_{k}$ as linear combination of
simple monomials

$\displaystyle f_{k}(x)=\sum_{n=0}^{k}c_{n}^{k}x^{n},\quad x\in[0,b].$ (14)

When $f(x)$ is a probability density function the integrals in Eq. (13) are
the statistical moments of the probability distribution. On $V({\cal{I}},w)$
we define the $i$-th generalised moment of a function as

$\displaystyle\mu_{i}=\int_{{\cal{I}}}x^{i}f(x)\,w(x)\ dx.$ (15)

#### Example 1: Polynomial approximation of noisy data using Legendre
polynomials (Type A)

We demonstrate the application of BP outlined above, using Legendre orthogonal
polynomials $\\{\tilde{P}_{n}(x)\\}_{n=0}^{k}$ defined on the interval
$[0,b]$, for calculating the coefficients of a polynomial regression problem.
In our example, we consider numerical data simulated using $N=501$ uniformly
distributed points on the interval ${\cal{I}}=[0,1]$,
$x_{i}=0.002i,\,i=0,\dots,500$, generated by the chirp function
$f_{i}^{o}=f(x_{i})=\cos(7\pi x_{i}^{2})+\epsilon_{i}$, where $\epsilon_{i}$
are independent Gaussian random variables with mean zero and variance equal to
$0.01$. The simulated data is shown in Fig. 1(a) by the green dots. Our goal
is to approximate the functional form of these data points by a polynomial. As
established by Theorems 1 and 3, the coefficients $c_{n}^{k}$, for
$,n=0,\dots,k$ in

$\displaystyle f_{k}(x)=\sum_{n=0}^{k}c_{n}^{k}x^{n},$

give rise to the polynomial that minimises the $L^{2}$-norm of the
approximation error. These can be computed using numerical integration
methods, in our case, the fourth-order Simpson’s 1/3 rule.

The green solid line in Fig. 1(b) is the resulting approximation obtained by
fixing $k=17$. Fig. 1(c) displays the downgraded version of that
approximation, where downgrading was achieved by removing three terms (one by
one), as prescribed in Remark III.2. It is worth noticing that the removal of
the three terms, corresponding to $x,x^{4}$ and $x^{17}$, does not
significantly change the quality of the approximation. The $L^{2}$-norm of the
error in representing the chirp with the approximation given in Fig. 1(b) is
$4.55\times 10^{-2}$ and with the downgraded approximation of Fig. 1(c) is
$5.18\times 10^{-2}$. In order to stress the suitability of the proposed
approach for reducing the number of terms in our polynomial approximation, in
Fig. 1(d), we instead present a polynomial approximation which is obtained
with the same number of terms, but fixing $k=14$. In this case the
$L^{2}$-norm error increases to $1.76\times 10^{-1}$. With these error values
we have additionally calculated the Bayesian information criterion (BIC) score
[18] which is used to assess the model selection for independent Gaussian
random data, taking into account the number of parameters $\gamma$ required in
the approximation,

$\displaystyle\text{BIC}=\gamma\log(N)+N\log\left(\frac{1}{N}\sum_{i=0}^{N}(f_{i}^{o}-f_{k}(x_{i}))^{2}\right).$

One can interpret the BIC score as a logarithm of the efficiency between the
model accuracy and its complexity, where a lower value represents a better
approximation. In the above definition we have, $\gamma=17$ for approximation
(a) and $\gamma=14$ for approximations (b) and (c). Even if the best BIC value
amongst the three approximations is produced by model (b), model (c) with
three fewer parameters would increase of the BIC value by $3.45\%$. Adopting
model (d) with also three fewer parameters would increase it by $45.30\%$.
Since approximations with smaller BIC are in general preferred, these scores
highlight the advantage of reducing the amount of terms in the approximation
by removing the three selected terms, rather than reducing the polynomial
order.

Figure 1: Polynomial approximations to the simulated noisy chirp data
represented as green points displayed in (a). The black dashed line in all
fourth plots represent the original (noiseless) chirp function $f(x)=\cos(7\pi
x^{2})$ on the interval $[0,1]$. (b) The green line corresponds to the
polynomial approximation using Legendre polynomials corresponding to $k=17$.
(c) The green line is the downgraded version of approximation displayed in
figure (b) by removing three terms $x,x^{4}$, and $x^{17}$ following our
methodology. (d) The green line is a polynomial approximation with the same
number of terms as in (c), but with fixing $k=14$.

####  III.1.2 Biorthogonal polynomials from Laguerre polynomials

In this case the orthogonal polynomials $\\{p_{n}(x)\\}_{n=0}^{k}$ generating
the biorthogonal polynomials Eq. (9) are the Laguerre polynomials
$\\{L_{n}(x)\\}_{n=0}^{k}$ in $V([0,\infty),w)$, with $w(x)=\exp(-x)$. These
polynomials can be obtained by formula [17]

$\displaystyle L_{j}(x)=\sum_{i=0}^{j}\binom{j}{i}\frac{(-1)^{i}}{i!}x^{i}.$
(16)

Consequently, the coefficients $a_{i}^{j}$ to obtain the polynomial
$\beta_{n}^{k}(x)$ through Eq. (9) in this particular case are

$\displaystyle a_{i}^{j}=\binom{j}{i}\frac{(-1)^{i}}{i!}.$

Thus, the functions $\beta_{n}^{k}(x)$ adopt the particular form

$\displaystyle\beta_{n}^{k}(x)$
$\displaystyle=\sum_{j=n}^{k}\binom{j}{n}\frac{(-1)^{n}}{n!}\sum_{i=0}^{j}\binom{j}{i}\frac{(-1)^{i}}{i!}x^{i}=\sum_{j=n}^{k}\sum_{i=0}^{j}(-1)^{n+i}\frac{1}{n!}\frac{1}{i!}\binom{j}{n}\binom{j}{i}x^{i},\quad
x\in[0,\infty)$

Subsequently, the coefficients $c_{n}^{k}$ producing the orthogonal projection
of $f$ onto $V_{k}$ as linear combination of simple monomials

$\displaystyle f_{k}(x)=\sum_{n=0}^{k}c_{n}^{k}x^{n},\quad x\in[0,\infty)$
(17)

are calculated as

$\displaystyle c_{n}^{k}$
$\displaystyle=\int_{0}^{\infty}f(x)\beta_{n}^{k}(x)\,e^{-x}\,dx=\sum_{j=n}^{k}\sum_{i=0}^{j}\frac{(-1)^{n+i}}{n!\,i!}\binom{j}{n}\binom{j}{i}\mu_{i}$
(18)

with

$\displaystyle\mu_{i}=\int_{0}^{\infty}x^{i}f(x)e^{-x}\,dx.$ (19)

#### Example 2: Analytical polynomial approximation and approximation
downgrading using Laguerre polynomials (Type A)

This example aims at deriving a closed polynomial form for approximating the
exponential decay $f(x)=\exp(-\alpha x)$, for $\alpha>0$. In the first
instance we use biorthogonal polynomials derived from Laguerre polynomials
$\\{L_{n}(x)\\}_{n=0}^{k}$ over the interval $[0,\infty)$. Consequently, the
generalised moments $\mu_{i}$ for these functions are analytically obtained as

$\displaystyle\mu_{i}=\int_{0}^{\infty}x^{i}e^{-\alpha
x}e^{-x}\,dx=\int_{0}^{\infty}x^{i}e^{-(\alpha+1)x}\,dx=\frac{i!}{(\alpha+1)^{i+1}}.$

On replacing these moments in Eq. (18) we obtain:

$\displaystyle c_{n}^{k}$
$\displaystyle=\sum_{j=n}^{k}\frac{(-1)^{n}}{n!}\binom{j}{n}\sum_{i=0}^{j}\binom{j}{i}\frac{(-1)^{i}}{(\alpha+1)^{i+1}}=\sum_{j=n}^{k}\frac{(-1)^{n}}{n!}\binom{j}{n}\frac{\alpha^{j}}{(\alpha+1)^{j+1}}.$

Thus, our Laguerre-based biorthogonal approximation to the exponential decay
has the closed form given by

$\displaystyle
f_{k}(x)=\sum_{n=0}^{k}\sum_{j=n}^{k}\frac{(-1)^{n}}{n!}\binom{j}{n}\frac{\alpha^{j}}{(\alpha+1)^{j+1}}\,x^{n},\quad
x\in[0,\infty).$ (20)

The approximation ${f}_{k}(x)$, for $\alpha=1$ and $k=14$, is shown in Fig.
2(a). We notice that for $\alpha=1$ and $x>10$ the exponential decay takes
vanishing small values. More precisely, the value of the exponential at $x=10$
is $f(10)=\exp(-10)<4.54\times 10^{-5}$. Consequently, we find it interesting
to compare the approximation given in Eq. (20) to the approximation resulting
from the Legendre-based BP defined on the closed interval $[0,10]$. The latter
approximation readily follows using in Eq (13) and the moments defined as

$\displaystyle\mu_{i}=\int_{0}^{b}x^{i}e^{-\alpha x}\,dx=e^{-\alpha
b}\sum_{j=0}^{i}(-1)^{i-j}\frac{i!}{j!(-\alpha)^{i-j+1}}b^{j}+\frac{i!}{\alpha^{i+1}}$

with $\alpha=1$ and $b=10$. The Legendre-based BP approximation defined on
$[0,10]$ is attained by placing the computed coefficients in Eq. (14). This
approximation for $k=9$ is shown in Fig. 2(a) (in the scale of the figure both
approximations coincide with the true function). The Laguerre-based
biorthogonal approximation produces the maximum point-wise absolute error
$E=\max_{x\in[0,10]}|f(x)-f_{14}(x)|=2.62\times 10^{-4}$ and the Legendre-
based biorthogonal approximation the value
$E=\max_{x\in[0,10]}|f(x)-{f}_{9}(x)|=2.20\times 10^{-4}$ involving fewer
terms. The pointwise absolute errors are shown in logarithmic scale in Fig.
2(b).

We now use both polynomial types to approximate the Gamma distribution
$f(x)=x\exp(-x)$ plotted in Fig. 2(c). In this case, to obtain similar
Laguerre-based and Legendre-based biorthogonal approximations we need $k=17$
and $k=11$ terms respectively. The corresponding maximum absolute errors are
equal to $E=3.90\times 10^{-4}$ for Laguerre-based and $E=8.52\times 10^{-4}$
for Legendre-based approximation. The pointwise absolute errors are shown in
logarithmic scale of Fig. 2(d).

Figure 2: (a) The exponential decay function $f(x)=\exp(-\alpha x)$ with
$\alpha=1$ is given by the black dashed curve defined on the interval
$[0,10]$. Approximations using our methodology outlined in the main text using
the Laguerre polynomials with $k=14$ (blue solid curve) and Legendre
polynomials with $k=9$ (red solid curve) are presented (both approximations
lie directly under each other). (b) Plot of the pointwise absolute error of
the approximations in logarithmic scale from (a). (c) The gamma distribution
$f(x)=x\exp(-\alpha x)$ (black dashed curve) is approximated by Laguerre
polynomials with $k=17$ (blue solid curve) and Legendre polynomials with
$k=11$ (red solid curve), again with both approximations lying directly on top
of each other. (d) Corresponding pointwise absolute error of the
approximations in logarithmic scale from (c).

###  III.2 Biorthogonal polynomials from Type B classic OPS

We consider now families of BP which are generated by classical OPS that we
refereed to as Type B (c.f.(2)). These are

  * (i)

Legendre polynomials $\\{P_{n}(x)\\}_{n=0}^{k}$ corresponding to
${\cal{I}}=[-1,1]$ and $w(x)=1$.

  * (ii)

Chebyshev polynomials $\\{T_{n}(x)\\}_{n=0}^{k}$ corresponding to
${\cal{I}}=[-1,1]$ and $w(x)=\frac{1}{\sqrt{1-x^{2}}}.$

#### Biorthogonal polynomials from Legendre and Chebyshev polynomials on
${\cal{I}}=[-1,1]$.

(i) Legendre polynomials $\\{P_{n}(x)\\}_{n=0}^{k}$ on ${\cal{I}}=[-1,1]$ are
generated as Ref. [17].

$\displaystyle P_{j}(x)=\frac{\sqrt{2j+1}}{2^{j}\sqrt{2}}\sum_{l=0}^{\lfloor
j/2\rfloor}(-1)^{l}\binom{j}{l}\binom{2j-2l}{j}x^{j-2l}.$

(ii) Chebyshev polynomials $\\{T_{n}(x)\\}_{n=0}^{k}$ on ${\cal{I}}=[-1,1]$
are generated as Ref. [17]

$\displaystyle T_{0}(x)$ $\displaystyle=\frac{1}{\sqrt{\pi}},\qquad
T_{j}(x)=\frac{\sqrt{2}}{\sqrt{\pi}}\frac{j}{2}\sum_{l=0}^{\lfloor
j/2\rfloor}(-1)^{l}2^{j-2l}\frac{(j-l-1)!}{l!(j-2l)!}x^{j-2l}\,\quad\text{for
$j=1,2,3,\dots$}.$

Then, for case (i)

$\displaystyle
P_{2i+n}(x)=\frac{\sqrt{4i+2n+1}}{2^{2i+n}\sqrt{2}}\sum_{l=0}^{\lfloor(2i+n)/2\rfloor}(-1)^{l}\binom{2i+n}{l}\binom{2(2i+n)-2l}{2i+n}x^{2i+n-2l}$
(21)

and

$\displaystyle
a_{l}^{2i+n}=(-1)^{l}\frac{\sqrt{4i+2n+1}}{2^{2i+n}\sqrt{2}}\binom{2i+n}{l}\binom{2(2i+n)-2l}{2i+n}.$
(22)

The BP arise using Eq. (21) and Eq. (22) in the equation below

$\displaystyle\beta_{n}^{k}(x)$
$\displaystyle=\sum_{i=0}^{\lfloor(k-n)/2\rfloor}a_{i}^{2i+n}P_{2i+n}(x),\quad
x\in[-1,1].$ (23)

For case (ii)

$\displaystyle T_{0}$ $\displaystyle=\frac{1}{\sqrt{\pi}},\qquad
T_{2i+n}(x)=\frac{\sqrt{2}}{\sqrt{\pi}}\frac{2i+n}{2}\sum_{l=0}^{\lfloor(2i+n)/2\rfloor}(-1)^{l}2^{2i+n-2l}\frac{(2i+n-l-1)!}{l!(2i+n-2l)!}x^{2i+n-2l}$
(24)

so that

$\displaystyle{a^{\prime}}_{0}^{0}$
$\displaystyle=\frac{1}{\sqrt{\pi}},\qquad{a^{\prime}}_{l}^{2i+n}=(-1)^{l}\frac{\sqrt{2}}{\sqrt{\pi}}\frac{2i+n}{2}2^{2i+n-2l}\frac{(2i+n-l-1)!}{l!(2i+n-2l)!}$
(25)

and the concomitant BP are obtained using in Eq. (23) the polynomials (24) and
coefficients (25), instead of (21) and (22).

#### Example 3: Analytical polynomial approximations using Legendre and
Chebyshev polynomials (Type B) with comparison to matrix inversion method

This example applies the Legendre and Chebyshev based classic BP to compute
the coefficients in the polynomial approximating of the continuous function:

$\displaystyle f(x)=(1-x^{2})e^{-x}\sin(8\pi x),\quad x\in[-1,1].$ (26)

Approximating this function up to an $L^{2}$-norm error of $10^{-4}$, using
either Legendre or Chebyshev polynomials, requires to consider an polynomial
approximation of order $k=36$. The function and the approximations are
represented in Fig. 3(a). Figure 3(b) shows the pointwise absolute value error
of the approximations in logarithmic scale.

Figure 3: (a) Approximations using our methodology outlined in the main text
using the Legendre polynomials with $k=36$ (red solid curve) and Chebyshev
polynomials with $k=36$ (grey solid curve) are presented (both approximations
lie directly under each other) to the function
$f(x)=(1-x^{2})\exp(-x)\sin(8\pi x)$ (Eq. (26)) given by the black dashed
curve on the interval $[-1,1]$. (b) Plot of the pointwise absolute error of
the approximations in logarithmic scale from (a). (c) Plot of the
approximation of the function (26) by inverting the Gram matrix using a
polynomial of order $k=36$ (orange solid curve). (d) The pointwise absolute
error arising from the polynomial produced by inverting the Gram matrix (solid
orange curve) compared to the approximation given in (a) by applying our
methodology using Legendre polynomials with $k=36$ (solid red line) which is
three orders of magnitude smaller on average.

######  Remark III.3.

It is worth commenting that the polynomial approximation of the functions (26)
can not be realised by the standard way of solving the least squares problem.
Indeed, the standard method would try to solve the normal equations:

$\displaystyle\left\langle
f,x^{j}\right\rangle=\sum_{n=0}^{k}c_{n}^{k}\left\langle
x^{n},x^{j}\right\rangle,\quad\text{for $j=0,\dots,k$},$

which introduces the need to invert the Gram matrix of elements

$\displaystyle G(n,j)=\left\langle
x^{n},x^{j}\right\rangle=\frac{1-(-1)^{n+j+1}}{n+j+1},\quad\text{for
$n=0,\dots,k$ and $j=0,\dots,k$}.$

This is a Hankel-type matrix which is very badly conditioned even for small
values of $k$ [9]. In particular for $k=36$ the determinant if close to zero
and the condition number is $10^{18}$. Figs. 3(c) and (d) show the functional
approximation and pointwise absolute error when computing using the matrix
inversion method respectively. Comparison to the error produced using our
biorthogonal polynomials via Legendre polynomials is presented in Fig. 3(d).

##  IV Conclusions

A recursive procedure for constructing dedicated sets of polynomials
$\\{\beta_{n}^{k}\\}_{n=0}^{k},\,k\in\mathbb{N}$ in inner product spaces
$V({\cal{I}},w)$ was proposed. Each set fulfils the following properties:

  * (i)

It is biorthogonal to the basis of monomials $\\{x^{n}\\}_{n=0}^{k}$ in
$V({\cal{I}},w)$.

  * (ii)

Span the same subspace $V_{k}$ as $\\{x^{n}\\}_{n=0}^{k},\,k\in\mathbb{N}$.

These two properties amount to a representation of the orthogonal projection
operation onto $V_{k}$, of an element $f\in V({\cal{I}},w)$, as given by

$\displaystyle\hat{\operatorname{P}}_{V_{k}}f=\sum_{n=0}^{k}x^{n}\left\langle
f,\beta_{n}^{k}\right\rangle.$

Since the proposed construction is based on classic orthogonal polynomials
there is no matrix inversion involved in the process. This feature is
particularly convenient to decide the order of a polynomial model by adjusting
the residual error. Additionally, it is straightforward to realise a posterior
reduction of the model by the removal of any monomial from the basis.

The procedure was illustrated by deriving biorthogonal polynomials based on
Legendre, Laguerre and Chebyshev polynomials. The consideration of other
orthogonal polynomials, such as Hermite polynomials, would readily follow.

The convenience of the proposal was highlighted by three numerical examples:
(i) Polynomial regression modelling of noisy data. (ii) Derivation of analytic
polynomial forms for the exponential decay function. (iii) The approximation
of a continuous function which requires a polynomial of degree at least $k=36$
to be well represented. Such an example, which cannot not be numerically
realised by standard techniques involving the inversion of a Hankel-type
matrix, gives an inside of the many applications that could benefit from the
proposed construction of biorthogonal polynomials.

###### Acknowledgements.

The authors thank Dr Juan Neirotti for fruitful discussions on this work. JL
was supported by the European Union’s Horizon 2020 research and innovation
program under the Marie Skłodowska-Curie Grant Agreement No. 823937 for the
RISE project HALT and the Leverhulme Trust Project Grant No. RPG-2021-014.
Appendix A Proof of Property 5. This proof follows the steps required to prove
Theorem 3.1 in Ref. [15]. Assuming properties (i) and (ii) stated in Property
5, we will now present proofs of properties (iii) and (iv). The proof of (iii)
in Property 5 is a direct consequence of the assumed biorthogonality property
$\left\langle x^{m},\beta_{n}^{k}\right\rangle=\delta_{n,m}$ for $m=0,\dots,k$
and $n=0,\dots,k$. To prove statement (iv) in Property 5, we use

$\displaystyle\beta_{n}^{k\backslash\ell_{j}}(x)=\beta_{n}^{k}(x)-\beta_{\ell_{j}}^{k}(x)\frac{\left\langle\beta_{\ell_{j}}^{k},\beta_{n}^{k}\right\rangle}{\left\langle\beta_{\ell_{j}}^{k},\beta_{\ell_{j}}^{k}\right\rangle},\quad
n=0,\ldots,\ell_{j-1},\ell_{j+1},\ldots,k,\quad x\in{\cal{I}},$ (A.1)

to write

$\displaystyle\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}=\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k\backslash\ell_{j}}\right\rangle=\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k}\right\rangle-\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{\ell_{j}}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle\frac{\left\langle\beta_{\ell_{j}}^{k},\beta_{n}^{k}\right\rangle}{\left\langle\beta_{\ell_{j}}^{k},\beta_{\ell_{j}}^{k}\right\rangle}.$
(A.2)

This will enable us to prove the following statements:

  * (a)

$\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}g=g$ for all $g\in
V_{k\backslash\ell_{j}}$.

  * (b)

$\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}g=0$ for all $g\in
V_{k\backslash\ell_{j}}^{\perp}$, with $V_{k\backslash\ell_{j}}^{\perp}$
denoting the orthogonal complement of $V_{k\backslash\ell_{j}}$ in $V_{k}$.

Since every $g\in V_{k\backslash\ell_{j}}$ can be expressed as a linear
combination $g=\displaystyle\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}a_{n}x^{n}$ for some coefficients $a_{n}$
for $n=0,\dots,j-1,j+1,\dots,N$, then using $\left\langle
x^{l},\beta_{n}^{k}\right\rangle=\delta_{n,l}$ for $l=0,\dots,k$, and
$n=0,\dots,k$ we have:

$\displaystyle\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}=\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{n}\sum_{\begin{subarray}{c}l=0\\\
l\neq{\ell_{j}}\end{subarray}}^{k}a_{l}\left\langle
x^{l},\beta_{n}^{k}\right\rangle+0=\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{n}a_{n}=g,$

which proves (a). Moreover, notice that since

$\displaystyle\sum_{\begin{subarray}{c}n=0\\\
n\neq{\ell_{j}}\end{subarray}}^{k}x^{n}\left\langle\cdot,\beta_{n}^{k}\right\rangle=\hat{\operatorname{P}}_{V_{k}}-x^{\ell_{j}}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle,$

Eq. (A.2) can be re-written as

$\displaystyle\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}$
$\displaystyle=\hat{\operatorname{P}}_{V_{k}}-x^{\ell_{j}}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle-\frac{\hat{\operatorname{P}}_{V_{k}}\beta_{\ell_{j}}^{k}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle}{\|\beta_{\ell_{j}}^{k}\|^{2}}+\frac{x^{\ell_{j}}\left\langle\beta_{\ell_{j}}^{k},\beta_{\ell_{j}}^{k}\right\rangle}{\|\beta_{\ell_{j}}^{k}\|^{2}}$
$\displaystyle=\hat{\operatorname{P}}_{V_{k}}-\frac{\hat{\operatorname{P}}_{V_{k}}\beta_{\ell_{j}}^{k}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle}{\|\beta_{\ell_{j}}^{k}\|^{2}}=\hat{\operatorname{P}}_{V_{k}}-\frac{\beta_{\ell_{j}}^{k}\left\langle\cdot,\beta_{\ell_{j}}^{k}\right\rangle}{\|\beta_{\ell_{j}}^{k}\|^{2}}.$
(A.3)

Since $V_{k\backslash\ell_{j}}^{\perp}$ is spanned by the single element
$\beta_{\ell_{j}}^{k}\in V_{k}$ for all $g\in
V_{k\backslash\ell_{j}}^{\perp}=\alpha\beta_{\ell_{j}}^{k},$ for some constant
$\alpha\in\mathbb{R}$. Then

$\displaystyle\hat{\operatorname{P}}_{V_{k\backslash\ell_{j}}}\alpha\beta_{\ell_{j}}^{k}=\alpha\beta_{\ell_{j}}^{k}-\alpha\beta_{\ell_{j}}^{k}=0,$

which demonstrates (b) and completes the proof of Property 5.

## References

  * [1] E. Ostertagová, “Modelling using Polynomial Regression”, Procedia Engineering, 48, 500–506 (2012).
  * [2] M. Ekum and A. Ogunsanya, “Application of hierarchical polynomial regression models to predict transmission of COVID-19 at global level”, Int J Clin Biostat Biom, 6 (2020) DOI: 10.23937/2469-5831/1510027.
  * [3] M. Ajona, P. Vasanthi and D.S. Vijayan, “Application of multiple linear and polynomial regression in the sustainable biodegradation process of crude oil”, Sustainable Energy Technologies and Assessments, 54, (2022), 102797.
  * [4] J. Munkhammar, L. Mattsson and J. Rydén, “Polynomial probability distribution estimation using the method of moments”, PLOS One 2017, https://doi.org/10.1371/journal.pone.0174573.
  * [5] G. Finlayson, M. Mackiewicz, A. Hurlbert, “Color Correction Using Root-Polynomial Regression” IEEE Transactions on Image Processing, 24, 1460–1470 (2015).
  * [6] M. Shi, W. Hu, M. Li, J. Zhang, X. Song and W. Sun, “Ensemble regression based on polynomial regression-based decision tree and its application in the in-situ data of tunnel boring machine”, Mechanical Systems and Signal Processing, 188 (2023) 110022
  * [7] J. Singh and S. Singh, “A review on Machine learning aspect in physics and mechanics of glasses”, Materials Science and Engineering: B, 284 (2022) 115858.
  * [8] G.G. Chrysos, S. Moschoglou, G. Bouritsas, J. Deng, Y. Panagakis and S. Zafeiriou, “Deep Polynomial Neural Networks”, IEEE Transactions on Pattern Analysis and Machine Intelligence 44, 4021–4034 (2022).
  * [9] N.L. Zamarashkin and E.E. Tyrtyshnikov, “Eigenvalue estimates for Hankel matrices”, Sbornik: Mathematics 192 537–550 (2001).
  * [10] C. Szegiö, Collected papers, Birkhaäuser, Boston, 1982.
  * [11] A. Iserles and S.P. Nørsett, “On the theory of bi-orthogonal polynomials”, Trans. Amer. Math. Soc. 306, 455–474 (1988).
  * [12] A. Iserles and S.P. Nørsett, Christoffel-Darboux-type formulae and a recurrence for bi-orthogonal polynomials. J. Constr. Approx. 5 437–453 (1989).
  * [13] A. Iserles and S.P. Nørsett, Explicit representations of biorthogonal polynomials, Numer. Algor. 10, 51–68 (1995).
  * [14] L. Rebollo-Neira and D. Lowe, “Optimised orthogonal matching pursuit approach”, IEEE Signal Process. Letters, 9, 137–140 (2002).
  * [15] L. Rebollo-Neira, “Backward adaptive biorthogonalisation”, Journal of Mathematics and Mathematical Sciences (2002).
  * [16] L. Rebollo-Neira, “Constructive updating/downdating of oblique projectors: a generalisation of the Gram-Schmidt process”, J. Phys. A: Math. Theor. 40, 6381–6394 (2007).
  * [17] E. Kreyszig, Introductory Functional Analysis with Applications, John Wiley & Sons, 1978.
  * [18] M. P. Priestley, Spectral Analysis and Time Series. Academic Press, 1981, ISBN 978-0-12-56492-3. .