Hermitian and skew-Hermitian splitting method on a progressive-iterative approximation for least squares fitting

1.   Introduction

Geometric iteration techniques (GI) are a subcategory of iterative linear algebraic methods for solving linear equations that have gained popularity as a data-fitting tool in recent years. Since its introduction, geometric iteration approaches in geometric design, including computer-aided design (CAD), differential equations, statistics, and approximation theory, have been widely used in many areas of research and practical applications. The progressive and iterative approximation (PIA) method presented by Lin et al. ^[1] is an adaptable and stable approach to data fitting that belongs to the geometric iteration methods category. The PIA method can be applied to interpolate and approximate curves and surfaces with normalized totally positive (NTP) bases.

Moreover, the EPIA method ^[2] is extended to approximate a given data points set. The EPIA organizes the data points into groups with equal to the amount of control points as there are control points. Each group represents a control point and then executes PIA-format iterative approximation. However, the limit of the generated curves and surfaces sequence is not the least square fitting (LSF) result to the data points set. Then, Deng and Lin proposed the LSPIA method ^[3], in which the least square fitting (LSF) result to the data set with the B-spline basis is the limit of the generated curves and surfaces sequence. Ebrahimi and Loghman ^[4] presented a composite iterative procedure for the progressive and iteration approximation for generates a sequence of matrices based on the Schulz iterative method which it uses in the adjusting vectors. In addition, Lin et al. ^[5] also gives an review of interpolation and approximation geometric iteration methods, as well as local characteristics and accelerating methods, as well as demonstrating convergence.

Recently, Wang ^[6] presented an extended the progressive and iterative approximation for least square fitting (ELSPIA) approach with global and local relaxation parameters, which converges faster than the LSPIA method ^[3] in terms of iterative steps and computing time. In this article, we propose the Hermitian and skew-hermitian splitting methods for least squares fitting on a progressive-iterative approximation, can theoretically reduce the approximation error to zero.

