# Tikhonov Regularized Variable Projection Algorithms for Separable Nonlinear Least Squares Problems

Authors: Zhengqing Fu; Lanlan Guo
Journal: Complexity
 (2019)

---

## Abstract

This paper considers the classical separable nonlinear least squares problem. Such problems can be expressed as a linear combination of nonlinear functions, and both linear and nonlinear parameters are to be estimated. Among the existing results, ill-conditioned problems are less often considered. Hence, this paper focuses on an algorithm for ill-conditioned problems. In the proposed linear parameter estimation process, the sensitivity of the model to disturbance is reduced using Tikhonov regularisation. The Levenberg–Marquardt algorithm is used to estimate the nonlinear parameters. The Jacobian matrix required by LM is calculated by the Golub and Pereyra, Kaufman, and Ruano methods. Combining the nonlinear and linear parameter estimation methods, three estimation models are obtained and the feasibility and stability of the model estimation are demonstrated. The model is validated by simulation data and real data. The experimental results also illustrate the feasibility and stability of the model.


---

## Body

## 1. Introduction

The separable nonlinear least squares problem is a special type of nonlinear least squares problem. It was proposed by Golub and Pereyra for estimating the parameters of the formula for atomic physical ion half-life [1]. In the separable nonlinear least squares model, the estimated model can be expressed as a linear combination of nonlinear problems.To utilise the special structure of the separable nonlinear least squares problem, in 1973, Golub and Pereyra proposed the variable projection (VP) algorithm [1], which eliminates linear parameters and simplifies the problem into a set of nonlinear parameter estimation problems. This method also reduces the computational complexity. Compared with the corresponding algorithm for a direct solution, the VP algorithm has fewer iteration steps and less initial point guessing. Moreover, when the original problem is ill posed, the degree of ill-posedness can be reduced [2–5]. As research has progressed, the mathematical models of many problems in actual engineering can be expressed as separable nonlinear least squares, for instance, in inverse problems and problems in signal processing, medical and biological imaging, neural networks, communication, electrical and electronic engineering, and differential equation dynamic systems [6–15].In 1980, Ruhe and Wedin analysed the separation and nonseparation of parameters and concluded that the parameter separation approach was simpler and more effective [16]. In 1990, Shen and Ypma used the variable projection method to transform a separable nonlinear least squares problem into a case containing only nonlinear parameters, which improved the computational efficiency of the function [17]. In 2003, Golub and Pereyra summarised the development and application of separable nonlinear least squares over the past 30 years [18]. In recent years, many studies and application of the separable nonlinear least squares problem have obtained various results. Chen studied the nonlinear four-parameter sine wave model and calculated it using the variable projection method. Experiments show that the convergence speed has a strong relationship with the frequency parameter [19]. Chung and Nagy applied separable nonlinear least squares to large-scale ill-posed problems and improved the distortion problem in image processing [20].There are a large number of research results for separable nonlinear least squares, but there are relatively few studies on potentially ill-conditioned problems with parameters. When a singular value becomes very small, the least squares model may be ill-conditioned. Regularisation is a commonly used method to solve ill-conditioned problems that causes regression coefficients to have smaller variance values, thus solving potentially ill-posed problems [21–28]. The Tikhonov regularisation (TR) method [21–23], truncated singular value method [24, 25], kernel function-based regularisation method [26, 27], and l1 norm regularisation method [28] are often used to solve ill-posed problems. When estimating nonlinear parameters, iterative search methods such as the Gauss–Newton method, steepest gradient method, and LM method are commonly used [29–34]. In these methods, the Jacobian matrix has a strong influence on the computational efficiency of the algorithm.This paper proposes a hybrid VP-based algorithm that combines the LM and TR methods and evaluates the effect of different Jacobian matrix algorithms on its accuracy and efficiency. The TR method is used to regularise the linear parameter estimation, and the LM algorithm is used to estimate the nonlinear parameters in a separable nonlinear least squares problem. The Jacobian matrix is calculated using three methods: Golub and Pereyra (GP), Kaufman (KAU), and Ruano (RJF). The algorithm proposed in this paper combines the advantages of the LM algorithm with those of a regularisation method and can improve the estimation efficiency. The model was validated using two examples: exponential fitting and the determination of waveform parameters for airborne radar sounding data.


