# The Convergence of Least-Squares Progressive Iterative Approximation with Singular Iterative Matrix

Developed in [Deng and Lin, 2014], Least-Squares Progressive Iterative Approximation (LSPIA) is an efficient iterative method for solving B-spline curve and surface least-squares fitting systems. In [Deng and Lin 2014], it was shown that LSPIA is convergent when the iterative matrix is nonsingular. In this paper, we will show that LSPIA is still convergent even the iterative matrix is singular.

## Authors

• 7 publications
• 15 publications
• 2 publications
08/18/2019

### On a progressive and iterative approximation method with memory for least square fitting

In this paper, we present a progressive and iterative approximation meth...
09/02/2019

### Implicit Progressive-Iterative Approximation for Curve and Surface Reconstruction

Implicit curve and surface reconstruction attracts the attention of many...
03/05/2020

### A generalized projection iterative methods for solving non-singular linear systems

In this paper, we propose and analyze iterative method based on projecti...
09/01/2020

### GMRES on singular systems revisited

In [Hayami K, Sugihara M. Numer Linear Algebra Appl. 2011; 18:449–469], ...
12/13/2021

### On using the complex step method for the approximation of Fréchet derivatives of matrix functions in automorphism groups

We show, that the complex step approximation Im(f(A+ihE))/h to the Fréch...
05/10/2020

### Partial least squares for function-on-function regression via Krylov subspaces

People employ the function-on-function regression (FoFR) to model the re...
07/28/2018

### Holographic Sensing

Holographic representations of data encode information in packets of equ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Least-squares fitting is a commonly employed approach in engineering applications and scientific research, including geometric modeling. With the advent of big data era, least-squares fitting systems with singular coefficient matrices often appear, when the number of the fitted data points is very large, or there are “holes” in the fitted data points. LSPIA deng2014progressive is an efficient iterative method for least-squares B-spline curve and surface fitting brandt2015optimal . In Ref. deng2014progressive , it was shown that LSPIA is convergent when the iterative matrix is nonsingular. In this paper, we will show that, when the iterative matrix is singular, LSPIA is still convergent. This property of LSPIA will promote its applications in large scale data fitting.

The motivation of this paper comes from our research practices, where some singular least-squares fitting systems emerge. For examples, in generating trivariate B-spline solids by fitting tetrahedral meshes lin2015constructing , and in fitting images with holes by T-spline surfaces lin2013efficient , coefficient matrices of least-squares fitting systems are singular. There, LSPIA was employed to solve the least-squares fitting systems, and converged to stable solutions. However, in Ref. lin2015constructing ; lin2013efficient , convergence of LSPIA for solving singular linear systems was not proved.

The progressive-iterative approximation (PIA) method was first developed in (lin2004constructing, ; lin2005totally, )

, which endows iterative methods with geometric meanings, so it is suitable to handle geometric problems appearing in the field of geometric design. It was proved that the PIA method is convergent for B-spline fitting

(lin2011extended, ; deng2014progressive, ), NURBS fitting (shi06iterative, ), T-spline fitting (lin2013efficient, ), subdivision surface fitting (cheng2009loop, ; fan2008subdivision, ; chen2008progressive, ), as well as curve and surface fitting with totally positive basis (lin2005totally, )

. The iterative format of geometric interpolation (GI)

maekawa2007interpolation is similar as that of PIA. While PIA depends on the parametric distance, the iterations of GI rely on the geometric distance. Moreover, the PIA and GI methods have been employed in some applications, such as reverse engineering (kineri2012b, ; yoshihara2012topologically, ), curve design (okaniwa2012uniform, ), surface-surface intersection (lin2014affine, ), and trivariate B-spline solid generation (lin2015constructing, ), etc.

The structure of this paper is as follows. In Section 2, we show the convergence of LSPIA with singular iterative matrix. In Section 3, an example is illustrated. Finally, Section 4 concludes the paper.

## 2 The iterative format and its convergence analysis

To integrate the LSPIA iterative formats for B-spline curves, B-spline patches, trivariate B-spline solids, and T-splines, their representations are rewritten as the following form,

 P(t)=n∑i=0PiBi(t). (1)

