Simultaneous Approximation of Measurement Values and Derivative Data using Discrete Orthogonal Polynomials

TL;DR

This paper introduces a stable, efficient polynomial approximation method that simultaneously fits measurement values and derivatives using a weighted Hermite approach with orthogonal polynomials, suitable for real-time applications.

Contribution

The paper develops a numerically stable, recurrence-based Hermite approximation method that integrates value and derivative data with covariance weighting, improving efficiency and stability over existing techniques.

Findings

The method achieves higher numerical stability in polynomial fitting.

It allows real-time approximation using inner products instead of matrix inverses.

Monte Carlo simulations validate the approach on synthetic data.

Abstract

This paper presents a novel method for polynomial approximation (Hermite approximation) using the fusion of value and derivative information. Therefore, the least-squares error in both domains is simultaneously minimized. A covariance weighting is used to introduce a metric between the value and derivative domain, to handle different noise behaviour. Based on a recurrence relation with full re-orthogonalization, a weighted polynomial basis function set is generated. This basis is numerically more stable compared to other algorithms, making it suitable for the approximation of data with high degree polynomials. With the new method, the fitting problem can be solved using inner products instead of matrix-inverses, yielding a computational more efficient method, e.g., for real-time approximation. A Monte Carlo simulation is performed on synthetic data, demonstrating the validity of the method. Additionally, various tests on the basis function set are presented, showing the improvement on the numerical stability.

Figures

Equations

$$\hat{\mathbf{y}}$$

$${\hat{\mathbf{y}}}^{\prime }$$

$$\mathbf{y}$$

$${\mathbf{y}}^{\prime }$$

$$\mathbf{s}$$

$$\mathbf{t}$$

$$E\left(\mathbf{\gamma }\right)$$

$$=$$