## 2. VP Model and Its Parameter Estimation Method

### 2.1. VP Model

The separable nonlinear least squares model can be expressed in the following form:(1)yt=∑j=1mθLϕjθN;tt,where ϕjθN;ttj=1,2,…,m,t=1,2,…,n are nonlinear functions; yt, t=1,2,…,n are observation data; tt is the relevant variable of yt; and θL=θL,1,θL,2,…,θL,m and θN=θN,1,θN,2,…,θN,k are the linear and nonlinear parameters, respectively, to be estimated.The optimal parametersθL and θN can be obtained by the following nonlinear function:(2)θL,θN=argminθL,θN∑t=1nyt−∑jmθLϕjθN;tt2.The above formula can be expressed using matrices as follows:(3)θL,θN=argminθL,θNy−ΦθNθL22,where the column vector of matrix ΦθN is a nonlinear function ϕjθN;tt, the element of vector y is yt, and ⋅22 is the Euclidean norm. For given nonlinear parameters θN, the linear parameters θL can be obtained by solving the following nonlinear least squares problem:(4)θL=ΦθNTΦθN−1ΦθNTy=ΦθN+y,where ΦθN+ is the pseudoinverse of ΦθN. If ΦθNTΦθN−1ΦθNT is computable, then (4) can be solved. Hence, the premise of (4) is that matrix ΦθNTΦθN−1 is nonsingular, but this premise is not true in some cases.


### 2.2. Determination of the Linear Parameters

Without loss of generality, the power system can be expressed as(5)y=ΦθNθL=ΦθL,where Φ∈Ra×b,a≥b, θL∈Rb, and y∈Ra.We can then perform singular value decomposition (SVD) onΦ, which can be expressed as(6)Φ=USVT=∑i=1luiσiviT,where U=u1,u2,…,ua∈Ra×a and V=v1,v2,…,vb∈Rb×b are the unitary matrices composed of column vectors. S=sij∈Ra×b, sij=σi,i=j0,i≠j, σ1≥σ2≥⋯≥σl>0,σl+1=σl+2=⋯=σb=0. In addition, σii=1l is the singular value of Φ. Using SVD to solve (5), we obtain(7)θL=∑i=1luiTyσivi.To reduce the sensitivity of the SVD method to disturbances, a filter factorfi can be added. The SVD form with the added filter factor is(8)θL=∑i=1lfiuiTyσivi.Different filtering factors will result in different regularisation methods, and the TR method is a widely applied method [35, 36]. The filtering factor obtained by the TR method can be expressed as fi=σi2/σi2+σ12, and (8) can be expressed as(9)θL=∑i=1lσi2σi2+σ12uiTyσivi.


### 2.3. Determination of the Nonlinear Parameters

