A weighted online regularization for a fully nonparametric model with heteroscedasticity

1.   Introduction

The underlying model of this article can be formulated quite easily: given some real valued variables $x$ and $y$ fulfill the heteroscedastic model:

where $f(x)$ belongs to Sobolev space $W^{2, 2}(0, 1)$, variance function $\sigma(x)$ is a positive function defined on $(0, +\infty)$ and the random error term $\varepsilon$ is independent of $x$ and satisfies $E(\varepsilon) = 0$ and $Var(\varepsilon) = 1$.

Such a problem is widely applied in various fields such as Computerized Tomography (CT) and the inverse problem of option pricing (IPOP) in ^[1,2,3]. Moreover, this is a classical ill-posed problem. There are plenty of regularization methods for this ill-posed problem in one dimension or higher dimensions ^[4,5]. Researchers ^[6,7,8,9] have proposed some statistical methods to solve this ill-posed problem. Zhang ^[10] employs a relatively small number of knots for interpolation to reduce computation cost. These classical approaches are based on the equal variance of the error term and, when $N$ is large enough, needs large amount of calculating. However, the variance $\sigma(x)$ in (1.1) is typically unequal, a phenomenon known as heteroscedasticity ^[11,12].

Heteroscedasticity, which is an econometric term, is the variances of perturbations $\sigma(x)$ in the model which are not exactly equal^[13]. Heteroscedasticity will lead to large errors in the model. In order to reduce the error caused by heteroscedasticity, the principal methods are mainly divided into three parts. The most common method is to take the logarithm of the data^[14]. However, not all data must be logarithmic. The second method is Robust Standard Error Regression, which is the most popular and effective treatment method at present^[15,16]. The main idea is to modify the objective function in the classical least squares regression, which is very sensitive to outliers. Robust method is only applicable to heteroscedasticity and independent observation. Another method is FGLS regression^[17]. For points with larger residual value, the smaller weight is given to solve the heteroscedasticity problem. It is still the least square method in essence. However, the convergence speed of this method is slow, that is, the deviation of limited samples will be large.

B-spline functions are crucial elements in various fields, especially in computer graphics, computer-aided design (CAD), and numerical analysis^[18,19,20]. They play a significant role in curve and surface representation, interpolation, approximation, and modeling^[21,22,23]. B-spline functions offer a versatile mathematical framework for representing complex shapes and data, enabling efficient and accurate solutions in various industries and scientific disciplines. Their ability to balance smoothness, flexibility, and local control makes them a cornerstone of computational modeling and design processes. Among them, cubic splines are the most representative. Cubic B-spline functions strike a balance between simplicity and expressiveness. They are more flexible than linear or quadratic B-splines, allowing for smoother curves and surfaces, while still being relatively straightforward to manipulate.

In this paper, based on B-splines function and FGLS regression, we propose a weighted online regularization method which apply weight coefficients and a relatively small number of knots for interpolation to reduce the effect of heteroscedasticity and computation cost, and exchange the amount of data for the accuracy of reconstruction. Moreover, this algorithm has the characteristics of smaller memory occupation and can process online data, which means that when new data is added, there is no need to re-process the processed data.

The paper is organized as follows. In Section 2, we introduce the issue which is to be studied and propose an online reconstruction algorithm. Error estimation is carried out in Section 3. Section 4 provides a principle of optimal selection for weight coefficients and the regularization parameter in objective function. The performance of the proposed algorithm is illustrated in Section 5. In Section 6, we present two applications of our reconstruction algorithm in inverse problems. Section 7 concludes main results.

2.   Problem statement and algorithm design

In this section, we consider the following problem: for a lager positive integer $N$, we try to reconstruct function $f(x)$ and its derivative $f'(x)$ in the model (1.1) from observation data $\{(x_i, y_i)\}^N_{i = 1}$ and function $\sigma(x)$.

2.1.   Regularized solution

First, we define the Tikhonov objective functional as follows:

where $g(x)\in W^{2, 2}(0, 1)$, $\{w_i\}^N_{i = 1}$ are weight coefficients and $\alpha$ is a regularization parameter. The minimization problem is reformulated as follows:

$f_N$ is defined as the approximated solution of $f(x)$ in (1.1).

In this paper, we reconstruct not only function $f(x)$ but also its derivative $f'(x)$. This is a classical ill-posed problem. It is computationally costly to solve the problem with local regression methods such as kernel regression and local polynomial regression. This problem is usually solved by regularization method and spline functions. Moreover, B-spline basis functions have compact support which makes it possible to speed up calculations. Thus, we choose B-spline function.

