Progressive Iterative Approximation with Preconditioners

Abstract

The progressive iterative approximation (PIA) plays an important role in curve and surface fitting. By using the diagonally compensated reduction of the collocation matrix, we propose the preconditioned progressive iterative approximation (PPIA) to improve the convergence rate of PIA. For most of the normalized totally positive bases, we show that the presented PPIA can accelerate the convergence rate significantly in comparison with the weighted progressive iteration approximation (WPIA) and the progressive iterative approximation with different weights (DWPIA). Furthermore, we propose an inexact variant of the PPIA (IPPIA) to reduce the computational complexity of the PPIA. We introduce the inexact solver of the preconditioning system by employing some state-of-the-art iterative methods. Numerical results show that both the PPIA and the IPPIA converge faster than the WPIA and DWPIA, while the elapsed CPU times of the PPIA and IPPIA are less than those of the WPIA and DWPIA.

1. Introduction

The progressive iterative approximation (PIA) is an important iterative method for fitting a given set of data points. Due to its clear geometric meaning, stable convergence and simple iterative scheme, the PIA has intrigued researchers for decades. In particular, said methodology has achieved great success in geometric design [1,2], data fitting [3,4], reverse engineering [5] and NURBS solid generation [6,7,8]. For more detailed research on this topic, we refer the reader to read a recent survey [9]. However, the PIA converges very slowly, especially when the number of interpolating points increases to some extent. This is because the collocation matrices resulting from some normalized totally positive (NTP) bases are usually ill-conditioned; see [10,11]. Therefore, much more attention was paid to accelerating the convergence rate of PIA; see the weighted progressive iteration approximation (WPIA) (also named modified Richardson method) in [12,13]; the QR-PIA in [11]; the Jacobi-PIA in [14]; the Chebyshev semi-iterative method in [15]; the progressive iterative approximation with different weights (DWPIA) in [16]; the GS-PIA in [17]; the PIA with mutually different weights in [18]; the Schulz iterative method in [19]; and a lot literature therein.

There are other remedy strategies, such as locally interpolating, i.e, interpolating partial points selected from the total given points, which results in the local progressive-iterative approximation format [20], the extended progressive-iterative approximation format [21] and the progressive iterative approximation format for least square fitting [22].

Let ${\left\{b_i\left(t\right)\right\}}_{i=0}^n$ be an NTP basis, such as the Bernstein basis [23] or the Said–Ball basis [24]. Given a set of data points ${\left\{{\mathit{v}}_i\right\}}_{i=0}^n$ in ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$, we assign a parameter value $t_i$ to the i-th point ${\mathit{v}}_i$. For PIA, the initial interpolating curve is given by

$${\gamma }^{\left(0\right)}\left(t\right)=\sum _{i=0}^n{\mathit{u}}_i^{\left(0\right)}{b}_i\left(t\right),{\mathit{u}}_i^{\left(0\right)}={\mathit{v}}_i,i=0,1,\cdots ,n,$$

and the $(k+1)$-th interpolating curve can be generated by

$${\gamma }^{(k+1)}\left(t\right)=\sum _{i=0}^n{\mathit{u}}_i^{(k+1)}{b}_i\left(t\right),k=0,1,2,\cdots ,$$

where

$${\mathit{u}}_i^{(k+1)}={\mathit{u}}_i^{\left(k\right)}+{\delta }_i^{(k+1)},{\delta }_i^{(k+1)}={\mathit{v}}_i-{\gamma }^{\left(k\right)}\left(t_i\right).$$

(1)

Equation (1) can be rewritten as the following form

$${\mathit{u}}_i^{(k+1)}={\mathit{u}}_i^{\left(k\right)}+{\mathit{v}}_i-\sum _{j=0}^n{\mathit{u}}_j^{\left(k\right)}{b}_j\left(t_i\right),i=0,1,\cdots ,n,$$

(2)

which generates a sequence of control polygons ${\mathit{u}}_i^{(k+1)}$, for $k=0,1,2,\cdots$.

Let ${\Delta }^{\left(k\right)}={\left[{\delta }_0^{\left(k\right)},{\delta }_1^{\left(k\right)},\cdots ,{\delta }_n^{\left(k\right)}\right]}^T$, ${\mathit{U}}^{\left(k\right)}={\left[{\mathit{u}}_0^{\left(k\right)},{\mathit{u}}_1^{\left(k\right)},\cdots ,{\mathit{u}}_n^{\left(k\right)}\right]}^T$, $\mathit{V}={[{\mathit{v}}_0,{\mathit{v}}_1,\cdots ,{\mathit{v}}_n]}^T$. Then the matrix form of Equations (1) and (2) is

$${\mathit{U}}^{(k+1)}={\mathit{U}}^{\left(k\right)}+{\Delta }^{\left(k\right)}=(I-B){\mathit{U}}^{\left(k\right)}+\mathit{V},$$

(3)

where ${\mathit{U}}^{\left(0\right)}=\mathit{V}$, ${\Delta }^{\left(k\right)}=\mathit{V}-B{\mathit{U}}^{\left(k\right)}$, I is the identity matrix and B is the collocation matrix resulting from the basis $\{b_0\left(t\right),b_1\left(t\right),\cdots ,b_n\left(t\right)\}$ at $t_i$, $i=0,1,\cdots ,n$. We note in [23] that the collocation matrix B resulting from any NTP basis is a stochastic totally positive matrix.

If $\underset{k\to \infty }{lim}{\gamma }^{\left(k\right)}\left(t_i\right)={\mathit{v}}_i$, then the sequence of curves ${\gamma }^{\left(k\right)}\left(t\right)$ is said to have the PIA property and Equation (2) is referred to as the PIA format. It is shown in [25] that the PIA property holds for any NTP basis.

In fact, the PIA is mathematically equivalent to the Richardson iteration for solving the collocation matrix equation

$$B\mathit{U}=\mathit{V};$$

(4)