For matrixFθNk,θLk=y−ΦθNθL, the Jacobian matrix JθNk is the first derivative of FθNk,θLk:(10)JθNk=∂FθNk,θLk∂θNkT.Assume(11)PΦθNk⊥=I−ΦθNkΦθNk+.The matrix can be expressed as(12)FθNk,θLk=I−ΦθNθLy.The Jacobian matrix can be calculated in the following three ways:(1)
The calculation method given by Golub and Pereyra is as follows:(13)JGP=DPΦ⊥y=−PΦ⊥DΦΦ−y−PΦ⊥DΦΦ−Ty,whereD is the Fréchet derivative of the map. In addition, Φ− is the symmetric generalised inverse of Φ that satisfies ΦΦ−Φ=Φ and ΦΦT=ΦΦ−(2)
Kaufman proposed the following simplified calculation method for the Jacobian matrix after studying (13):(14)JKAU=−PΦ⊥DΦΦ−y.(3)
Ruano proposed an even simpler calculation method for the Jacobian matrix based on Kaufman’s method:(15)JRJF=−DΦΦ−y.Ruano et al. [37] proved that (15) is valid and can be obtained from Kaufman’s Jacobian matrix (14). It can be proved that the same gradient vector can be obtained by the above three Jacobian matrices [38].The following equation can be used to iterate nonlinear parameters:(16)θNk+1=θNk+βkdk,where βk is a scalar step size that ensures that the objective function y−ΦθNθ˜L22 is decreasing and dk is the search direction. In the LM algorithm, βk can be determined by the following equation [39]:(17)rθk+ρmkdk≤rθk+ρmkgkTdk,where mk is the smallest nonnegative integer satisfying (17), ρ=0.5, gk=Jθkrθk, and βk=ρmk.d
k can be determined by the following equation:(18)JθNkTJθNk+γkIdk=−JθNkTFθNk,θLk,where γk is the damping factor and affects both the gradient and dk. If γk tends to infinity, dk tends to the fastest direction. If γk tends to zero, dk tends to the Gauss–Newton direction. There are many methods of selecting damping factor γk [40]. After obtaining γk, the search direction dk can be obtained by solving (17). Step parameter βk can be obtained by a linear search algorithm with mixed quadratic interpolation and cubic interpolation [41–43].Combining the above methods for the linear and nonlinear parameters, three algorithms can be obtained. The LM algorithm is used for the nonlinear parameters, and the three methods of GP, KAU, and RJF are used to calculate the Jacobian matrix in this algorithm. For the linear parameters, the TR method is used. The formula is as follows:(19)θL,θN=argminθNy−ΦθNθL22,θL=∑i=1lσi2σi2+σ12uiTyσivi.Moreover, the specific algorithm is as follows:Step 1: give the initial valueθN0 of the nonlinear parameter, and predefine the maximum number of iteration steps kmax and error εStep 2: calculate the initial valueθL0 of the linear parameter using the TR method of (5) and (9)Step 3: alternately update the nonlinear parameters and linear parametersStep 4: repeat the third step until either the maximum number of iteration stepskmax is reached or θNk+1−θNk<ε


## 2.1. VP Model

The separable nonlinear least squares model can be expressed in the following form:(1)yt=∑j=1mθLϕjθN;tt,where ϕjθN;ttj=1,2,…,m,t=1,2,…,n are nonlinear functions; yt, t=1,2,…,n are observation data; tt is the relevant variable of yt; and θL=θL,1,θL,2,…,θL,m and θN=θN,1,θN,2,…,θN,k are the linear and nonlinear parameters, respectively, to be estimated.The optimal parametersθL and θN can be obtained by the following nonlinear function:(2)θL,θN=argminθL,θN∑t=1nyt−∑jmθLϕjθN;tt2.The above formula can be expressed using matrices as follows:(3)θL,θN=argminθL,θNy−ΦθNθL22,where the column vector of matrix ΦθN is a nonlinear function ϕjθN;tt, the element of vector y is yt, and ⋅22 is the Euclidean norm. For given nonlinear parameters θN, the linear parameters θL can be obtained by solving the following nonlinear least squares problem:(4)θL=ΦθNTΦθN−1ΦθNTy=ΦθN+y,where ΦθN+ is the pseudoinverse of ΦθN. If ΦθNTΦθN−1ΦθNT is computable, then (4) can be solved. Hence, the premise of (4) is that matrix ΦθNTΦθN−1 is nonsingular, but this premise is not true in some cases.


## 2.2. Determination of the Linear Parameters

