# Sequential Quadratic Programming Method for Nonlinear Least Squares Estimation and Its Application

Authors: Zhengqing Fu; Goulin Liu; Lanlan Guo
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this study, we propose a direction-controlled nonlinear least squares estimation model that combines the penalty function and sequential quadratic programming. The least squares model is transformed into a sequential quadratic programming model, allowing for the iteration direction to be controlled. An ill-conditioned matrix is processed by our model; the least squares estimate, the ridge estimate, and the results are compared based on a combination of qualitative and quantitative analyses. For comparison, we use two equality indicators: estimated residual fluctuation of different methods and the deviation between estimated and true values. The root-mean-squared error and standard deviation are used for quantitative analysis. The results demonstrate that our proposed model has a smaller error than other methods. Our proposed model is thereby found to be effective and has high precision. It can obtain more precise results compared with other classical unwrapping algorithms, as shown by unwrapping using both simulated and real data from the Jining area in China.


---

## Body

## 1. Introduction

Most problems faced during the course of processing measurement data are nonlinear in nature. In the classical least squares approach, the nonlinear approximation function is expanded at an approximate value. This transformation can linearize the data, allowing the problem to be solved through the linear least squares approach. The precondition of this model is that initial parameter values must reach the adjusted value; otherwise, the model error is significant.Nonlinear least squares is an optimization method used to solve nonlinear problems [1–4]. Sequential quadratic programming is another important and effective optimization method [5–7]. A nonlinear least squares problem can be transformed into a sequential quadratic programming model and then solved [8–11].A direction-controlled estimation model is proposed in this paper. The nonlinear least squares problem is transformed into a sequential quadratic programming model. The iterative point enters the feasible region via the penalty function, after which a problem’s optimal solution can be obtained through sequential quadratic programming [12–14]. In the model, the penalty function is used to reduce the constraints of the initial values [15–19]. Sequential quadratic programming possesses quadratic convergence. The model is contrasted with least squares, and ridge estimate and simulated data are used to check model performance [20–23]. The root-mean-squared error (RMSE) and standard deviation between the estimated and true values, and the ε value representing the quality of unwrapping, are taken as indices of performance [24, 25]. We conducted experiments that showed the feasibility and effectiveness of the model. Unwrapping experiments using real data from the Jining area in China show that our proposed algorithm achieves more precise results than the least squares unwrapping algorithms. Quantitative indexes include differences in RMSE between rewrapped results and the original wrapped phase, computation time, and ε~ values [26–28]. The sequential quadratic programming method achieves better results with respect to three indices.


## 2. Sequential Quadratic Programming Model

Quadratic optimization with equality constraints is defined as follows:(1)minfx=12xTGx+gTxs.t.Ax=b,x≥0where G∈Rn×n is a symmetric matrix, g∈Rn×1,A∈Rm×n,b∈Rm×1.The Lagrange function of (1) is(2)Lx,u=12xTGx+gTx+vTAx-bwhere v is the Lagrange multiplier. At x^, the Kuhn–Tucker (K–T) condition of (2) is(3)Gx^+g+ATv=0Ax^=bIfG is a positive definite matrix, then the unique solution of (3) is(4)x^=-G-1g+G-1ATAG-1ATAG-1g+bIn the sequential quadratic programming model, an optimal subproblem is set up to search for the next feasible point in the current iterative point.dk is the solution of the optimal subproblem, that is(5)xk+1=xk+αkdkwhere αk∈(0,1]. The process for more iterations is repeated, and the optimization solution can then be obtained.The optimal subproblem of (1) can be defined as follows:(6)minx∇fTxkx-xk+12x-xkTGxkx-xks.t.sxk+∇sTxkx-xk=0where s(x)=Ax-b, ∇f(xk)=∂f/∂x1,∂f/∂x2,⋯,∂f/∂xnTx=xk, ∇s(xk)=∂s/∂x1,∂s/∂x2,⋯,∂s/∂xnTx=xkG ( x k ) is the Hessian matrix of f(x):(7)Gxk=∇2fxk=∂2f∂x2=∂2fx∂x12∂2fx∂x1∂x2⋯∂2fx∂x1∂xt∂2fx∂x1∂x2∂2fx∂x22⋯∂2fx∂x2∂xt⋯⋯⋯⋯∂2fx∂x1∂xt∂2fx∂x2∂xt⋯∂2fx∂xt2x=xkThe penalty function of model (1) is defined as follows:(8)P~x,π=12xTGx+gTx+πAx-bwhere π>0 is the penalty coefficient.If the objective and constraint function satisfy certain conditions, then the penalty method is convergent; i.e., the algorithm converges to the unique, optimal result [15, 29, 30].The penalty function method is simple to compute. As the penalty coefficient increases, the initial point is iterated to the optimal solution. However, when the penalty coefficient is very large, computation time increases.