Specifically, T-spline patches sederberg2004t and trivariate T-spline solids zhang2012solid can be naturally written as the form (1). Moreover,

• If  (1) is a B-spline curve, then, is a scalar , and , where is a B-spline basis function.

• If  (1) is a B-spline patch with control points, then, , and , where and are B-spline basis functions. In the control net of the B-spline patch, the original index of is , and the original index of is , where represents the maximum integer not exceeding , and is the module of by .

• If is a trivariate B-spline solid with control points, then , and . In the control net of the trivariate B-spline solid, the original index of is , the original index of is , and the original index of is .

Suppose we are given a data point set

 {Qj=(xi,yi,zi),j=0,1,⋯,m}, (2)

each of which is assigned a parameter . Let the initial form be,

 P(0)(t)=n∑i=0P(0)iBi(t), n≤m. (3)

It should be noted that, though the initial control points

are usually chosen from the given data points, the initial control points are unrelated to the convergence of LSPIA. To perform LSPIA iterations, data points are classified into groups. All of data points with parameters

satisfying are classified into the group, corresponding to the control point (3).

After the iteration of the LSPIA, the form is generated,

 P(k)(t)=n∑i=0P(k)iBi(t).

To produce the form , we first calculate the difference vectors for data points (DVD) (Fig. 1),

 δ(k)j=Qj−P(k)(tj), j=0,1,⋯,m.

And then, two procedures are performed, i.e., vector distribution and vector gathering (Fig. 1). In the vector distribution procedure, all of DVDs corresponding to data points in the group are distributed to the control point ; in the vector gathering procedure, all of DVDs distributed to the control point are weighted averaged to generate the difference vector for control point (DVC) (Fig. 1),

 Δ(k)i=∑j∈IiBi(tj)δj∑j∈IiBi(tj), i=0,1,⋯,n,

where is the index set of the data points in the group. Then, the new control point is produced by adding the DVC to , i.e.,

 P(k+1)i=P(k)i+Δ(k)i, i=0,1,⋯,n, (4)

 P(k+1)(t)=n∑i=0P(k+1)iBi(t). (5)

In this way, we get a sequence of iterative forms . Let,

 P(k) =[P(k)0,P(k)1,⋯,P(k)n]T, (6) Q =[Q0,Q1,⋯,Qm]T. (7)

From Eq. (4), it follows,

 P(k+1)i=P(k)i+1∑j∈IiBi(tj)∑j∈IiBi(tj)(Qj−P(k)(tj))=P(k)i+1∑j∈IiBi(tj)∑j∈IiBi(tj)(Qj−n∑l=0P(k)lBl(tj))

Therefore, we get the LSPIA iterative format in matrix form,

 P(k+1)=P(k)+ΛAT(Q−AP(k)), k=0,1,⋯ (8)

where, is a diagonal matrix, and,

 A=⎡⎢ ⎢ ⎢ ⎢ ⎢⎣B0(t0)B1(t0)⋯Bn(t0)B0(t1)B1(t1)⋯Bn(t1)⋮⋮⋮B0(tm)B1(tm)⋯Bn(tm)⎤⎥ ⎥ ⎥ ⎥ ⎥⎦(m+1)×(n+1).
###### Remark 0

The iterative format (8) is slightly different from that developed in Ref. deng2014progressive , where diagonal elements of the diagonal matrix are equal to each other. Although the difference of their iterative formats is slight, the convergence analysis of the iterative format (8) is a bit more difficult lin2013efficient .

###### Remark 0

Because diagonal elements of the diagonal matrix in the iterative format (8) are all positive, the diagonal matrix is nonsingular.

To show the convergence of the LSPIA iterative format (8), it is rewritten as,

 P(k+1)=(I−ΛATA)P(k)+ΛATQ. (9)

In Ref. deng2014progressive , it was shown that, when the iterative matrix is nonsingular, the LSPIA iterative format is convergent. In the following, we will show that, even the matrix is not of full rank, and then is singular, the iterative format (8) is still convergent.

We first show some lemmas.

###### Lemma 0

The eigenvalues

of the matrix are all real, and satisfy .

Proof: On one hand, suppose is an arbitrary eigenvalue of the matrix