Then we solve suitable $g(x)$ in the finite dimensional function space of cubic B-spline function instead of the infinite dimensional space. Inspired by the method in ^[25], we just select some equidistant points not sample points as interpolation node.

Let $M$ be a positive integer and mesh size $d = M^{-1}$. Equidistant knots $\{p_j\}^M_{j = 0}$ are defined as

The vector space of all cubic B-spline functions with knots $\{p_j\}^M_{j = 0}$ is called space $\mathcal{V}_m$. Assume function $g(x)$ belongs to space $\mathcal{V}_m$, $g(x)$ can be written as

where $\{\lambda_j\}^{M+1}_{j = -1}$ are constants and $\phi_j(x) = \phi(\frac{x-p_j}{d})$. $\phi(x)$ is standard cubic B-spline function defined in ^[24].

Denote the weight coefficient $W = diag(\sqrt{w_1}, \sqrt{w_2}, \cdots, \sqrt{w_N})$, the noisy sample $\mathit{\boldsymbol{y}} = (y_1, y_2, \cdots, y_N)^T$, the function parameter $\mathit{\boldsymbol{\lambda}} = (\lambda_{-1}, \lambda_0, \cdots, \lambda_{M+1})^T$ and the row vector $H_x = (\phi_{-1}(x), \phi_{0}(x), \cdots, \phi_{M+1}(x))$. Through the Eq (2.3) and the definition above, the Eq (2.2) can be rewritten as

where $H: = (H_{x_1}, H_{x_2}, \cdots, H_{x_N})^T$ and $P\in \mathbb{R}^{(M+3)\times(M+3)}$ is given by

When $\mathit{\boldsymbol{\lambda}}_N$ is known, the approximated solution $f_N(x) = H_x\mathit{\boldsymbol{\lambda}}_N$.

Theorem 1. Suppose $N > 2$, and the observation points $\{x_i\}^N_{i = 1}$ are not identical. Then, the minimization problem (2.4) has a unique minimizer $\mathit{\boldsymbol{\lambda}}_N$ which satisfies:

Proof. Since $J_1(\mathit{\boldsymbol{\lambda}})$ in (2.4) is a quadratic form with respect to $\mathit{\boldsymbol{\lambda}}$, the derivative of $J_1(\mathit{\boldsymbol{\lambda}})$ takes zero only at $\mathit{\boldsymbol{\lambda}} = \mathit{\boldsymbol{\lambda}}_N$, that is

Moreover, $P$ is positive definite and $H^TW^TWH$ is positive semidefinite. Thus, $\mathit{\boldsymbol{\lambda}}_N$ is unique. □

2.2.   Algorithm

Through the theorem above, we propose a online Tikhonov regularization algorithm in Algorithm 1. In "Require" sentence, we input the observation data $\{(x_i, y_i)\}^N_{i = 1}$ and the necessary parameters–the number of knots $M$, weight coefficients $W$ in (4.6) and regularization parameter $\alpha$ in Table 1. The algorithm proceeds from Line 1 to Line 8. In Line 1–2, we initialize matrix $A_0 = 0 \in \mathbb{R}^{(M+3)\times(M+3)}$ and $b_0 = 0 \in \mathbb{R}$ and generate the matrix $P\in \mathbb{R}^{(M+3)\times(M+3)}$. Line 3–6 provides a method to update $A$ and $b$ in (2.5) according to the newly entered data. Through solving linear system in Line 7, we have $\mathit{\boldsymbol{\lambda}}_N$ and the approximated solution $f_N$ satisfies $f_N(x) = H_x\mathit{\boldsymbol{\lambda}}_N$.

Remark 1. When a new data is inputed, one just run the algorithm from Line 4 to Line 6. The computational complexity at Line 4, 5 and 6 are $\mathcal{O}(1)$ and the computational complexity at Line 7 is $\mathcal{O}(M)$. Moreover, the total data storage of Algorithm 1 is $\mathcal{O}(M)$.

Weight coefficients W are introduced to reduce the effect of heteroscedasticity. The selection of weight coefficients W will be thoroughly discussed in Section 3.

Note that, $N$ continues to increase with the influx of new data, when $N$ becomes so large that $M\ll N^{\frac{1}{5}}$, we need to increase $M$ to $M\approx N^{\frac{1}{5}}$ and restart the algorithm.

3.   Error estimation

In this chapter, we analyze the reconstruction error of function $f(x)$. Some assumptions are necessary before proving. Suppose that the matrix $P$ is positive definited. Without losing generality, it is advisable to set the domain of variable $x$ to $[0, 1]$. The weighting coefficient matrix $W$ can be redefined as