## 3. Sequential Quadratic Programming for Nonlinear Least Squares

Synthesizing both methods, the penalty function is simple to compute in earlier iterations, and the sequential quadratic programming method allows for quadratic convergence. A direction-controlled nonlinear least squares estimation algorithm using the penalty function and sequential quadratic programming is proposed herein. The least squares model is transformed into the quadratic optimization model. The iterative point enters the feasible region by the penalty function, and the optimal solution can then be obtained by sequential quadratic programming. Efficiency in problem solving can be improved by combining the two methods.The nonlinear model is defined as follows:(9)L+Δ=fxwhere L=(L1,L2,⋯,Ln)T is an n×1 vector; f is a nonlinear function;X = ( X 1 , X 2 , ⋯ , X t ) T is a t×1 vector to be estimated, and Δ=(Δ1,Δ2,⋯,Δn)T is an n×1 observational error vector. The error equation of (9) is(10)Vx=fx-LThe least squares principle must be satisfied, that is,(11)min⁡Fx=VTxPVx=fx-LTPfx-LEquation (11) can be converted to quadratic optimization with equality constraints as follows:(12)minFΔxk=VTPVs.t.V=fxk-LThe subproblem of (12) can be obtained by linearizing the equality constraint at xk. The expression is as follows:(13)minFΔxk=VTPVs.t.V=BxkΔxk+Cxk-Lwhere B(xk)=∂f(x)/∂xx=xk, C(xk)=f(xk), and the positive definite matrix P is the weight matrix. Equation (12) shows convex programming.The penalty function of (12) is(14)P~x,π=VTPV+πfx-LThe steps of the algorithm using the penalty function and sequential quadratic programming are as follows.( 1 ) The initial values are defined as x(0) and π0; the control error is e>0,ε>0; and the number of iterations is k(k=0).( 2 ) The initial value x(0) is substituted into (14) to calculate the iterative value f(x(k))-L. If f(x(k))-L<e, the iteration terminates.( 3 ) If x(k+1)-x(k)≤ε, the result is x∗=x(k+1); otherwise, set x^=x(k) and move to the next step.( 4 ) By taking x(0)=x^ as the new initial value, (12) is solved using sequential quadratic programming. The correction x(k) value can then be obtained.( 5 ) If x(k+1)-x(k)≤ε, the result is x∗=x(k+1); otherwise, return to step (4) until the condition is met.In the course of running the algorithm, the penalty coefficient of (6) is as follows [31]:(15)π=max⁡VTΔxkTP000VΔxkΔxk,2π-+βwhere β≥0 is the constant, and π- is the penalty coefficient of the previous iteration.


## 4. Experimental Results and Analysis

### 4.1. Simulated Data

