On smoothing of data using Sobolev polynomials

1.   Introduction

Data smoothing is a method commonly used to obtain a smooth approximation of a crude data set. If $\{ f_1, f_2, \dots, f_n\}$ is a sequence of $n$ data points, the goal of data smoothing is to find a corresponding sequence of graduated data points $\{u_1, u_2, \dots, u_n\}$ that provides a better representation of the underlying unknown true values ^[7,33]. Data smoothing is also referred to as data graduation in actuarial science and is commonly used in actuarial studies to compute smooth mortality rates from crude data.

Goodness-of-fit and smoothness are the two most important criteria in data smoothing ^[7]. The Whittaker-Henderson method (WHM), a common non-parametric technique often presented as an alternative to the moving average method ^[18,20] in data smoothing, generates the smoothed data points while balancing both criteria. In WHM, the sequence of smoothed data points is the minimizer $u = \{u_i\}_{i = 1}^n$ of the function

where

The first term is a measure of fidelity to the original data, while the second term is a measure of smoothness. The weights $\omega_i$'s are pre-determined and data-dependent. Note that for high values of $\lambda$, the smoothing is favored, while lower values result in smoothed data that is closer to actual values. A WHM where the order of the difference operator $\Delta$ is set to $2$ is known as the Hodrick-Prescott (HP) filter ^[15].

Several techniques for data graduation can be found in ^{[12,16,33,43,44]}. These methods treat data graduation as a minimization problem in $\mathbb{R}^n$. In ^[26], the WHM in (1.1) is generalized as a minimization problem in the function space $L^2(\Omega)$. By treating the data $\{u_1, u_2, \ldots, u_n\}$ as a piecewise linear function, the data graduation problem in (1.1) is solved by minimizing the functional $\tilde{J}(u): L^2(\Omega) \longrightarrow \mathbb{R}$ given by

where $f$ is an interpolation of the data points $\{ f_1, f_2, \dots, f_n\}$, $\omega$ is an interpolation of the weights $\{ \omega_1, \omega_2, \dots, \omega_n\}$, and ${\textbf D}$ is the derivative operator.

The functional (1.2) can be viewed as the infinite-dimensional generalization of (1.1). Note that the second term in (1.2) is a penalty term for data smoothing. This technique of adding a term for smoothing is called regularization, which has gained popularity in various applications ^{[2,30,31,34,39]}. In this work, we modify the smoothing term of the functional $\tilde{J}$ in (1.2) and consider

where $m$ is the order of the differential operator ${\textbf D}$. By doing this, the solution gains more regularity. To minimize $J$ in (1.3), we project the problem to a finite-dimensional polynomial space. Hence, the solution that we obtain is a polynomial. This approach is motivated by ^[23], where it is argued that the WHM is useful in obtaining smooth data points when these graduated data points approximate a polynomial. One advantage of having a polynomial solution is that it is easier to analyze the solution. With a polynomial solution, one can easily obtain extreme values, the relevant points at which the graduated data are increasing or decreasing, and concavity of the smooth data. Furthermore, having a continuous solution instead of discrete points makes interpolation and extrapolation much easier.

Our proposed method can be used to understand time series trends by analyzing the resulting polynomial solution. In ^[40], the use of polynomials is demonstrated in the analysis of COVID-19 cases. By partitioning the daily new COVID-19 cases into sections and fitting these sections into either a logarithmic, exponential, or linear function, the authors characterized the trend of the disease and identified periods when the incidence of the disease is rapidly or slowly increasing. However, it is not clear how the data set is partitioned. Identifying all possible partitions of data points and evaluating the fit of these points to each of the three functions may be computationally expensive. With polynomials, partitioning can be set at inflection points.

The function $Q$ in (1.1) uses the Newton's advancing operator $\Delta u_i$, which assumes that data points of independent variables must be evenly spaced. Hence, the WHM cannot be used for data with missing terms. Our proposed scheme does not have this limitation because one can still compute the polynomial interpolation of the data, regardless if it is not evenly spaced or some of its terms are missing. Another disadvantage of WHM is that it solves a linear system whose dimension is equal to the number of variables in (1.1). Hence, smoothing high frequency data using WHM requires solving a high-dimensional linear system. In our proposed method, the number of variables depends only on the degree of the polynomial. Thus, regardless of having incomplete or high frequency data, the resulting minimization problem has the same dimension.