Without loss of generality, the power system can be expressed as(5)y=ΦθNθL=ΦθL,where Φ∈Ra×b,a≥b, θL∈Rb, and y∈Ra.We can then perform singular value decomposition (SVD) onΦ, which can be expressed as(6)Φ=USVT=∑i=1luiσiviT,where U=u1,u2,…,ua∈Ra×a and V=v1,v2,…,vb∈Rb×b are the unitary matrices composed of column vectors. S=sij∈Ra×b, sij=σi,i=j0,i≠j, σ1≥σ2≥⋯≥σl>0,σl+1=σl+2=⋯=σb=0. In addition, σii=1l is the singular value of Φ. Using SVD to solve (5), we obtain(7)θL=∑i=1luiTyσivi.To reduce the sensitivity of the SVD method to disturbances, a filter factorfi can be added. The SVD form with the added filter factor is(8)θL=∑i=1lfiuiTyσivi.Different filtering factors will result in different regularisation methods, and the TR method is a widely applied method [35, 36]. The filtering factor obtained by the TR method can be expressed as fi=σi2/σi2+σ12, and (8) can be expressed as(9)θL=∑i=1lσi2σi2+σ12uiTyσivi.


## 2.3. Determination of the Nonlinear Parameters

For matrixFθNk,θLk=y−ΦθNθL, the Jacobian matrix JθNk is the first derivative of FθNk,θLk:(10)JθNk=∂FθNk,θLk∂θNkT.Assume(11)PΦθNk⊥=I−ΦθNkΦθNk+.The matrix can be expressed as(12)FθNk,θLk=I−ΦθNθLy.The Jacobian matrix can be calculated in the following three ways:(1)
The calculation method given by Golub and Pereyra is as follows:(13)JGP=DPΦ⊥y=−PΦ⊥DΦΦ−y−PΦ⊥DΦΦ−Ty,whereD is the Fréchet derivative of the map. In addition, Φ− is the symmetric generalised inverse of Φ that satisfies ΦΦ−Φ=Φ and ΦΦT=ΦΦ−(2)
Kaufman proposed the following simplified calculation method for the Jacobian matrix after studying (13):(14)JKAU=−PΦ⊥DΦΦ−y.(3)
Ruano proposed an even simpler calculation method for the Jacobian matrix based on Kaufman’s method:(15)JRJF=−DΦΦ−y.Ruano et al. [37] proved that (15) is valid and can be obtained from Kaufman’s Jacobian matrix (14). It can be proved that the same gradient vector can be obtained by the above three Jacobian matrices [38].The following equation can be used to iterate nonlinear parameters:(16)θNk+1=θNk+βkdk,where βk is a scalar step size that ensures that the objective function y−ΦθNθ˜L22 is decreasing and dk is the search direction. In the LM algorithm, βk can be determined by the following equation [39]:(17)rθk+ρmkdk≤rθk+ρmkgkTdk,where mk is the smallest nonnegative integer satisfying (17), ρ=0.5, gk=Jθkrθk, and βk=ρmk.d
k can be determined by the following equation:(18)JθNkTJθNk+γkIdk=−JθNkTFθNk,θLk,where γk is the damping factor and affects both the gradient and dk. If γk tends to infinity, dk tends to the fastest direction. If γk tends to zero, dk tends to the Gauss–Newton direction. There are many methods of selecting damping factor γk [40]. After obtaining γk, the search direction dk can be obtained by solving (17). Step parameter βk can be obtained by a linear search algorithm with mixed quadratic interpolation and cubic interpolation [41–43].Combining the above methods for the linear and nonlinear parameters, three algorithms can be obtained. The LM algorithm is used for the nonlinear parameters, and the three methods of GP, KAU, and RJF are used to calculate the Jacobian matrix in this algorithm. For the linear parameters, the TR method is used. The formula is as follows:(19)θL,θN=argminθNy−ΦθNθL22,θL=∑i=1lσi2σi2+σ12uiTyσivi.Moreover, the specific algorithm is as follows:Step 1: give the initial valueθN0 of the nonlinear parameter, and predefine the maximum number of iteration steps kmax and error εStep 2: calculate the initial valueθL0 of the linear parameter using the TR method of (5) and (9)Step 3: alternately update the nonlinear parameters and linear parametersStep 4: repeat the third step until either the maximum number of iteration stepskmax is reached or θNk+1−θNk<ε