In this case,AX=L, where(16)A~=2.0000-5.00001.00001.0000-9.5000-2.00004.00001.0000-1.05008.5000-2.00001.00001.0000-1.00002.4000-1.00002.50004.0000-0.50007.0000-1.00003.20004.0000-0.50008.40001.00001.0000-3.00000.40000.49003.00007.0000-3.00001.500012.70005.0000-1.0000-2.00002.5000-3.00004.00002.0000-2.00002.01003.00004.00003.0000-2.00002.00005.0000.The observations areL~=[-10.510.451.41214.1-0.1121.21.59.0112]T.X = [ x 1   x 2   x 3   x 4   x 5 ] Tare unknown parameters. The true values are X^=[11111]T.The white noise isΔ~N(0,σ2I), where I denotes the identity matrix and σ=0.1, e=10-2ε=10-2.The condition number ofA is 41397, so this is an ill-conditioned problem. Using sequential quadratic programming, least squares, and the ridge estimate, the estimated results can be compared, as shown in Table 1.Table 1
Estimated values by sequential quadratic programming, least squares, and the ridge estimate.
True value Least squaresmethod Ridge method Sequential quadratic programming method x ^ 1 2.8552 0.8654 1.2205 1 -1.3035 1.3872 1.2598 1 1.2584 0.8283 1.1896 1 2.3954 0.5180 0.9064 1 -1.9306 1.1375 1.1034  RMSE  1.9672 0.2996 0.1852  The standard deviation  4.3989 0.6699 0.4141In Table1, the ridge parameter of the ridge estimation is 1.4568. Figure 1 shows the ridge parameter from the L-curve [32].Figure 1
Ridge parameter from the L-curve.In this ill-conditioned problem, all methods obtain estimated values. The distance between the initial and true values is large. Sequential quadratic programming obtains the closest solution (Table1). The RMSE and standard deviation show that the sequential quadratic programming model is the best, and the ridge estimate method is better than the least squares method. In sequential quadratic programming, the iterative point enters the feasible region via the penalty function and the optimal solution can then be obtained via sequential quadratic programming. The constraints on the initial values of sequential quadratic programming are reduced; however, the computational cost is increased. Ridge estimation performs much better than least squares, but determining the ridge parameters is difficult.


### 4.2. Phase Unwrapping Experiments

#### 4.2.1. Simulated Data

To verify the effectiveness of the algorithm, simulated data were used to conduct experiments. In a simulated interferogram with a size of 100 pixels × 100 pixels, a coherence coefficient of 0.85 and Gaussian noise with a standard deviation of 0.8160 radians were added to the image,e=ε=10-3. The least squares and sequential quadratic programming methods were used to unwrap images. Figure 2 shows the unwrapping result of each method.Figure 2
Unwrapping results obtained using the above-mentioned two methods.
(a)
Noise-free wrapped phase (b)
Wrapped phase after noise (c)
Least squares method (d)
Sequential quadratic programming methodBoth methods obtain good unwrapping results, without lines or islands. The image obtained by sequential quadratic programming is smoother than that obtained by least squares and is more similar to the noise-free image.Figure3 shows a cross-sectional view of the unwrapping phase.Figure 3
Fiftieth line cross-sectional view of the unwrapping phase using least squares and sequential quadratic programming.
(a)
Noise-free (b)
Least squares (c)
Sequential quadratic programmingAs shown in Figure3, the abscissa represents the horizontal position (x-coordinate) of the pixel in the 50th line of the image. The ordinate denotes the phase. The maximum value of the cross-sectional view of sequential quadratic programming is close to the maximum value of the simulated surface. The maximum value of the cross-sectional view of the least squares method is far from the maximum value of the simulated surface. The curve obtained for of the least squares method is also much rougher than that of sequential quadratic programming.Table2 shows the running time, RMSE, and ε~ value of sequential quadratic programming and least squares obtained by phase unwrapping. Both methods result in relatively reliable unwrapping. The running time is similar, indicating that the efficiency of the two algorithms is almost the same. The RMSE of sequential quadratic programming is smaller than that of least squares, indicating that sequential quadratic programming is more precise. ε~ indicates the quality of the wrap; the smaller the ε~, the better the unwrapping quality. Thus, the sequential quadratic programming method provides better unwrapping than the least squares method.Table 2
Unwrapping time, RMSE, andε~ value.
Method Computation time (s) RMSE (rad) ε ~ (rad) Least squares 12.142 0.7594 0.6214  Sequential quadratic programming 11.866 0.5393 0.4228The experiment with simulated data shows that sequential quadratic programming is a feasible method: the model is effective and has high precision.


#### 4.2.2. Observational Data