see [15,26]. Therefore, the state-of-the-art iterative method, the generalized minimal residual (GMRES) method for solving general linear systems of equations, was suggested in [26] for solving (4). It is known that the GMRES method works efficiently when the coefficient matrix B is positive definite. It is true that the GMRES has better convergence behavior than PIA in the same required precision when fitting 19 points on the helix. However, we found that the collocation matrix B resulting from some NTP bases is not always positive definite, even if B is a totally positive matrix. This led us to further study the PIA. Obviously, the convergence of PIA depends on the spectral radius of the iteration matrix $I-B$. The smaller the spectral radius is, the faster the PIA converges. This motivated us to use the preconditioning techniques to reduce the spectral radius of the iteration matrix $I-B$. The PIA consists of the interpolatory PIA and the approximating PIA [9]. In the first case, the number of the control points is the same as that of the given data points. In the second case, the number of the selected control points is less than that of the given data points. In this paper, we focus on preconditioning the interpolatory PIA. The proposed preconditioning technique can be extended straightforwardly to the approximating PIA.

The paper is organized as follows. In Section 2, we first exploit the diagonally compensated reduction of the collocation matrix B to construct a class of preconditioners for B; then we establish the PPIA (preconditioned PIA) format and analyze its convergence. In order to improve the computational efficiency of the PPIA, we further propose an IPPIA (inexact PPIA) format and analyze its convergence in Section 3. Finally, some numerical examples are given to illustrate the effectiveness of our proposed methods in Section 4.

2. The Preconditioned Progressive Iterative Approximation (PPIA)

As mentioned earlier, the PIA does always converge although slowly. In this section, we therefore consider the preconditioning technique to accelerate the convergence rate of PIA. Preconditioning means replacing the initial system (4) with the system

$$P^{-1}B\mathit{U}=P^{-1}\mathit{V},$$

(5)

where P is an approximation of B with the following properties: $P^{-1}B$ is well conditioned or has few extreme eigenvalues, and $P\mathit{z}=\mathit{d}$ is easy to solve or $P^{-1}$ is easy to obtain.

A careful and problem-dependent choice of P can often make the condition number of $P^{-1}B$ much smaller than that of B and thus accelerate convergence dramatically.

2.1. Preconditioners

In this subsection, we focus on preconditioner P for the collocation matrix B resulting from an NTP basis. Very often, the preconditioner P should be chosen to preserve certain structures of the original matrix B. There are many applications that generate structured matrices, and by exploiting the structure, one may be able to design faster and/or more accurate algorithms; see for example [27,28,29]; furthermore, structure may also help in producing solutions which have more precise physical meanings. Structure comes in many forms, including Hamiltonian, Toeplitz, Vandermonde matrices and so on. Note that B is a stochastic matrix, so we hope that the preconditioner P is also a stochastic matrix. Moreover, in order to get $P^{-1}$ easily, we restrict P to be a banded matrix with small bandwidth. Based on those two requirements, we construct P as follows.

We first split B into $B=B_q+R$, where

$$R=\left(\begin{array}{cccccc}0& \cdots & 0& b_{q+1}\left(t_0\right)& \cdots & b_n\left(t_0\right)\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \ddots & \ddots & \ddots & \ddots & b_n\left(t_{q+1}\right)\\ b_0\left(t_{q+1}\right)& \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ b_0\left(t_n\right)& \cdots & b_{q+1}\left(t_n\right)& 0& \cdots & 0\end{array}\right)\ge 0,$$

(6)

and $B_q=B-R$ is a banded matrix with bandwidth $2q+1(0\le q\le n)$. Then we define a diagonal matrix

$$D=diag(d_0,d_1,\cdots ,d_n)$$

(7)

such that $D\mathit{e}=R\mathit{e}$ with $\mathit{e}={(1,1,\cdots ,1)}^T$. Finally, let

$$P=B_q+D,$$

(8)

be the diagonally compensated reduction matrix of B. The matrix R is called the reduced matrix and D is called the compensation matrix for R; see for instance [30,31]. It is easy to check that P is a stochastic banded matrix and its construction time is small enough and can be neglected, because it only involves a matrix-vector multiplication. If P in (8) is invertible, it can serve as a preconditioner for the system (4). Unfortunately, we cannot give a theoretical proof about the nonsingularity of P, but numerous numerical experiments demonstrate P is invertible for different order n with various q. Therefore, we have

$$P^{-1}B\mathit{U}=P^{-1}\mathit{V}.$$

(9)

Accordingly, the PIA for (9) becomes the following format

$${\mathit{U}}^{(k+1)}=(I-P^{-1}B){\mathit{U}}^{\left(k\right)}+P^{-1}\mathit{V},k=0,1,\cdots ,$$

where ${\mathit{U}}^{\left(0\right)}=\mathit{V}$, or equivalently

$${\mathit{U}}^{(k+1)}={\mathit{U}}^{\left(k\right)}+P^{-1}{\Delta }^{\left(k\right)},$$

(10)

where ${\Delta }^{\left(k\right)}=\mathit{V}-B{\mathit{U}}^{\left(k\right)}$.

For convenience, we refer to (10) as the PPIA format.

2.2. Convergence Analysis of PPIA

Before analyzing the convergence, we first investigate some properties of the collocation matrix B resulting from the Bernstein basis.

Proposition 1.

For the Bernstein basis, if $t_i=\frac{i}{n}$, $i=0,1,\cdots ,n$, then the entries of the collocation matrix B are

$$b_j\left(t_i\right)=\left(\genfrac{}{}{0pt}{}{n}{j}\right){(1-t_i)}^{n-j}{t}_i^j,i=0,1,\cdots ,n;j=0,1,\cdots ,n,$$

and possess the following properties:

(a)

The diagonal entry $b_j\left(t_j\right)$ is the maximum of all entries which locate at the j-th row and j-th column; i.e., $b_j\left(t_j\right)=\underset{0\le i\le n}{max}{b}_j\left(t_i\right)$ and $b_j\left(t_j\right)=\underset{0\le i\le n}{max}{b}_i\left(t_j\right)$.

(b)

$b_j\left(t_j\right)>b_j\left(t_{j-1}\right)>\cdots >b_j\left(t_0\right)\ge 0$, and $b_j\left(t_j\right)>b_j\left(t_{j+1}\right)>\cdots >b_j\left(t_n\right)\ge 0$.

(c)

$b_j\left(t_j\right)>b_{j-1}\left(t_j\right)>\cdots >b_0\left(t_j\right)\ge 0$, and $b_j\left(t_j\right)>b_{j+1}\left(t_j\right)>\cdots >b_n\left(t_j\right)\ge 0$.

Proof. 