## 3. Numerical Experiments

The experimental environment used in this study was MATLAB 2016b, running on a PC equipped with a 2.30 GHz CPU, 4 GB of memory, and a Windows 10 operating system.


### 3.1. Index-Fitting Model

Two experiments were used to evaluate the three algorithms proposed in this paper. The first of these experiments was derived from the observation of the decay of radioactive materials. The model describes the sum of the exponential terms of two unknown attenuation factors and is written as follows:(20)yβ,λ;t≈β1e−λ1t+β2e−λ2t,where β=β1,β2T and λ=λ1,λ2T are the linear and nonlinear unknown parameters to be estimated, respectively. Moreover, there are 21 ti and corresponding observations yi, i=1,2,…,21. We minimise the following functions according to the principle of least squares:(21)12∑i=121yi−yβ,λ;t2=12∑i=121yi−β1e−λ1t+β2e−λ2t2=12∑i=121yi−∑j=12βjϕjλ,ti2.According to the composition of the function matrix in the second part, the nonlinear function part can be expressed as a matrix form asΦλ=∑j=12φjλ,ti. After SVD, for Φλ, the estimated values of the linear parameters θLk are obtained from (9). For the nonlinear parameters, the objective function FθNk,θLk is first constructed by (12), the Jacobian matrix is then obtained by the three different methods, and the LM algorithm is substituted into (16) to determine the parameters. When iterating, we set the initial value to λ0=7,3T, the maximum number of iteration steps to 100, and ε=10−6. Performing SVD on Φ yields singular values of 1.8480 and 0.3586.Figure1 shows the observations and their fitting curves.Figure 1
Observation data and fitting curves. (a) VPGP. (b) VPKAU. (c) VPRJF.
(a)
(b)
(c)The curves of the observation value fit the results of the three methods very well. The curves obtained by the three methods are basically the same, demonstrating that all three methods are feasible and the results are reliable.Figure2 shows the change in parameters λ1 and λ2 with respect to the number of iterations.Figure 2
Change in parameters with respect to the number of iterations. (a) VPGP. (b) VPKAU. (c) VPRJF.
(a)
(b)
(c)As Figure2 shows, for λ1, the change in parameter curves for the two methods VPGP+TR and VPRJF+TR are basically the same, and the specified accuracy can be achieved after nine iterations. Visually, it is clear that after five iterations, the solution is very close to the final result. In contrast, the VPKAU+TR method reached the maximum number of iterations (100). However, after 30 iterations, even though the result changes a little, it is close to the final result.Forλ2, the change in parameter curves for the two methods VPGP+TR and VPRJF+TR are also basically the same, and the specified accuracy is also achieved after nine iterations. The graph shows that, after four iterations, the solution is very close to the final calculation result. The VPKAU+TR method again reached the maximum number of calculations, but the solution is close to the final result after 20 iterations.Table1 lists the results of the parameters obtained by the three methods.Table 1
Parameter results.
Method
λ
1
λ
2
β
1
β
2
VPGP+TR
3.0067
2.8893
10.5868
1.4004
VPKAU+TR
3.0068
2.8891
10.5868
1.4004
VPRJF+TR
3.0067
2.8893
10.5868
1.4004The table shows that, as for the curves in Figure1, the results obtained by the three methods are basically the same.Table2 gives the runtime (in seconds), number of iterations, residual sum, and root mean square error (RMSE) for the three methods.Table 2
Quantitative comparison of performance for each method.
Method
Runtime (s)
Number of iterations
Residual sum
RMSE
VPGP+TR
40.596687
10
0.1477
0.0859
VPKAU+TR
17.186718
10
0.1477
0.0859
VPRJF+TR
31.073001
100
0.1477
0.0859The quantitative indicators (residual sums and RMSE) are equal and very small, indicating that the three algorithms achieve the same calculation accuracy, which is relatively high. The difference between them is only the number of iterations and runtime. The VPGP+TR method takes the longest to calculate because this method uses the most complicated way to calculate the Jacobian matrix. In contrast, VPRJF+TR is a simple method for calculating the Jacobian matrix. However, to meet the accuracy requirements, the number of iterations is also large, so the calculation time is also long. The VPKAU+TR method somewhat simplifies the calculation of the Jacobian matrix, and to some extent, the accuracy of the calculation is guaranteed, so the runtime is the shortest. For this simulation experiment, the VPKAU+TR method performs the best of all methods.