Satellite data collected by ENVISAT ASAR, covering the Jining area of Shandong province, China, were used to test the model’s performance.Table3 shows the basic parameters of the data.Table 3
Selected parameters of ASAR image data.
Number Sensor Date Type Latitude and longitudeof center Polarization Ascending Wavelength 022403 ENVISAT 20050306 SLC 116.314 6, 35.305 1 VV Ascending C  022404 ENVISAT 20050130 SLC 116.314 4, 35.305 1 VV Ascending CThe images are 1024 × 1024 pixels, with a resolution of 30 m. Figure4 shows the unwrapping result of each method, e=ε=10-3.Figure 4
Unwrapping of ASAR images, using least squares and sequential quadratic programming.
(a)
Wrapped phase (b)
Least squares method (c)
Sequential quadratic programming methodBoth unwrapping methods were run successfully. In the interferogram, the stripe noise is large and the fringes are unclear. The phase ripple of sequential quadratic programming is smaller than that of least squares, demonstrating that sequential quadratic programming experiences less interference.Table4 shows the computation time, RMSE, and ε~ value.Table 4
Unwrapping time, RMSE, andε values.
Method computation time (s) RMSE (rad) ε ~ (rad) Least squares 725.940 1.7594 0.6214  Sequential quadratic programming 1336.948 0.5393 0.4228The RMSE andε~ values of sequential quadratic programming are smaller than those of the least squares method. This illustrates that the difference between the rewrapped results and the original wrapped phase is small. The unwrapping phases of the sequential quadratic programming method are smoother and more reliable. However, the computation time of sequential quadratic programming is greater than that of the least squares method. This is due to high levels of noise and the penalty function increasing the computation time of sequential quadratic programming.


## 4.1. Simulated Data

In this case,AX=L, where(16)A~=2.0000-5.00001.00001.0000-9.5000-2.00004.00001.0000-1.05008.5000-2.00001.00001.0000-1.00002.4000-1.00002.50004.0000-0.50007.0000-1.00003.20004.0000-0.50008.40001.00001.0000-3.00000.40000.49003.00007.0000-3.00001.500012.70005.0000-1.0000-2.00002.5000-3.00004.00002.0000-2.00002.01003.00004.00003.0000-2.00002.00005.0000.The observations areL~=[-10.510.451.41214.1-0.1121.21.59.0112]T.X = [ x 1   x 2   x 3   x 4   x 5 ] Tare unknown parameters. The true values are X^=[11111]T.The white noise isΔ~N(0,σ2I), where I denotes the identity matrix and σ=0.1, e=10-2ε=10-2.The condition number ofA is 41397, so this is an ill-conditioned problem. Using sequential quadratic programming, least squares, and the ridge estimate, the estimated results can be compared, as shown in Table 1.Table 1
Estimated values by sequential quadratic programming, least squares, and the ridge estimate.
True value Least squaresmethod Ridge method Sequential quadratic programming method x ^ 1 2.8552 0.8654 1.2205 1 -1.3035 1.3872 1.2598 1 1.2584 0.8283 1.1896 1 2.3954 0.5180 0.9064 1 -1.9306 1.1375 1.1034  RMSE  1.9672 0.2996 0.1852  The standard deviation  4.3989 0.6699 0.4141In Table1, the ridge parameter of the ridge estimation is 1.4568. Figure 1 shows the ridge parameter from the L-curve [32].Figure 1
Ridge parameter from the L-curve.In this ill-conditioned problem, all methods obtain estimated values. The distance between the initial and true values is large. Sequential quadratic programming obtains the closest solution (Table1). The RMSE and standard deviation show that the sequential quadratic programming model is the best, and the ridge estimate method is better than the least squares method. In sequential quadratic programming, the iterative point enters the feasible region via the penalty function and the optimal solution can then be obtained via sequential quadratic programming. The constraints on the initial values of sequential quadratic programming are reduced; however, the computational cost is increased. Ridge estimation performs much better than least squares, but determining the ridge parameters is difficult.


## 4.2. Phase Unwrapping Experiments

### 4.2.1. Simulated Data

