The weighted GS-PIA algorithm for cubic B-spline curve interpolations and convergence analysis

THE WEIGHTED GS-PIA ALGORITHM FOR CUBIC B-SPLINE CURVE INTERPOLATIONS AND CONVERGENCE ANALYSIS ZHONGYUN LIU∗, JIAN YANG∗, XIAOFEI XU∗, MENGZHU LIN∗,AND YULIN ZHANG† Abstract. The weighted Gauss-Seidel-progressive iterative approximation (WGS-PIA) algorithm for cubic B-spline curve interpolations is considered in this paper. The convergence of the WGS-PIA algorithm is analyzed, and an upper bound which is strictly smaller than one for the contraction factor of this WGS-PIA algorithm is derived. It is shown that for cubic B-spline curve interpolations, the GS-PIA algorithm converges faster than the Jacobi-PIA (J-PIA) algorithm, and that there always exists a positive weight ωsuch that the WGS-PIA converges faster than GS-PIA. Particularly, we derive a formula for the effective weight ω⋆and the “theoretically optimal” weight ωm, which significantly improves the performance of the WGS-PIA algorithm with minimal additional cost. The numerical experiments are shown that for a given termination tolerance, the number of iterations and the CPU time required by the WGS-PIA algorithm are less than those required by the GS-PIA algorithm. Key words. curve interpolations; WGS-PIA algorithm; convergence; cubic B-spline basis; optimal weight AMS subject classifications. 15A48, 15A51, 65D17, 65F10 1. Introduction. The progressive iterative approximation (PIA) method [10] (also called the geometric iterative method [11]), plays an important role in curve and surface-fitting for a given set of data points. This method is of clear geometric meaning, stable convergence and simple iterative scheme. Above all, the PIA avoids solving a system of linear equations directly, which may cause numerical instability when the number of given data points becomes larger. For those reasons, the PIA method and its variants have intrigued researchers for decades. In particular, the PIA method has been widely applied in CAGD, data fitting and reverse engineering, mesh generation and so on [9, 11]. Given an ordered set of data points {ql}n l=0 in R2or R3, each point qlis assigned a parameter value tl for all l= 0,1,· · · , n, and they follow t0= 0, tl= l−1 X i=0 ∥qi+1 −qi∥2, l = 1,· · · , n. It is clear that t0< t1<· · · < tn. We then define t−3=t−2=t−1=t0, tn=tn+1 =tn+2 =tn+3 to get {tl}n+3 l=−3, generating the cubic B-spline basis {bl(t)}n l=0 by the de Boor-Cox formula [21]. These functions are a set of locally effective piecewise polynomial functions, used to construct complex smooth curves and surfaces [10, 12, 14, 18]. Throughout this paper, all PIA-type algorithms are assumed to be cubic B-spline curve interpolation algorithms if not specifically mentioned. The classical PIA process for cubic B-spline curve interpolations is as follows. The PIA scheme. Setting {p(0) l}n l=0 ={ql}n l=0 (the initial control polygon) and constructing c(0)(t) = Pn l=0 p(0) lbl(t) (the initial interpolating curve), for k= 1,2,· · · ,until {p(k) l}n l=0 converges, ∗School of Mathematics and Statistics, Changsha University of Science and Technology, Changsha 410076, P. R. China (Liu: [email protected]; Yang: [email protected]; Xu: [email protected]; Lin: [email protected]) †Centro de Matem´atica, Universidade do Minho, 4710-057 Braga, Portugal (Yulin Zhang: [email protected]). 1 2Z. Y. Liu, J. Yang, X. F. Xu, M. Z. Lin, and Y. L. Zhang •to compute δ(k) l=ql−c(k−1)(tl),for l= 0,1,· · · , n; •to update p(k) l=p(k−1) l+δ(k) l,for l= 0,1,· · · , n; (1.1) •to construct the k-th interpolating curve c(k)(t) = Pn l=0 p(k) lbl(t). Eq. (1.1) generates the k-th control polygon {p(k) l}n l=0 and satisfies the following relation p(k) l=p(k−1) l+ql− n X j=0 p(k−1) jbj(tl),for l= 0,1,· · · , n. (1.2) Denote P(k)= [p(k) 0,p(k) 1,· · · ,p(k) n]T,Q= [q0,q1,· · · ,qn]T; then, in matrix-matrix form, Eq. (1.2) can be expressed as follows, P(k)=P(k−1) + (Q−BP(k−1)),for k= 1,2,· · · ,(1.3) where B=     b0(t0)b1(t0)· · · bn(t0) b0(t1)b1(t1)· · · b0(t0) . . .. . .. . . b0(tn)b1(tn)· · · bn(tn)      (1.4) is the collocation matrix resulting from the basis functions {bl(t)}n l=0 at tj, j = 0,1,· · · , n. Very clearly, the PIA (1.3) is a special Richardson iteration for matrix equation BP=Q, see, for instance, references [3, 16]. To accelerate the convergence rate of PIA, some variants of PIA were proposed. Those variants can be written as the following form P(k)=P(k−1) +M−1(Q−BP(k−1)).(1.5) From the viewpoint of numerical linear algebra, Eq. (1.5) can be viewed as a classical splitting iteration for matrix equation BP=Q, with the splitting defined as B=M−N. Usually, the matrices M,M−1Band N are called preconditioning, preconditioned and residual matrices, respectively, see, for example, the reference [2]. So we refer to Eq. (1.5) as the preconditioned PIA (PPIA). For example, the PPIA (1.5) becomes •the PIA, if M=I; •the WPIA [15], if M=1 ωI; •the J-PIA [14], if M=D; •the GS-PIA [18], if M=D+L; •the PIA with different weights (DWPIA) [20], if M=diag(ω1,· · · , ωn+1), where B=D+L+U(1.6) with L,D,Udenoting strictly lower triangular part, diagonal part, strictly upper triangular part of the collocation matrix B(1.4), respectively. In this paper, we focus on the acceleration of the GS-PIA algorithm for cubic B-spline curve interpolations [18]. Motivated by the work for the WPIA algorithm in [15], we consider the WGS-PIA algorithm for cubic B-spline curve interpolations. Additionally, we offer a more succinct proof of convergence for both GS-PIA The Weighted GS-PIA algorithm 3 and J-PIA compared to the analyses presented in [14, 18]. Notably, we demonstrate that GS-PIA exhibits faster convergence rates compared to J-PIA. The organization of this paper is as follows: after recalling some basic definitions and known results in next section, we introduce the WGS-PIA algorithm and analyze its convergence in Section 3, a highly efficient weight selection strategy is introduced and some numerical examples are shown in Section 4, and a brief conclusion is followed in last section. Notation. The notations |K|,∥K∥,ρ(K) denote the absolute value, the 2-norm, and the spectral radius of any matrix K, respectively. The absolute value |K|is defined as the matrix obtained by taking the absolute value of all elements of K. For any matrices A= (aij) and B= (bij), A > B (A≥B) means that aij > bij (aij ≥bij) for all indices i, j. 2. Preliminaries. In this section we recall some basic definitions and known results which will be used later on. Definition 2.1. [1, 7] Let Zn×ndenote the set of all real n×nmatrices which have all non-positive off-diagonal entries. A nonsingular matrix A∈Zn×nis called M-matrix if Ais a nonsingular matrix and A−1≥0. Definition 2.2. [1] For an n×nmatrix A= (alj ), we define a new matrix ⟨A⟩= (αlj ), where αlj =(|alj |, j =l −|alj|, j =l, and we called this matrix ⟨A⟩the comparison matrix of A. Definition 2.3. [1] If ⟨A⟩is an M-matrix, then Ais said to be an H-matrix. Lemma 2.4. [1] If Ais a real n×n H-matrix, then |A−1| ≤ ⟨A⟩−1. Definition 2.5. [5, 7, 8] Let A=M−Nbe a splitting. Then this splitting is called: •regular, if M−1≥0and N≥0; •weak regular, if M−1≥0and M−1N≥0; •H-splitting if ⟨M⟩−|N|is an M-matrix; •H-compatible splitting if ⟨A⟩=⟨M⟩−|N|. Lemma 2.6. [4] Let A−1≥0and A=M1−N1=M2−N2 be weak regular splittings of A. In either of the following cases a) N1≤N2; b) M−1 1≥M−1 2,N1≥0; c) M−1 1≥M−1 2,N2≥0; the inequality ρ(M−1 1N1)≤ρ(M−1 2N2) holds. Lemma 2.7. [5] Let Abe an n×n H-matrix. •If the splitting A=M−Nis an H-splitting, then ρ(M−1N)<1. •If the splitting A=M−Nis an H-compatible splitting, then it is an H-splitting. 4Z. Y. Liu, J. Yang, X. F. Xu, M. Z. Lin, and Y. L. Zhang 3. The WGS-PIA algorithm. Let B= (D+L) + Ube a splitting defined as in (1.6). Similar to the WPIA algorithm [15], by taking M=1 ω(D+L) in Eq. (1.5), we then get the following WGS-PIA algorithm P(k)=P(k−1) +ω(D+L)−1(Q−BP(k−1)),(3.1) where the iteration matrix is Gωgs =I−ω(D+L)−1B. Obviously, when ω= 1, the WGS-PIA reduces to the GS-PIA. Therefore, we can expect that the WGS-PIA has a better convergence behavior than the GS-PIA for some ω. 3.1. The convergence theorems for the J-PIA and GS-PIA. It is shown in [12, 14, 18] that collocation matrix Bresulting from cubic B-spline basis functions is a totally nonnegative H-matrix. Thus, we can give more concise proofs for Theorem 1 in [14] and Theorem 2 in [18]. Theorem 3.1. ([14, Theorem 1], [18, Theorem 2]) Let Bin (1.4) be the collocation matrix resulting from a cubic B-spline basis and B=D+L+Ube defined as in (1.6). Then, the iteration matrix of the J-PIA method [14] is GJ=−D−1(L+U)and the iteration matrix of the GS-PIA method [18] is GGS =−(D+L)−1U. Furthermore, ρ(GJ)<1and ρ(GGS )<1. Proof. From the hypothesis, we have that B=D+L+Uis a totally nonnegative H-matrix. Then, from Definition 2.5, we have the splitting B=D+ (L+U) and B= (D+L) + Uare H-compatible splittings, respectively. Hence, by Lemma 2.7, we obtain ρ(GJ)<1 and ρ(GGS)<1. Thus, the proof is complete. In particular, we observe from the numerical tests in [18] that for cubic B-spline curve interpolations, the GS-PIA converges faster than the Jacobi-PIA. Now we give a theoretical proof. Theorem 3.2. Under the hypothesis of Theorem 3.1, ρ(GGS )≤ρ(GJ). Proof. From the hypothesis, we have that B=D+L+Uis a totally nonnegative H-matrix, which implies D≥0, L≥0, U≥0 and D+Lbeing an H-matrix. Note that the iteration matrix of the GS-PIA is GGS =−(D+L)−1Uand the iteration matrix of the Jacobi-PIA is GJ=−D−1(L+U). Then we have ρ(GGS) = ρ((D+L)−1U) ≤ρ[|(D+L)−1||U|] ≤ρ[⟨D+L⟩−1U] (by Lemma 2.4) =ρ[(D−L)−1U]. Since D−1≥0, (D−L)−1≥0, D−1(L+U)≥0 and (D−L)−1U≥0, so the splittings ⟨B⟩=D−(L+U) and ⟨B⟩= (D−L)−Uare both weak regular. Due to ⟨B⟩−1≥0 and U≤L+U, by Lemma 2.6, we then have ρ((D−L)−1U)≤ρ(D−1(L+U)) = ρ(−D−1(L+U)) = ρ(GJ). This means ρ(GGS)≤ρ(GJ). The proof is complete. 3.2. The convergence theorem. Firstly, we review Theorem 3.3 ([19, Theorems 1-2]), the convergence theorem of the extrapolation scheme of the general stationary iteration method, as follows. Lemma 3.3. ([19, Theorems 1-2]) Let S={νj=γj+iηj,j= 1,· · · , n}be the set of all eigenvalues of the iteration matrix Gof the stationary iteration scheme xm+1 =Gxm+c, m = 0,1,· · · .(3.2) Let γM, γmand ηMbe the largest and the smallest real part, and the largest magnitude imaginary part of eigenvalues of G, respectively, and γm≤γj≤γM<1holds. Then The Weighted GS-PIA algorithm 5 (1) The extrapolation scheme of (3.2) xm+1 =Gωxm+ωc, m = 0,1,· · · ,(3.3) where Gω= (1 −ω)I+ωG, converges if and only if 0< ω < ζ with ζ= minj{2(1−γj) (1−γj)2+η2 j}. (2) Let ωm=ω1if ϑ ≤ψ ω∗if ϑ ≥ψ,(3.4) where ω1=1−γM (γM−1)2+γ2 j,ω∗=2 2−(γm+γM),ϑ= (1 −γM)(γM−γm)and ψ= 2η2 M, we have min ωρ(Gω)≤ρ(Gωm) =    ηM ((γm−1)2+η2 M)1/2<1, if ωm=ω1, [(γM−γm)2+4β2 M]1/2 2−γm−γM<1, if ωm=ω∗, (3.5) with equality holding iff (1) γM+iηM∈S, when ϑ≤ψ, or (2) γm+iηM,γM+iηM∈S, when ϑ≥ψ. Now, we can give the convergence theorem of the WGS-PIA algorithm (3.1) for cubic B-spline curve interpolations. Theorem 3.4. Let µj=αj+iβj,j= 1,2,· · · , n, be all eigenvalues of the matrix H= (D+L)−1B. Let αm=min 1≤j≤n{αj},αM=max 1≤j≤n{αj},βM=max 1≤j≤n|βj|,ζ=minj{2αj α2 j+β2 j},ϑ=αm(αM−αm),ψ= 2β2 M, ω1=αm α2 m+α2 M,ω∗=2 αm+αM. Then, under the hypothesis of Theorem 3.1, we have (1) when 0< ω < ζ, the WGS-PIA (3.1) converges. (2) when ωm=ω1if ϑ ≤ψ ω∗if ϑ ≥ψ,(3.6) an upper bound of the spectral radius of the iteration matrix Gωgs of the WGS-PIA (3.1) reaches its minimum, and the following relation min ωρ(Gωgs)≤   βM (α2 m+β2 M)1/2<1, if ωm=ω1, [(αM−αm)2+4β2 M]1/2 αm+αM<1, if ωm=ω∗, holds. Proof. We notice that H= (D+L)−1B=I+ (D+L)−1U=I−GGS. According to Theorem 3.1, we have that ρ(GGS )<1, i.e., |νj|<1, which implies γm< γj< γM<1. Thus, by Lemma 3.3, we have derived this conclusions of (1) and (2). Theorem 3.4 shows that the performance of the WGS-PIA algorithm is seriously dependent on the choice of weights. We remark that the ωmin (3.6) minimizes an upper bound of the spectral radius of the iteration matrix Gωgs of WGS-PIA, and this ωmdoes not minimize ρ(Gωgs) itself. Thus, it is a crucial and challenging task to find an approximation ω⋆of the optimal weight ωopt that minimizes ρ(Gωgs). 4. Numerical experiments. In this section, we first present an inexpensive approach to obtain ω⋆, which is a numerical approximation of ωopt. This approach needs to compute the maximum and minimum eigenvalues of Hin Theorem 3.4, which can be computed using the power and inverse iteration methods [17]. However, the computational complexity of these algorithms are rather large. Therefore, exploiting the 6Z. Y. Liu, J. Yang, X. F. Xu, M. Z. Lin, and Y. L. Zhang structure of collocation matrix, we first give a strategy to obtain the approximation of those eigenvalues in Section 4.1. Next, we establish the relationship between ω⋆and those eigenvalues through a large number of numerical examples in Section 4.2. Finally, we give some numerical examples to show the effectiveness of WGS-PIA in Section 4.3. 4.1. Approximation of the eigenvalues of preconditioned matrix. By Theorem 3.4, we guess that there is a correlation between ω⋆and ωmin (3.6). The latter requires the pre-computation of the maximum and minimum eigenvalues of the preconditioned matrix H. Fortunately, based on the experimental data so far, ω⋆changes only slightly when the eigenvalues of H are slightly perturbed. This enables us to approximately compute the maximum and minimum eigenvalues of H. Since D+Lis an H-matrix and ⟨D+L⟩=⟨D⟩−|L|, according to Lemma 2.7, we can derive that ρ(D−1L) = ρ(˜ L)<1. [16, Theorem 1.11] ensures that (I+˜ L)−1= ∞ X n=0 (−˜ L)n,(4.1) Thus, the following is an approximate form of H: H=(D+L)−1B =I+ (D+L)−1U =I+ (I+˜ L)−1˜ U ≈I+˜ U−˜ L˜ U(4.2) with ˜ U=D−1Uand first two items of (4.1) are retained to yield (4.2). It is clear that only O(n) flops are needed to compute (4.2), which is much less than the O(n2) flops needed to compute H. Moreover, (4.2) is an upper bidiagonal matrix, so its eigenvalues lie on the diagonal. Thus, it is straightforward to obtain approximate maximal and minimal eigenvalues of H. In contrast, it is extremely expensive to compute the eigenvalues of lower Hessenberg matrix Hdirectly. To sum up, the above method can approximately derive eigenvalues of H, which consumes much less CPU time than the original matrix. 4.2. Computation of ω⋆.Usually, ω⋆should be in the interval ωm+1 2, ωm. Through numerous numerical experiments, we observe that ω⋆is linearly related to ωm. Furthermore, by linear regression, we establish the following mathmatical model ω⋆=2ωm+ 1 3.(4.3) The following will show the relationship between the linear regression model and the sample sets. On the Archimedes spiral of Example 4.4, we select interpolation points of different sizes and pseudo-random environments as sample sets. Concretely, we use the MATLAB function rng to initialize the Mersenne Twister generator, selecting random number seeds from 1 to 5. We obtain a total of 100 experimental sample sets {(ω(i) m, ω(i))}i=1,··· ,100 from 20 groups of data points (n= 10,20,...,200), each conforming to a uniform distribution over [0,1]. Each ω(i) mis computed by combining (4.2) and (3.6), and the ω(i)for each sample set are determined by traversing ω(i)= 0.1:0.01 : 2, where a MATLAB colon is used to denote that ωare taken every 0.01 between the interval [0.1,2]. In Fig. 4.1, Sample points and Prediction curve denote all 100 experimental sample sets {(ω(i) m, ω(i))}i=1,··· ,100 and model curve (4.3), respectively. Notice that all of these sample points are distributed around the curve. The Weighted GS-PIA algorithm 7 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Sample points Prediction curve Fig. 4.1 Comparison of sample sets {(ω(i) m, ω(i))}i=1,··· ,100 on the Archimedes spiral of Example 4.4 To demonstrate the effectiveness of the ω⋆in (4.3), we test ω= 0.1 : 0.01 : 2 and ω=ω⋆in Examples 4.2 and 4.5 which are given in next subsection. The results are displayed in Fig.4.2, where ϵ(k)= max 0≤l≤n∥ck(tl)−ql∥denotes the interpolation errors at the k-th step of WGS-PIA Algorithm (3.1), and IT denotes the corresponding number of iterations. We observe that all gray lines are clustered to the right of the red line, which implys that ω⋆is a good weight for the WGS-PIA algorithm. 0 200 400 600 800 1000 1200 IT 10-10 10-8 10-6 10-4 10-2 100 (k) Fig. 4.2 ϵ(k)for WGSPIA with different weights in Example 4.2 and 4.5 In the next subsection, we use the WGS-PIA method with ω⋆in (4.3) for computing Examples 4.1-4.5. All numerical experiments show that our approach is effective. For more details, see subsection 4.3. 4.3. Numerical examples. We use the non-uniform cubic B-spine basis to verify the effectiveness of the WGS-PIA algorithm. For comparison, we also test the GS-PIA algorithm. The performance of the algorithms is evaluated by the number of iterations (IT) and the elapsed CPU time/s (CPU/s), and CPU time is calculated by averaging 100 repeated experiments. All numerical experiments were performed on computer DESKTOP-9TC1VPG with 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz by MATLAB(R2022a). We give five numerical examples as follows. These numerical examples are taken from [18] and 8Z. Y. Liu, J. Yang, X. F. Xu, M. Z. Lin, and Y. L. Zhang http://paulbourke.net/geometry/. Example 4.1. Consider a parametric function given by          x= cos u2−cos 2u 2k+1  y= sin u2−cos 2u 2k+1  z=−sin 2u 2k+1  where 0 ≤u≤(4k+ 2)π. It is called Cinquefoil Knot if k = 2. The 500 interpolation points on it were selected non-uniformly. Example 4.2. Consider the helix given by      x=rcos(τ(θ2−θ1) + θ1) y=rsin(τ(θ2−θ1) + θ1) z=hτ where r= 30, h = 50, θ1=π 6, θ2= 40π. The value of τat 2000 in [0, 1] is uniformly selected to obtain the interpolation points. Example 4.3. Consider 400 interpolation points on the Freeths Nephroid ρ= 1 + 2sin(θ/2)(0 ≤θ≤ 6π). Example 4.4. Consider the data interpolation of 3000 points sampled from Archimedes spiral ρ=θ(0 ≤ θ≤30π). Example 4.5. Given 987 scattered interpolation points. We select Examples 4.1-4.2 for the 3D curve fits, Examples 4.3-4.5 for the 2D curve fits. Among them, interpolation points of the Examples 4.1 and 4.3 are selected at a smaller scale than Example 4.2 and Example 4.4. Example 4.5 was given as an example of a set of scatter points. In Table 4.1-4.5, we record the number of iterations (the elapsed CPU time) required by GS-PIA algorithm and WGS-PIA algorithm when ϵ(k)<10−6,10−7,· · · ,10−11. The weight ω⋆are computed via (4.3), and its values(CPU time) are listed independently in the second column of the tables. Table 4.1 The number of iteration steps (CPU time) for some given tolerances in Example 4.1 Algorithm ω⋆(CPU/s) IT(CPU/s) ϵ(k)<10−6ϵ(k)<10−7ϵ(k)<10−8ϵ(k)<10−9ϵ(k)<10−10 ϵ(k)<10−11 GS-PIA - 322(0.0296) 432(0.0387) 542(0.0492) 652(0.0594) 762(0.0679) 872(0.0796) WGS-PIA 1.6396(0.0015) 195(0.0212) 262(0.0293) 328(0.0340) 395(0.0400) 462(0.0476) 528(0.0541) Table 4.2 The number of iteration steps (CPU time) for some given tolerances in Example 4.2 Algorithm ω⋆(CPU/s) IT(CPU/s) ϵ(k)<10−6ϵ(k)<10−7ϵ(k)<10−8ϵ(k)<10−9ϵ(k)<10−10 ϵ(k)<10−11 GS-PIA - 774(1.6800) 980(2.0701) 1186(2.7445) 1391(3.3275) 1600(3.5826) 1825(3.6609) WGS-PIA 1.4380(0.0345) 467(1.0994) 592(1.4877) 715(1.7682) 839(1.9915) 965(2.2573) 1101(2.6753) The Weighted GS-PIA algorithm 9 Table 4.3 The number of iteration steps (CPU time) for some given tolerances in Example 4.3 Algorithm ω⋆(CPU/s) IT(CPU/s) ϵ(k)<10−6ϵ(k)<10−7ϵ(k)<10−8ϵ(k)<10−9ϵ(k)<10−10 ϵ(k)<10−11 GS-PIA - 1113(0.0870) 1575(0.1237) 2038(0.1621) 2500(0.2039) 2962(0.2390) 3424(0.2750) WGS-PIA 1.6601(0.0012) 670(0.0469) 948(0.0738) 1226(0.1137) 1504(0.1353) 1782(0.1589) 2060(0.1858) Table 4.4 The number of iteration steps (CPU time) for some given tolerances in Example 4.4 Algorithm ω⋆(CPU/s) IT(CPU/s) ϵ(k)<10−6ϵ(k)<10−7ϵ(k)<10−8ϵ(k)<10−9ϵ(k)<10−10 ϵ(k)<10−11 GS-PIA - 1903(5.9570) 2540(7.9414) 3176(9.8952) 3812(11.8881) 4449(13.8767) 5085(15.8941) WGS-PIA 1.6740(0.0648) 1144(4.2880) 1527(5.7635) 1909(7.0314) 2292(7.9557) 2674(9.2660) 3056(10.5765) Table 4.5 The number of iteration steps (CPU time) for some given tolerances in Example 4.5 Algorithm ω⋆(CPU/s) IT(CPU/s) ϵ(k)<10−6ϵ(k)<10−7ϵ(k)<10−8ϵ(k)<10−9ϵ(k)<10−10 ϵ(k)<10−11 GS-PIA - 50(0.0282) 65(0.0340) 80(0.0402) 96(0.0474) 111(0.0549) 126(0.0647) WGS-PIA 1.5094(0.0076) 37(0.0208) 45(0.0269) 56(0.0275) 64(0.0354) 71(0.0389) 81(0.0443) In Figs. 4.3-4.7, we display the initial control node and their WGS-PIA interpolation figures for the 1-st iteraion and the 10-th iteraions, as well as the relation diagrams of iteration steps and interpolation errors of the GS-PIA and WGS-PIA. It can be seen with Tables 4.1-4.5 that, at least from the interpolation error less than 10−6, the number of iterations of WGS-PIA is only about 60% of that of GS-PIA. This conclusion can be intuitively perceived by the graphs in the lower right corner of Figs. 4.3-4.7. Initially, the interpolation errors of WGS-PIA are reduced more slowly compared to GS-PIA. When ϵis roughly less than 10−4, the errors of WGS-PIA decay faster than that of GS-PIA. On the other hand, WGS-PIA takes less CPU time compared to GS-PIA with the same interpolation error. And in the WGS-PIA algorithm, the computation time of ω⋆is almost negligible compared to the iteration time. Since the weight ω⋆only needs to be solved once in each experiment, as the required convergence accuracy increases, the percentage of computation time for this part becomes increasingly smaller. 5. Conclusions. In this paper, we have developed the WGS-PIA algorithm by combining the GS-PIA [18] with the WPIA [15]. In Theorem 3.4, we have shown that there always exists a positive weight ωsuch that WGS-PIA algorithm converges for cubic B-spline curve interpolations. Moreover, we have given an ωm that minimizes an upper bound of the spectral radius of the iteration matrix Gωgs of the WGS-PIA (3.1). Meanwhile, we have given a more concisely proof in Theorem 3.1, which shows the convergence of the GS-PIA [18] and the J-PIA [14] by splitting theories [5]. In the numerical experiments of [18], we observed GS-PIA has a better performance than J-PIA for cubic B-spline curve interpolations. We then proved it theoretically, and this conclusion has been shown in Theorem 3.2. Theoretically, we can chose ωto be any constant in the interval (0, ζ) with ζdefined in Theorem 3.4.