$$\left({\mathsf{B}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\mathsf{B}+{\mathsf{B}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathsf{B}}^{\prime }\right)\mathbf{\gamma }={\mathsf{B}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\hat{\mathbf{y}}+{\mathsf{B}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\hat{\mathbf{y}}}^{\prime }.$$

$${\mathsf{P}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\mathsf{P}+{\mathsf{P}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathsf{P}}^{\prime }=\mathsf{I},$$

$$\mathbf{\gamma }={\mathsf{P}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\hat{\mathbf{y}}+{\mathsf{P}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\hat{\mathbf{y}}}^{\prime },$$

$${\gamma }_i={\mathbf{p}}_{\text{i}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\hat{\mathbf{y}}+{\mathbf{p}}_{\text{i}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\hat{\mathbf{y}}}^{\prime },$$

$$\tilde{\mathbf{y}}$$

$${\tilde{\mathbf{y}}}^{\prime }$$

$${\mathbf{p}}_{\text{k+1}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{p}}_{\text{k+1}}+{\mathbf{p}}_{\text{k+1}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{p}}_{\text{k+1}}^{\prime }=1$$

$${\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{p}}_{\text{k+1}}+{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{p}}_{\text{k+1}}^{\prime }=0,$$

$$p_{\text{k+1}}\left(x\right)=\alpha x{p}_{\text{k}}\left(x\right)-\beta p_{\text{k}}\left(x\right)-\gamma p_{\text{k-1}}\left(x\right).$$

$$p_{\text{k+1}}^{\prime }\left(x\right)=\alpha x{p}_{\text{k}}^{\prime }\left(x\right)+\alpha x^{\prime }{p}_{\text{k}}\left(x\right)-\beta p_{\text{k}}^{\prime }\left(x\right)-\gamma p_{\text{k-1}}^{\prime }\left(x\right).$$

$${\mathbf{p}}_{\text{k+1}}=\alpha \mathbf{x}\circ {\mathbf{p}}_{\text{k}}-\beta {\mathbf{p}}_{\text{k}}-\gamma {\mathbf{p}}_{\text{k-1}}$$

$${\mathbf{p}}_{\text{k+1}}^{\prime }=\alpha \mathbf{x}\circ {\mathbf{p}}_{\text{k}}^{\prime }+\alpha {\mathbf{x}}^{\prime }\circ {\mathbf{p}}_{\text{k}}-\beta {\mathbf{p}}_{\text{k}}^{\prime }-\gamma {\mathbf{p}}_{\text{k-1}}^{\prime }.$$

$${\mathbf{p}}_{\text{k+1}}=\alpha \mathbf{x}\circ {\mathbf{p}}_{\text{k}}-{\mathsf{P}}_{\text{k}}\mathbf{\beta }$$

$${\mathbf{p}}_{\text{k+1}}^{\prime }=\alpha \mathbf{x}\circ {\mathbf{p}}_{\text{k}}^{\prime }+\alpha {\mathbf{x}}^{\prime }\circ {\mathbf{p}}_{\text{k}}-{\mathsf{P}}_{\text{k}}^{\prime }\mathbf{\beta }.$$

$${\mathbf{u}}_{\text{k}}=\mathbf{x}\circ {\mathbf{p}}_{\text{k}}$$

$${\mathbf{v}}_{\text{k}}=\mathbf{x}\circ {\mathbf{p}}_{\text{k}}^{\prime }+{\mathbf{x}}^{\prime }\circ {\mathbf{p}}_{\text{k}},$$

$${\mathbf{p}}_{\text{k+1}}=\alpha {\mathbf{u}}_{\text{k}}-{\mathsf{P}}_{\text{k}}\mathbf{\beta }$$

$${\mathbf{p}}_{\text{k+1}}^{\prime }=\alpha {\mathbf{v}}_{\text{k}}-{\mathsf{P}}_{\text{k}}^{\prime }\mathbf{\beta }.$$

$$\left({\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathsf{P}}_{\text{k}}+{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathsf{P}}_{\text{k}}^{\prime }\right)\mathbf{\beta }=\alpha \left({\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{u}}_{\text{k}}+{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{v}}_{\text{k}}^{\prime }\right).$$

$${\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathsf{P}}_{\text{k}}+{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathsf{P}}_{\text{k}}^{\prime }=\mathsf{I},$$

$$\mathbf{\beta }=\alpha \left({\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{u}}_{\text{k}}+{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{v}}_{\text{k}}\right).$$

$${\mathbf{p}}_{\text{k+1}}=\alpha \left({\mathbf{u}}_{\text{k}}-{\mathsf{P}}_{\text{k}}{\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{u}}_{\text{k}}-{\mathsf{P}}_{\text{k}}{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{v}}_{\text{k}}\right)=\alpha {\mathbf{c}}_{\text{k+1}}$$

$${\mathbf{p}}_{\text{k+1}}^{\prime }=\alpha \left({\mathbf{v}}_{\text{k}}-{\mathsf{P}}_{\text{k}}^{\prime }{\mathsf{P}}_{\text{k}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{u}}_{\text{k}}-{\mathsf{P}}_{\text{k}}^{\prime }{\mathsf{P}}_{\text{k}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{v}}_{\text{k}}\right)=\alpha {\mathbf{c}}_{\text{k+1}}^{\prime }.$$

$$\alpha =\sqrt{\frac{1}{{\mathbf{c}}_{\text{k+1}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}{\mathbf{c}}_{\text{k+1}}+{\mathbf{c}}_{\text{k+1}}^{\prime \text{T}}{\mathsf{W}}_{\mathrm{d}\mathbf{y}}{\mathbf{c}}_{\text{k+1}}^{\prime }}}.$$

$${\mathbf{p}}_{\text{0}}=\frac{\mathbf{e}}{\sqrt{{\mathbf{e}}^{\text{T}}{\mathsf{W}}_{\mathbf{y}}\mathbf{e}}},$$

Full text

A Monte Carlo simulation is performed on synthetic data, demonstrating the validity of the method. Additionally, various tests on the basis function set are presented, showing the improvement on the numerical stability.

I Motivation

Measurements with a fusion of value and derivative information are common in geotechnical monitoring, e.g. Fig. 1. In this case, inclinometer measurements and reference measurements from total stations are used to reconstruct the measurement line and detect unwanted changes or deflections. These measurements are perturbed by noise of different characteristics, i.e., the reference measurements act as constraints and are normally more precise than inclinometer measurements. Especially in engineering problems, polynomials are the model of choice for the approximation of such data, due to their properties and their relation to the underlying physical model.

Within this work, we develop a framework for the approximation with polynomials of given values and derivatives with different noises characteristics using covariance weighted discrete orthogonal polynomials. The developed method proves to be suitable for high degree polynomial approximation, since it is numerically more stable than other methods.

The main contribution of this paper span:

The derivation of a novel methodology for the generation of a discrete orthogonal polynomial basis function set which can be used to approximate values and derivatives which are subjected to noise of different levels. It uses a three-term recurrence relation with full re-orthogonalization (to increase numerical stability) together with covariance weighting of the residual (for introducing a metric between the value and derivative domain). 2. 2.

The derivation of covariance propagation for the coefficients and the resulting approximation. 3. 3.

The introduction of various measures for the evaluation of the numerical stability.

The paper is structured as follows: A review of literature is found in Sec. II. In Sec. III all the derivations for the generation of a covariance weighted discrete orthogonal polynomial basis function set are introduced. The derivation of the covariance propagation is presented in Sec. III-D. To demonstrate the validity of the novel method, a Monte Carlo experiment for high degree polynomial approximation is performed on a synthetic data set in Sec. IV. It is shown that the residual between the data and the approximated values is reduced to noise with the same parameters as used in the generation of the synthetic noise. The numerical stability of the generated basis is validated in Sec. V demonstrating the novel method to be advantageous compared to using Vandermonde type basis functions, especially for higher degrees. Various quality measures are tested and presented for both, complete and incomplete basis function sets.

II Review of Literature

The use of value and derivative information for polynomial approximation (Hermite approximation) is not common in literature, whereas the use value and derivative information for interpolation is well-known in literature.

Based on the idea of a Newton type interpolation, Hermite introduced a similar methodology for the interpolation of values and derivatives. The original idea can be found in [1]. Based on this, several books, e.g. [2], introduce this idea to a broader audience.

A generalized version for arbitrary orders of derivatives at given points can be found in [3]. This idea is extended to bivariate functions in [4]. An extensive study on multivariate Hermite interpolation can be found in [5].

Since the calculation of those interpolating polynomials is straight forward, Hermite type interpolating polynomials gained popularity for the approximation of functions using a two point approach rather than Taylor expansion [6, 7]. They are well studied in terms of error bounds [8, 9] and compared to other standard methods [10]. They are used to solving ordinary differential equation [11, 12] and partial differential equation [13, 14, 15] by approximating the equations using Hermite type interpolation. It is pointed out, that using this method lead to higher order of approximation, which improves step-size for iterative methods compared to standard methods.

Another use of the Hermite interpolation is presented in [16]. In this work it is used to perform a moving-least-squares approximation based on local basis functions which use value and derivative information.

The idea of including derivative information within the approximation of data is mostly found in constrained polynomial approximation, e.g., [17, 18, 19, 20]. There it is assumed, that some information (e.g. reference points) is $100\text{\,}\mathrm{\char 37\relax}$ certain. Approximation of data with uncertain value and derivatives is introduced in [21].

The idea of including a covariance weighting on the residuals to get a valid metric between the value and derivative domain is introduced in [18] and extended to higher order derivatives in [22].

Within this paper we use this idea of covariance weighting within the generation of a discrete orthogonal polynomial basis function set to improve the numerical stability, making it suitable for high degree polynomial fitting.

III Theoretical Framework

The herein presented framework addresses the approximation of a polynomial model given noisy values and collocated derivatives using discrete orthogonal polynomials, i.e., systems of polynomials that satisfy a discrete orthogonality constraint.

III-A Modelling of Measured Values

The measured noisy values and derivatives are collected in the vectors $\hat{\boldsymbol{y}}=\begin{bmatrix}\hat{y}_{1}&\hat{y}_{2}&\dots&\hat{y}_{n}\end{bmatrix}^{\textrm{T}}$ and $\hat{\boldsymbol{y}}^{\prime}=\begin{bmatrix}\hat{y}_{1}^{\prime}&\hat{y}_{2}^{\prime}&\dots&\hat{y}_{n}^{\prime}\end{bmatrix}^{\textrm{T}}$. They are sampled at the positions $\boldsymbol{x}=\begin{bmatrix}x_{1}&x_{2}&\dots&x_{n}\end{bmatrix}^{\textrm{T}}$. For modelling the measurement, we assume that we measure the true value plus covariant noise. This covariant noise is generated from a vector of gaussian random variables with zero mean and unit variance (i.i.d. noise) together with the according covariance matrix. Mathematically this is described as

[TABLE]

where $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime}$ are the true values, $\mathsf{\Lambda}_{\boldsymbol{y}}$ and $\mathsf{\Lambda}_{\,\mathrm{d}\boldsymbol{y}}$ are the associated covariance matrices and $\boldsymbol{s}$ and $\boldsymbol{t}$ are vectors of i.i.d. noise.

III-B Approximation of Values and Derivatives

The goal is now to approximate the given data in a least-squares sense, using a set of discrete basis functions collected in the columns of the matrix $\mathsf{B}$ and their derivatives $\mathsf{B}$’. The true values $\boldsymbol{y}$ and derivates $\boldsymbol{y}^{\prime}$ are modelled as

[TABLE]

where $\mathsf{B}^{\prime}$ denotes the first derivative of the basis functions $\mathsf{B}$ with respect to $x$ and $\boldsymbol{\gamma}$ is the coefficient vector. To apply Gauss’s least-squares theorem, the involved errors must be i.i.d., so we solve for the gaussian random variables $\boldsymbol{s}$ and $\boldsymbol{t}$ using Eqn. (1)-(4), yielding

[TABLE]

This results in the following functional to be minimized,

[TABLE]

Clearly, this introduces covariance weighting on the residual and thus a metric between value and derivative domain is established as stated in [18] and[22]. Substituting $\mathsf{W}_{\boldsymbol{y}}^{\frac{1}{2}}\triangleq\mathsf{\Lambda}_{\boldsymbol{y}}^{-\frac{1}{2}}$ and $\mathsf{W}_{\,\mathrm{d}\boldsymbol{y}}^{\frac{1}{2}}\triangleq\mathsf{\Lambda}_{\,\mathrm{d}\boldsymbol{y}}^{-\frac{1}{2}}$ and minimizing Eqn. (7) leads to the normal equations for weighted regression [23],

[TABLE]

To solve this equation for $\boldsymbol{\gamma}$ one can use standard methods, e.g, inverting $\left(\mathsf{B}^{\textrm{T}}\mathsf{W}_{\boldsymbol{y}}\mathsf{B}+\mathsf{B}^{\prime\textrm{T}}\mathsf{W}_{\,\mathrm{d}\boldsymbol{y}}\mathsf{B}^{\prime}\right)$ which is known to be computational costly.

To overcome this problem, we developed a method to find a polynomial basis function set $\mathsf{P}$ and its derivative $\mathsf{P}^{\prime}$, which fulfil

[TABLE]

i.e., a discrete orthogonality condition. The matrices $\mathsf{P}=\begin{bmatrix}\boldsymbol{p}_{\text{0}}&\boldsymbol{p}_{\text{1}}&\dots&\boldsymbol{p}_{\text{i}}&\dots&\boldsymbol{p}_{\text{d}}\end{bmatrix}$ and $\mathsf{P}^{\prime}=\begin{bmatrix}\boldsymbol{p}_{\text{0}}^{\prime}&\boldsymbol{p}_{\text{1}}^{\prime}&\dots&\boldsymbol{p}_{\text{i}}^{\prime}&\dots\boldsymbol{p}_{\text{d}}^{\prime}\end{bmatrix}$ are a collection of discrete polynomial basis functions and their derivatives. Each column $\boldsymbol{p}_{\text{i}}$ represents a polynomial of degree $i$ and $\boldsymbol{p}_{\text{i}}^{\prime}$ is the first derivative of that discrete polynomial. They are sorted on increasing degrees yielding a basis function set of degree $d$. Respectively a linear combination of those basis functions yield a polynomial of degree $d$ at most.

The coefficients $\boldsymbol{\gamma}$ for the approximating polynomial can then be easily calculated from Eqn. (9) using,

[TABLE]

which are inner products of the basis functions and the covariance weighted measurements. That is, the coefficient $\gamma_{i}$ for a certain basis function $\boldsymbol{p}_{\text{i}}$ of degree $i$ can be directly calculated as

[TABLE]

which is similar to the calculation of discrete Fourier series. Thus, the computational efficiency is improved compared to solving the fitting problem using standard algorithms including matrix inverses. To calculate the approximated values $\tilde{\boldsymbol{y}}$ and $\tilde{\boldsymbol{y}}^{\prime}$, we use the estimated parameters $\boldsymbol{\gamma}$ within the model equation (3) and (4), yielding

[TABLE]

III-C Synthesis of Weighted Discrete Orthogonal Basis

In this section, a novel method for the synthesis of a set of weighted discrete orthogonal polynomials fulfilling Eqn. (10) is presented.

Two important conditions can be directly derived form the identity in Eqn. (10). These are the normal condition

[TABLE]

and the orthogonality condition

[TABLE]

which requires, that a basis function of degree $k+1$ is orthogonal to all basis functions of lower degree. The matrices $\mathsf{P}_{\text{k}}=\begin{bmatrix}\boldsymbol{p}_{\text{0}}&\boldsymbol{p}_{\text{1}}&\dots&\boldsymbol{p}_{\text{k}}\end{bmatrix}$ and $\mathsf{P}_{\text{k}}^{\prime}=\begin{bmatrix}\boldsymbol{p}_{\text{0}}^{\prime}&\boldsymbol{p}_{\text{1}}^{\prime}&\dots&\boldsymbol{p}_{\text{k}}^{\prime}\end{bmatrix}$ collect the basis vectors up to degree $k$.

For the generation of a set of suitable polynomial basis functions $\mathsf{P}$ and their derivatives $\mathsf{P}^{\prime}$ up to a certain degree $d$, we use a recurrence relation with full re-orthogonalisation (as studied in [24, 17]) together with covariance weighting.

Starting from the classical three-term recurrence relation from functional analysis for orthogonal polynomials, e.g. [25], also known as Gram-Schmidt process [23], a polynomial $p_{\text{k+1}}\left(x\right)$ of degree $k+1$ can be generated from lower degree polynomials using

[TABLE]

The derivative of this polynomial with respect to $x$ is calculated as

[TABLE]

This is the definition for the continuous case.

The discrete formulations of Eqn. (17) and (18) for a given vector $\boldsymbol{x}=\begin{bmatrix}x_{1}&x_{2}&\dots&x_{n}\end{bmatrix}^{\textrm{T}}$ are111The operation with the symbol $\circ$ denotes the Hadamard product, i.e., the element-wise product.

[TABLE]

and

[TABLE]

The vectors $\boldsymbol{p}_{\text{i}}=\begin{bmatrix}p_{\text{i}}\left(x_{1}\right)&p_{\text{i}}\left(x_{2}\right)&\dots&p_{\text{i}}\left(x_{n}\right)\end{bmatrix}^{\textrm{T}}$ and $\boldsymbol{p}_{\text{i}}^{\prime}=\begin{bmatrix}p_{\text{i}}^{\prime}\left(x_{1}\right)&p_{\text{i}}^{\prime}\left(x_{2}\right)&\dots&p_{\text{i}}^{\prime}\left(x_{n}\right)\end{bmatrix}^{\textrm{T}}$ are the polynomial of degree $i$ and its derivative evaluated at the points $\boldsymbol{x}$. As shown in [26], the normal three-term recurrence relation is numerically unstable, so we use a complete re-orthogonalization as suggested in [24]. Using this improvement, Eqn. (19) and (20) read as

[TABLE]

and

[TABLE]

Using the substitutions

[TABLE]

and

[TABLE]

Eqn. (21) and (22) read as

[TABLE]

and

[TABLE]

Using the orthogonality condition (16) together with (25) and (26) yields

[TABLE]

Since, per definition, the previous generated basis functions have to fulfil the orthogonality and normal condition

[TABLE]

$\boldsymbol{\beta}$ can be expressed as

[TABLE]

Using this, Eqn. (25) and (26) can be rewritten as

[TABLE]

and

[TABLE]

$\boldsymbol{c}_{\text{k+1}}$ and $\boldsymbol{c}_{\text{k+1}}^{\prime}$ represent a basis function and (its derivative) which is orthogonal to all previous basis functions but not yet normed. To fulfil the norm condition in Eqn. (15), the scaling factor $\alpha$ is calculated as

[TABLE]

This scaling factor and the coefficient vector $\boldsymbol{\beta}$ have to be calculated for each newly generated basis function.

To start the recurrence, the first basis functions $\boldsymbol{p}_{\text{0}}$, $\boldsymbol{p}_{\text{1}}$ and their derivatives have to be defined in an initial step, to meet the above mentioned conditions. The first basis function of degree $d=0$ is defined as

[TABLE]

which is a normalized constant vector. $\boldsymbol{e}$ denotes a vector of ones. The first derivative of this basis function is the zero vector

[TABLE]

For the second basis function $\boldsymbol{p}_{\text{1}}$ we first generate the vector

[TABLE]

which we project onto the orthogonal complement of $\boldsymbol{p}_{\text{0}}$ to meet the orthogonality condition, yielding

[TABLE]

This is a scaled version of $\boldsymbol{p}_{\text{1}}$. To get the slope for the derivatives of the basis functions right, the vector $\boldsymbol{x}^{\prime}$ is calculated as

[TABLE]

Using this, we calculate a scaled version of $\boldsymbol{p}_{\text{1}}^{\prime}$ based on Eqn. (24) and (26), yielding

[TABLE]

To generate the basis functions fulfilling the norm conditions, we calculate the scaling factor

[TABLE]

From this we calculate the second pair of basis functions

[TABLE]

and

[TABLE]

This is the prerequisite to start the synthesis of higher order basis functions. The final set of basis functions $\mathsf{P}\triangleq\mathsf{P}_{\text{d}}$ and $\mathsf{P}^{\prime}\triangleq\mathsf{P}_{\text{d}}^{\prime}$ of a certain degree $d$ can now be used in Eqn. (11), (13) and (14) to calculate the coefficients and to finally approximate perturbed values and its derivatives with a polynomial of degree $d$. A set of basis functions is shown in Fig. 2.

The generation of a basis function set using the presented method is a generalized Gram-Schmidt process. As it can be seen, the generation of the basis depends only on the relative locations of the $x$ values222To improve numerical stability, $\boldsymbol{x}$ is transformed to be centered at the origin and scaled to unit norm. and not on the values itself. If the abscissa values $\boldsymbol{x}$ do not change, the basis function set can be calculated a priori and the solution to the weighted fitting problem reduces to a simple matrix-vector multiplication (see Eqn. (11)). This is advantageous when implemented in smart sensors or low power controllers. A further advantage is, that the weighting matrices $\mathsf{W}_{\boldsymbol{y}}$ and $\mathsf{W}_{\,\mathrm{d}\boldsymbol{y}}$ can be rank-deficient, e.g., points can be weighted with [math] if they should not be considered. This can be helpful to suppress outliers.

III-D Covariance Propagation

Since we are dealing with noisy data, covariance propagation is a prerequisite for making assumptions about the quality of the approximated values. Based on Eqn. (11) the covariance \mathsf{\Lambda}_{\text{\boldsymbol{\gamma}}} for the coefficients $\boldsymbol{\gamma}$ is calculated as

[TABLE]

with the block matrices

[TABLE]

Similarly, the covariance matrices for the approximated values and derivatives can be propagated as

[TABLE]

and

[TABLE]

Since $\mathsf{W}_{\text{c}}\triangleq\mathsf{\Lambda}_{\text{c}}^{-1}$ the middle-term results in $\mathsf{W}_{c}\mathsf{\Lambda}_{\text{c}}\mathsf{W}_{c}^{\textrm{T}}=\mathsf{W}_{c}$. Together with the identity from Eqn. (10), the above equations simplify to

[TABLE]

[TABLE]

and

[TABLE]

These covariance matrices for the approximated coefficients and values, can be used to calculate confidence intervals or predictions intervals. Note: As it can be seen in Eqn. (46), the presented method decorrelates the noise to i.i.d. noise, to accord with Gauss’s theorem.

IV Numerical Example

To test the validity of the herein presented method to approximate a polynomial given perturbed values and derivatives, a synthetic dataset is generated. The underlying function is defined as

[TABLE]

with its analytical first derivative

[TABLE]

The function and its derivative are evaluated in the range $\left[-2\pi,2\pi\right]$ at $n=500$ equally spaced nodes, yielding the vectors of values and derivatives $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime}$. Gaussian noise with different gains $\sigma_{y_{i}}=\sigma_{y}=0.1$ and $\sigma_{\,\mathrm{d}y_{i}}=\sigma_{\,\mathrm{d}y}=2$ is added to those vectors yielding the noisy measurement vectors $\hat{\boldsymbol{y}}$ and $\hat{\boldsymbol{y}}^{\prime}$. A polynomial of degree $d=35$ is used for approximating the noisy data set. To test the developed method, a Monte Carlo simulation is performed with $n_{\text{iter}}=1000$ iterations. As a measure, the standard deviation of the residuals $\operatorname{std}\{\boldsymbol{r}_{\text{y}}\}=\operatorname{std}\{\boldsymbol{y}-\hat{\boldsymbol{y}}\}$ and \operatorname{std}\{\boldsymbol{r}_{\text{\,\mathrm{d}y}}\}=\operatorname{std}\{\boldsymbol{y}^{\prime}-\hat{\boldsymbol{y}}^{\prime}\} are calculated in each run. Since the presented method uses covariance weighting, the standard deviation of the result should be the same as $\sigma_{y}\approx\operatorname{std}\{\boldsymbol{y}-\hat{\boldsymbol{y}}\}$ and $\sigma_{\,\mathrm{d}y}\approx\operatorname{std}\{\boldsymbol{y}^{\prime}-\hat{\boldsymbol{y}}^{\prime}\}$. The mean value of the standard deviations over all runs is shown in Fig. 3.

As it can be seen, although the noise gains are very different, the presented method which uses covariance weighted approximation delivers the correct results for both, values and derivatives demonstrating the method to be valid.

As expected, the method using the Vandermonde basis as presented in [22], is not stable for such a high degree. As one can inspect, the approximation does not follow the signal.

V Numerical Quality of Basis

To verify the numerical quality of herein presented method, a meaningful measure has to be found. As Wilkinson [27] pointed out, a posteriori estimation of error bounds is preferred to a priori error predictions in such cases. The identity in Eqn. (10) can be written in terms of the block matrices as

[TABLE]

Rewriting in terms of a unitary matrix $\mathsf{U}$ yields

[TABLE]

with

[TABLE]

The residual matrix

[TABLE]

should contain zeros in a perfect reconstruction. Due to numerical errors, this does not hold (see Fig. 4).

To summarize the numerical quality of the generated basis, the following error measures are tested in order to find the appropriate measure:

Maximum norm. The maximum norm is the largest single element within the residual matrix, i.e., $\epsilon_{\text{max}}=\lVert\mathsf{R}\rVert_{\text{max}}=\max\{|r_{ij}|\}$. Since this norm depends only on one specific entry, this may lead to wrong conclusions. 2. 2.

Frobenius norm. The Frobenius norm is the square root of the sum of the squares of all entries in the residual matrix, i.e., $\epsilon_{\text{F}}=\lVert\mathsf{R}\rVert_{\text{F}}=\sqrt{\sum_{i}\sum_{j}r_{ij}^{2}}$. This norm is a measure for the total error. 3. 3.

Determinant. The determinant of a matrix is a theoretical quality measure and equals $1$ for an ideal unitary matrix. The variation from this, i.e., $\epsilon_{\text{det}}=1-\det\{U\}$ is a measure for the quality of the tested basis function set. 4. 4.

Condition number. The condition number of a matrix is connected to the error propagation. For a unitary matrix the condition number is $1$. Thus, the measure tested is $\epsilon_{\text{cond}}=1-\operatorname{cond}\{U\}$. 5. 5.

Rank. Since the presented method generates a weighted orthogonal basis function set, the rank of $\mathsf{U}$ should be full rank, i.e., no linear dependencies. Due to round off errors this can be used to find linear dependencies, i.e., $\epsilon_{\text{rank}}=n-\,\mathrm{rank}\left\{\,\mathsf{U}\,\right\}$ .

Based on those measures the approximate number of significant digits $\eta$ is calculated as

[TABLE]

These measures are calculated for a complete basis function set where the number of basis functions equals the number of data points (values and derivatives). Therefore, the degree of the resulting polynomial is $d=2n-1$. In Fig. 5 the different measures are presented for varying degrees.

As it can be seen, the Frobenius norm and the condition number are the most meaningful measures, since they show the highest dependency on the degree of the resulting polynomial. The rank measure is not visible, since the proposed algorithm generates full-rank basis function sets, so there is no error visible. Since the Frobenius norm is a measure for the total error, this measure is chosen to compare to other algorithms in the following.

The new method is compared to the one presented in [22]. It uses a Vandermonde basis function set to solve the same problem. As [17] pointed out, the Vandermonde basis for normal polynomial regression gets degenerate at high degrees. This behaviour can also be inspected within this paper. As it can be seen in Fig. 6, the new method generates a more stable result also for high degrees.

Since this new method can also be used for approximation (overdetermined system of equations), the quality of the basis is determined for an incomplete basis function set for $n=1000$. The result for various degrees of polynomial is visualized in Fig. 7, showing that the new method performs also better for an incomplete basis.

VI Conclusion

The herein presented method introduces a novel polynomial fitting framework for the approximation of value and derivative data. Including both sources of information within the fitting procedure improves the quality of the fitted polynomial improving both, reconstruction of values and derivatives. The method uses a recurrence relation with full re-orthogonalization together with covariance weighting for introducing a metric between value and derivative domain, yielding a set of discrete orthogonal polynomials. As it is shown, the generated basis function set is numerically more stable compared to other methods, especially for high degree polynomials. Using this basis, the fitting problem is reduced to inner products, which is beneficial in terms of computational efficiency. Due to the covariance weighting, the noise associated with the channels is decorrelated to i.i.d. noise, to accord with Gauss’s theorem. Thus, there is no bias based on different noise parameters. The validity of the method is tested and presented on a numerical example, where a polynomial of degree $d=35$ is fitted to a periodic function. In future research, this method will be adapted for fitting data given noisy constraints, based on discrete orthogonal polynomials.

Acknowledgement

This work was partially funded under the auspices of the EIT - KIC Raw materials program within the project “Maintained Mine and Machine” (MaMMa) with the grant agreement number: [EIT/RAW MATERIALS/SGA2018/1].