It was proven in [23] that the Bernstein basis has the unimodality property in the interval $[0,1]$; i.e., $b_j\left(t\right)$ has a unique extremum at $t=j/n$ in the interval $[0,1]$. Moreover, for any fixed $t^{*}\in [0,1]$, there exists a corresponding index j such that

$$b_0\left(t^{*}\right)\le \cdots \le b_{j-1}\left(t^{*}\right)\le b_j\left(t^{*}\right)\mathrm{and}{b}_j\left(t^{*}\right)\ge b_{j+1}\left(t^{*}\right)\ge \cdots \ge b_n\left(t^{*}\right).$$

In other words, $b_j\left(t\right)$ reaches its maximum at the value $t^{*}$ which is the closest to $j/n$; thus Property (a) holds. Properties (b) and (c) can be proved in a similar way, by using the unimodality property of the Bernstein basis. □

On one hand, according to the definition of D in (7), we have

$${\parallel D\parallel }_{\infty }={\parallel R\parallel }_{\infty }.$$

(11)

From (8) and (11), we further have

$${\parallel P-B\parallel }_{\infty }={\parallel D-R\parallel }_{\infty }\le 2{\parallel R\parallel }_{\infty }.$$

(12)

On the other hand, from (6), we have

$$\begin{array}{cc}\hfill {\parallel R\parallel }_{\infty }& =\underset{0\le j\le n}{max}\sum _{|i-j|>q}{b}_i\left(t_j\right)\hfill \\ & =\underset{0\le j\le n}{max}\left\{b_0\left(t_j\right)+\cdots +b_{j-q-1}\left(t_j\right)+b_{j+q+1}\left(t_j\right)+\cdots +b_n\left(t_j\right)\right\},\hfill \end{array}$$

where the term $b_k\left(t_j\right)=0$, if $|j-k|\le q$. By Proposition 1, we have

$$b_{j+q+1}\left(t_j\right)>\cdots >b_n\left(t_j\right)=0,b_{j-q-1}\left(t_j\right)>\cdots >b_0\left(t_j\right)=0.$$

For each fixed j, it is easy to verify that both $b_{j+q+1}\left(t_j\right)$ and $b_{j-q-1}\left(t_j\right)$ are monotonically decreasing and finally decay to zero, as q increases. Hence, any of sum of each row’s entries outside of the q-th diagonal decays to 0. This means that ${\parallel R\parallel }_{\infty }$ approaches zero if $q\to n$. Combined with (12), we conclude that the matrix P is a good approximation of B. Thus, the spectrum of $P^{-1}B$ approaches 1 for large q.

Remark 1.

We remark here that the spectral radius $\rho (I-P^{-1}B)$ may decrease as q increases; i.e., the bigger the q is, the faster the PPIA format converges. However, according to (10), one needs to compute $P^{-1}{\Delta }^{\left(k\right)}$ or equivalently solve an additional matrix equation $P{\mathit{X}}^{\left(k\right)}={\Delta }^{\left(k\right)}$ at the $k(k=0,1,\cdots )$-th iteration. This will result in a large computational complexity, especially for large q. Therefore, we need to seek a tradeoff between the convergence rate and computational complexity. Experimentally, we found that $q\approx \frac{n}{2}$ is a suitable choice.

Remark 2.

We emphasize that most of the NTP bases have the unimodality property, so Proposition 1 also holds for other NTP bases with unimodality. This means that our PPIA will work well for such a class of bases, such as the Said–Ball basis. This was experimentally verified in our numerical tests.

3. The Inexact Preconditioned Progressive Iterative Approximation (IPPIA)

As mentioned in Remark 1, we need to compute $P^{-1}{\Delta }^{\left(k\right)}$ or equivalently solve $P{\mathit{X}}^{\left(k\right)}={\Delta }^{\left(k\right)}$ at the k-th iteration. This is very costly and impractical in actual implementations for large n. To reduce the computational complexity, we propose an IPPIA format in which we inexactly solve $P{\mathit{X}}^{\left(k\right)}={\Delta }^{\left(k\right)}$ by employing a state-of-the-art iterative method, that is, the conjugate gradient (CG) method for solving the related system of normal equations.

3.1. The Inexact Solver