Make $w_ {min}$ and $w_ {max}$ are the minimum and maximum values of parameters $\{w_j\}_{j = 1}^N$ respectively. Then for the observation data $\{(x_i, y_i)\}^N_{i = 1}$, $y_i$ can be expressed as

where $\varepsilon_j$ are independent with expectation 0 and variance 1.

Recording error function as

Next, consider the influence of regularization function on function itself and random noise.

When the data is divided into deterministic part and random part, the reconstruction result can also be divided into two parts

where the deterministic part $f_ {abs}$ satisfies

here $\mathit{\boldsymbol{\varepsilon}}$ is a column vector of the observed noise

let $\sigma^2_* = \max_{x_j}\sigma^2(x_j)$. And the random part $f_{noise}$ satisfies

Then, the error function can be rewritten as

Two important lemmas are introduced before error estimation.

Firstly, the indicator function is defined to describe the distribution of observation points.

Definition 1. Divide the interval $[0, 1]$ into $M$ disjoint subspaces $I_i$, $1\leq i\leq M$,

where $d = 1/M$ is the grid spacing. There are $N_i$ observation data on each subspace $I_i,$ 1 $\leq i\leq M$. Define the indicator function $\rho_M(x)$ as

Lemma 1. (^[10], Lemma 3.3) Assume that the indicator function $\rho_M(x)$ has an upper bound $\rho_{u}$ in the interval $[0, 1]$,

Then, for the first order continuous differentiable function $f (x), \; x\in [0, 1]$, the mean square sum of values at all observation points satisfies:

Next, the error estimation of cubic spline interpolation is introduced.

Lemma 2. (^[24], Theorem 1.56) Assume the objective function $f(x)\in H^2(a, b), \ x_0 = a, \ x_M = b$, the grid spacing is $d$. Let $s_f$ be a cubic spline interpolation function of f with natural or fixed boundary conditions, then the interpolation error can be estimated as

In addition, if $s_f$ has natural boundary conditions, then the relation between the second derivatives of $f$ and $s_f$ is

Now, the error estimate of the deterministic part is analyzed.

Theorem 2. Assume that $\rho_M$ is an indicator function on the interval $[0, 1]$ and has an upper bound $\rho_ u$. $f(x) \in H^2(0, 1)$. Then, the mean square sum of $e_{abs}(x)$ at the observation points is estimated as follows:

Proof. The deterministic part. Note

where $s_f(x)$ is a cubic spline interpolation function satisfies natural or fixed boundary conditions on the interval $[0, 1]$. By Lemmas 1 and 2, we have

And

For $f_{abs}$ is the smallest element of the weighted regular functional $J(\cdot; \alpha, \mathit{\boldsymbol{y}}-\mathit{\boldsymbol{\varepsilon}})$,

Then,

From Eqs (3.3)–(3.5), we have

On the other hand, by Lemma 2,

Combining (3.2) and (3.6),

□

Remark 2. According to formula (3.7), in order to control the error effectively, it is necessary to select the same order of regularization parameter $\alpha$ as $d^4$ to control the interpolation error.

Next, consider the error estimation of the random part. The reconstruction result of the random part is

By the conclusion in chapter 4, we have

where $\mathit{\boldsymbol{\lambda}}^{\varepsilon}$ satisfies

Theorem 3. Under the above assumptions, the mean square sum of the random part $f_{noise}$ on the observation points is satisfied with the probability of $1-\beta$

Proof. The proof is divided into three steps.

Step 1: Rewrite $\mathit{\boldsymbol{\lambda}}^{\varepsilon}$ as

where $I$ is an identity matrix of $N\times N$.

Step 2: Feature decomposition $S$. For $S$ satisfying:

Obviously, it is semi-positive definite, thus $S$ can be rewritten as

where $U$ is an orthogonal matrix and $T$ is a diagonal matrix composed of the eigenvalues of $S$. For

$T$ can be rewritten as

where $\{t_j\}_{j = 1}^{M+3}$ is a monotonically decreasing sequence,

Step 3: Obtaining estimates in the form of confidence intervals.

For the mean square sum of $f_{noise}$ on observation points satisfies:

And

Combined with Markov inequality: for non-negative variable $X$, for any given $a > 0$, the following formula holds

Therefore, the estimation of noise mean square sum error of noise under the probability of $1-\beta$ satisfies:

□

Finally, we estimate the error $||e_N||^2_{L^2(0, 1)}$. Through Theorems 2 and 3, the errors of the deterministic and random parts of the reconstruction function are analyzed. Moreover, we can obtain the mean square error of $f(x)$ under the probability of $1-\beta$.