This paper is organized as follows: Section 2 discusses the derivation of the minimizer of (1.3) in a finite-dimensional polynomial space. Section 3 presents our proposed algorithm. In Section 4, we present the applications of our method on mortality rates and COVID-19 data. We also apply our method to data with missing terms and high frequency data. Finally, we give our conclusions and recommendations in Section 5.

2.   Theoretical framework

Because we expect the minimizer of (1.3) to be at least $m$-times differentiable, we consider the Sobolev space $H^m(\Omega)$ ^[9]. We present a first-order optimality condition for (1.3). We show that solving this optimality condition is equivalent to minimizing $J$ in (1.3).

Theorem 1. Suppose $f$ and $\omega$ are sufficiently smooth functions. Because $\omega$ is an interpolation of the weights, then we can assume that $0 < \omega(x)\leq \bar{\omega}$, $\forall x\in\Omega$, for some $\bar{\omega} > 0$. Then, $u\in H^m(\Omega)$ is a minimizer of (1.3), if and only if, $u$ satisfies

for all $v\in H^m (\Omega)$ with $m\in\mathbb{N}$.

Proof. Let $v\in H^m (\Omega)$ with $m\in\mathbb{N}$. Define $r:\mathbb{R}\rightarrow \mathbb{R}$ by $r(t): = J(u+tv)-J(u)$. Suppose $u$ is the minimizer of (1.3). Then $J(u)\leq J(u+tv)$ for any $t\in\mathbb{R}$. Hence, $r(t)\geq 0$, $\forall t\in\mathbb{R}$ and $r(t) = 0$ at $t = 0$. By the definition of $J$, we can simplify $r(t)$ as

Hence,

Since $r$ is optimal at $t = 0$, then

Rearranging the above equation gives us (2.1).

Now, suppose $u$ satisfies (2.1). Let $z\in H^m (\Omega)$ and define $v: = z-u$. Then,

From the assumption, the first two terms of the above equation is 0. Thus,

Therefore, for any $z\in H^m(\Omega)$, $J(z)\geq J(u)$, which means that $u$ is a minimizer of $J$.

Theorem 1 tells us that it is sufficient to solve (2.1) to minimize $J$ in (1.3). To show that (2.1) has a unique minimizer, we first need the following results.

Lemma 1. Let $\Omega\subset\mathbb{R}^d$ be bounded with Lipschitz boundary. Suppose $g$ satisfies the following conditions:

(1) $g:H^m(\Omega)\rightarrow [0, +\infty)$ is a seminorm.

(2) There exists a positive constant $C$ such that $0\leq g(v)\leq C\|v\|_{H^m(\Omega)}$, for all $v\in H^m(\Omega)$.

(3) If $v\in\mathcal{P}_{k-1}: = \{\mathit{{space\ of\ polynomials\ of\ degree}}k-1\}$ and $g(v) = 0$, then $v\equiv0$. Then the norm

is equivalent to the norm

Proof. This is a special case of the result proven in ^[38], where the equivalence was shown on a system of functionals and on a more general Sobolev space $W^{m, p}(\Omega)$.

Lemma 2. For $\Omega\subset\mathbb{R}$, we define

Then the norm in (2.4) is equivalent to the norm in (2.2), or equivalently, $\exists \rho_1, \rho_2 > 0$ such that

Proof. The norm in (2.4) is equivalent to

We can set

Note that $g$ is a weighted $L^2(\Omega)$-norm, and hence the first condition of Lemma 1 is satisfied. Because $\omega\leq\bar{\omega}$, for all $x\in\Omega$ and from the definition of norm in (2.2), we get

