Abstract: Progressive-iterative approximation(PIA) is an intuitive and effective method for data fitting. Classical PIA method requires that the number of control points is equal to the number of the data points. It is not suitable for fitting mass data. In order to improve the classical PIA method, the algorithm for fitting data points with triangular surfaces based on PIA method is studied, especially for the low-degree case usually used in practice. It is proved that the quadratic, cubic and quartic non-uniform triangular Bézier surfaces have the property of progressive-iterative approximation for least square fitting(LSPIA). And the limit of the sequence of triangular Bézier surfaces obtained by iteration is just the least square fitting of the data points. Meanwhile, a method is provided to show how to choose the value of the weight so that the iteration has the fastest convergence speed. A numerical example is presented to validate the effectiveness of the LSPIA method.