Theorem 4. Assume that $\rho_M$ is an indicator function on the interval $[0, 1]$ and has an upper bound $\rho_ u$. Let $f(x)\in H^2(0, 1)$, for any $\beta\in [0, 1]$, such that the error $e_N$ under the probability of $1-\beta$ satisfies:

Proof. By Theorems 2 and 3, the conclusion is easy to get. □

4.   Parameters selection

From inequality (3.10), the estimation error is determined by the model parameters. Thus, it is a very important problem to select appropriate model parameters. In this chapter, we are going to determine the regularization parameter $\alpha$, the number of knots $M$, and weight coefficients $W$ in Algorithm 1.

First, we consider the regularization parameter $\alpha$ and the number of knots $M$. Through Theorem 4, the three formulas

must be in the same order of magnitude. Otherwise, it will result in oversize error $||e_N||^2_{L^2(0, 1)}$ or over-fitting. Assume the value of $\rho_u$ and $||f''||^2_{L^2(0, 1)}$ do not effect the order of magnitude. Thus, we have

where $d = 1/M$. Through the equation above, the regularization parameter $\alpha$ and the number of knots $M$ satisfy:

Next, we consider the optimal choice of weight cofficients $W$. For weighted linear regression, researchers provide a method to choose weight cofficients. Thus, based on the method of weighted linear regression and the linear form of the reconstruction results, we provide Theorem 3.1 to choose weight coefficients.

Let $e_N = f_N - f$ be the error function of the proposed regularization algorithm. From error function $e_N$ and Eq (2.5), the error comes from two aspects: the interpolation error caused by pre selected interpolation nodes and the random error caused by observation noise. The error $e_N$ can not be reduced easily. But we can select suitable weight coefficients $W$ to reduce the variance of the error $e_N$. Denote $f^*(x)$ as a best approximated solution of $f(x)$ in space $V_m$ as follows:

where the error $e^* = f^*-f$ just comes from the interpolation error. Similar the process of Theorem 1, $f^*$ can also be written as $f^*(x) = H_x\mathit{\boldsymbol{\lambda}}^*_N$ and $\mathit{\boldsymbol{\lambda}}^*_N$ is given by

where $\mathit{\boldsymbol{f}}: = (f(x_1), f(x_2), \cdots, f(x_N))^T$.

Theorem 5. Assume regularization parameter $\alpha\ll 1$, $M$ is large enough and variance function $\sigma(x)$ satisfies $\sigma:[0, 1]\rightarrow (0, +\infty)$. If $W^*$ is given by

then, $\mathit{\boldsymbol{\lambda}}_N(W^*)$ is the approximated minimum empirical variance unbiased estimation of $\mathit{\boldsymbol{\lambda}}^*_N$, that is, $\mathit{\boldsymbol{\lambda}}^*_N\approx \mathbb{E}[\mathit{\boldsymbol{\lambda}}_N(W^*)]$ and, for any $W\in \mathbb{R}^{N\times N}$,

Proof. $\mathit{\boldsymbol{f}}$ can be written as

where $\mathit{\boldsymbol{\varepsilon_1}}$ is the interpolation error which tends to zero, when $M$ large enough. From the equation above, Eqs (2.5) and (1.1), we have

Denote $\mathit{\boldsymbol{\varepsilon}} = (\sigma(x_1)\varepsilon, \sigma(x_2)\varepsilon, \cdots, \sigma(x_N)\varepsilon)$, diagonal matrix $D: = W^TW$ and diagonal matrix $V: = diag(\sigma^2(x_1), \sigma^2(x_2), \cdots, \sigma^2(x_N))$. Through the definition above and Eqs (2.5), (1.1) and (4.5), for any $W\in \mathbb{R}^{N\times N}_+$, if $\varepsilon_1\rightarrow 0$ and $\alpha\rightarrow 0$, we have

Moreover, when $W = W^*$,

□

From Theorem 2, when $\alpha \rightarrow 0$ and $M$ is large enough, we can choose $W^*$ to reduce the variance of estimator $\mathit{\boldsymbol{\lambda}}_N(W)$. Moreover, the value $\{1/\sigma(x_i)\}^N_{i = 1}$ in $W^*$ can not be infinity. When some weight coefficients $1/\sigma(x_i)$ are too large, the reconstruction results are only affected by the data corresponding to these large weight coefficients. Thus, we suggest to choose the suitable weight coefficients $W$ as follows:

where $k$ is a small positive constant such as 0.1 and 0.01 and $w_{i}$ satisfies:

Combining Eqs (4.1) and (4.7), we have the parameter selection Table 1.

5.   An illustrative example