Therefore, the second condition of Lemma 1 is satisfied. The last condition easily follows because $\dfrac{\omega(x)}{\lambda} > 0$ for all $x\in\Omega$. Therefore, by Lemma 1, our assertion holds.

Theorem 2. The variational formulation in (2.1) has a unique solution in $H^m(\Omega)$.

Proof. We define

and

We use Lax-Milgram's Lemma to prove the existence of the unique solution. For the discussion of this lemma, we refer the readers to ^[9]. To do this, we need to show that $a$ is bilinear, bounded, and coercive, and $b$ is linear and bounded. The bilinearity and linearity of $a$ and $b$, respectively, follows directly from their respective definitions.

We now show that $a$ is bounded, that is, $\exists C_1 > 0$ such that $|a(u, v)|\leq C_1\Vert u\Vert_{H^m(\Omega)} \Vert v\Vert_{H^m(\Omega)}$. We use triangle inequality, Cauchy-Schwarz inequality, and the definition of the norm in (2.2). Thus,

We use (2.5) from Lemma 2 to show that $a$ is coercive, that is, $\exists C_2 > 0$ such that $|a(u, u)|\geq C_2\Vert u \Vert^2_{H^m(\Omega)}$ for all $u\in H^m(\Omega)$. Indeed,

Finally, we show that $b$ is bounded, that is, $\exists C_3 > 0$ such that $|b(v)|\leq C_3\Vert v \Vert_{L^2(\Omega)}$. Using Cauchy-Schwarz inequality, we get

We have shown in Theorem 1 that we can minimize $J$ in (1.3) by solving an equivalent variational formulation, which we have shown in Theorem 2 to have a unique solution in $H^m(\Omega)$. Because of the equivalence of $\Vert \cdot \Vert_{H^m(\Omega)}$ and $\vert \cdot \vert_{H^m(\Omega)}'$, we can use the following inner product for $H^m(\Omega)$:

where $\omega$ is a continuous function in $\Omega$, and $\lambda > 0$. With the above inner product, we can rewrite the optimality condition in (2.1) as

In this work, we pose the above variational formulation in a finite-dimensional subspace $S$ of $H^m(\Omega)$, that is, we solve

We consider the subspace $S$ which is spanned by a set of orthonormal polynomials $\{p_1, p_2, \dots, p_l, \}$ for some $l\in\mathbb{N}$. We refer to $\{p_1, p_2, \dots, p_l\}$ as Sobolev polynomials because they are constructed using the inner product in Eq (2.6), which involves derivatives, such that each polynomial in this subspace are also in the Sobolev space $W^{m, p}(\Omega)$. We were motivated to use Sobolev polynomials because they are the best polynomial approximation to a function $f$ with respect to the $L^2(\Omega)$ norm ^[28]. A comprehensive review of the history and recent development in the study of Sobolev polynomials can be found in ^[27].

Since $S = {\rm{span}}\{p_1, p_2, \dots, p_l\}$, it is sufficient to find $u\in S$ such that

If $u\in S$, then $u = \sum\limits_{i = 1}^{l} u_ip_i$ and so, we solve

Equivalently,

Because $\langle p_i, p_j \rangle _H = \delta_{ij}$, we obtain

Therefore, the solution $u\in S$ of

is given by

We construct the orthonormal polynomials using Gram-Schmidt orthonormalization process ^[24]. We will also use polynomial interpolations for $\omega$ and $f$ so that the solution $u$ is a linear combination of polynomials. The Gram-Schmidt process is started using Chebyshev polynomials as the initial basis functions ^[36].

3.   Proposed algorithm

To compute $u$ in (2.8), we first obtain the polynomial interpolations of $\omega$ and $f$. Then we get a set of orthonormal polynomials from an independent set of polynomials. In our case, we use the Chebyshev polynomials $\{c_1, c_2, \ldots, c_l\}$, for some user-defined degree $l-1$. Then, using Gram-Schmidt orthonormalization process, we construct a set of orthonormal polynomials $\{p_1, p_2, \ldots, p_l\}$. The order of differentiation $m$ is also user-defined. The construction of the Sobolev polynomials is summarized in Algorithm 1.