To verify the effectiveness of the algorithm, simulated data were used to conduct experiments. In a simulated interferogram with a size of 100 pixels × 100 pixels, a coherence coefficient of 0.85 and Gaussian noise with a standard deviation of 0.8160 radians were added to the image,e=ε=10-3. The least squares and sequential quadratic programming methods were used to unwrap images. Figure 2 shows the unwrapping result of each method.Figure 2
Unwrapping results obtained using the above-mentioned two methods.
(a)
Noise-free wrapped phase (b)
Wrapped phase after noise (c)
Least squares method (d)
Sequential quadratic programming methodBoth methods obtain good unwrapping results, without lines or islands. The image obtained by sequential quadratic programming is smoother than that obtained by least squares and is more similar to the noise-free image.Figure3 shows a cross-sectional view of the unwrapping phase.Figure 3
Fiftieth line cross-sectional view of the unwrapping phase using least squares and sequential quadratic programming.
(a)
Noise-free (b)
Least squares (c)
Sequential quadratic programmingAs shown in Figure3, the abscissa represents the horizontal position (x-coordinate) of the pixel in the 50th line of the image. The ordinate denotes the phase. The maximum value of the cross-sectional view of sequential quadratic programming is close to the maximum value of the simulated surface. The maximum value of the cross-sectional view of the least squares method is far from the maximum value of the simulated surface. The curve obtained for of the least squares method is also much rougher than that of sequential quadratic programming.Table2 shows the running time, RMSE, and ε~ value of sequential quadratic programming and least squares obtained by phase unwrapping. Both methods result in relatively reliable unwrapping. The running time is similar, indicating that the efficiency of the two algorithms is almost the same. The RMSE of sequential quadratic programming is smaller than that of least squares, indicating that sequential quadratic programming is more precise. ε~ indicates the quality of the wrap; the smaller the ε~, the better the unwrapping quality. Thus, the sequential quadratic programming method provides better unwrapping than the least squares method.Table 2
Unwrapping time, RMSE, andε~ value.
Method Computation time (s) RMSE (rad) ε ~ (rad) Least squares 12.142 0.7594 0.6214  Sequential quadratic programming 11.866 0.5393 0.4228The experiment with simulated data shows that sequential quadratic programming is a feasible method: the model is effective and has high precision.


### 4.2.2. Observational Data

Satellite data collected by ENVISAT ASAR, covering the Jining area of Shandong province, China, were used to test the model’s performance.Table3 shows the basic parameters of the data.Table 3
Selected parameters of ASAR image data.
Number Sensor Date Type Latitude and longitudeof center Polarization Ascending Wavelength 022403 ENVISAT 20050306 SLC 116.314 6, 35.305 1 VV Ascending C  022404 ENVISAT 20050130 SLC 116.314 4, 35.305 1 VV Ascending CThe images are 1024 × 1024 pixels, with a resolution of 30 m. Figure4 shows the unwrapping result of each method, e=ε=10-3.Figure 4
Unwrapping of ASAR images, using least squares and sequential quadratic programming.
(a)
Wrapped phase (b)
Least squares method (c)
Sequential quadratic programming methodBoth unwrapping methods were run successfully. In the interferogram, the stripe noise is large and the fringes are unclear. The phase ripple of sequential quadratic programming is smaller than that of least squares, demonstrating that sequential quadratic programming experiences less interference.Table4 shows the computation time, RMSE, and ε~ value.Table 4
Unwrapping time, RMSE, andε values.
Method computation time (s) RMSE (rad) ε ~ (rad) Least squares 725.940 1.7594 0.6214  Sequential quadratic programming 1336.948 0.5393 0.4228The RMSE andε~ values of sequential quadratic programming are smaller than those of the least squares method. This illustrates that the difference between the rewrapped results and the original wrapped phase is small. The unwrapping phases of the sequential quadratic programming method are smoother and more reliable. However, the computation time of sequential quadratic programming is greater than that of the least squares method. This is due to high levels of noise and the penalty function increasing the computation time of sequential quadratic programming.


## 4.2.1. Simulated Data