with eigenvector

, i.e.,

 ΛATAv=λv. (10)

By multiplying at both sides of Eq. (10), we have,

 AΛAT(Av)=λ(Av).

It means that is also an eigenvalue of the matrix with eigenvector . Moreover, , because,

 xTAΛATx=xTAΛ12(Λ12)TATx=(xTAΛ12)(xTAΛ12)T≥0,

the matrix is a positive semidefinite matrix. Eigenvalues of a semidefinite matrix are all nonnegative, so is real, and .

On the other hand, because the B-spline basis functions are nonnegative and form a partition of unity, it holds, . Together with , we have,

 ∥∥ΛATA∥∥∞≤∥∥ΛAT∥∥∞∥A∥∞=1.

Therefore, the eigenvalue of matrix satisfies,

 λ≤∥∥ΛATA∥∥∞≤1.

In conclusion, eigenvalues of the matrix are all real, and satisfy .

Because  (9) is singular, is also singular, and then is its eigenvalue. The following lemma deals with the relationship between the algebraic multiplicity and geometric multiplicity of the zero eigenvalue of .

###### Remark 0

In this paper, we assume that the dimension of the zero eigenspace of

is . So, the rank of the matrix is

 rank(ATA)=n−n0+1.

Because is nonsingular (refer to Remark 2), we have

 rank(ΛATA)=n−n0+1.
###### Lemma 0

The algebraic multiplicity of the zero eigenvalue of matrix is equal to its geometric multiplicity.

Proof: The proof consists of three parts.

(1) The algebraic multiplicity of the zero eigenvalue of matrix is equal to its geometric multiplicity. Because

is a positive semidefinite matrix, it is a diagonalizable matrix. Then, for any eigenvalue of

(including the zero eigenvalue), its algebraic multiplicity is equal to its geometric multiplicity. In Remark 4, we assume that the dimension of the zero eigenspace of , i.e., the geometric multiplicity of zero eigenvalue of , is . So, the algebraic multiplicity and geometric multiplicity of zero eigenvalue of are both .

(2) The geometric multiplicity of the zero eigenvalue of matrix is equal to that of matrix . Denote the eigenspaces of matrices and associated with the zero eigenvalue as and , respectively. The geometric multiplicities of the zero eigenvalue of matrices and are dimensions of and , respectively.

Note that the matrix is nonsingular (Remark 2). On one hand, , leading to . So, . On the other hand, , resulting in . So, . In conclusion, . Therefore, the geometric multiplicity of the zero eigenvalue of matrix is equal to that of matrix .

(3) The algebraic multiplicity of the zero eigenvalue of matrix is equal to that of matrix . Denote as an identity matrix,

 ATA=(bij)(n+1)×(n+1), and, Λ=diag(d0,d1,⋯,dn),

where .

The characteristic polynomial of and can be written as (horn1985matrix, , pp.42),

 pATA(λ)=det(λI−ATA)=λn+1−E1(ATA)λn+E2(ATA)λn−1+⋯+(−1)n+1En+1(ATA), (11)

and,

 pΛATA(λ)=det(λI−ΛATA)=λn+1−E1(ΛATA)λn+E2(ΛATA)λn−1+⋯+(−1)n+1En+1(ΛATA), (12)

where are the sums of the principal minors of , and are the sums of the principal minors of .

On one hand, because the algebraic multiplicity of zero eigenvalue of is (see Part (1)), its characteristic polynomial (11) can be represented as,

 pATA=λn0(λn−n0+1−E1(ATA)λn−n0+⋯+(−1)n−n0+1En−n0+1(ATA)).

where . Moreover, because is positive semi-definite, all of its principal minors are nonnegative. Therefore, we have . Consequently, all of principal minors of are nonnegative, and there is at least one principal minor of is positive.

On the other hand, because (Remark 4), all of () principal minors of are zero. Therefore,

 El(ΛATA)=0, l>n−n0+1. (13)

Denote and are the principal minors of and , respectively. Now, consider a principal minor of .

 MΛATA(i1,i2,⋯,ik)=det⎛⎜⎝di1bi1,i1⋯di1bi1,ik⋯⋯⋯dikbik,i1⋯dikbik,ik⎞⎟⎠=(k∏j=1dij)MATA(i1,i2,⋯,ik),