The user can also specify the value for the smoothing parameter $\lambda$. However, by default, the algorithm uses the approach presented in ^[12,42]. The smoothing parameter $\lambda$ is computed as the minimizer of the generalized cross validation (GCV) score, which is expressed as

where

$n$ is the number of data points, $\gamma_i's$ are the eigenvalues of $D^T D$, and $D$ is a tridiagonal matrix with entries given by

with $h_i$ representing the step between $\hat{f}_i$ and $\hat{f}_{i+1}$.

For a detailed discussion on how to use GCV to identify the smoothing parameter, we refer the readers to ^{[11,12,14,41,42]}. To compute the GCV value numerically, we used the Matlab code provided in the Appendix of ^[12]. The $GCV$ function in (3.1) may have multiple local minima so we have to use a global minimizer. In this study, we use Genetic Algorithm (GA) which has growing applications in science and engineering because of its capability to estimate the global minimum and its non-reliance on the derivative of the objective function ^[21,22,37]. We utilize the Matlab built-in function $\texttt{ga}$, which only requires the upper and lower bounds for the smoothing parameter. Although GA can converge to the global minimum, it can still sometimes get stuck at a local minimum. To guarantee global convergence, we run $\texttt{ga}$ 10 times and store the best solution $\lambda_{GA}$. To further improve accuracy, we hybridize GA with $\texttt{fmincon}$, a Matlab built-in code for interior point algorithm ^[6], which is a local search technique. We use $\lambda_{GA}$ as the initial guess of $\texttt{fmincon}$. We denote $\lambda^\star$ as the minimizer of the hybrid GA-interior point method. Upon obtaining $\lambda^\star$, we can use the orthogonal polynomial calculated in Algorithm 1 and obtain the smooth polynomial approximation $u$ in (2.8). We summarize the proposed method in Algorithm 2.

Given a data $\{f_1, f_2, \ldots, f_n\}$, a unique smooth approximation $u$ is calculated using (2.8). Note that the orthogonal polynomials $p_i$ do not rely on the data. This means that if the data has a particular error structure, uncertainty might occur in the coefficient terms $u_i$ in (2.7). In this study, we rely on a general bootstrap approach presented in ^[8,10] to quantify the uncertainty of the coefficients and construct confidence intervals. The step-by-step procedure is stated in Algorithm 3.

4.   Results and discussion

In this section, we test the proposed algorithm to determine the graduated values of different datasets. All numerical simulations were performed in Matlab 2021a on a computer with Intel(R) Core(TM) i7-8550U CPU clocked at 1.80GHz and 1.99 GHz with 8 GB of RAM that runs Windows 10 OS. By default, we set the order of derivative to 10 and the degree of the polynomial to 8. In all simulations, the bounds for the smoothing parameter is set to $\left[\lambda_{\min}, \lambda_{\max} \right]: = [0.1, 5]$.

For the first application, we test our method to calculate the graduated value of male mortality rates. We obtained the crude mortality rates from the 2017 Philippine Intercompany Mortality Study of the Actuarial Society of the Philippines ^[1]. The illustration of our proposed scheme is shown in the left panel of Figure 1. The red curve represents the linearly interpolated data and the black dashed curve shows our proposed method. The smooth approximation using Whittaker-Henderson graduation is shown in green. It can be seen how our proposed method obtained similar results as the Whittaker-Henderson method. Both methods produced mortality rates that corrected the fluctuations observed in the crude data. Moreover, plots produced from both methods also continuously increase across all ages. While graduated data produced by the discrete Whittaker-Henderson method may appear smooth because of interpolation, the graduated data produced by our method is guaranteed to be smooth because it is represented by a polynomial. Note that the crude data for the mortality rates for ages 80 to 85 are not available. The Whittaker-Henderson method cannot be used to obtain the graduated values outside the range of the data set. In ^[1], extrapolation was used to calculate the graduated values for ages 80–85. The authors used the Gompertz-Makeham model (blue curve) to do this. An advantage of our proposed method is that extrapolation is straightforward. One simply needs to evaluate the polynomial at the points outside the range of the data set. The computed graduated value using the Gompertz-Makeham model and our method are both shown in Figure 1. The right panel shows the plot of the $GCV$ function in (3.1) for $\lambda\in[0.1, 5]$. Observe that the plot (blue curve) is nonlinear and multi-modal which justifies the use of the hybrid GA and interior point method as a global minimization algorithm. The red dot shows the global minimize ($\lambda^\star = 0.5285$), which we set as the smoothing parameter for this problem. Figure 2 presents our uncertainty analysis for this problem using Algorithm 3. For this simulation, we assume a Gamma error structure ^[5]. The uncertainty in our calculated coefficient values translates into the confidence bounds (red dashed lines) around the smooth approximation of the data. Since the degree of the polynomial is 8, the dimension of the subspace $S$ is 9, which means that we calculate 9 coefficients $u_i$'s. The histograms and confidence intervals are presented in the bottom panels of Figure 2. Observe that all the computed coefficients fall within their corresponding 95% confidence intervals.