In this section, we provide a illustrative example to show the feasibility and advantages of the weighted regularization method. Assume that the heteroscedastic model satisfies:

where $f(x): = 40(x-0.5)^2$ and $\sigma(x): = 3\sin(150\pi x)+3.5$. We consider following problem: for a lager positive integer $N = 300$, given observation data $\{(x_i, y_i)\}^N_{i = 1}$ and function $\sigma(x)$, we try to reconstruct function $f(x)$ and its derivative $f'(x)$. We select the number of knots $M = 10$, mesh size $d = 1/M = 0.1$, regularization parameter $\alpha = 10^{-4}$ and weight coefficients $W$ satisfying (4.6) with $k = 0.1$.

Figure 1 compares the reconstruction results-"improve" in this paper with the results-"tradition" in ^[10]. When function $\sigma(x)$ changes sharply, our results are closer to exact function than that those in ^[10].

6.   Applications in inverse problems

In this section, we will give two applications of our numerical differentiation method in two kinds of inverse problems.

6.1.   The inverse problems of identifying coefficients in the boundary value problem

We consider a nonlinear ill-posed problem–the identification of the coefficient $c$ in the boundary value problem

from the solution $f(x)$. We can get $c(x)$ directly from

The example we take here is that $f(x) = 10e^{-\sin(\pi x)}$, $p(x) = 0$ and $h_0 = 10$. In our computation, $\{(x_i, y_i)\}^N_{i = 1}$ is denoted as the observation data where $y$ satisfies:

where $\sigma(x) = 5 \sin(150\pi x) + 5.1$. From the observation data above and the improved method in this paper, we reconstruct $f(x)$ and $f'(x)$. The results are present in Figure 2.

Through the reconstruction results of $f(x)$ and $f'(x)$ and Eq (6.1), we have the estimator of $c(x)$ in Figure 3, where our method is better than the traditional.

6.2.   Determination of the interface in computerized tomography

In this subsection, we apply the method in this paper to a simple but interesting problem in Computerized Tomography (CT).

We assume that the attenuation coefficient of an object with respect to an X-ray at $x$ is $p(x)$. We scan the cross-section by an X-ray beam L of unit intensity. The intensity past the object is $e^{-\int_Lp(x)dx}$. We denote function $f$ as follows:

The main problem in CT is to recover the function $p$ from $f$. However, in many cases, it will be enough if one can reconstruct the interface of the different mediums, that is, the discontinuous points of $p$ which is related directly with the nondifferentiable points of $f$. Thus we just need to determine the nondifferentiable points in $f(x)$.

$D$ is denoted as the triangular cross section of a object. The attenuation coefficient is $0$ outside the object and $1$ inside the object, that is:

see Figure 4.

From Figure 4 and Eq (6.2), the function $f(x)$ can be calculated directly:

In our computation, $\{(x_i, y_i)\}^N_{i = 1}$ is denoted as the observation data where $y$ satisfies Eq (1.1) and $\sigma(x) = 5 \sin(150\pi x) + 5$. Let the number of data $N = 300$, the number of knots $M = 50$, mesh size $d = 1/M = 0.02$, regularization parameter $\alpha\approx 10^{-4}$ and the parameter $k = 0.01$ in (4.6).

In Figure 5, we present observation data, noise data and the reconstructions of $f(x)$ and $f''(x)$. Compared with the traditional method in ^[10], through the improved approach in this paper, we can reconstruct $f(x)$ better and find three nondifferentiable points of $f(x)$: $0.2$, $0.5$ and $0.8$.

7.   Conclusions

In this paper, we propose a weighted online regularization method–Algorithm 1, which updates reconstruction results through new input data without old data again. The error estimation in Section 3 illustrates that the algorithm can reach the best convergence order when the relationship between the amount of observed data $N$ and the dimension $M$ of linear space is $M \approx N^{\frac{1}{5}}$. Compared with the classical B-spline method, weight coefficients and a relatively small number of knots for interpolation are used to reduce the effect of heteroscedasticity and computation cost. Moreover, a feasible selection method of the regularization parameter, the number of knots, and weight coefficients are provided in Table 1.

The biggest advantage of this algorithm is to solve the problem of low data accuracy caused by heteroscedasticity by processing a large amount of data. It can significantly reduce the error of observation data while ensuring the calculation efficiency.

In the next work, we will consider the higher dimensional heteroscedasticity problem, which requires higher regularity of the function and analyze other radial basis functions, such as Sobolev radial basis functions, to solve high-dimensional problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by Doctoral Research Initiation Fund of Shenzhen Institute of Information Technology with grant (No. SZIIT2023KJ003).

Conflict of interest

The author declare that have no competing interest.