### 3.2. Waveform Parameter Calculation for Airborne Laser Radar Sounding Data

The second experiment used laser radar sounding waveform data. The experimental data were flight experimental data obtained from a flight over the sea near Sanya City, Hainan Province, China.The waveform sample data are stored as an array, and the amplitude of each sample point is recorded separately. The sample information is presented as a ‘waveform’ form on a two-dimensional plane. The onboard laser system records discrete echoes that simulate the “waveform” of the raw sampled data using the most appropriate functional model. To improve the sounding accuracy, it is highly desirable to choose a high-precision theoretical model.This experiment has a total of 320 waveform data, and the function model is as follows:(22)y=a+be−x−μ2/2σ2.Selecting reasonable data can effectively reduce the influence of noise. In this study, we selected some of the data that have no gross errors. We extracted 42 data from the waveform data to determine the parameters. When iterating, we set the initial value toλ0=15,3T, the maximum number of iteration steps to 100, and the error to ε=10−3. SVD for Φ yielded singular values of 6.5941 and 1.9585.Table3 lists the calculation results of the model parameters.Table 3
Parameter results.
Method
a
b
μ
σ
VPGP+TR
207.1076
498.6737
14.0070
3.4534
VPKAU+TR
207.0941
498.6461
14.0068
3.4540
VPRJF+TR
207.0518
498.5602
14.0085
3.4560The results given in the table show that the results obtained by the three methods are basically the same, indicating that all three methods are feasible approaches. The result obtained by VPRJF+TR is slightly different from those of the other two methods because this method simplifies the calculation of the Jacobian matrix the most, resulting in a gap.Figure3 shows the fitting curves of the model.Figure 3
Waveform data and fitting curves. (a) VPGP+TR. (b) VPKAU+TR. (c) VPRJF+TR.
(a)
(b)
(c)In Figure3, all three methods fit the curve well, which also intuitively illustrates the feasibility of these three methods.Figure4 shows the change in parameters μ and σ with respect to the number of iterations.Figure 4
Change in parameters with respect to the number of iterations. (a) VPGP+TR. (b) VPKAU+TR. (c) VPRJF+TR.
(a)
(b)
(c)The results in Figure4 show that for parameters μ and σ, when error ε=10−3, the parameter change is basically stable after five iterations. This shows that the convergence of these three methods is fast and accurate.Table4 lists the runtime (in seconds), number of iterations, and root mean square error (RMSE).Table 4
Quantitative comparison of performance for each method.
Method
Runtime (s)
Number of iterations
RMSE
VPGP+TR
263.797378
10
38.0979
VPKAU+TR
91.929637
10
38.0980
VPRJF+TR
9.817583
10
38.0983Forε=10−3, number of iterations and RMSE (i.e., the convergence speed and accuracy) of these three methods are basically the same. With respect to runtime, the VPGP+TR method again takes the longest because of its calculation of the Jacobian matrix. Hence, the same number of iterations takes longer. Similarly, the Jacobian matrix of the VPRJF+TR method is the simplest and requires the shortest time to calculate. For the VPKAU+TR method, the computational complexity is between those of the other two methods, and it takes moderate amount of time to calculate.


## 3.1. Index-Fitting Model