In the proposed algorithm, the user can choose the degree of polynomial. The left panel of Figure 3 demonstrates our next example to determine the graduation of male mortality rates for ages 0–30 using various degrees of polynomial while fixing the order of derivative to 10. For this example, we calculated $\lambda^\star = 1.3798$. As presented, the graduation of mortality rates can produce a range of polynomial functions that practitioners can choose from depending on their purpose. However, it is also possible to set the degree of polynomial based on the fit of the smoothing curve. One approach is to iterate over a set of positive integers and select the degree that produces the polynomial that has least distance to the piecewise linear interpolation of data. The right panel of Figure 3 demonstrates this and shows that a polynomial with degree 8 ($l = 9$) has the best fit. This polynomial is shown as the red curve in the left panel.

We also applied our proposed method to obtain a smooth time series data for daily new COVID-19 cases in Germany ^[17]. This dataset spans new COVID-19 cases from 1 September 2020 to 31 May 2021. Here, $\lambda^\star = 0.7227$. The left panel in Figure 4 presents the crude daily new COVID-19 cases in Germany as blue line, and the graduated data as red line. As shown, the plot of the graduated data points exhibits peaks and troughs that follow the crude data. The underlying trend in the data can also be observed after smoothing.

The right panel in Figure 4 compares the results of smoothing using different methods, implemented using the Matlab built-in function $\texttt{smoothdata}$. Simple moving average, gaussian or kernel smoothing, and smoothing using locally weighted scatter plot smoothing (LOWESS) all use a window size of 7 because this allows coverage of both the incubation period and the time from the first appearance of symptoms to diagnosis ^[19].

Daily new COVID-19 cases can be viewed as the change in the total number of cases for every small change in time. Hence, by integrating the interpolating polynomial, the cumulative number of COVID-19 cases can be approximated. As an illustration, we integrate the polynomial obtained in Figure 4 to compute the cumulative number of cases. In Figure 5 (left panel), we illustrate how the integral of the polynomial closely approximates the actual cumulative data. The resulting approximation of the cumulative data is also a polynomial, which makes finding the inflection points and concavity of the curve easier. These changes in concavity are shown in Figure 5 (right panel). The segments of the plot in red indicate a deceleration in the number of cases while the blue segments indicate acceleration. Plots like this are useful in evaluating the effectiveness of policies in containing the pandemic, or in identifying events that may have contributed to a more rapid spread of the disease.

We also applied our method to obtain the trend in the total COVID-19 deaths in China from 23 January 2020 to 22 February 2020 ^[29]. For this example, we computed $\lambda^\star = 0.8018$. It was shown in ^[3] that the growth in the number of deaths related to COVID-19 in the early stages of the pandemic follows a quadratic trend. By setting the degree of polynomial to 2, we obtain the graduation in Figure 6, which captures the trend of the data of COVID-19 deaths.