where (Remark 2). In other words, the principal minor of is the product of a principal minor of and a positive value . Together with that all of principal minors of are nonnegative, and there is at least one principal minor of is positive, the sum of all principal minors of , namely, , is positive. That is,

 En−n0+1(ΛATA)>0. (14)

By Eqs. (13) and (14), the characteristic polynomial of  (12) can be transformed as,

 pΛATA=λn0(λn−n0+1−E1(ΛATA)λn−n0+⋯+(−1)n−n0+1En−n0+1(ΛATA)),

where . It means that the algebraic multiplicity of zero eigenvalue of is , equal to the algebraic multiplicity of zero eigenvalue of .

Combing results of part (1)-(3), we have shown that the algebraic multiplicity of the zero eigenvalue of matrix is equal to its geometric multiplicity.

Denote as a matrix block,

 Jr(a,b)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝abab\huge{0}⋱⋱\huge{0}⋱ba⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠r×r. (15)

Specifically, is a Jordan block. Lemma 3 and 5 result in Lemma 6 as follows.

###### Lemma 0

The Jordan canonical form of matrix  (9) can be written as,

 J=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Jn1(λ1,1)Jn2(λ2,1)\huge{0}⋱Jnk(λk,1)0\huge{0}⋱0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠(n+1)×(n+1), (16)

where are nonzero eigenvalues of , which need not be distinct, and  (15) is an Jordan block, .

Proof: Based on Lemma 3, eigenvalues of are all real and lie in , so the Jordan canonical form of can be written as,

 J=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Jn1(λ1,1)⋱\huge{0}Jnk(λk,1)Jm1(0,1)\huge{0}⋱Jml(0,1)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠(n+1)×(n+1),

where are nonzero eigenvalues of , which need not be distinct, and  (15) is an Jordan block, ; is an Jordan block corresponding to the zero eigenvalue of , .

According to the theory on Jordan canonical form (horn1985matrix, , p.129), the number of Jordan blocks corresponding to an eigenvalue is the geometric multiplicity of the eigenvalue, and the sum of orders of all Jordan blocks corresponding to an eigenvalue equals its algebraic multiplicity. Based on Lemma 5, the algebraic multiplicity of the zero eigenvalue of matrix is equal to its geometric multiplicity, so the Jordan blocks corresponding to the zero eigenvalue of are all matrix . This proves Lemma 6.

Denote as the Moore-Penrose (M-P) pseudo-inverse of the matrix . We have the following lemma.

###### Lemma 0

There exists an orthogonal matrix

, such that,

 VT(ATA)+(ATA)V=diag(1,1,⋯,1n−n0+1,0,0,⋯,0n0). (17)

Proof: Because (Remark 4), and

is a positive semidefinite matrix, it has singular value decomposition (SVD),

 ATA=Vdiag(δ1,δ2,⋯,δn−n0+1,0,⋯,0n0)VT, (18)

where is an orthogonal matrix, are singular values of . Then, the M-P pseudo-inverse of is,

 (ATA)+=Vdiag(1δ1,1δ2,⋯,1δn−n0+1,0,⋯,0n0)VT.

Therefore,

 (ATA)+(ATA)=Vdiag(1,1,⋯,1n−n0+1,0,0,⋯,0n0)VT,

where is an orthogonal matrix.

Based on the Lemmas above, we can show the convergence of the iterative format (9) when is singular.

###### Theorem 8

When  (9) is singular, the iterative format (9) is convergent.

Proof: By Lemma 6, the Jordan canonical form of matrix  (9) is  (16

). Then, there exists a invertible matrix

, such that,

 ΛATA=W−1JW.

Therefore (refer to Eq. (16)),

 \footnotesizeI−ΛATA=W−1⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Jn1(1−λ1,−1)Jn2(1−λ2,−1)\Large{0}⋱Jnk(1−λk,−1)1\Large{0}⋱1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠W,