Two experiments were used to evaluate the three algorithms proposed in this paper. The first of these experiments was derived from the observation of the decay of radioactive materials. The model describes the sum of the exponential terms of two unknown attenuation factors and is written as follows:(20)yβ,λ;t≈β1e−λ1t+β2e−λ2t,where β=β1,β2T and λ=λ1,λ2T are the linear and nonlinear unknown parameters to be estimated, respectively. Moreover, there are 21 ti and corresponding observations yi, i=1,2,…,21. We minimise the following functions according to the principle of least squares:(21)12∑i=121yi−yβ,λ;t2=12∑i=121yi−β1e−λ1t+β2e−λ2t2=12∑i=121yi−∑j=12βjϕjλ,ti2.According to the composition of the function matrix in the second part, the nonlinear function part can be expressed as a matrix form asΦλ=∑j=12φjλ,ti. After SVD, for Φλ, the estimated values of the linear parameters θLk are obtained from (9). For the nonlinear parameters, the objective function FθNk,θLk is first constructed by (12), the Jacobian matrix is then obtained by the three different methods, and the LM algorithm is substituted into (16) to determine the parameters. When iterating, we set the initial value to λ0=7,3T, the maximum number of iteration steps to 100, and ε=10−6. Performing SVD on Φ yields singular values of 1.8480 and 0.3586.Figure1 shows the observations and their fitting curves.Figure 1
Observation data and fitting curves. (a) VPGP. (b) VPKAU. (c) VPRJF.
(a)
(b)
(c)The curves of the observation value fit the results of the three methods very well. The curves obtained by the three methods are basically the same, demonstrating that all three methods are feasible and the results are reliable.Figure2 shows the change in parameters λ1 and λ2 with respect to the number of iterations.Figure 2
Change in parameters with respect to the number of iterations. (a) VPGP. (b) VPKAU. (c) VPRJF.
(a)
(b)
(c)As Figure2 shows, for λ1, the change in parameter curves for the two methods VPGP+TR and VPRJF+TR are basically the same, and the specified accuracy can be achieved after nine iterations. Visually, it is clear that after five iterations, the solution is very close to the final result. In contrast, the VPKAU+TR method reached the maximum number of iterations (100). However, after 30 iterations, even though the result changes a little, it is close to the final result.Forλ2, the change in parameter curves for the two methods VPGP+TR and VPRJF+TR are also basically the same, and the specified accuracy is also achieved after nine iterations. The graph shows that, after four iterations, the solution is very close to the final calculation result. The VPKAU+TR method again reached the maximum number of calculations, but the solution is close to the final result after 20 iterations.Table1 lists the results of the parameters obtained by the three methods.Table 1
Parameter results.
Method
λ
1
λ
2
β
1
β
2
VPGP+TR
3.0067
2.8893
10.5868
1.4004
VPKAU+TR
3.0068
2.8891
10.5868
1.4004
VPRJF+TR
3.0067
2.8893
10.5868
1.4004The table shows that, as for the curves in Figure1, the results obtained by the three methods are basically the same.Table2 gives the runtime (in seconds), number of iterations, residual sum, and root mean square error (RMSE) for the three methods.Table 2
Quantitative comparison of performance for each method.
Method
Runtime (s)
Number of iterations
Residual sum
RMSE
VPGP+TR
40.596687
10
0.1477
0.0859
VPKAU+TR
17.186718
10
0.1477
0.0859
VPRJF+TR
31.073001
100
0.1477
0.0859The quantitative indicators (residual sums and RMSE) are equal and very small, indicating that the three algorithms achieve the same calculation accuracy, which is relatively high. The difference between them is only the number of iterations and runtime. The VPGP+TR method takes the longest to calculate because this method uses the most complicated way to calculate the Jacobian matrix. In contrast, VPRJF+TR is a simple method for calculating the Jacobian matrix. However, to meet the accuracy requirements, the number of iterations is also large, so the calculation time is also long. The VPKAU+TR method somewhat simplifies the calculation of the Jacobian matrix, and to some extent, the accuracy of the calculation is guaranteed, so the runtime is the shortest. For this simulation experiment, the VPKAU+TR method performs the best of all methods.