To verify the effectiveness of the algorithm, simulated data were used to conduct experiments. In a simulated interferogram with a size of 100 pixels × 100 pixels, a coherence coefficient of 0.85 and Gaussian noise with a standard deviation of 0.8160 radians were added to the image,e=ε=10-3. The least squares and sequential quadratic programming methods were used to unwrap images. Figure 2 shows the unwrapping result of each method.Figure 2
Unwrapping results obtained using the above-mentioned two methods.
(a)
Noise-free wrapped phase (b)
Wrapped phase after noise (c)
Least squares method (d)
Sequential quadratic programming methodBoth methods obtain good unwrapping results, without lines or islands. The image obtained by sequential quadratic programming is smoother than that obtained by least squares and is more similar to the noise-free image.Figure3 shows a cross-sectional view of the unwrapping phase.Figure 3
Fiftieth line cross-sectional view of the unwrapping phase using least squares and sequential quadratic programming.
(a)
Noise-free (b)
Least squares (c)
Sequential quadratic programmingAs shown in Figure3, the abscissa represents the horizontal position (x-coordinate) of the pixel in the 50th line of the image. The ordinate denotes the phase. The maximum value of the cross-sectional view of sequential quadratic programming is close to the maximum value of the simulated surface. The maximum value of the cross-sectional view of the least squares method is far from the maximum value of the simulated surface. The curve obtained for of the least squares method is also much rougher than that of sequential quadratic programming.Table2 shows the running time, RMSE, and ε~ value of sequential quadratic programming and least squares obtained by phase unwrapping. Both methods result in relatively reliable unwrapping. The running time is similar, indicating that the efficiency of the two algorithms is almost the same. The RMSE of sequential quadratic programming is smaller than that of least squares, indicating that sequential quadratic programming is more precise. ε~ indicates the quality of the wrap; the smaller the ε~, the better the unwrapping quality. Thus, the sequential quadratic programming method provides better unwrapping than the least squares method.Table 2
Unwrapping time, RMSE, andε~ value.
Method Computation time (s) RMSE (rad) ε ~ (rad) Least squares 12.142 0.7594 0.6214  Sequential quadratic programming 11.866 0.5393 0.4228The experiment with simulated data shows that sequential quadratic programming is a feasible method: the model is effective and has high precision.


## 4.2.2. Observational Data

Satellite data collected by ENVISAT ASAR, covering the Jining area of Shandong province, China, were used to test the model’s performance.Table3 shows the basic parameters of the data.Table 3
Selected parameters of ASAR image data.
Number Sensor Date Type Latitude and longitudeof center Polarization Ascending Wavelength 022403 ENVISAT 20050306 SLC 116.314 6, 35.305 1 VV Ascending C  022404 ENVISAT 20050130 SLC 116.314 4, 35.305 1 VV Ascending CThe images are 1024 × 1024 pixels, with a resolution of 30 m. Figure4 shows the unwrapping result of each method, e=ε=10-3.Figure 4
Unwrapping of ASAR images, using least squares and sequential quadratic programming.
(a)
Wrapped phase (b)
Least squares method (c)
Sequential quadratic programming methodBoth unwrapping methods were run successfully. In the interferogram, the stripe noise is large and the fringes are unclear. The phase ripple of sequential quadratic programming is smaller than that of least squares, demonstrating that sequential quadratic programming experiences less interference.Table4 shows the computation time, RMSE, and ε~ value.Table 4
Unwrapping time, RMSE, andε values.
Method computation time (s) RMSE (rad) ε ~ (rad) Least squares 725.940 1.7594 0.6214  Sequential quadratic programming 1336.948 0.5393 0.4228The RMSE andε~ values of sequential quadratic programming are smaller than those of the least squares method. This illustrates that the difference between the rewrapped results and the original wrapped phase is small. The unwrapping phases of the sequential quadratic programming method are smoother and more reliable. However, the computation time of sequential quadratic programming is greater than that of the least squares method. This is due to high levels of noise and the penalty function increasing the computation time of sequential quadratic programming.


## 5. Conclusions

A direction-controlled nonlinear least squares estimation algorithm using the sequential quadratic programming method is proposed herein. The least squares model is transformed into the nonlinear programming model, and the iteration direction can then be controlled by sequential quadratic programming. The ill-conditioned matrix and phase unwrapping are used to verify the method. The proposed algorithm achieves better results than alternative approaches in terms of RMSE andε~ values, and its effectiveness is demonstrated using observational satellite data. In future research, the role change if the inequality constrains are taken into account will be considered.



---

*Source: 3087949-2019-06-18.xml*