The PIA ^[1] is an iterative approach for solving a system of linear equations $BP = Q$, as we know, where $B = \{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ is a collocation matrix with $0 = t_1 < t_2 < t_3 < \dots < t_m = 1$ and $n = m$, $\{Q_{j}\}_{j = 1}^{m}$ refers to a collection of data points and $\{P_{i}\}_{i = 1}^{n}$ refers to a collection of control points. When all of the entries in each row of the blending basis matrix $B$ add up to 1, $B$ is normalized. Consider that if all of a matrix's minors are non-negative, it is totally positive. The basis is considered to be totally positive if all of the collocation matrices are totally positive. As a result, the blending basis matrix $B$ is normalized totally positive (NTP), and the NTP basis is the corresponding basis. The NTP bases have traditionally been related to shape-preserving representations ^[7]. The majority of curve and surface representations used in computer-aided design, such as the Bernstein basis and the B-spline basis, are based on NTP bases.

Using the Hermitian and skew-Hermitian splits, Bai et al. ^[8] presented an iteration approach (HSS). Hu et al. ^[9] developed a new PIA iterative algorithm and its weighted form for progressive iterative interpolation of data points based on the HSS-iteration approach, using NTP bases, namely HPIA and WHPIA, respectively. These two approaches' iterative processes formally consist of two iterations, with the iterative distinction vectors in the two iterations differing.

As we know, the LSPIA ^[3,10] and ELSPIA ^[6] are the iteration approaches for solving a system of linear equations $\mu B^{T}BP = B^{T}Q$, where $B = \{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ is a collocation matrix with $0 = t_1 < t_2 < t_3 < \dots < t_m = 1$ and $m \geq n$, $\mu$ is constant satisfying the condition $0 < \mu < \frac{2}{\lambda_{0}}$ with $\lambda_{0}$ is the largest eigenvalue of $B^{T}B$, $\{Q_{j}\}_{j = 1}^{m}$ refers to a collection of data points and $\{P_{i}\}_{i = 1}^{n}$ refers to control points.

According to theoretical analysis ^[8], the HSS-iteration consistently converges unconditionally to the exact solution of the system of linear equations $Ax = b$, where $A \in \mathbb{C}^{n \times n}$ and non-singular. Using the HSS-iteration approach, we propose a new procedure, called the HSS-LSPIA method, for solving a system of linear equations of LSPIA such that $A = \mu B^{T}B$, $x = P$ and $b = B^{T}Q$, where $B = \{B_{i}(t_{j})\}_{i, j = 1}^{n, m} \in \mathbb{R}^{m \times n}$ is a collocation matrix with $m \geq n$.

The paper is organized as the following: Section 2 presents preliminaries. Section 3 shows the iterative format of the HSS-LSPIA method schemes for curves and surfaces with NTP bases respectively. And proves the convergence of HSS-LSPIA in the case of curves and surfaces in this section. In Section 4, we present the experimental results of five examples show that the HSS-LSPIA method is more effectively efficient than the LSPIA ^[3], ELSPIA ^[6] and WHPIA^[9]. Finally, in Section 5, some conclusions are summarized.

2.   Preliminaries

2.1.   LSPIA and ELSPIA

For any blending curves and surfaces fitting, we'll go through the progressive iterative approximation property. For the ordered point set $\{Q_{j}\}_{j = 1}^{m}$, select control point set $\{P_{i}\}_{i = 1}^{n}$ from $\{Q_{j}\}_{j = 1}^{m}$ with $m \geq n$, let $\{Q_{j}\}_{j = 1}^{m}$ be the assigned parameters to $\{Q_{j}\}_{j = 1}^{m}$, choose $B = \{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ as an NTP basis. For the ELSPIA method^[6], let $\gamma^{0}_{j} = 0, \quad j = 1, 2, \dots, m,$ and $P^{0}_{i} = P_{i}, \quad i = 1, 2, \dots, n$ be weighted sum and the initial control points, respectively, we can obtain the $(k + 1)$th blending curve,

We update the control points

by adjusting vector

with the inner iteration

and the difference vector

The parameters $\omega$ and $\mu$ both satisfy the condition in this case

where $\sigma_{0}$ stands for the collocation matrix B's largest singular value.

When the parameters $\omega$ and $\mu$ are satisfy as to the following

the LSPIA method ^[3] is a special case of the ELSPIA method for blending curve fitting.

In the case of the surface, two NTP bases are required. In both directions $t$ and $v$, a Kronecker product is formed by $\{B_i(t)\}_{ i = 1}^{n_{1}}$ and $\{N_j(v)\}_{ j = 1}^{n_{2}}$. We given the initial control points $\{P_{ij}^{0}\}_{i, j = 1}^{n_{1}, n_{2}}$, with a real increasing parameters set $\{t_{h}^{0}\}_{h = 1}^{m_{1}}$, i.e. $0 = t_1 < t_2 < t_3 < \dots < t_{m_1} = 1$ and $\{v_{l}^{0}\}_{l = 1}^{m_{2}}$, i.e. $0 = v_1 < v_2 < v_3 < \dots < v_{m_2} = 1$. We select $\{P_{ij}\}_{i, j = 1}^{n_{1}, n_{2}}$ from the ordered point set $\{Q_{ij}\}_{i, j = 1}^{m_{1}, m_{2}}$. The sequence of surface iterations $C^{(k+1)}(t, v)$ is related to that of the curve case with $k \geq 0$;

We update the control points

by adjusting vector

with the inner iteration

and the difference vector

with $i = 1, 2, \dots, m_{1}$ and $\ j = 1, 2, \dots, m_{2},$ where $\mu$ is the global relaxation parameter and $\omega$ is the local relaxation parameter satisfy

where $\sigma_{B0}$ and $\sigma_{N0}$ are the largest singular values of the collocation matrix B and N, which correspond to the bases $\{B_{i}(t)\}_{i = 1}^{n_{1}}$ and $\{N_{j}(v)\}_{j = 1}^{n_{2}}$, respectively.

When the parameters $\omega$ and $\mu$ are satisfy as to the following

the LSPIA method ^[3] is a special case of the ELSPIA method for blending surface fitting.

Let

and

According to updating the control points in the case of curve (2.2), we revise the iterative form with (2.3)–(2.5) as the matrix iterative form in the same matrix form as updating the control points in the case of surface (2.9). We then rewrite with (2.10)–(2.12) as follows

2.2.   HSS-iteration

Many problems in scientific computing lead to a system of linear equations

with $B$ a large sparse non-Hermitian and positive definite matrix.

Iterative methods for the system of linear equations (2.16) require efficient splittings of the coefficient matrix $B$. For example, the Jacobi and the Gauss-Seidel iterations ^{[11,12,13,14,15,16]} split the matrix $B$ into its diagonal and off-diagonal part. Because the matrix $B$ naturally possesses a Hermitian/skew-Hermitian splitting (HSS) ^[17,18]

where

Bai et al. presented the Hermitian/skew-hermitian splitting (HSS) iteration approach ^[9] as a new way to solve a system of linear equations (2.16). We'll pretend that we've recieved the iteration sequence until ${P^{k+1}}$ ($k \geq 1$) converges. Then, the HSS iteration can compute from the format as

where $I$ is the identity matrix and $\alpha$ is a positive constant.

They have also proved this method converges unconditionally to the exact solution of the system of linear equations (2.16). The upper bound of the contraction factor of the HSS iteration is dependent on the spectrum of the Hermitian part $H$ but is independent of the spectrum of the skew-Hermitian part $S$ as well as the eigenvectors of the matrices $H$, $S$, and $B$.

2.3.   WHPIA iteration

Hu et al. ^[9] presented a new PIA iterative technique and its weighted edition for progressive iterative interpolation of data points based on the HSS-iteration approach, using NTP bases called HPIA and WHPIA, respectively. The WHPIA method typically converges faster than the HPIA method in the majority of cases. Following that, we briefly review the WHPIA' iterative processes, which formally consist of two iterations with different iterative distinction vectors in the two iterations, and set up a function based on NTP bases like a perturbation term in the iteration process.

The case of a curve with an NTP basis $\{B_{i}(t)\}_{i = 1}^{m}$ and a collection of data points to be interpolated $\{Q_{i}\}_{i = 1}^{m}$. Then, with a real increasing parameters set $\{t_{i}\}_{i = 1}^{m}$ i.e. $0 = t_1 < t_2 < t_3 < \dots < t_{m} = 1$, use these data points to be the initial control points $\{P_{i}^{0}\}_{i = 1}^{m}$ and $\{P_{i}^{\frac{1}{2}}\}_{i = 1}^{m}$. The initial blending curve can be formed as follows:

and

Then, the sequence of curves $C^{k}(t)$ for $k \geq 1$ as the following

where

In (2.23), $C^{k}(B_{i})$ is defined as follows

which is a function that receives the blending bases as variables.

Surface fitting with two NTP bases $\{B_{i}(t)\}_{i = 1}^{n_{1}}$, $\{N_{j}(v)\}_{j = 1}^{n_{2}}$, and a collection of data points $\{Q_{ij}\}_{i, j = 1}^{n_{1}, n_{2}}$ to be interpolated. Then, using these data points to be the initial control points $\{P_{ij}^{0}\}_{i, j = 1}^{n_{1}, n_{2}}$, with a real increasing parameters set $\{t_{i}^{0}\}_{i = 1}^{n_{1}}$, i.e. $0 = t_1 < t_2 < t_3 < \dots < t_{n_1} = 1$ and $\{v_{j}^{0}\}_{j = 1}^{n_{2}}$, i.e. $0 = v_1 < v_2 < v_3 < \dots < v_{n_2} = 1$. The initial blending surface can be created as follows:

Then, the sequence of surfaces $C^{k+1}(t, v)$ for $k \geq 1$, can be compute as follows

where

We can use the following expression in (2.26) as a perturbation term in the iteration process, even as we appear to have done with a curve;

The WHPIA format of curves and surface, which involves two iterations and substituting the sequential iterative process length of each control point in PIA, is referred to as (2.23) and (2.26), respectively.

3.   Iterative format of HSS-LSPIA

As we know, the PIA ^[1] is an iterative approach for solving a system of linear equations $BP = Q$, $B \in \mathbb{C}^{m \times n}$ with $n = m$. The matrix form of PIA method as

and the LSPIA ^[3,10] and ELSPIA ^[6] are the iteration approaches for solving a system of linear equations $\mu B^{T}BP = B^{T}Q$, $B \in \mathbb{C}^{m \times n}$ with $m \geq n$. The matrix form of LSPIA method as

For an NTP basis $\{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ and a data set points $\{Q_{j}\}_{j = 1}^{m}$ to be fitted in either $\mathbb{R}^{m \times2}$ or $\mathbb{R}^{m \times 3}$. We use all these data points to be the initial control points $\{P_{i}^{0}\}_{i = 1}^{n}$ in $\mathbb{R}^{n \times2}$ or $\mathbb{R}^{n \times 3}$ and combine the HSS-iteration with the LSPIA method for solving a system of linear equations

The matrix form of HSS-LSPIA method as the following

where

and

We can rewrite (3.4) as

Then, we have

and

By (3.8) and (3.9), we get

So,

Let

and

then

Remark 3.1. Since the skew-Hermitian part $S = 0$, HSS-LSPIA can called shifted splitting LSPIA ^[19].

3.1.   The curve case

Considering an NTP basis $\{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ with $m \geq n$ and a data set points $\{Q_{j}\}_{j = 1}^{m}$ to be fitted in either $\mathbb{R}^{m \times 2}$ or $\mathbb{R}^{m \times 3}$. The initial blending curve ${C^{0}(t)}$ and ${C^{\frac{1}{2}}(t)}$ can then be generated using all these data points to be the initial control points $\{P_{i}^{0}\}_{i = 1}^{n}$ in $\mathbb{R}^{n \times2}$ or $\mathbb{R}^{n \times 3}$ and $\{P_{i}^{\frac{1}{2}}\}_{i = 1}^{n}$, with a real increasing parameters set $\{t_{i}\}_{i = 1}^{m}$ i.e. $0 = t_1 < t_2 < t_3 < \dots < t_{m} = 1$. The initial blending curve can be formed as follows:

and

Then, the remaining curve of the sequence $C^{k+1}(t)$ for $k \geq 0$, can be compute as following:

where

In (3.16), $C^{k+\frac{1}{2}} (t_{j})$ is defined as the following:

Both global and local relaxation parameters, $\mu$ and $\alpha$ are satisfied as

where $\lambda_{max}$ denotes the largest eigenvalue of $B^{T}B$. If (3.16) converge, then

the HSS-LSPIA format of curves, including of two iterations and substituting the sequential iterative process from each control point in LSPIA, is recognized as (3.16).

Remark 3.2. From (3.16), we will get $\lim_{k \to \infty} C^ {k +1}(t_{i})$ exist, if for any $\epsilon > 0$ there is a natural number $A$, when $k > A$, $\| P^ {k +1}_{i} - P^{k}_{i} \| < \epsilon$. Furthermore, whether the spectral radius of the iterative matrix $P^ {k +1}$ in (3.11) is less than $1$ affects convergence.

3.2.   Convergence analysis of the HSS-LSPIA curve fitting

Lemma 3.3. ^[1] Given an NTP basis $\{B_{i}(t_{j})\}_{i, j = 1}^{n, m}$ with $m \geq n$ and $B^TB$ is nonsingular collocation matrix. And assuming that $\lambda_{i}(B^TB), \ i = 1, 2, \dots, n$ is eigenvalue of $B^TB$, then $0 < \lambda_{i}(B^TB) \leq 1$.

Let $B \in \mathbb{C}^{m \times n}, \ B = M_{i} - N_{i} \ (i = 1, 2)$ be two splitting of the matrix $B$. In the matrix form, Eq (3.16)'s two iterative processes can be written as

where $M_{1} = (\alpha I + \mu B^{T}B), \ N_{1} = \alpha I, \ M_{2} = \alpha I, N_{2} = (\alpha I - \mu B^{T}B)$, $I$ is identity matrix. Then we can rewrite (3.20) as a unified matrix iteration format as follows

or

where $M = M_{2}^{-1} N_{2} M_{1}^{-1} N_{1}$ and $M_{0} = M_{2}^{-1} (I +N_{2} M_{1}^{-1})$.

Now we have to show that the iterative control point sequence $\{P^{k+1}\}$ leads to the unique solution $P^{*}$, i.e. $\rho (M) < 1$.

Theorem 3.4. Let $B \in \mathbb{C}^{m \times n}$ with $m \geq n$ be a positive definite matrix, let $H = \frac{\mu}{2} \left(B^{T}B + \left(B^{T}B \right) \right) = \mu B^{T}B$ and $S = \frac{\mu}{2} \left(B^{T}B - \left(B^{T}B \right) \right) = 0$ be its Hermitian ans skew-Hermitian part, and given $\alpha$ be a positive constant. The HHS-LSPIA iteration's iteration matrix $M(\alpha)$ is then provided by

or

and the spectral radius of $\hat{M}(\alpha)$ is $\rho(\hat{M}(\alpha))$, bounded by

where the spectral set of the matrix $H$ is $\lambda(H)$.

As a consequence, it maintains that

i.e., the HHS-LSPIA iteration converges to the unique solution $P^{*}$ of the linear equation (3.3).

Proof. By the resemblance invariance to the matrix spectrum, we have

We know that $(\alpha I)(\alpha I)^{-1} = I$ which is unitary matrix i.e., $|| (\alpha I)(\alpha I)^{-1} || _{2} = 1$. Then, we get

Since $\lambda_{i} > 0 \ (i = 1, 2, 3, \dots, n)$ and $\alpha$ is positive constant, it is easy to see that $\rho(M(\alpha)) \leq \rho(\alpha) < 1$.

Corollary 3.5. Under the same assumption of Theorem 3.4. Let $\lambda_{max}$ and $\lambda_{min}$ be the maximum and the minimum eigenvalues of $B^{T}B$, respectively. Then, we have the approximate fastest convergence rate of HHS-LSPIA methed for the curve fitting when

where $\mu = \frac{2}{ \lambda_{max} + \lambda_{min} }$ presented in ^[3].

Proof. In ^[8], it is shown that the optimal spectral radius. Now,

To compute an approximate optimal $\alpha > 0$ such that the convergence factor $\rho(M(\alpha))$ of the HSS-LSPIA iteration is minimized, we can minimize the upper bound $\sigma(\alpha)$ of $\rho(M(\alpha))$ instead. If $\alpha^{*}$ is such a minimum point, then it must satisfy $\alpha^{*} - \lambda_{min} > 0, \alpha^{*} - \lambda_{max} < 0$, and

Therefore,

or

Remark 3.6. From (3.16) and Corollary 3.5, we will get $\frac{\mu}{\alpha} = \frac{\mu}{ \mu \sqrt{ \lambda_{max}(B^{T}B) \lambda_{min}(B^{T}B)}} = \frac{1}{ \sqrt{ \lambda_{max}(B^{T}B) \lambda_{min}(B^{T}B)}} = \omega$. Then, we use the $\omega = \frac{1}{ \sqrt{ \lambda_{max}(B^{T}B) \lambda_{min}(B^{T}B)}}$ for compute in (3.16) instead of $\frac{\mu}{\alpha}$.

3.3.   The surface case

Let two NTP bases $\{B_{i}(t_{i'})\}_{i, i' = 1}^{n_{1}, m_{1}}$ $\in$ $\mathbb{C}^{m_{1} \times n_{1} }$, $\{N_{j}(v_{j'})\}_{j, j' = 1}^{n_{2}, m_{2}}$ $\in$ $\mathbb{C}^{m_{2} \times n_{2}}$, and a data set points $\{Q_{i'j'}\}_{i', j' = 1}^{m_{1}, m_{2}}$ $\in$ $\mathbb{R}^{(m_{1} \times m_{2}) \times 3}$ to be fitted, where $m_{1} \geq m_{2} \geq n_{1} \geq n_{2}$. The initial blending surface ${C^{0}(t, v)}$ and ${C^{\frac{1}{2}}(t, v)}$ can then be generated using all these data points to be the initial control points $\{P_{ij}^{0}\}_{i, j = 1}^{n_{1}, n_{2}}$ $\in$ $\mathbb{R}^{(n_{1} \times n_{2}) \times 3}$ and $\{P_{ij}^{\frac{1}{2}}\}_{i, j = 1}^{n_{1}, n_{2}}$, with a real increasing parameters set $\{t_{i'}\}_{i' = 1}^{m}$ i.e. $0 = t_1 < t_2 < t_3 < \dots < t_{m} = 1$. The initial blending surface can be formed as follows:

and

Let $D = B^{T}B \otimes N^{T}N$ where $\otimes$ is Kronecker product. Then, the remaining surfaces of the sequence $C^{k+1}(t, v)$ for $k \geq 0$, can be compute as follows:

where

We can use the following expression in (3.28) as in the case of curves:

If (3.28) converge, then

the HSS-LSPIA format of surfaces, including two iterations and substituting the sequential iterative process length from each control point in LSPIA, is recognized as (3.28).

3.4.   Convergence analysis of the HSS-LSPIA surface fitting

Lemma 3.7. ^[1] Given an NTP bases $\{B_{i}(t_{i'})\}_{i, i' = 1}^{n_{1}, m_{1}}$ $\in$ $\mathbb{C}^{m_{1} \times n_{1} }$ and $\{N_{j}(v_{j'})\}_{j, j' = 1}^{n_{2}, m_{2}}$ $\in$ $\mathbb{C}^{m_{2} \times n_{2}}$, $B^{T}B$ and $N^{T}N$ are nonsingular collocation matrices, where $m_{1} \geq m_{2} \geq n_{1} \geq n_{2}$. Assuming that $\lambda_{i}(B^{T}B), \ i = 1, 2, \dots, n_{1}$ and $\lambda_{i}(N^{T}N), \ i = 1, 2, \dots, n_{2}$ are their eigenvalues respectively. Then,

$(1) \quad 0 < \lambda_{i}(B^{T}B) \leq 1, 0 < \lambda_{i}(N^{T}N) \leq 1,$

$(2) \quad 0 < \lambda_{i}(D) \leq 1, \quad D = B^{T}B \otimes N^{T}N$ and here $\otimes$ is Kronecker product.

Let $D = \hat{M}_{i} - \hat{N}_{i} \ (i = 1, 2)$ be two splitting of the matrix $D$. The two iterative processes of (3.28) can be written in the matrix form

where $\hat{M}_{1} = (\alpha I + \mu D^{T}D), \ \hat{N}_{1} = \alpha I, \ \hat{M}_{2} = \alpha I, \hat{N}_{2} = (\alpha I - \mu D^{T}D)$, $I$ is identity matrix. Then we can rewrite (3.31) as a unified matrix iteration format as follows

or

where $\hat{M} = \hat{M}_{2}^{-1} \hat{N}_{2} \hat{M}_{1}^{-1} \hat{N}_{1}$ and $\hat{M}_{0} = \hat{M}_{2}^{-1} (I +\hat{N}_{2} \hat{M}_{1}^{-1})$.

Now we need to show that the iterative control point sequence $\{P^{k+1}\}$ leads to the unique solution $P^{*}$, i.e. $\rho (M) < 1$.

Theorem 3.8. Let $D$ be a positive definite matrix, let $\hat{H} = \frac{\mu}{2} \left(D^{T}D + \left(D^{T}D \right) \right) = \mu D^{T}D$ and $\hat{S} = \frac{\mu}{2} \left(D^{T}D - \left(D^{T}D \right) \right) = 0$ be its Hermitian ans skew-Hermitian part, and let $\alpha$ be a positive constant. The HHS-LSPIA iteration's matrix $M(\alpha)$ is then provided by

or

and the spectral radius of $\hat{M}(\alpha)$ is $\rho(\hat{M}(\alpha))$, bounded by

where $\lambda(\hat{H})$ is the spectral set of the matrix $\hat{H}$.

As a consequence, it maintains that

i.e., the HHS-LSPIA iteration converges to the unique solution $P^{*}$ of the linear equation (3.3).

Proof. By the resemblance invariance to the matrix spectrum, we get

We know that $(\alpha I)(\alpha I)^{-1} = I$ which is unitary matrix i.e., $|| (\alpha I)(\alpha I)^{-1} || _{2} = 1$. Then, we get

Since $\lambda_{i} > 0 \ (i = 1, 2, \dots, n)$ and $\alpha$ is positive constant, it is simpler to see that $\rho(\hat{M}(\alpha)) \leq\rho(\alpha) < 1$.

Corollary 3.9. Under the same assumption of Theorem 3.8. Let $\lambda_{max}$ and $\lambda_{min}$ be the maximum and the minimum eigenvalues of $D^{T}D$, respectively. Then, we have the fastest convergence rate of HHS-LSPIA methed for the surface fitting when

where $\mu = \frac{2}{ \lambda_{max} + \lambda_{min} }$ presented in ^[3].

The proof to the Corollary 3.9 is omitted since it is similar to the case of Corollary 3.5.

Remark 3.10. From (3.28) and Corollary 3.9, we will get $\frac{\mu}{\alpha} = \frac{\mu}{ \mu \sqrt{ \lambda_{max}(D^{T}D) \lambda_{min}(D^{T}D)}} = \frac{1}{ \sqrt{ \lambda_{max}(D^{T}D) \lambda_{min}(D^{T}D)}} \\ = \hat{\omega}$. Then, we use the $\hat{\omega} = \frac{1}{ \sqrt{ \lambda_{max}(D^{T}D) \lambda_{min}(D^{T}D)}}$ for compute in (3.28) instead of $\frac{\mu}{\alpha}$.

4.   Implementation and examples

In this section, we demonstrated how the HHS-LSPIA effectiveness for three B-spline curve fitting examples and three examples of B-spline tensor product surface fitting as follows:

● Example 1: A subdivision curve generated using the incenter subdivision technique yielded 100 points (music notation).

● Example 2: $y = 3 cos(\theta) \ (\theta \in [0, 6.28])$ is a sample of 68 points taken evenly from the polar coordinate equation.

● Example 3: G-shape font characteristics measured and smoothed with 315 points.

● Example 4: $z = y^2 cos(x)$ the surface model with $21 \times 21 \ (441)$ points.

● Example 5: $z = y sin(x^2) - x cos(y^2)$ the surface model with $21 \times 21 \ (441)$ points.

● Example 6: $z = x^2 sin(y^2)$ the surface model with $33 \times 33 \ (1089)$ points.

For all the examples, the number of control points is equal to the number of data points and the number of control points is the same in each method. We show the details of the implementation for B-spline curve fitting with the HSS-LSPIA method in the following sections, and we deduce the specifics surface for B-spline tensor product surface fitting in just the same way. Given an ordered point set $\{ Q_j\}_{j = 1}^{m}$, we assign the parameters $\{ t_{j} \}_{j = 1}^{m}$ for $\{ Q_{j} \}_{j = 1}^{m}$ with the normalized accumulated chord parameterization method ^[3,20], that is, $t_1 = 0, \ t_m = 1,$

where $D = \sum_{j = 2}^{m} \|Q_j- Q_{j-1}\|$ is the sum of chord length. In addition, the knots for the cubic B-spline fitting curve $C(t) = \sum_{i = 1}^{n} B_{i, 3}(t)P_{i}$, is defined as $\{ 0, 0, 0, 0, \bar t_4, \bar t_5, \dots, \bar t_n, 1, 1, 1, 1\}$, where

where $d = \frac{m}{n-2}, \; \alpha = id-i, \; j = \lfloor{id} \rfloor$ $i \in \{1, 2, 3 \dots, n-3\}$.

In our implementation, we choose the parameter $\alpha^{*}$ defined by (3.24) for the HSS-LSPIA method. The iteration process is stopped for the comparison of the LSPIA ^[3], ELSPIA ^[6], WHPIA ^[9] and HSS-LSPIA methods if $\| E_{k+1} - E_{k} \| < 10 ^{-7}$, where the fitting error of the $k-$th iteration is computed by

and we compare $E_{k}$ of each iteration to show the efficiency and validity of the HSS-LSPIA.

All the examples are operated on a laptop with a 2.4 GHz Quad-Core Intel Core i5 processor and 8 GB memory via MATLAB R2019a.

In Examples 1–3, the cubic B-spline curves created by the LSPIA ^[3], ELSPIA ^[6], WHPIA ^[9] and HSS-LSPIA methods with the fitting result are shown in Figures 1–3 respectively. In Examples 4–6, the B-spline tensor product surface fitting created by the LSPIA ^[3], ELSPIA ^[6], WHPIA ^[9] and HSS-LSPIA methods with the fitting result are shown in Figures 4–6 respectively. The blue points in each of the three instances represent known points that form a dotted limit curve to be fitted, green points represent control points gained at the appropriate step, and the red line shows the curve of each iteration.

In Figure 7a–7f, we can see that the HSS-LSPIA method has the fitting inaccuracy of the $k-$th iteration less than the ELSPIA method, LSPIA method and the WHPIA method respectively.

In Table 1 columns "IT" (iteration numbers), "CPU time(s)" (ten times the average amount of CPU times(s)) and "$\| E_{k+1} \|$" (fitting error of the $(k+1)$th iteration) for the LSPIA ^[3], ELSPIA ^[6], WHPIA ^[9] and HSS-LSPIA methods are listed. We can observe that the HSS-LSPIA methods requires fewer iteration numbers, uses less CPU time, and has a smallest fitting error at $(k+1)$th iteration. We can now conclude from the experiments that, the HSS-LSPIA methods is faster than the ELSPIA method, LSPIA method and the WHPIA method respectively, in terms of convergence speed. The examples provided in this section are just intended to show how effective the methods propose by us. Because the spectral radius will vary depends on the NTP bases, we cannot guarantee that our methods will converge faster than ELSPIA, LSPIA, and WHPIA in other examples.

5.   Conclusions

In this research, based on the HSS iterative approach for solving system of linear equations of LSPIA, a new progressive and iterative approximation for least square fitting method called HSS-LSPIA is developed to solve the problems of curves and surfaces approximation with NTP bases. The HSS-LSPIA format including two iterations with the iterative distinction vectors in which the two iterations are different from each other. And then, we provided the approximate value of the fastest convergence positive constant. In terms of iterative steps and computing time, numerical results indicate that the HSS-LSPIA method is more effectively efficient than the ELSPIA, LSPIA, and WHPIA methods, respectively. It has demonstrated that the HSS-LSPIA method can be faster than the ELSPIA, LSPIA, and WHPIA methods in terms of convergence speed.

Acknowledgments

The authors acknowledge the financial support provided by the Center of Excellence in Theoretical and Computational Science (TaCS-CoE), KMUTT. The first author was supported by the Petchra Pra Jom Klao Ph.D. Research Scholarship from King Mongkut's University of Technology Thonburi (Grant No. 54/2561). Moreover, this research project is supported by Thailand Science Research and Innovation (TSRI) Basic Research Fund: Fiscal year 2022 under project number FRB650048/0164.

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this article.