## 3.2. Waveform Parameter Calculation for Airborne Laser Radar Sounding Data

The second experiment used laser radar sounding waveform data. The experimental data were flight experimental data obtained from a flight over the sea near Sanya City, Hainan Province, China.The waveform sample data are stored as an array, and the amplitude of each sample point is recorded separately. The sample information is presented as a ‘waveform’ form on a two-dimensional plane. The onboard laser system records discrete echoes that simulate the “waveform” of the raw sampled data using the most appropriate functional model. To improve the sounding accuracy, it is highly desirable to choose a high-precision theoretical model.This experiment has a total of 320 waveform data, and the function model is as follows:(22)y=a+be−x−μ2/2σ2.Selecting reasonable data can effectively reduce the influence of noise. In this study, we selected some of the data that have no gross errors. We extracted 42 data from the waveform data to determine the parameters. When iterating, we set the initial value toλ0=15,3T, the maximum number of iteration steps to 100, and the error to ε=10−3. SVD for Φ yielded singular values of 6.5941 and 1.9585.Table3 lists the calculation results of the model parameters.Table 3
Parameter results.
Method
a
b
μ
σ
VPGP+TR
207.1076
498.6737
14.0070
3.4534
VPKAU+TR
207.0941
498.6461
14.0068
3.4540
VPRJF+TR
207.0518
498.5602
14.0085
3.4560The results given in the table show that the results obtained by the three methods are basically the same, indicating that all three methods are feasible approaches. The result obtained by VPRJF+TR is slightly different from those of the other two methods because this method simplifies the calculation of the Jacobian matrix the most, resulting in a gap.Figure3 shows the fitting curves of the model.Figure 3
Waveform data and fitting curves. (a) VPGP+TR. (b) VPKAU+TR. (c) VPRJF+TR.
(a)
(b)
(c)In Figure3, all three methods fit the curve well, which also intuitively illustrates the feasibility of these three methods.Figure4 shows the change in parameters μ and σ with respect to the number of iterations.Figure 4
Change in parameters with respect to the number of iterations. (a) VPGP+TR. (b) VPKAU+TR. (c) VPRJF+TR.
(a)
(b)
(c)The results in Figure4 show that for parameters μ and σ, when error ε=10−3, the parameter change is basically stable after five iterations. This shows that the convergence of these three methods is fast and accurate.Table4 lists the runtime (in seconds), number of iterations, and root mean square error (RMSE).Table 4
Quantitative comparison of performance for each method.
Method
Runtime (s)
Number of iterations
RMSE
VPGP+TR
263.797378
10
38.0979
VPKAU+TR
91.929637
10
38.0980
VPRJF+TR
9.817583
10
38.0983Forε=10−3, number of iterations and RMSE (i.e., the convergence speed and accuracy) of these three methods are basically the same. With respect to runtime, the VPGP+TR method again takes the longest because of its calculation of the Jacobian matrix. Hence, the same number of iterations takes longer. Similarly, the Jacobian matrix of the VPRJF+TR method is the simplest and requires the shortest time to calculate. For the VPKAU+TR method, the computational complexity is between those of the other two methods, and it takes moderate amount of time to calculate.


## 4. Conclusion

For the separable nonlinear least squares problem, this paper establishes a solution model consisting of an LM algorithm combined with TR. For linear parameter estimation, the sensitivity of the model to disturbance is reduced by the TR method. For the nonlinear parameter estimation, the LM algorithm is adopted and the Jacobian matrix is calculated by the Golub and Pereyra, Kaufman, and Ruano methods. Combining nonlinear and linear parameter estimation methods, three estimation models are obtained and the stability and accuracy of model estimation are demonstrated. The model was verified by simulation data in the experiment in Section3.1 and real data in the experiment in Section 3.2. The experimental results also illustrate the feasibility of the model.In future research, we will apply the Jacobian matrix decomposition to the VP algorithm. This decomposition can simplify the model and improve the efficiency of the algorithm.



---

*Source: 4861708-2019-11-25.xml*