Denote by ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$ the iteration sequence generated from IPPIA format defined below. Since the sequence ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$ is an approximation of ${\left\{{\mathit{U}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$, ${\Delta }^{\left(k\right)}$ in $P{\mathit{X}}^{\left(k\right)}={\Delta }^{\left(k\right)}$ should be replaced by ${\tilde{\Delta }}^{\left(k\right)}=\mathit{V}-B{\tilde{\mathit{U}}}^{\left(k\right)}$.

Now, let ${\mathit{X}}^{(l,k)}$ be the l-th approximation of the solution of

$$P{\mathit{X}}^{\left(k\right)}={\tilde{\Delta }}^{\left(k\right)},$$

(13)

which can be obtained by employing the CG method to solve the system of normal equations

$$P^{T}P{\mathit{X}}^{\left(k\right)}=P^T{\tilde{\Delta }}^{\left(k\right)}.$$

(14)

The CG method used here can be thought of as an inner iterative format and terminates if

$$\parallel P^T{\tilde{\Delta }}^{\left(k\right)}-P^{T}P{\mathit{X}}^{(l,k)}{\parallel }_F/\parallel P^T{\tilde{\Delta }}^{\left(k\right)}{\parallel }_F<{\tau }_k,$$

(15)

where ${\tau }_k$ is the prescribed stopping tolerance and $\parallel ·\parallel$ is the Frobenius norm.

Once the approximation ${\mathit{X}}^{(l,k)}$ to (13) satisfies (15), we can form ${\tilde{\mathit{U}}}^{(k+1)}$ by the following outer iteration

$${\tilde{\mathit{U}}}^{(k+1)}={\tilde{\mathit{U}}}^{\left(k\right)}+{\mathit{X}}^{(l,k)}.$$

(16)

We refer to (16) as the IPPIA format.

Very often, to ensure the computational efficiency, the iteration (16) can also be terminated if the stopping criterion,

$${\epsilon }^{(k+1)}\le \theta {\epsilon }^{\left(k\right)},\mathrm{and}\theta \in [0,1],$$

(17)

is satisfied, where ${\epsilon }^{\left(k\right)}$ defined by

$${\epsilon }^{\left(k\right)}=\underset{0\le i\le n}{max}∥{\mathit{v}}_i-{\gamma }^{\left(k\right)}\left(t_i\right)∥$$

(18)

is the interpolation error of the k-th interpolating curve ${\gamma }^{\left(k\right)}\left(t\right)$.

Since the term $P^{-1}{\Delta }^{\left(k\right)}$ in (10) is replaced by an inexact solution ${\mathit{X}}^{(l,k)}$, the iterative sequence ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$ does not necessarily converge to the exact solution ${\mathit{U}}^{*}$. Clearly, the selection of ${\tau }_k$ will affect the convergence of iterative sequence ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$. If ${\tau }_k$ equals to zero, the solution obtained by the CG method is the exact one for the system (14). In this case, the IPPIA format is reduced into the PPIA format. Therefore, one can expect that the iterative sequence has a good convergence rate if ${\tau }_k$ is chosen to be small enough. The selection of ${\tau }_k$ is discussed in the following subsection.

3.2. Convergence Analysis of IPPIA

Theorem 1.

Let P be the preconditioner defined as in (8) for B and ${\mathit{U}}^{*}$ be the exact solution to (4). Suppose that the CG method starts from an initial guess ${\mathit{X}}^{(0,k)}$, terminates if (15) is satisfied and produces an inner iterative sequence $\left\{{\mathit{X}}^{(l,k)}\right\}$ for each outer iterative index k, where ${\mathit{X}}^{(l,k)}$ is an approximation of $P^{-1}{\Delta }^{\left(k\right)}$. Then, for the iterative sequence ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$ generated by IPPIA, we have

$${∥{\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}∥}_F\le {\pi }_k{∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F,$$

where ${\pi }_k=\parallel I-P^{-1}{B\parallel }_F+{\tau }_k\parallel P^{-1}{\parallel }_F\parallel P^{-T}{\parallel }_F{\parallel B\parallel }_F$. In particular, if

$${\pi }_{max}=\parallel I-P^{-1}{B\parallel }_F+{\tau }_{max}\parallel P^{-1}{\parallel }_F\parallel P^{-T}{\parallel }_F{\parallel B\parallel }_F<1,$$

(19)

then the iterative sequence ${\left\{{\tilde{\mathit{U}}}^{\left(k\right)}\right\}}_{k=0}^{\infty }$ converges to ${\mathit{U}}^{*}$, where ${\tau }_{max}=\underset{k}{max}\left\{{\tau }_k\right\}$.

Proof. 

Note that ${\Delta }^{\left(k\right)}=\mathit{V}-B{\mathit{U}}^{\left(k\right)}$, ${\tilde{\Delta }}^{\left(k\right)}=\mathit{V}-B{\tilde{\mathit{U}}}^{\left(k\right)}$, from (10), we have

$$\begin{array}{cc}\hfill {\mathit{U}}^{(k+1)}& ={\mathit{U}}^{\left(k\right)}+P^{-1}\left(\mathit{V}-B{\mathit{U}}^{\left(k\right)}\right)\hfill \\ \hfill & ={\mathit{U}}^{\left(k\right)}+P^{-1}\left(\mathit{V}-B{\tilde{\mathit{U}}}^{\left(k\right)}+B{\tilde{\mathit{U}}}^{\left(k\right)}-B{\mathit{U}}^{\left(k\right)}\right)\hfill \\ \hfill & ={\mathit{U}}^{\left(k\right)}+P^{-1}{\tilde{\Delta }}^{\left(k\right)}+P^{-1}B\left({\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{\left(k\right)}\right).\hfill \end{array}$$

(20)

Let ${\mathit{U}}^{*}$ be the exact solution of (4). Then

$${\mathit{U}}^{(k+1)}-{\mathit{U}}^{*}=\left(I-P^{-1}B\right)\left({\mathit{U}}^{\left(k\right)}-{\mathit{U}}^{*}\right).$$

(21)

According to (16), (20) and (21), we have

$$\begin{array}{cc}\hfill P^{T}P\left({\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}\right)=& P^{T}P\left({\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{(k+1)}+{\mathit{U}}^{(k+1)}-{\mathit{U}}^{*}\right)\hfill \\ \hfill =& P^{T}P\left({\tilde{\mathit{U}}}^{\left(k\right)}+{\mathit{X}}^{(l,k)}\right)-P^{T}P\left[{\mathit{U}}^{\left(k\right)}+P^{-1}{\tilde{\Delta }}^{\left(k\right)}+P^{-1}B\left({\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{\left(k\right)}\right)\right]\hfill \\ \hfill & +P^{T}P\left(I-P^{-1}B\right)\left({\mathit{U}}^{\left(k\right)}-{\mathit{U}}^{*}\right)\hfill \\ \hfill =& P^T(P-B)\left({\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}\right)+P^{T}P{\mathit{X}}^{(l,k)}-P^T{\tilde{\Delta }}^{\left(k\right)}.\hfill \end{array}$$

(22)

Due to the nonsingularity of P, (22) can be transformed into

$${\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}=\left(I-P^{-1}B\right)\left({\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}\right)+P^{-1}{P}^{-T}\left(P^{T}P{\mathit{X}}^{(l,k)}-P^T{\tilde{\Delta }}^{\left(k\right)}\right).$$

(23)

Since $B{\mathit{U}}^{*}=\mathit{V}$, from the hypothesis (15), we have

$${∥P^T{\tilde{\Delta }}^{\left(k\right)}-P^{T}P{\mathit{X}}^{(l,k)}∥}_F<{\tau }_k{∥P^T{\tilde{\Delta }}^{\left(k\right)}∥}_F\le {\tau }_k{∥P^T∥}_F{∥B∥}_F{∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F.$$

(24)

Combining (23) with (24), we have

$$\begin{array}{cc}\hfill {∥{\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}∥}_F& \le {∥\left(I-P^{-1}B\right)\left({\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}\right)∥}_F+{∥P^{-1}∥}_F{∥P^{-T}∥}_F{∥P^{T}P{\mathit{X}}^{(l,k)}-P^T{\tilde{\Delta }}^{\left(k\right)}∥}_F\hfill \\ \hfill & \le \left({∥I-P^{-1}B∥}_F+{\tau }_k{∥P^{-1}∥}_F{∥P^{-T}∥}_F{∥B∥}_F\right){∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F.\hfill \end{array}$$

(25)

Since ${∥I-P^{-1}B∥}_F<1$, we have from (25) that ${\pi }_{max}<1$ for small enough ${\tau }_{max}$. The proof is thus complete. □

Theorem 1 tells us that the error of the IPPIA iteration consists of two parts: ${∥I-P^{-1}B∥}_F$ and ${\tau }_k{∥P^{-1}∥}_F{∥P^{-T}∥}_F{∥B∥}_F$. The first part is the error upper bound for PPIA iteration. The second part is the error upper bound resulting from the CG method, depending on the choice of the stopping tolerance ${\tau }_k$. This ${\tau }_k$ is hard to be optimally determined between the convergence rate and computational costs in practice.

The following theorem presents one possible way of choosing the tolerance ${\tau }_k$.

Theorem 2.

Let the assumptions in Theorem 1 be satisfied. Suppose that $\left\{{\zeta }_k\right\}$ is a nondecreasing and nonnegative sequence satisfying ${\zeta }^{\left(k\right)}\ge 1$, and $\underset{k\to \infty }{lim}sup{\zeta }^{\left(k\right)}=+\infty$, and that χ is a real constant in the interval $(0,1)$ satisfying

$${\tau }_k\le c{\chi }^{{\zeta }^{\left(k\right)}},k=0,1,\cdots ,$$

(26)

where c is a nonnegative constant. Then we have

$$\underset{k\to \infty }{lim}sup\frac{{∥{\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}∥}_F}{{∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F}={∥I-P^{-1}B∥}_F,$$

i.e., the convergence rate of IPPIA is asymptotically the same as that of PPIA.

Proof. 

From (25) and (26), we have

$$\begin{array}{cc}\hfill {∥{\tilde{\mathit{U}}}^{(k+1)}-{\mathit{U}}^{*}∥}_F\le & \left({∥I-P^{-1}B∥}_F+{\tau }_k{∥P^{-1}∥}_F{∥P^{-T}∥}_F{∥B∥}_F\right){∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F\hfill \\ \hfill \le & \left({∥I-P^{-1}B∥}_F+c{\chi }^{{\zeta }^{\left(k\right)}}{∥P^{-1}∥}_F{∥P^{-T}∥}_F{∥B∥}_F\right){∥{\tilde{\mathit{U}}}^{\left(k\right)}-{\mathit{U}}^{*}∥}_F.\hfill \end{array}$$

Since $\underset{k\to \infty }{lim}c{\chi }^{{\zeta }^{\left(k\right)}}=0$, the result follows straightforwardly. □

Remark 3.

Theorem 2 tells us that the convergence rate of IPPIA is asymptotically the same as that of PPIA if the tolerance sequence $\left\{{\tau }_k\right\}$ approaches 0 as k increases. To reduce the computational complexity, it is not necessary for ${\tau }_k$ to approach zero in actual implementations, which will be shown numerically in the next section.

4. Numerical Experiments

In this section, some numerical examples are given to illustrate the effectiveness of our methods. All the numerical experiments were done on a PC with Intel(R) Core(TM) i5-5200U CPU @2.20 GHz by Matlab R2012b.

In our numerical experiments, we used the Bézier curves and the Said–Ball curves to test the following three examples [12,17,26].

For comparison, we also tested the WPIA and the DWPIA. The parameter values $t_i$ and the optimal weight $\omega$ used in our numerical experiments were the same as those used in [12]; i.e., $t_i=i/n$, for $i=0,1,\cdots ,n$ and $\omega =2/(1+{\lambda }_n)$ with ${\lambda }_n$ being the smallest eigenvalue of B. Moreover, we set the parameter $a=1.9$ when we tested the DWPIA.

Example 1.

Consider the lemniscate of Gerono given by

$$(x\left(t\right),y\left(t\right))=\left(cost,sintcost\right),t\in \left[-\frac{\pi }{2},\frac{3\pi }{2}\right].$$

A sequence of $n+1$ points ${\left\{{\mathit{v}}_i\right\}}_{i=0}^n$ is taken from the lemniscate of Gerono as follows:

$${\mathit{v}}_i=\left(x\left(-\frac{\pi }{2}+i\frac{2\pi }{n}\right),y\left(-\frac{\pi }{2}+i\frac{2\pi }{n}\right)\right),i=0,1,\cdots ,n.$$

Example 2.

Consider the helix of radius 5 given by

$$(x\left(t\right),y\left(t\right),z\left(t\right))=\left(5cost,5sint,t\right),t\in [0,6\pi ].$$

A sequence of $n+1$ points ${\left\{{\mathit{v}}_i\right\}}_{i=0}^n$ is taken from the helix as below

$${\mathit{v}}_i=\left(5cos\left(\frac{6\pi }{n}i\right),5sin\left(\frac{6\pi }{n}i\right),\frac{6\pi }{n}i\right),i=0,1,\cdots ,n.$$

Example 3.

Consider the four-leaf clover given by

$$(x\left(t\right),y\left(t\right))=\left(4sin\left(t\right)sin\left(4t\right),4cos\left(t\right)sin\left(4t\right)\right),t\in [0,2\pi ].$$

A sequence of $n+1$ points ${\left\{{\mathit{v}}_i\right\}}_{i=0}^n$ is taken from the four-leaf clover in the following way

$${\mathit{v}}_i=\left(4sin\left(i\frac{2\pi }{n}\right)sin\left(i\frac{8\pi }{n}\right),4cos\left(i\frac{2\pi }{n}\right)sin\left(i\frac{8\pi }{n}\right)\right),i=0,1,\cdots ,n.$$

4.1. Tests for PPIA

Numerically, we found that the half-bandwidth of the preconditioner P in (8) $q\approx n/2$ is a good choice. Therefore, when we used Bézier curves and Said–Ball curves to test Example 1 with $n=10$, we took $q=5$ and 6, respectively. Similarly, when we used Bézier curves and Said–Ball curves to test Example 2 with $n=18$, we took $q=9$ and 12, respectively; when we used Bézier curves and Said–Ball curves to test Example 3 with $n=21$, we took $q=11$ and 14, respectively.

The spectral distributions of $I-B$ and $I-P^{-1}B$ resulting from Bézier curves and Said–Ball curves in Examples 1–3 are illustrated in Figure 1, Figure 2 and Figure 3, respectively. From those figures, we can see that most of the eigenvalues of the preconditioned iteration matrices $I-P^{-1}B$ are far away from 1 and clustered around 0 in each of six cases.

To further illustrate the effectiveness of the PPIA, we list in

Table 1. The spectral radii of iteration matrices of WPIA, DWPIA and PPIA in Examples 1–3.

	Bézier Curve			Said–Ball Curve
	WPIA	DWPIA	PPIA	WPIA	DWPIA	PPIA
Example 1	0.99927450	0.99924155	0.16473699	0.99964097	0.99947135	0.53425173
Example 2	0.99999967	0.99999967	0.80078821	0.99999992	0.99999988	0.82319850
Example 3	1.00000000	0.99999998	0.72221739	1.00000000	1.00000000	0.98124671

the spectral radii of the iteration matrices of the WPIA, DWPIA and PPIA formats resulting from the considered Bézier and Said–Ball curves. Clearly, the spectral radii of the iteration matrices of the PPIA are much less than those of the corresponding WPIA, DWPIA. Additionally, we see that the spectral radii of PPIA decline sharply when the half-bandwidth of the preconditioner P is near $n/2$. Those results entirely coincide with the aforementioned conclusions in Section 2. Therefore, we can expect that the PPIA method has a good convergence behavior.

We denote the number of iterations by k, the k-th interpolation error by ${\epsilon }^{\left(k\right)}$ and the elapsed CPU time after k iterations by $t^{\left(k\right)}$ (in seconds), respectively. We emphasize here that the elapsed CPU time for the same routine in Matlab was different every time. Therefore, we took the average of all runtimes of 100 independent runs as the elapsed CPU time. The detailed numerical results of the PPIA, DWPIA and WPIA resulting from Bézier curves and Said–Ball curves are listed in

Table 2. The approximation errors and runtimes of WPIA, DWPIA and PPIA in Example 1.

Curve	k	WPIA		DWPIA		PPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$4.620\times {10}^{-1}$	–	$4.620\times {10}^{-1}$	–	$4.620\times {10}^{-1}$	–
	1	$2.632\times {10}^{-1}$	$8.29\times {10}^{-4}$	$2.601\times {10}^{-1}$	$4.83\times {10}^{-5}$	$7.500\times {10}^{-4}$	$8.67\times {10}^{-5}$
	10	$5.447\times {10}^{-2}$	$8.96\times {10}^{-4}$	$5.208\times {10}^{-2}$	$8.94\times {10}^{-5}$	$5.017\times {10}^{-14}$	$1.42\times {10}^{-4}$
	411	$8.611\times {10}^{-4}$	$2.65\times {10}^{-3}$	$7.510\times {10}^{-4}$	$1.83\times {10}^{-3}$	–	–
	430	$7.499\times {10}^{-4}$	$2.80\times {10}^{-3}$	–	–	–	–
Said–Ball	0	$5.470\times {10}^{-1}$	–	$5.470\times {10}^{-1}$	–	$5.470\times {10}^{-1}$	–
	1	$4.143\times {10}^{-1}$	$8.42\times {10}^{-4}$	$3.414\times {10}^{-1}$	$4.79\times {10}^{-5}$	$4.048\times {10}^{-3}$	$8.62\times {10}^{-5}$
	10	$8.979\times {10}^{-2}$	$9.13\times {10}^{-4}$	$5.524\times {10}^{-2}$	$8.79\times {10}^{-5}$	$1.733\times {10}^{-8}$	$1.56\times {10}^{-4}$
	20	$5.252\times {10}^{-2}$	$9.27\times {10}^{-4}$	$2.909\times {10}^{-2}$	$1.32\times {10}^{-4}$	$2.257\times {10}^{-14}$	$2.11\times {10}^{-4}$
	310	$6.058\times {10}^{-3}$	$2.20\times {10}^{-3}$	$4.045\times {10}^{-3}$	$1.43\times {10}^{-3}$	–	–
	478	$4.051\times {10}^{-3}$	$2.99\times {10}^{-3}$	–	–	–	–

,

Table 3. The approximation errors and runtime of WPIA, DWPIA and PPIA in Example 2.

Curve	k	WPIA		DWPIA		PPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$4.625\times {10}^0$	–	$4.625\times {10}^0$	–	$4.625\times {10}^0$	–
	1	$3.868\times {10}^0$	$1.08\times {10}^{-3}$	$3.866\times {10}^0$	$6.19\times {10}^{-5}$	$8.396\times {10}^{-5}$	$1.20\times {10}^{-4}$
	10	$9.821\times {10}^{-1}$	$1.26\times {10}^{-3}$	$1.001\times {10}^0$	$1.15\times {10}^{-4}$	$1.943\times {10}^{-6}$	$2.12\times {10}^{-4}$
	40	$2.014\times {10}^{-1}$	$1.38\times {10}^{-3}$	$2.012\times {10}^{-1}$	$3.55\times {10}^{-4}$	$2.858\times {10}^{-9}$	$5.47\times {10}^{-4}$
	163,447	$8.417\times {10}^{-5}$	$1.22\times {10}^0$	$8.396\times {10}^{-5}$	$1.24\times {10}^0$	–	–
	163,902	$8.396\times {10}^{-5}$	$1.24\times {10}^0$	–	–	–	–
Said–Ball	0	$5.698\times {10}^0$	–	$5.698\times {10}^0$	–	$5.698\times {10}^0$	–
	1	$6.614\times {10}^0$	$1.19\times {10}^{-3}$	$5.585\times {10}^0$	$6.15\times {10}^{-5}$	$3.281\times {10}^{-4}$	$1.23\times {10}^{-4}$
	10	$1.812\times {10}^0$	$1.30\times {10}^{-3}$	$1.822\times {10}^0$	$1.19\times {10}^{-4}$	$7.809\times {10}^{-6}$	$2.21\times {10}^{-4}$
	40	$3.319\times {10}^{-1}$	$1.37\times {10}^{-3}$	$2.429\times {10}^{-1}$	$3.24\times {10}^{-4}$	$2.078\times {10}^{-8}$	$5.60\times {10}^{-4}$
	72,224	$3.281\times {10}^{-4}$	$5.36\times {10}^{-1}$	$3.358\times {10}^{-4}$	$5.31\times {10}^{-1}$	–	–
	74,606	–	–	$3.281\times {10}^{-4}$	$5.97\times {10}^{-1}$	–	–

and

Table 4. The approximation errors and runtime of WPIA, DWPIA and PPIA in Example 3.

Curve	k	WPIA		DWPIA		PPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$4.625\times {10}^0$	–	$4.625\times {10}^0$	–	$4.625\times {10}^0$	–
	1	$9.493\times {10}^{-1}$	$9.54\times {10}^{-4}$	$9.487\times {10}^{-1}$	$5.39\times {10}^{-5}$	$6.363\times {10}^{-7}$	$1.20\times {10}^{-4}$
	10	$6.657\times {10}^{-2}$	$1.14\times {10}^{-3}$	$7.524\times {10}^{-2}$	$9.74\times {10}^{-5}$	$2.104\times {10}^{-10}$	$1.95\times {10}^{-4}$
	40	$6.615\times {10}^{-3}$	$1.30\times {10}^{-3}$	$6.903\times {10}^{-3}$	$2.67\times {10}^{-4}$	$1.245\times {10}^{-14}$	$4.21\times {10}^{-4}$
	49,109	$6.363\times {10}^{-7}$	$2.34\times {10}^{-1}$	$6.405\times {10}^{-7}$	$2.39\times {10}^{-1}$	–	–
	49,369	–	–	$6.363\times {10}^{-7}$	$2.44\times {10}^{-1}$	–	–
Said–Ball	0	$4.049\times {10}^0$	–	$4.049\times {10}^0$	–	$4.049\times {10}^0$	–
	1	$2.375\times {10}^0$	$9.30\times {10}^{-4}$	$3.176\times {10}^0$	$7.08\times {10}^{-5}$	$2.127\times {10}^{-5}$	$1.25\times {10}^{-4}$
	10	$6.678\times {10}^{-2}$	$1.11\times {10}^{-3}$	$4.639\times {10}^{-1}$	$1.06\times {10}^{-4}$	$1.638\times {10}^{-7}$	$1.88\times {10}^{-4}$
	40	$2.519\times {10}^{-2}$	$1.25\times {10}^{-3}$	$3.202\times {10}^{-2}$	$2.57\times {10}^{-4}$	$7.217\times {10}^{-8}$	$4.24\times {10}^{-4}$
	2776	$4.285\times {10}^{-5}$	$1.44\times {10}^{-2}$	$2.128\times {10}^{-5}$	$1.39\times {10}^{-2}$	–	–
	5354	$2.127\times {10}^{-5}$	$2.60\times {10}^{-2}$	–	–	–	–

, respectively. Under the requirement of the same precision, the number of iterations by PPIA is much less than those of WPIA and DWPIA, and the runtime by PPIA is much less than the runtimes of WPIA and DWPIA. This means that our PPIA can accelerate the convergence rate of PIA significantly.

In

Table 5. The detailed runtime of WPIA, DWPIA and PPIA in Example 3.

Curve	${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$
Bézier	$1.12\times {10}^{-3}$	$1.83\times {10}^{-5}$	$2.75\times {10}^{-5}$
Said–Ball	$1.06\times {10}^{-3}$	$1.70\times {10}^{-5}$	$2.48\times {10}^{-5}$

, we show the detailed runtimes of WPIA, DWPIA and PPIA for when we tested Example 3 by using Bézier curves and Said–Ball curves. We denote the runtime of computing the optimal weight for WPIA by $T_1$, the runtime of computing different weights for DWPIA by $T_2$ and the runtime of constructing the preconditioner P for PPIA by $T_3$, respectively. We can see that the cost of constructing P is small enough and can be neglected.

We remark here that all curves in Figure 4, Figure 5 and Figure 6 are those obtained at 10th iteration step by employing different iterations with different bases, respectively. In Figure 4, Figure 5 and Figure 6, the curves generated by PPIA and WPIA formats are denoted by the solid and dashed lines respectively. We can see from Figure 4, Figure 5 and Figure 6 that the approximations obtained by PPIA are better than those obtained by WPIA, no matter whether Bézier curves or Said–Ball curves are used. Clearly, the PPIA converges faster than the WPIA.

4.2. Tests for IPPIA

From the above numerical experiments, we know that the PPIA works much better than the WPIA and DWPIA if the number of the interpolating points is moderate. However, when we fit 53 points sampled from the lemniscate of Gerono, the helix and the four-leaf clover in Examples 1–3 by using Bézier curves and Said–Ball curves, we found that the WPIA and DWPIA converge very slowly. At the same time, we found that the PPIA does not work, because the condition number of P is very large and the solution to $P{\Delta }^{\left(k\right)}={\widehat{\Delta }}^{\left(k\right)}$ cannot be directly obtained. In this case, we can employ the IPPIA to get the corresponding fitting curves. The numerical results (including the number of iterations, the approximation errors and the runtime) obtained by the WPIA, DWPIA and IPPIA are listed in

Table 6. Numerical results when fitting 53 points sampled from the lemniscate of Gerono.

Curve	k	WPIA		DWPIA		IPPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$1.606\times {10}^{-1}$	–	$1.606\times {10}^{-1}$	–	$1.606\times {10}^{-1}$	–
	1	$7.987\times {10}^{-2}$	–	$7.289\times {10}^{-2}$	–	$1.507\times {10}^{-3}$	–
	3	$4.125\times {10}^{-2}$	–	$3.313\times {10}^{-2}$	–	$4.537\times {10}^{-7}$	–
	5	$3.018\times {10}^{-2}$	–	$2.126\times {10}^{-2}$	–	$3.378\times {10}^{-9}$	$2.39\times {10}^{-3}$
	137	$3.451\times {10}^{-9}$	$2.66\times {10}^{-3}$	$1.192\times {10}^{-8}$	–	–	–
	381	–	–	$3.456\times {10}^{-9}$	$2.94\times {10}^{-3}$	–	–
Said–Ball	0	$4.308\times {10}^{-1}$	–	$4.308\times {10}^{-1}$	–	$4.308\times {10}^{-1}$	–
	1	$3.789\times {10}^{-1}$	–	$4.225\times {10}^{-1}$	–	$1.265\times {10}^{-2}$	–
	2	$3.602\times {10}^{-1}$	–	$3.139\times {10}^{-1}$	–	$2.369\times {10}^{-4}$	–
	3	$9.127\times {10}^{-2}$	–	$1.305\times {10}^{-1}$	–	$9.892\times {10}^{-6}$	–
	4	$8.820\times {10}^{-2}$	–	$1.360\times {10}^{-1}$	–	$3.396\times {10}^{-6}$	$2.42\times {10}^{-3}$
	291	$3.382\times {10}^{-6}$	$3.43\times {10}^{-3}$	$6.918\times {10}^{-6}$	–	–	–
	358	–	–	$3.382\times {10}^{-6}$	$2.82\times {10}^{-3}$	–	–

,

Table 7. Numerical results when fitting 53 points sampled from the helix.

Curve	k	WPIA		DWPIA		IPPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$2.882\times {10}^0$	–	$2.882\times {10}^0$	–	$2.882\times {10}^0$	–
	1	$7.021\times {10}^{-1}$	–	$6.469\times {10}^{-1}$	–	$5.007\times {10}^{-2}$	–
	2	$3.359\times {10}^{-1}$	–	$2.827\times {10}^{-1}$	–	$5.533\times {10}^{-4}$	–
	4	$1.997\times {10}^{-1}$	–	$1.506\times {10}^{-1}$	–	$6.341\times {10}^{-7}$	$3.03\times {10}^{-3}$
	233	$1.536\times {10}^{-5}$	–	$6.349\times {10}^{-7}$	$2.56\times {10}^{-3}$	–	–
	316	$6.416\times {10}^{-7}$	$4.23\times {10}^{-3}$	–	–	–	–
Said–Ball	0	$5.773\times {10}^0$	–	$5.773\times {10}^0$	–	$5.773\times {10}^0$	–
	1	$6.681\times {10}^0$	–	$5.763\times {10}^0$	–	$2.508\times {10}^{-1}$	–
	3	$4.316\times {10}^0$	–	$5.213\times {10}^0$	–	$2.913\times {10}^{-4}$	–
	5	$2.086\times {10}^0$	–	$2.062\times {10}^0$	–	$9.780\times {10}^{-5}$	$3.26\times {10}^{-3}$
	359	$2.110\times {10}^{-4}$	–	$9.750\times {10}^{-5}$	$3.84\times {10}^{-3}$	–	–
	438	$9.765\times {10}^{-5}$	$5.40\times {10}^{-3}$	–	–	–	–

and

Table 8. Numerical results when fitting 53 points sampled from the four-leaf clover.

Curve	k	WPIA		DWPIA		IPPIA
		${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$	${\mathbf{\epsilon }}^{\left(\mathit{k}\right)}$	${\mathit{t}}^{\left(\mathit{k}\right)}$
Bézier	0	$1.219\times {10}^0$	–	$1.219\times {10}^0$	–	$1.219\times {10}^0$	–
	1	$3.815\times {10}^{-1}$	–	$3.460\times {10}^{-1}$	–	$1.967\times {10}^{-2}$	–
	2	$2.270\times {10}^{-1}$	–	$1.915\times {10}^{-1}$	–	$1.526\times {10}^{-4}$	–
	4	$1.226\times {10}^{-1}$	–	$9.344\times {10}^{-2}$	–	$2.802\times {10}^{-7}$	$2.46\times {10}^{-3}$
	125	$1.378\times {10}^{-4}$	–	$2.791\times {10}^{-7}$	$1.18\times {10}^{-3}$	–	–
	316	$2.799\times {10}^{-7}$	$3.60\times {10}^{-3}$			–	–
Said–Ball	0	$3.338\times {10}^0$	–	$3.338\times {10}^0$	–	$3.338\times {10}^0$	–
	1	$2.227\times {10}^0$	–	$3.025\times {10}^0$	–	$9.302\times {10}^{-2}$	–
	2	$2.271\times {10}^0$	–	$3.543\times {10}^0$	–	$4.527\times {10}^{-3}$	–
	4	$4.344\times {10}^{-1}$	–	$1.495\times {10}^0$	–	$3.743\times {10}^{-5}$	$2.58\times {10}^{-3}$
	218	$3.741\times {10}^{-5}$	$2.99\times {10}^{-3}$	$1.505\times {10}^{-4}$	–	–	–
	379	–	–	$3.712\times {10}^{-5}$	$3.02\times {10}^{-3}$	–	–

respectively, where the parameters in IPPIA are taken as $q=26$, $\theta =0.98$ and ${\tau }_k={10}^{-2}$.

From Table 6, Table 7 and Table 8, we can see that the IPPIA converges much faster than the WPIA and DWPIA. For example, when we fit the points sampled from the lemniscate of Gerono by using Bézier curves, we observed that the errors caused by performing 5 IPPIA iterations are about the same as those caused by performing 137 WPIA iterations or 381 DWPIA iterations. Again, when we fit 53 points sampled from the lemniscate of Gerono by using Said–Ball curves, we found that the errors caused by performing 4 IPPIA iterations are about the same as those caused by performing 291 WPIA iterations or 358 DWPIA iterations. Furthermore, under the requirement of the same accuracy, the runtime of IPPIA iterations was less than that of WPIA iterations and it was comparable to that of DWPIA iterations.

When fitting 53 points sampled from the lemniscates of Gerono, helix and four-leaf clover, we sketched the Bézier curves and the Said–Ball curves generated by performing five IPPIA and WPIA iterations in Figure 7, Figure 8 and Figure 9, respectively. The solid and dashed lines in Figure 7, Figure 8 and Figure 9 represent the curves generated by IPPIA and WPIA respectively. Clearly, the IPPIA outperforms the WPIA in the sense of the convergence rate.

5. Conclusions

In this paper, we have exploited the properties of collocation matrix B, resulting from an NTP basis, to design an efficient preconditioner P for PIA. The main tool to derive our preconditioner is the diagonally compensated reduction. By applying that P to PIA, we have developed a novel version of PIA, called the PPIA. By both convergence theory and numerical experiments, we have shown that the PPIA outperforms the WPIA and DWPIA, in the sense of convergence rate and elapsed CPU time. In particular, to reduce the computational complexity of PPIA, we have also proposed an inexact variant of PPIA (i.e., IPPIA) and shown its effectiveness, especially when the number of interpolating points is large.

As mentioned in Section 1, the PIA has extensive successful applications. Therefore, we can expect that our methods have a wide range of potential applications.