One of the important applications of our proposed method is in handling time series data with missing terms. Figure 7 demonstrates this in the case of missing time series data points for Schistosomiasis cases in years $2009$–$2011$ and $2013$ ^[35]. Here, $\lambda^\star = 1.5577$. As shown in the left panel, our method produced graduated data points that are also consistent with the general trend in the dataset. Note that both WHM and moving average technique require that the time series data be evenly spaced. We also applied other methods in extrapolating the missing data points. The right panel of Figure 7 compares the results of our smoothing algorithm when paired with either linear, spline, piecewise cubic hermite interpolating polynomial (PCHIP), and the modified Akima piecewise cubic hermite interpolation (MAkima) to fill the missing data. We used the Matlab built-in function $\texttt{fillmissing}$ using different methods in our numerical experiments. Our method was applied after extrapolating the data. For this example, smoothing using the points generated by the piecewise linear interpolation produced the least $L^2$-norm error. The user can choose which approach they prefer when filling in missing data. This is a pre-processing step that should be implemented before applying our method to data with missing values.

Our method offers a useful alternative in smoothing high-frequency datasets. WHM entails high computational cost in smoothing large datasets because it requires solving a linear system whose dimension is equal to the number of data points. On the contrary, our method is more suitable because the number of unknowns depends only on the degree of the polynomial used in the graduation. We demonstrate this in Figure 8 using temperature data from ^[32] collected every $5$ minutes from 24 April 2019 to 1 October 2020. The data set considered includes $140$ thousand data points. For this example, we calculated $\lambda^\star = 0.9586$. Applying WHM to this dataset is not possible because our computer can only handle a linear system of a dimension of at most 30 thousand. Moreover, observe that there are missing data points shown in Figure 8 as a break in the blue line. Based on the results of Figure 7, we used linear interpolation to fill the missing data. Our method was able to provide predicted graduated values for the missing sequence in this high-frequency dataset. Moreover, the smooth approximation also exhibits the trend in the data.

5.   Conclusions

In this paper, we presented an alternative method in data graduation that uses Sobolev polynomials. We have proven that the resulting minimization problem has a unique solution in a suitable function space. Furthermore, we formulated an approach using the Gram-Schmidt orthogonalization process to find the solution in an approximate polynomial space. We applied this method in obtaining the graduated values of male mortality rates in the Philippines, and COVID-19 data. We also demonstrated the usefulness of our method in addressing missing data points and smoothing high frequency data. Because the obtained results are polynomials, interpolation and extrapolation become straightforward. Furthermore, inflection points and concavity of the the graduated values are easier to identify. The regularity of polynomials makes the analysis simpler.

Our study has some limitations that can be explored in future research. First, although one can calculate the smoothing parameter by minimizing the GCV score, the order of the derivative and the degree of polynomial are entirely user-defined. If the user does not want to use the default values, the user can fix the order of the derivative and choose the degree of the polynomial that will yield the least $L^2$-error. This can be time-consuming especially if the user sets the order of the derivative to a high value. As a future study, we can explore how both of these relevant parameters can be chosen based solely on data. Second, our study is posed in one-dimensional time series data. For future work, one can extend our method to solve multi-dimensional data smoothing problems. For example, one can consider mortality rates that depend on age and policy duration. Third, one can also consider using other regularization terms (e.g., total variation) and other non-polynomial basis functions. Fourth, in the case of high-frequency data, one can cut the computational cost if the interval is subdivided into segments to reduce the data size ^[15]. Once the data is segmented, the smooth approximation can be implemented using parallel computing. This means that the solution in each segment represents a spline of the entire smooth approximation. However, this will require continuity conditions between segments, which needs a rigorous theoretical analysis. Fifth, another area that can be explored is the relation of our method with space-state models. Since we are treating the smoothing as a minimization in an infinite dimensional Lebesgue space, an extensive study on how to formulate the corresponding difference equations needs to be carried out. These are exciting research directions but demand a thorough investigation and would require other numerical optimization algorithms.

Acknowledgments

This research was supported by a grant from the Computational Research Laboratory of the Institute of Mathematics, University of the Philippines Diliman.

Conflict of interest

The authors declare that they have no conflict of interest.