where . Then, together with Lemma 7, it holds,

 liml→∞(I−ΛATA)l=W−1diag(0,⋯,0n−n0+1,1,⋯,1n0)W=I−W−1diag(1,⋯,1n−n0+1,0,⋯,0n0)W=I−W−1VT(ATA)+(ATA)VW=I−(VW)−1(ATA)+(ATA)(VW) (19)

Now, consider the linear system (refer to Eq. (7)). It has solutions if and only if james1978generalised ,

 (ATA)(ATA)+(ATQ)=ATQ. (20)

Subtracting from both sides of the iterative format (9), together with Eq. (20), we have,

 P(k+1)−(ATA)+ATQ=(I−ΛATA)P(k)+ΛATQ−(ATA)+ATQ=(I−ΛATA)P(k)+Λ(ATA)(ATA)+(ATQ)−(ATA)+ATQ=(I−ΛATA)P(k)−(I−ΛATA)(ATA)+ATQ=(I−ΛATA)(P(k)−(ATA)+ATQ)=(I−ΛATA)k+1(P(0)−(ATA)+ATQ). (21)

Owing to Eq. (19), it follows,

 P(∞)−(ATA)+ATQ=limk→∞(P(k+1)−(ATA)+ATQ)=limk→∞(I−ΛATA)k+1(P(0)−(ATA)+ATQ)=(I−(VW)−1(ATA)+(ATA)(VW))(P(0)−(ATA)+ATQ) (22)

By simple computation, Eq. (22) changes to,

 P(∞)=(VW)−1(ATA)+(ATA)VW(ATA)+ATQ+(I−(VW)−1(ATA)+(ATA)(VW))P(0). (23)

Therefore, the iterative format (9) is convergent when is singular. Theorem 8 is proved.

###### Remark 0

Returning to Eq. (23), if is the inverse matrix of , i.e., , it becomes,

 P(∞)=(ATA)+(ATA)(ATA)+ATQ+(I−(ATA)+(ATA))P(0)=(ATA)+ATQ+(I−(ATA)+(ATA))P(0), (24)

where, is an arbitrary initial value. Eq. (24) is the M-P pseudo-inverse solution of the linear system , which is the normal equation of the least-squares fitting to the data points (2). Because is an arbitrary value, there are infinite solutions to the normal equation . Within these solutions, is the one with minimum Euclidean norm horn1985matrix .

Actually, if diagonal elements of matrix  (9) are equal to each other, denoting as , iterative format (9) can be written as,

 P(k+1)=(I−αATA)P(k)+αATQ. (25)

In this case, we have the following theorem.

###### Theorem 10

If is singular, and the spectral radius , the iterative format (25) converges to the M-P pseudo-inverse solution of the linear system . Moreover, if the initial value , the iterative format  (25) converges to , i.e., the M-P pseudo-inverse solution of the linear system with the minimum Euclidean norm.

Proof: Because is both a normal matrix and a positive semidefinite matrix, its eigen decomposition is the same as its singular value decomposition horn1985matrix , with the form presented in Eq. (18). So, we have,

 αATA=Vdiag(αδ1,αδ2,⋯,αδn−n0+1,0,⋯,0n0)VT,

where is an orthogonal matrix, and are both the nonzero eigenvalues and nonzero singular values of . Because , it holds

 0<αδi≤1,i=1,2,⋯,n−n0+1.

Then, based on Lemma 7, we have,

 liml→∞(I−αATA)l=Vdiag(0,⋯,0n−n0+1,1,⋯,1n0)VT=I−Vdiag(1,⋯,1n−n0+1,0,⋯,0n0)VT=I−VVT(ATA)+(ATA)VVT=I−(ATA)+(ATA).

Same as the deduction in the proof of Theorem 8 (Eqs. (21) (22)), we have,

 P(k+1)−(ATA)+ATQ=(I−αATA)k+1(P(0)−(ATA)+ATQ),

and,

 P(∞)−(ATA)+ATQ=limk→∞(P(k+1)−(ATA)+ATQ)=limk→∞(I−αATA)k+1(P(0)−(ATA)+ATQ)=(I−(ATA)+(ATA))(P(0)−(ATA)+ATQ).

Therefore,

 P(∞)=(ATA)+ATQ+(I−(ATA)+(ATA))P(0),

where