Data fitting is an extensively employed modeling tool in geometric design. With the advent of the big data era, the data sets to be fitted are made larger and larger, leading to more and more leastsquares fitting systems with singular coefficient matrices. LSPIA (least-squares progressive iterative approximation) is an efficient iterative method for the least-squares fitting. However, the convergence of LSPIA for the singular least-squares fitting systems remains as an open problem. In this paper, the authors showed that LSPIA for the singular least-squares fitting systems is convergent. Moreover, in a special case, LSPIA converges to the Moore-Penrose (M-P) pseudo-inverse solution to the leastsquares fitting result of the data set. This property makes LSPIA, an iterative method with clear geometric meanings, robust in geometric modeling applications. In addition, the authors discussed some implementation detail of LSPIA, and presented an example to validate the convergence of LSPIA for the singular least-squares fitting systems.