# Efficient Parameters Estimation Method for the Separable Nonlinear Least Squares Problem

Authors: Ke Wang; Guolin Liu; Qiuxiang Tao; Min Zhai
Journal: Complexity
 (2020)

---

## Abstract

In this work, we combine the special structure of the separable nonlinear least squares problem with a variable projection algorithm based on singular value decomposition to separate linear and nonlinear parameters. Then, we propose finding the nonlinear parameters using the Levenberg–Marquart (LM) algorithm and either solve the linear parameters using the least squares method directly or by using an iteration method that corrects the characteristic values based on the L-curve, according to whether or not the nonlinear function coefficient matrix is ill posed. To prove the feasibility of the proposed method, we compared its performance on three examples with that of the LM method without parameter separation. The results show that (1) the parameter separation method reduces the number of iterations and improves computational efficiency by reducing the parameter dimensions and (2) when the coefficient matrix of the linear parameters is well-posed, using the least squares method to solve the fitting problem provides the highest fitting accuracy. When the coefficient matrix is ill posed, the method of correcting characteristic values based on the L-curve provides the most accurate solution to the fitting problem.


---

## Body

## 1. Introduction

The separable nonlinear least squares problem is a special case of nonlinear least squares, and its model can be expressed as a linear combination of nonlinear functions. More generally, one can also consider that there are two sets of unknown parameters, where one set is dependent on the other and can be explicitly eliminated. This problem was first proposed in parameter estimation of the atomic physics particle half-life formula developed by Golub and Pereyra in 1973 [1]. In real life, models of this type are very common. For example, in the machine learning community, neural networks, their numerous variants [2, 3], and some neuro fuzzy systems [4, 5] take the form of a linear ensemble of some nonlinear basis functions. In the field of signal processing, Prony’s method [6, 7], which takes the sum of complex exponentials, is frequently used to analyze the frequency components of a signal. In the field of algorithm application, waveform data decomposition is one of the key steps in processing based on airborne full-waveform light detection and ranging (LiDAR) data. The full waveform can be decomposed into a linear combination of multiple Gaussian functions. Through waveform data decomposition, discrete point cloud and waveform parameter information can be obtained [8]. Furthermore, this approach has many applications in areas such as mechanical systems [9], telecommunications [10], robotics [11], and environmental sciences [12]. These applications can be viewed as nonlinear data-fitting problems in terms of numerical expression. However, data-fitting problems are often quite challenging numerically. Fortunately, by exploiting the special structure of separable nonlinear models, efficient algorithms can be obtained.Based on the special structure of the separable nonlinear least squares problem, Golub and Pereyra [1] proposed the variable projection (VP) algorithm to eliminate linear parameters and obtain simplified functions involving only nonlinear parameters and also used the Levenberg–Marquart (LM) algorithm for its solution. Dimensionality reduction of parameter space improves the possibility of obtaining a global optimal solution [13]. Therefore, improvement and application [14] of VP have been conducted on that basis. Kaufman [15] proposed a modified VP algorithm based on trapezoidal decomposition to simplify the calculation, which reduced its computational complexity and improved computational efficiency. Subsequently, Ruano et al. [16] proposed an improved VP algorithm based on QR decomposition for the sparse case of nonlinear function matrix, which also effectively improved operation efficiency. Further, Gan and Li [17] proposed a VP algorithm based on Gram–Schmidt matrix decomposition for the case in which the number of observations is much larger than the number of linear parameters, which reduces the amount of calculation required.In this study, in view of the ill-posed condition of a nonlinear function matrix, singular value decomposition (SVD) [18] is adopted to simplify the VP algorithm and improve the stability of the calculation. Then, after linear parameters are eliminated by the improved VP algorithm, the separable least squares problem is transformed into an optimization problem containing only nonlinear parameters [19]. As regards an estimation method for nonlinear parameters, Liu et al. [20] combined this with sequential quadratic programming (SQP), developing a gradient-based optimization algorithm to determine the optimal time-delays and system parameters in a novel dynamic optimization problem for nonlinear multistage systems with time delays. However, this approach was restricted to parameter identification problems.The common nonlinear least squares iterative solutions [21] are the gradient descent method [22], Gauss–Newton method [23], and Levenberg–Marquart (LM) method [24, 25]. For example, in [26], the nonlinear least squares problem of the distributional robust parameter identification model for time-delay systems is transformed into a single-level optimization problem and a gradient-based optimization method is developed to solve the transformed problem. The method only involves the first-order moment information and is simple to calculate. To obtain robust estimates against the noise in measurements, Liu et al. [27] proposed a robust estimation formulation, in which the cost function was the variance of the error function and an additional constraint indicates an allowable sacrifice from the optimal expectation value of the classical estimation problem. On this basis, a gradient-based optimization algorithm to numerically solve the classical and robust parameter estimation problems was developed. This is more efficient than the existing methods used for problems, where optimization parameters outnumber constraints. This method involves simple calculations, but its convergence speed is generally slower than that of the Gauss–Newton method.Torres et al. [28] used the sequentially semiseparable matrix to calculate the Jacobian matrix [29] and Hessian matrix [30], employing the Gauss–Newton method to optimize the output error of the global system. The effectiveness of the algorithm was verified by numerical examples. Bellavia et al. [31] improved the approximation function by controlling the accuracy level when the accuracy was too low to be optimized, and then proposed the LM method based on the dynamic precision relationship between the evaluation function and gradient for solving large-scale nonlinear least squares problems. They proved the global and local convergence and complexity of this method. The Gauss–Newton method has the advantages of fast convergence and high precision. However, the Jacobian matrix is required to be of full rank in the iterative process. If the problem has high nonlinearity or the residual is large, the method may not produce convergence. Therefore, the LM algorithm overcomes this shortcoming and adjusts the damping parameters according to the idea of a trust region, to effectively control the direction of iterative descent. Once the nonlinear parameters were determined, the least squares (LS) method [32] is used to estimate the linear parameters.In view of a general nonlinear multistage system with time-delay and system parameters, Liu et al. [33] proposed a new parameter estimation formulation, in which the cost function is the variance of the error function and the constraint indicates an allowable sacrifice from the optimal expectation value of the classical parameter estimation problem. This parameter estimation approach is capable of solving parameter estimation problems with multiple stages and multiple time delays and, compared with classical parameter estimation, it is able to withstand the uncertainty in the distribution of measurement data. Nevertheless, this method has the limitation of relying on the statistical distribution of the noisy measurement output.To enhance the estimation accuracy, Ding et al. [34, 35] presented a filtering-based gradient iterative algorithm and a filtering-based least squares iterative algorithm, which resulted in improved convergence speed. However, when nonlinear parameter estimation led to the ill condition of the linear least squares coefficient matrix, the estimation results of the linear parameters obtained by the LS method were unstable and even sometimes significantly different from the true values. For this problem, the regularization method is often used to solve for ill-posed problems in linear parameter estimation [36], such as the Tikhonov regularization method [37], truncated singular value method [38], and iteration by correcting characteristic values [39, 40].In addition, Chen et al. [41] proposed a weighted generalized crossvalidation method to determine Tikhonov regularization parameters for the regularization of separable nonlinear least squares ill-posed problems based on the VP algorithm and verified its effectiveness experimentally. Aiming at the randomness of parameter selection in the iteration by correcting characteristic values in the process of linear least squares solution, Zhai et al. [42] constructed the L-curve [43, 44] of the relationship between the norm of the parameter solution and the residual. They selected the maximum curvature point as the regularization parameter and verified the correctness of the method through numerical experiments.The existing literature presents several effective solutions for the parameter estimation problem, but only a few studies have been conducted on the structural transformation of separable nonlinear models. In this study, we consider the special structure of the separable nonlinear least squares problem, separating two types of parameters using a VP algorithm based on SVD. We then use the LM algorithm to estimate the nonlinear parameters while using the LS method or iteration by correcting characteristic values based on the L-curve method. By comparing the experimental results of parameter estimation of the Gauss function fitting model and fractional fitting model with those of the LM method without parameter separation, the validity of the VP algorithm based on SVD is evaluated and the accuracy of different linear parameter estimation methods is analyzed.The remainder of this paper is outlined as follows. In Section2, based on an improved VP algorithm derived from SVD, the methods of nonlinear parameter estimation and linear parameter estimation are explained and the algorithm for solving separable nonlinear least squares problems is provided. In Section 3, the Gaussian function fitting model and the fractional fitting model experiments are used to compare and analyze the proposed method with the traditional LM algorithm without unseparated parameters. Finally, we present our conclusions in Section 4.


## 2. Solution of Parameters for Separable Nonlinear Least Squares

### 2.1. Variable Projection Algorithm Based on SVD

Consider a set of observationsti,yii=1,2,…,m. The problem of parameter estimation is to find the optimal parameters a^=a1,a2,…,apΤ and b^=b1,b2,…,bqΤ when the following formula reaches the minimum value:(1)a^,b^=argmina,b∑i=1myti−∑j=1pajϕjb;ti2,where φjb;tii=1,2,…,m;j=1,2,…,p is a nonlinear function and p and q are the number of linear and nonlinear parameters to be estimated, respectively. The above formula can be written in a matrix form as(2)a^,b^=argmina,b∑i=1my−Φba2=argmina,by−Φba22,where the column of Φb corresponds to the nonlinear function φjb;t associated with the parameter b and ⋅2 denotes the Euclidean norm.Let(3)fa,b=y−Φba22.For the given nonlinear parametersb, the linear parameters a^ can be estimated by solving the following linear least squares problem:(4)a^=ΦΤbΦb−1Φby=Φ+by,where Φ+b is the pseudoinverse of Φb. Inserting (4) into (3),(5)b^=argmina,by−ΦbΦ+by22=argmina,bI−ΦbΦ+by22=argmina,bI−PΦby22=argmina,bPΦb⊥y22,where PΦb is the orthogonal projector on the linear space spanned by the columns of Φb and PΦb⊥ is the projector on the orthogonal complement of the column space.To simplify the calculation, the matrixΦb, which is composed of nonlinear functions, is decomposed by SVD:(6)Φb=USVΤ,where U is an m×m orthogonal matrix, S is an m×p diagonal matrix, and V is a p×p orthogonal matrix. We obtain Φ+b as(7)Φ+b=VS−1UΤ.Then,(8)PΦb=ΦbΦ+b=USVΤVS−1UΤ=UIp000UΤ.Suppose the rank of the matrixΦb is r, the first r elements on the diagonal of S are not zero, i.e., r≤p. U can be divided into m×r matrix U1 and m×m−r matrix U2. V can be divided into p×r matrix V1 and p×p−r matrix V2. Then,(9)Φb=U1︸m×rU2︸m×m−rSr000V1Τ︸r×pV2Τ︸p−r×p,PΦb⊥=I−U1︸m×rU2︸m×m−rIp000U1Τ︸r×mU2Τ︸m−r×m=−U1︸m×rU2︸m×m−r000Ip−rU1Τ︸r×mU2Τ︸m−r×m=U2U2Τ.The matrix composed ofm residual functions is simplified to the following equation by VP based on SVD: (10)Fb=PΦb⊥y=U2U2Τy.Then, the objective function of the separable nonlinear least squares problem is simplified to(11)fb=Fb22=PΦb⊥y22=U2U2Τy22.Theorem 1.
We assume that in the open setΩ∈ℜk, Φb has constant rank r=minp,q.(a)
Ifb^ is a critical point (or a global minimizer in Ω) of fb, and combined with (4), then a^,b^ is a critical point of fa,b (or a global minimizer for b∈Ω) and fa,b=fb(b)
Ifa^,b^ is a global minimizer of fa,b for b∈Ω, then b^ is a global minimizer of fb in Ω and fb=fa,b. Furthermore, if there is a unique b^ among the minimizing pairs of fa,b, then a^ must satisfy (4).Proof of Theorem 1.
From (3) we see that fa,b=y−Φba22. Therefore,(12)12∇fa^,b^=−PΦb⊥yΤΦbΦ+by⊕Φb,where ⊕ signifies a direct sum.
Assume thatb^ is a global minimizer of fb in Ω and a is defined by (4); then,(13)12∇fa^,b^=−PΦb⊥yΤΦbΦ+by⊕Φb=12∇fb^⊕0.
BecauseyΤPΦb⊥Φb=0, then a^,b^ is a critical point of fa,b.
Assume thatb^ is a global minimizer of fb in Ω and a^ satisfies (4). Then, clearly, fa^,b^=fb^. If we assume that there exists a∗,b∗,a∗∈Ω, such that fa∗,b∗<fa^,b^. Because for any b we have fb≤fa,b, then it follows that fb∗≤fa∗,b∗<fa^,b^=fb^, which is a contradiction to the fact that b^ is a global minimizer of fb in Ω. Therefore, a^,b^ is a global minimizer of fa,b in Ω, and part (a) of the theorem is proved.
Conversely, suppose thata^,b^ is a global minimizer of fa,b in Ω. Then, as mentioned above fb^≤fa^,b^. Now, let a∗=Φ+b^y. Then, we have fb^=fa∗,b^≤fa^,b^; however, because a^,b^ is a global minimizer, we must have equality. If there is a unique a among the minimizers of fa,b, then a∗≡a^. We still have that b^ is a global minimizer of fb^. Assume that it is not. Thus, there will be b¯∈Ω, such that fb¯<fb^. Let a¯ be equal to Φ+b¯y. Then, fb¯=fa¯,b¯<fb^=fa^,b^, which is a contradiction to the fact that a^,b^ is a global minimizer of fa,b. This completes the proof.


### 2.2. Nonlinear Parameter Estimation Using the LM Algorithm

For a separable least squares problem with only nonlinear parameters, the LM algorithm is adopted for the solution and the nonlinear parametersb is updated as(14)bk+1=bk+αkdk,where bk is current nonlinear parameter vector. αk is a small step length that ensures the decrease of the objection function (11), which is calculated by an imprecise search method, such as line search, in which we let mk be the smallest nonnegative integer m satisfying fxk+ρmdk≤fxk+σρmgkΤdk in the kth iteration process. Then,(15)αk=ρmk,where dk is search direction, which can be determined by the following equation:(16)dk=JkΤJk+μkI−1JkFk,where Jk∈ℜm×p is the Jacobian matrix of Fb. Kaufman [15] formulated an explicit analytic expression for the Jacobian. This ensures the efficiency and reliability of the VP algorithms. In addition, combined with the simplified objective function by SVD, the Jacobian matrix Jk of the residual vector F is expressed as(17)Jk=Fb=PΦb⊥y=−PΦb⊥ΦbΦ+by=−U2U2ΤΦbV1Sr−1U1Τy,where μk is damping parameter. It is adjusted with a strategy similar to the trust region radius [45], and a quadratic function is defined at the current iteration point as follows:(18)vd=fxk+JkΤFkΤd+12dΤJkΤJkΤd.Then, the incremental ratio ofvxk+dk to the objective function is calculated:(19)ηk=ΔfdkΔvdk=fxk+1−fxkJkΤFkΤdk+1/2dkΤJkΤJkΤdk.Whenηk is close to one, the fitting between the quadratic function vd and objective function is good at point xk. LM is used to solve the nonlinear least squares problem. The parameter μ should be smaller; that is, in this case, the Gauss–Newton method is more effective. When it is close to zero, the fitting between the quadratic function vd and objective function is poor at point xk. LM is used to solve the nonlinear least squares problem, and the parameter μ should be larger. When ηk is neither close to zero nor one, then μk is suitable and does not need to be adjusted. The critical values of η are usually 0.25 and 0.75, and the adjustment rules of μ are as follows.Whenηk is close to one, the fitting between quadratic function vd and the objective function is good at point xk. The parameter μ should be smaller when the LM algorithm is used to solve the nonlinear least squares problem. When ηk is close to zero, the fitting between the quadratic function vd and the objective function is poor at point xk. The parameter μ should be larger when LM is used to solve the nonlinear least squares problem. When ηk is neither close to zero nor one, μ is suitable and does not need to be adjusted. The critical values are usually 0.25 and 0.75, and the adjustment rules of μ are as follows:(20)μk+1=0.1μk,ηk>0.75,μk,0.25≤ηk≤0.75,10μk,ηk<0.25.


### 2.3. Linear Parameter Estimation by Correcting Characteristic Value Based on L-Curve

The nonlinear parameter estimationsb^ are solved in Section (2); the linear parameters a^ can then be calculated by (4). When ΦΤbΦb is a nonsingular matrix, the least squares method is used to calculate it directly. However, when ΦΤbΦb is singular or the condition number of Φb exceeds 100, solution of the linear least squares method is unstable or cannot be solved. Therefore, it needs to be solved by the regularization method. In this study, an iteration method that corrects characteristic values based on the L-curve is used to solve the problem.The least squares normal equation for the linear parameter estimation is(21)ΦΤbΦba^=ΦΤby.Addingλa to both sides of (21), we obtain(22)ΦΤbΦb+λIa^=ΦΤby+λa^.LetN=ΦΤbΦb and W=Φby, the step of the iteration method that corrects the characteristic value is(23)a^k+1=N+λI−1W+λa^k=N+λI−1+N+λI−12+⋯+N+λI−1k+1W+N+λI−1k+1a^k,where a^k is the current nonlinear parameter vector and λ is a regularization parameter selected according to the L-curve method. The L-curve describes the relationship between the norm of the regularized solution aλ22 and the residual norm y−Φba^λ22, which corresponds to each set of regularized parameter values.The convergence of the spectral correction iterative method based on the L-curve is proved as follows.Let the initial estimated value ofa be a^0, NLP=N+λI−1. The iterative calculation process of spectral correction can then be written as(24)a^1=N+λI−1W+N+λI−1a^0=NLPW+NLPa^0,a^2=NLPW+NLPa^1=NLPW+NLP⋅NLPW+NLP⋅NLPa^0,⋮a^i=NLPW+NLPa^i−1=NLPW+NLP⋅NLPW+⋯+NLP⋯NLP︸iW+NLP⋯NLP︸ia^0.Wheni=m, m is the number of inadequate iterations of the spectral modified iteration method; therefore, the relationship between estimated value a^m and iteration initial value a^0 is(25)a^m=NLP+NLP⋅NLP+⋯+NLP⋯NLP︸mW+NLP⋯NLP︸ma^0.Wheni=j, j is the number of sufficient iterations of the spectral modified iterative method, and m<j; therefore, the relationship between estimated value a^m and full iteration result a^0 is(26)a^m=NLP⋯NLP︸m−ja^j−NLP⋯NLP︸m−jNLP+NLP⋯NLP+⋯+NLP⋯NLP︸j−mW.According to (23), the relation between the estimated result a^j of spectral correction and the initial iteration value a^0 is as follows:(27)a^i=NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW+NLP⋯NLP︸ia^0=I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW−W+NLP⋯NLP︸ia^0.In the case ofRankN=t, N is the tth order positive definite matrix. An orthogonal matrix of order t is Q that follows the expression(28)N+λI=QDQΤ+λI=QD+λIQΤ,where D=diagσ1,σ2,…,σt,σ1≥σ2≥⋯≥σt>0, is the eigenvalue of N.LetD=diagσ1+λ,σ2+λ,…,σt+λ, then(29)NLP=N+λI−1=QD+λI−1QΤ=QD1QΤ,NLP22=maxσNLPΤNLP1/2=max1≤l≤tσl+1−1=σt+1−1<1.According to the convergence of the matrix sequence, it can be known that, in the case ofNLP<1, matrix NLP⋯NLP︸i is convergent and limi⟶∞NLP⋯NLP︸i=0. Matrix power series ∑i=0∞NLP⋯NLP︸i is absolutely convergent, and the sum is I−NLP−1.Therefore,(30)limi⟶∞a^i=limi⟶∞I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW−W+limi⟶∞NLP⋯NLP︸ia^0=I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸i+⋯W−W=I−NLP−1−IW=I−NLP−1I−I−NLPW=NLP−1I−NLP−1W=NLP−1−I−1W=N−1W.That is, in the case ofRankN=t, the spectral correction iterative method based on L-curve converges to the least squares solution for any initial value. Therefore, the relationship between the estimated result a^m and the iteration termination value a^j can be further rewritten as(31)a^m=NLP⋯NLP︸m−jN−1W−NLP⋯NLP︸m−j+1+NLP⋯NLP︸m−j+2+⋯+NLP−1+IW.In the case ofRankN<t, the initial value of a^0 is equal to 0 in the spectral correction iterative method. Further, limi⟶∞a^i=N−1W.Therefore, this method is convergent with the least squares method and is also feasible.


### 2.4. Summarization of the Algorithm

According to the estimation methods of linear parameters and nonlinear parameters proposed in Sections2.1–2.3, the separable nonlinear least squares solution method—in which nonlinear parameters in the nonlinear least squares problem with SVD variable projection separation are solved with the LM algorithm and linear parameters are solved directly with the least squares method—is referred to as LMVP + LS. When the linear parameters are solved using the iteration method that corrects characteristic values based on the L-curve, the method is referred to as LMVP + CCVL. The traditional LM method without separation of parameters is referred to as LMunSep. The steps for solving the separable nonlinear least squares problem are summarized as follows.A flowchart of Algorithm1 is presented in Figure 1.Algorithm 1: Parameter estimation algorithm for the separable nonlinear least squares problem.
Step 1: given initial valuesa0∈ℜp,b0∈ℜq, set μ0=0.1, ρ,σ∈0,1, k=0. ε is a very small positive number.
Step 2: calculate the matrixΦbk of the nonlinear function and Jacobian matrix Jk of the objective function after VP, using (17).
Step 3: calculate the iteration step length and iteration direction using (15) and (16), and then update the nonlinear parameters using (14).
Step 4: ifgk2=JkFk2≤ε, terminate the algorithm; otherwise, k=k+1, go to Step 3.
Step 5: calculate linear parameters whenΦΤbΦb is a nonsingular matrix. The least squares method is used to calculate this directly; when ΦΤbΦb is singular or the condition number of Φb exceeds 100, the solution is obtained using (23).Figure 1
Flowchart for the parameter estimation process.


## 2.1. Variable Projection Algorithm Based on SVD

Consider a set of observationsti,yii=1,2,…,m. The problem of parameter estimation is to find the optimal parameters a^=a1,a2,…,apΤ and b^=b1,b2,…,bqΤ when the following formula reaches the minimum value:(1)a^,b^=argmina,b∑i=1myti−∑j=1pajϕjb;ti2,where φjb;tii=1,2,…,m;j=1,2,…,p is a nonlinear function and p and q are the number of linear and nonlinear parameters to be estimated, respectively. The above formula can be written in a matrix form as(2)a^,b^=argmina,b∑i=1my−Φba2=argmina,by−Φba22,where the column of Φb corresponds to the nonlinear function φjb;t associated with the parameter b and ⋅2 denotes the Euclidean norm.Let(3)fa,b=y−Φba22.For the given nonlinear parametersb, the linear parameters a^ can be estimated by solving the following linear least squares problem:(4)a^=ΦΤbΦb−1Φby=Φ+by,where Φ+b is the pseudoinverse of Φb. Inserting (4) into (3),(5)b^=argmina,by−ΦbΦ+by22=argmina,bI−ΦbΦ+by22=argmina,bI−PΦby22=argmina,bPΦb⊥y22,where PΦb is the orthogonal projector on the linear space spanned by the columns of Φb and PΦb⊥ is the projector on the orthogonal complement of the column space.To simplify the calculation, the matrixΦb, which is composed of nonlinear functions, is decomposed by SVD:(6)Φb=USVΤ,where U is an m×m orthogonal matrix, S is an m×p diagonal matrix, and V is a p×p orthogonal matrix. We obtain Φ+b as(7)Φ+b=VS−1UΤ.Then,(8)PΦb=ΦbΦ+b=USVΤVS−1UΤ=UIp000UΤ.Suppose the rank of the matrixΦb is r, the first r elements on the diagonal of S are not zero, i.e., r≤p. U can be divided into m×r matrix U1 and m×m−r matrix U2. V can be divided into p×r matrix V1 and p×p−r matrix V2. Then,(9)Φb=U1︸m×rU2︸m×m−rSr000V1Τ︸r×pV2Τ︸p−r×p,PΦb⊥=I−U1︸m×rU2︸m×m−rIp000U1Τ︸r×mU2Τ︸m−r×m=−U1︸m×rU2︸m×m−r000Ip−rU1Τ︸r×mU2Τ︸m−r×m=U2U2Τ.The matrix composed ofm residual functions is simplified to the following equation by VP based on SVD: (10)Fb=PΦb⊥y=U2U2Τy.Then, the objective function of the separable nonlinear least squares problem is simplified to(11)fb=Fb22=PΦb⊥y22=U2U2Τy22.Theorem 1.
We assume that in the open setΩ∈ℜk, Φb has constant rank r=minp,q.(a)
Ifb^ is a critical point (or a global minimizer in Ω) of fb, and combined with (4), then a^,b^ is a critical point of fa,b (or a global minimizer for b∈Ω) and fa,b=fb(b)
Ifa^,b^ is a global minimizer of fa,b for b∈Ω, then b^ is a global minimizer of fb in Ω and fb=fa,b. Furthermore, if there is a unique b^ among the minimizing pairs of fa,b, then a^ must satisfy (4).Proof of Theorem 1.
From (3) we see that fa,b=y−Φba22. Therefore,(12)12∇fa^,b^=−PΦb⊥yΤΦbΦ+by⊕Φb,where ⊕ signifies a direct sum.
Assume thatb^ is a global minimizer of fb in Ω and a is defined by (4); then,(13)12∇fa^,b^=−PΦb⊥yΤΦbΦ+by⊕Φb=12∇fb^⊕0.
BecauseyΤPΦb⊥Φb=0, then a^,b^ is a critical point of fa,b.
Assume thatb^ is a global minimizer of fb in Ω and a^ satisfies (4). Then, clearly, fa^,b^=fb^. If we assume that there exists a∗,b∗,a∗∈Ω, such that fa∗,b∗<fa^,b^. Because for any b we have fb≤fa,b, then it follows that fb∗≤fa∗,b∗<fa^,b^=fb^, which is a contradiction to the fact that b^ is a global minimizer of fb in Ω. Therefore, a^,b^ is a global minimizer of fa,b in Ω, and part (a) of the theorem is proved.
Conversely, suppose thata^,b^ is a global minimizer of fa,b in Ω. Then, as mentioned above fb^≤fa^,b^. Now, let a∗=Φ+b^y. Then, we have fb^=fa∗,b^≤fa^,b^; however, because a^,b^ is a global minimizer, we must have equality. If there is a unique a among the minimizers of fa,b, then a∗≡a^. We still have that b^ is a global minimizer of fb^. Assume that it is not. Thus, there will be b¯∈Ω, such that fb¯<fb^. Let a¯ be equal to Φ+b¯y. Then, fb¯=fa¯,b¯<fb^=fa^,b^, which is a contradiction to the fact that a^,b^ is a global minimizer of fa,b. This completes the proof.


## 2.2. Nonlinear Parameter Estimation Using the LM Algorithm

For a separable least squares problem with only nonlinear parameters, the LM algorithm is adopted for the solution and the nonlinear parametersb is updated as(14)bk+1=bk+αkdk,where bk is current nonlinear parameter vector. αk is a small step length that ensures the decrease of the objection function (11), which is calculated by an imprecise search method, such as line search, in which we let mk be the smallest nonnegative integer m satisfying fxk+ρmdk≤fxk+σρmgkΤdk in the kth iteration process. Then,(15)αk=ρmk,where dk is search direction, which can be determined by the following equation:(16)dk=JkΤJk+μkI−1JkFk,where Jk∈ℜm×p is the Jacobian matrix of Fb. Kaufman [15] formulated an explicit analytic expression for the Jacobian. This ensures the efficiency and reliability of the VP algorithms. In addition, combined with the simplified objective function by SVD, the Jacobian matrix Jk of the residual vector F is expressed as(17)Jk=Fb=PΦb⊥y=−PΦb⊥ΦbΦ+by=−U2U2ΤΦbV1Sr−1U1Τy,where μk is damping parameter. It is adjusted with a strategy similar to the trust region radius [45], and a quadratic function is defined at the current iteration point as follows:(18)vd=fxk+JkΤFkΤd+12dΤJkΤJkΤd.Then, the incremental ratio ofvxk+dk to the objective function is calculated:(19)ηk=ΔfdkΔvdk=fxk+1−fxkJkΤFkΤdk+1/2dkΤJkΤJkΤdk.Whenηk is close to one, the fitting between the quadratic function vd and objective function is good at point xk. LM is used to solve the nonlinear least squares problem. The parameter μ should be smaller; that is, in this case, the Gauss–Newton method is more effective. When it is close to zero, the fitting between the quadratic function vd and objective function is poor at point xk. LM is used to solve the nonlinear least squares problem, and the parameter μ should be larger. When ηk is neither close to zero nor one, then μk is suitable and does not need to be adjusted. The critical values of η are usually 0.25 and 0.75, and the adjustment rules of μ are as follows.Whenηk is close to one, the fitting between quadratic function vd and the objective function is good at point xk. The parameter μ should be smaller when the LM algorithm is used to solve the nonlinear least squares problem. When ηk is close to zero, the fitting between the quadratic function vd and the objective function is poor at point xk. The parameter μ should be larger when LM is used to solve the nonlinear least squares problem. When ηk is neither close to zero nor one, μ is suitable and does not need to be adjusted. The critical values are usually 0.25 and 0.75, and the adjustment rules of μ are as follows:(20)μk+1=0.1μk,ηk>0.75,μk,0.25≤ηk≤0.75,10μk,ηk<0.25.


## 2.3. Linear Parameter Estimation by Correcting Characteristic Value Based on L-Curve

The nonlinear parameter estimationsb^ are solved in Section (2); the linear parameters a^ can then be calculated by (4). When ΦΤbΦb is a nonsingular matrix, the least squares method is used to calculate it directly. However, when ΦΤbΦb is singular or the condition number of Φb exceeds 100, solution of the linear least squares method is unstable or cannot be solved. Therefore, it needs to be solved by the regularization method. In this study, an iteration method that corrects characteristic values based on the L-curve is used to solve the problem.The least squares normal equation for the linear parameter estimation is(21)ΦΤbΦba^=ΦΤby.Addingλa to both sides of (21), we obtain(22)ΦΤbΦb+λIa^=ΦΤby+λa^.LetN=ΦΤbΦb and W=Φby, the step of the iteration method that corrects the characteristic value is(23)a^k+1=N+λI−1W+λa^k=N+λI−1+N+λI−12+⋯+N+λI−1k+1W+N+λI−1k+1a^k,where a^k is the current nonlinear parameter vector and λ is a regularization parameter selected according to the L-curve method. The L-curve describes the relationship between the norm of the regularized solution aλ22 and the residual norm y−Φba^λ22, which corresponds to each set of regularized parameter values.The convergence of the spectral correction iterative method based on the L-curve is proved as follows.Let the initial estimated value ofa be a^0, NLP=N+λI−1. The iterative calculation process of spectral correction can then be written as(24)a^1=N+λI−1W+N+λI−1a^0=NLPW+NLPa^0,a^2=NLPW+NLPa^1=NLPW+NLP⋅NLPW+NLP⋅NLPa^0,⋮a^i=NLPW+NLPa^i−1=NLPW+NLP⋅NLPW+⋯+NLP⋯NLP︸iW+NLP⋯NLP︸ia^0.Wheni=m, m is the number of inadequate iterations of the spectral modified iteration method; therefore, the relationship between estimated value a^m and iteration initial value a^0 is(25)a^m=NLP+NLP⋅NLP+⋯+NLP⋯NLP︸mW+NLP⋯NLP︸ma^0.Wheni=j, j is the number of sufficient iterations of the spectral modified iterative method, and m<j; therefore, the relationship between estimated value a^m and full iteration result a^0 is(26)a^m=NLP⋯NLP︸m−ja^j−NLP⋯NLP︸m−jNLP+NLP⋯NLP+⋯+NLP⋯NLP︸j−mW.According to (23), the relation between the estimated result a^j of spectral correction and the initial iteration value a^0 is as follows:(27)a^i=NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW+NLP⋯NLP︸ia^0=I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW−W+NLP⋯NLP︸ia^0.In the case ofRankN=t, N is the tth order positive definite matrix. An orthogonal matrix of order t is Q that follows the expression(28)N+λI=QDQΤ+λI=QD+λIQΤ,where D=diagσ1,σ2,…,σt,σ1≥σ2≥⋯≥σt>0, is the eigenvalue of N.LetD=diagσ1+λ,σ2+λ,…,σt+λ, then(29)NLP=N+λI−1=QD+λI−1QΤ=QD1QΤ,NLP22=maxσNLPΤNLP1/2=max1≤l≤tσl+1−1=σt+1−1<1.According to the convergence of the matrix sequence, it can be known that, in the case ofNLP<1, matrix NLP⋯NLP︸i is convergent and limi⟶∞NLP⋯NLP︸i=0. Matrix power series ∑i=0∞NLP⋯NLP︸i is absolutely convergent, and the sum is I−NLP−1.Therefore,(30)limi⟶∞a^i=limi⟶∞I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸iW−W+limi⟶∞NLP⋯NLP︸ia^0=I+NLP+NLP⋅NLP+⋯+NLP⋯NLP︸i+⋯W−W=I−NLP−1−IW=I−NLP−1I−I−NLPW=NLP−1I−NLP−1W=NLP−1−I−1W=N−1W.That is, in the case ofRankN=t, the spectral correction iterative method based on L-curve converges to the least squares solution for any initial value. Therefore, the relationship between the estimated result a^m and the iteration termination value a^j can be further rewritten as(31)a^m=NLP⋯NLP︸m−jN−1W−NLP⋯NLP︸m−j+1+NLP⋯NLP︸m−j+2+⋯+NLP−1+IW.In the case ofRankN<t, the initial value of a^0 is equal to 0 in the spectral correction iterative method. Further, limi⟶∞a^i=N−1W.Therefore, this method is convergent with the least squares method and is also feasible.


## 2.4. Summarization of the Algorithm

According to the estimation methods of linear parameters and nonlinear parameters proposed in Sections2.1–2.3, the separable nonlinear least squares solution method—in which nonlinear parameters in the nonlinear least squares problem with SVD variable projection separation are solved with the LM algorithm and linear parameters are solved directly with the least squares method—is referred to as LMVP + LS. When the linear parameters are solved using the iteration method that corrects characteristic values based on the L-curve, the method is referred to as LMVP + CCVL. The traditional LM method without separation of parameters is referred to as LMunSep. The steps for solving the separable nonlinear least squares problem are summarized as follows.A flowchart of Algorithm1 is presented in Figure 1.Algorithm 1: Parameter estimation algorithm for the separable nonlinear least squares problem.
Step 1: given initial valuesa0∈ℜp,b0∈ℜq, set μ0=0.1, ρ,σ∈0,1, k=0. ε is a very small positive number.
Step 2: calculate the matrixΦbk of the nonlinear function and Jacobian matrix Jk of the objective function after VP, using (17).
Step 3: calculate the iteration step length and iteration direction using (15) and (16), and then update the nonlinear parameters using (14).
Step 4: ifgk2=JkFk2≤ε, terminate the algorithm; otherwise, k=k+1, go to Step 3.
Step 5: calculate linear parameters whenΦΤbΦb is a nonsingular matrix. The least squares method is used to calculate this directly; when ΦΤbΦb is singular or the condition number of Φb exceeds 100, the solution is obtained using (23).Figure 1
Flowchart for the parameter estimation process.


## 3. Numerical Examples

In this study, we used three examples, the Gaussian function fitting model, the fractional fitting model, and decomposition of full-waveform LiDAR data, to verify the method proposed in Section2 for solving separable nonlinear least squares problems. The results were then compared with the LM algorithm with unseparated parameters with respect to iteration times, function calculation times, and fitting accuracy. The experiments were performed using MATLAB 2016b on a 2.3-GHz desktop PC running Windows 10.


### 3.1. Example 1: Gaussian Function Fitting Model

The Gaussian function fitting model often exists in the parameter estimation of signal processing. For the generalized Gaussian model, there is a linear combination of nonlinear functions that can be solved by separable nonlinear least squares based on the VP algorithm. The form of the model is as follows:(32)ya,b;ti=∑j=1pajexp−ti−uj22θj2+ξti,i=1,2,…,m,where m is the number of observations, p is the number of linear parameters, and ξti is the Gaussian noise. a=a1,a2,a3,a4,a5,a6Τ is the linear parameter vector to be estimated. b=u1,θ1,u2,θ2,u3,θ3,u4,θ4,u5,θ5,u6,θ6Τ is the nonlinear parameter vector to be estimated. Equation (32) is written in the form of a matrix as(33)ya,b;t=Φm×pu,θ;tap×1+ξm×1,where Φu,θ;t=ϕ1,ϕ2,…,ϕp and φi=exp−t−ui2/2θi2. The true parameter values of the model in (32) are a∗=50,40,60,60,40,50Τ, u∗=50,100,150,200,250,300Τ, and θ∗=10,5,15,12,20,15Τ. The observation error ξ is a normal distribution with a standard deviation of 0.5. According to the above parameters, ten groups of observations are randomly generated, where each group has 400 data points. The initial values of the nonlinear and linear parameters are randomly selected according to the following uniform distribution:(34)a10∼U45,55,u10∼U40,60,θ10∼U8,12,a20∼U35,45,u20∼U90,110,θ20∼U3,8,a30∼U55,65,u30∼U140,160,θ30∼U10,20,a40∼U55,65,u40∼U180,220,θ40∼U10,15,a50∼U35,45,u50∼U220,280,θ50∼U15,25,a60∼U55,65,u60∼U280,330,θ60∼U15,20.The error equations were listed based on the model in (32) and 400 pairs of observed values ti,yi, and the parameters are solved using the LMunSep, LMVP + LS, and LMVP + CCVL methods, respectively. The parameter estimates are shown in Table 1, which are obtained by taking the average of 10 times.Table 1
Parameter estimation and true values of Example 1.
Parameter
True value
LM
unSep
LM
VP + LS
LM
VP + CCVL
a
1
50
49.9938
49.9938
49.9435
a
2
40
40.1286
40.1286
40.1981
a
3
60
59.9751
59.9751
59.9469
a
4
60
60.0776
60.0776
60.0413
a
5
40
40.0097
40.0097
39.9811
a
6
50
50.0243
50.0243
50.3550
u
1
50
49.9994
49.9994
49.9994
θ
1
10
9.9845
9.9845
9.9845
u
2
100
99.9979
99.9979
99.9979
θ
2
5
4.9961
4.9961
4.9961
u
3
150
150.0055
150.0055
150.0055
θ
3
15
15.0132
15.0132
15.0132
u
4
200
200.0086
200.0086
200.0086
θ
4
12
11.9807
11.9807
11.9807
u
5
250
249.9897
249.9897
249.9897
u
6
300
300.0146
300.0146
300.0146
θ
6
15
14.9772
14.9772
14.9772The fitting curve is shown in Figure2, and the difference between the parameter estimation and true value is shown in Figure 3. Because the estimated parameters obtained by the three methods do not differ significantly, and the fitting curves are basically coincident, only one fitting curve is drawn.Figure 2
Observation points and fitting curves obtained for Example 1 for the three methods.Figure 3
Differences between true value and estimation of parameters obtained by the three methods for Example 1.We can see from Table1 that the results of the three methods are completely consistent in their estimation of the nonlinear parameters (μ^ and θ^). For the estimation of linear parameters, the results obtained by the LMunSep and LMVP + LS methods are identical and the results obtained by the LMVP + CCVL method deviate slightly from the true values compared to those obtained by the other methods. From Figure 2, we can see that the parameters estimated by the three methods can fit the observations very well and that their fitting curves are basically the same. Figure 3 shows that the difference value is closer to zero, and the parameter valuation is closer to the truth value. Therefore, the parameter estimation results of the LMunSep and LMVP + LS methods may be taken to be identical and the difference value is small.To better understand the comparison of the parameter estimation results obtained using the three methods, Table2 lists the maximum, minimum, and average values of the sum of squares of residuals between all the parameter estimation values and the true values (all_SSR), which were calculated ten times. Because the latter two methods eliminate linear parameters using the VP algorithm, the sum of squares of residuals between the nonlinear estimated parameters and the true value (nonpar_SSR) are also listed in Table 2.Table 2
Sum of the squared residuals between the estimated parameters obtained using the three methods and the true values in Example 1.
Method
all_SSR (max)
all_SSR (min)
all_SSR (mean)
nonpar_SSR (max)
nonpar_SSR (min)
nonpar_SSR (mean)
LM
unSep
0.2930
0.0371
0.1715
0.1699
0.0113
0.0424
LM
VP + LS
0.2930
0.0371
0.1715
0.1699
0.0113
0.0424
LM
VP + CCVL
0.9767
0.1864
0.5134
0.1699
0.0113
0.0424We compare the fitting residual sum of squares and calculation process of the three methods in Table3.Table 3
Sum of the squared residuals between predicted values obtained using the three methods and the true observations in Example 1.
Method
y_SSR (max)
y_SSR (min)
y_SSR (mean)
MSE
Iteration times (mean)
Function count (mean)
LM
unSep
5.0763
3.0089
4.2680
0.0107
11.9
251.5
LM
VP + LS
5.0763
3.0089
4.2680
0.0107
9.8
144.9
LM
VP + CVCL
14.3523
7.3535
9.9158
0.0248
9.8
144.9Table2 indicates that, in terms of numerical value, the all_SSR values obtained by LMVP + LS and LMunSep methods are equal. The results of nonpar_SSR using LMVP + CCVL are equal to those of the other two methods; however, the results of all_SSR are higher, indicating that the estimation of linear parameters is worse than that of the other two methods.Table3 indicates that the sum of squares of residual errors between the predicted and true values obtained using LMVP + CCVL is the largest, and those of the other two methods are equal. The mean square errors of the LMunSep, LMVP + LS, and LMVP + CCVL methods are 0.0107, 0.0107, and 0.0248, respectively, which are relatively small, indicating that the results obtained by all three methods are reliable. Compared with the number of iterations and average number of calculations of the function, the method based on VP exhibits a considerable reduction compared with the method in which the parameters are not separated.The residual change of the objective function in the iterative process is shown in Figure4.Figure 4
Residual changes of the three methods in the iteration process for Example 1.
(a)
(b)
(c)According to Figure4, in the case of Gaussian function fitting, on the whole, the residual changes of the three methods are all reduced and the estimated results of the final parameters are very close to the truth value, suggesting that the three methods are convergent.In addition, according to the structural characteristics of the model, when the nonlinear parameters are fixed in theLMVP + LS and LMVP + CCVL methods, the nonlinear function matrix is full rank and the number of conditions is less than five for all ten calculations; thus, there is no ill-condition problem.Based on the experimental results, we draw two conclusions: (1) Because there is no rank deficiency or ill posing in the matrix of the nonlinear function, the LS method helps to attain more accurate linear parameters directly than the regularization method. (2) TheLMVP + LS method eliminates the linear parameters using the VP algorithm based on SVD, reduces the dimension of the parameters to be estimated, and improves the efficiency of the iteration process. The number of iterations required in the LMVP + LS method is less than that for the LMunSep method and produces the best effect under the same precision.


### 3.2. Example 2: Fractional Fitting Model

The fractional fitting model is also a commonly used curve-fitting method. The model is described by the following expression:(35)ya,b;t=a1+a2t+a3t2−a411+b1+0.5b2−t/b32+11+b1−0.5b2−t/b32−a511+b4+0.5b5−t/b62+11+b4−0.5b5−t/b62−a611+b7−t/b82+ξt,where a=a1,a2,a3,a4,a5,a6Τ is a linear parameter vector to be estimated and b=b1,b2,b3,b4,b5,b6,b7,b8Τ is the nonlinear parameter vector. The model in (35) is written in the form of a matrix as follows:(36)ya,b;t=Φ200×6b;ta6×1+ξ200×1,where Φu,θ;t=φ1,φ2,…,φ6 and φii=1,2,…,6 is the coefficient term corresponding to the linear parameter ai. The true values of the parameters in the model in (35) are a∗=5.0,0.5,0.3,0.4,0.5,0.7Τ and b∗=0.5,1,2,0.8,0.8,2.5,0.8,3Τ. Observation error ξ is a normal distribution with a standard deviation of 0.001. According to the above parameters, ten groups of observations are randomly generated, where each group has 200 data points. The initial values of the nonlinear and linear parameters are randomly selected according to the following uniform distribution:(37)a0=a∗+ξa,ξa∼N0,0.01,b0=b∗+ξb,ξb∼N0,0.01.The error equations listed are based on the model in (32) and 200 pairs of observed values ti,yi, and the parameters are solved using the LMunSep, LMVP + LS, and LMVP + CCVL methods. The parameter estimates obtained by taking the average of the ten calculations are listed in Table 4.Table 4
Parameters estimation and true values for Example 2.
Order
Parameter
True value
LM
unSep
LM
VP + LS
LM
VP + CCVL
1
a
^
1
5
4.9631
−14618.7461
4.9955
2
a
^
2
0.5
0.7870
−724.1217
0.5343
3
a
^
3
0.3
0.2143
672.7237
0.2999
4
a
^
4
0.4
0.4850
284.4777
0.3869
5
a
^
5
0.5
0.6116
13585.6419
0.5080
6
a
^
6
0.7
0.6835
−41575.4314
0.7272
7
b
^
1
0.5
0.5356
0.5784
0.5784
8
b
^
2
1
0.9094
0.9502
0.9502
9
b
^
3
2
1.2188
2.0133
2.0133
10
b
^
4
0.8
0.9804
0.8459
0.8459
11
b
^
5
0.8
0.6996
0.7919
0.7919
12
b
^
6
2.5
2.2681
2.5049
2.5049
13
b
^
7
0.8
0.8597
0.8080
0.8080
14
b
^
8
3
3.0277
3.0042
3.0042Because the estimated parameters obtained by the three methods do not differ significantly, and the fitting curves are basically coincident, only one fitting curve is drawn. The fitting curve is shown in Figure5.Figure 5
Observation points and fitting curves obtained using the three methods for Example 2.Table4 and Figure 5 indicate that the three methods produce similar true values for the nonlinear parameters. The calculation results in the LMVP + CCVL and LMVP + LS methods are the same because the separation of the objective function is equal, and the nonlinear parameters are all estimated by the LM algorithm. When the LMVP + LS method is used, significant deviation is observed between the true values of the linear parameters. Combined with Figure 5, the fitting curves of the three methods are similar, which shows that the LMVP + LS method cannot obtain the optimal parameters even though it satisfies the least squares principle in solving the linear parameters. This is because the condition number of the nonlinear function matrix is considerably greater than 100 in the process of solving the separation problem after VP. The number of conditions is above 107 for all ten calculations, and these were seriously ill-posed problems.Because the differences in the linear parameters solved by theLMVP + LS method are much higher than those of the nonlinear parameters, the original difference was reduced by 104 and then plotted. The differences between the parameter estimation solutions of the three methods and the true value are shown in Figure 6.Figure 6
Difference between the true value and the estimation of parameters obtained using the three methods in Example 2.Figure6 indicates that the difference between the parameter estimation and the true value obtained by the LMVP + CCVL method is the smallest among the three methods, that is, the result is closest to the true values. The difference obtained by the LMunSep method is next. The difference obtained by the LMVP + LS method is the largest, even if that is reduced by 104.To better understand the comparison of the parameter estimation results of the three methods, Table5 lists the maximum, minimum, and average values of the sum of squares of residuals between all the parameter estimation values and the true values (all_SSR) obtained from the ten calculations. Because the latter two methods eliminate linear parameters using the VP algorithm, the sums of squares of residuals between nonlinear estimated parameters and the true value (nonpar_SSR) are also listed in Table 5.Table 5
Sums of the squared residuals between the estimated parameters obtained using the three methods and the true values in Example 2.
Method
all_SSR (max)
all_SSR (min)
all_SSR (mean)
nonpar_SSR (max)
nonpar_SSR (min)
nonpar_SSR (mean)
LM
unSep
5.4669
0.9927
2.4807
0.0227
0.0159
0.01915
LM
VP + LS
1.44E + 11
6.59E + 08
2.58E + 10
0.0226
0.0155
0.01896
LM
VP + CCVL
0.307
0.0168
0.1232
0.0226
0.0155
0.01896Table5 indicates that the sum of the residuals of all the parameters using the LMVP + LS method is much larger than the other two methods. The main reason is that the matrix of the nonlinear function is seriously ill posed, and the results of direct solution using the least squares method are not the optimal parameters. Therefore, the LMVP + CCVL method using the iteration method by correcting characteristic values based on the L-curve obtained the optimal estimates that are similar to the true values. The sum of the residual squares of parameter estimates and true values is the smallest. The overall result is also superior to the LMunSep method when the parameters are not separated. For the sum of residual squares of nonlinear parameters, the LMVP + CCVL and LMVP + LS methods exhibit the same results, which are lower than the values obtained by the LMunSep method.A comparison of the three methods in terms of model fitting and calculation processes is presented in Table6.Table 6
Sums of the squared residuals between predicted value obtained by the three methods and the true observations in Example 2.
Method
y_SSR (max)
y_SSR (min)
y_SSR (mean)
MSE
Iteration times (mean)
Function count (mean)
LM
unSep
7.30E − 04
2.60E − 04
3.78E − 04
1.89E − 06
161.4
2516.2
LM
VP + LS
8.37E − 04
3.60E − 04
5.84E − 04
2.92E − 06
20.8
207.7
LM
VP + CCVL
4.12E − 04
3.17E − 05
1.90E − 04
9.52E − 07
20.8
207.7Table6 indicates that in terms of the sum of squared residues between the predicted value and the true value, irrespective of the maximum, minimum, or average value, the lowest values are obtained when the LMVP + CCVL method is used, followed by the LMunSep method; highest values are obtained using the LMVP + LS method. In terms of the mean square error, the LMVP + CCVL method provides the lowest values. The number of iterations and average number of functions considerably decrease for the VP algorithm; the average number of iterations was reduced from 161.4 to 20.8 and the average calculation number of functions was reduced from 2516.2 to 207.7.The residual change in the objective function of the iterative process is shown in Figure7.Figure 7
Residual changes of the three methods in the iteration process of Example 2.
(a)
(b)
(c)According to Figure7, in the case of Gaussian function fitting, on the whole, the residual changes of the three methods are all reduced and the estimated results of the final parameters are very close to the truth value, so it can be seen that the three methods are convergent.The experimental analysis indicates that the problem is seriously ill posed, and when considering separable nonlinear least squares, the results obtained by theLMVP + CCVL method afford the highest fitting precision. The computational efficiency is also much better than that for the LM method without separation.


### 3.3. Example 3: Decomposition of Full-Waveform LiDAR Data

This section of the experimental data originated from the full-waveform LiDAR data of a region collected in 2016. The full-waveform lidar system record the backscattered energy at different elevation points in a certain elevation range in the form of waveforms. The survey area mainly contains buildings, trees, and roads. The point cloud and waveform information are stored in the LAS file in the LAS1.3 standard format. The waveform data sampling interval is 1 ns and the number of samples is 287. The echo waveform is regarded as the superposition of several Gaussian functions. The model form is as follows:(38)ya,b;ti=∑j=14ajexp−ti−uj22θj2+ξti,i=1,2,…,287,where ξti is the Gaussian noise, a=a1,a2,a3,a4Τ is the linear parameter vector to be estimated, and b=μ1,θ1,μ2,θ2,μ3,θ3,μ4,θ4Τ is the nonlinear parameter vector to be estimated.The parameters were solved by theLMunSep, LMVP + LS, and LMVP + CCVL methods. The fitting curves are shown in Figures 8(a)–8(c).Figure 8
Fitting curves of theLMunSep, LMVP + LS, and LMVP + CCVL methods.
(a)
(b)
(c)It can be seen from Figure8 that, in the case of the LMunSep method, only the first decomposition waveform can be fitted to the observation value and the latter three decomposition waveforms are not well fitted. LMVP + LS and LMVP + CCVL provide an accurate fitting result.As there is no truth value in the actual measurement, the estimated result is no longer compared with the truth value. To better understand the comparison of the parameter estimation results of the three methods, the mean square error between the fitting results and the observed values, the maximum fit difference (Diff-max), the minimum fit difference (Diff-min), the number of iterations, and function count calculations of the three methods are presented in Table7.Table 7
Comparison of the solution results of the three methods.
Method
Diff-max
Diff-min
MSE
Iteration times
Function count
LM
unSep
48.0734
0.0191
285.8365
20
290
LM
VP + LS
1.4111
0.0076
0.2577
10
106
LM
VP + CCVL
1.8342
0.0004
0.2577
10
106It can be seen from Table7 that, like the simulation experiment results, owing to parameter separation, the reduced dimension of the parameter estimation improves the chance of convergence. LMVP + LS is the best of the three methods, with the least number of iterations and function calculations. LMVP + CCVL is second. However, the LMunSep method, in which parameters are not separated, does not produce the correct result.The residual changes of the objective function in the iterative process are shown in Figure9.Figure 9
Residual changes of theLMunSep, LMVP + LS, and LMVP + CCVL methods in the iteration process of Example 3.
(a)
(b)
(c)From Figure9, we can know that all three methods tend to converge, but LMunSep, whose parameters are not separated, does not converge to the optimal result, and the sum of squared residual convolutions is still large. The sum of the residuals of the remaining two methods tends to zero. The number of iterations is reduced from 20 to 10, and function count is also decreased from 290 to 106. The calculation is greatly simplified.In short, these experiments indicate that the estimation method of parameter separation shows a great improvement in terms of the number of iterations, number of function calculations, and calculation result.


## 3.1. Example 1: Gaussian Function Fitting Model

The Gaussian function fitting model often exists in the parameter estimation of signal processing. For the generalized Gaussian model, there is a linear combination of nonlinear functions that can be solved by separable nonlinear least squares based on the VP algorithm. The form of the model is as follows:(32)ya,b;ti=∑j=1pajexp−ti−uj22θj2+ξti,i=1,2,…,m,where m is the number of observations, p is the number of linear parameters, and ξti is the Gaussian noise. a=a1,a2,a3,a4,a5,a6Τ is the linear parameter vector to be estimated. b=u1,θ1,u2,θ2,u3,θ3,u4,θ4,u5,θ5,u6,θ6Τ is the nonlinear parameter vector to be estimated. Equation (32) is written in the form of a matrix as(33)ya,b;t=Φm×pu,θ;tap×1+ξm×1,where Φu,θ;t=ϕ1,ϕ2,…,ϕp and φi=exp−t−ui2/2θi2. The true parameter values of the model in (32) are a∗=50,40,60,60,40,50Τ, u∗=50,100,150,200,250,300Τ, and θ∗=10,5,15,12,20,15Τ. The observation error ξ is a normal distribution with a standard deviation of 0.5. According to the above parameters, ten groups of observations are randomly generated, where each group has 400 data points. The initial values of the nonlinear and linear parameters are randomly selected according to the following uniform distribution:(34)a10∼U45,55,u10∼U40,60,θ10∼U8,12,a20∼U35,45,u20∼U90,110,θ20∼U3,8,a30∼U55,65,u30∼U140,160,θ30∼U10,20,a40∼U55,65,u40∼U180,220,θ40∼U10,15,a50∼U35,45,u50∼U220,280,θ50∼U15,25,a60∼U55,65,u60∼U280,330,θ60∼U15,20.The error equations were listed based on the model in (32) and 400 pairs of observed values ti,yi, and the parameters are solved using the LMunSep, LMVP + LS, and LMVP + CCVL methods, respectively. The parameter estimates are shown in Table 1, which are obtained by taking the average of 10 times.Table 1
Parameter estimation and true values of Example 1.
Parameter
True value
LM
unSep
LM
VP + LS
LM
VP + CCVL
a
1
50
49.9938
49.9938
49.9435
a
2
40
40.1286
40.1286
40.1981
a
3
60
59.9751
59.9751
59.9469
a
4
60
60.0776
60.0776
60.0413
a
5
40
40.0097
40.0097
39.9811
a
6
50
50.0243
50.0243
50.3550
u
1
50
49.9994
49.9994
49.9994
θ
1
10
9.9845
9.9845
9.9845
u
2
100
99.9979
99.9979
99.9979
θ
2
5
4.9961
4.9961
4.9961
u
3
150
150.0055
150.0055
150.0055
θ
3
15
15.0132
15.0132
15.0132
u
4
200
200.0086
200.0086
200.0086
θ
4
12
11.9807
11.9807
11.9807
u
5
250
249.9897
249.9897
249.9897
u
6
300
300.0146
300.0146
300.0146
θ
6
15
14.9772
14.9772
14.9772The fitting curve is shown in Figure2, and the difference between the parameter estimation and true value is shown in Figure 3. Because the estimated parameters obtained by the three methods do not differ significantly, and the fitting curves are basically coincident, only one fitting curve is drawn.Figure 2
Observation points and fitting curves obtained for Example 1 for the three methods.Figure 3
Differences between true value and estimation of parameters obtained by the three methods for Example 1.We can see from Table1 that the results of the three methods are completely consistent in their estimation of the nonlinear parameters (μ^ and θ^). For the estimation of linear parameters, the results obtained by the LMunSep and LMVP + LS methods are identical and the results obtained by the LMVP + CCVL method deviate slightly from the true values compared to those obtained by the other methods. From Figure 2, we can see that the parameters estimated by the three methods can fit the observations very well and that their fitting curves are basically the same. Figure 3 shows that the difference value is closer to zero, and the parameter valuation is closer to the truth value. Therefore, the parameter estimation results of the LMunSep and LMVP + LS methods may be taken to be identical and the difference value is small.To better understand the comparison of the parameter estimation results obtained using the three methods, Table2 lists the maximum, minimum, and average values of the sum of squares of residuals between all the parameter estimation values and the true values (all_SSR), which were calculated ten times. Because the latter two methods eliminate linear parameters using the VP algorithm, the sum of squares of residuals between the nonlinear estimated parameters and the true value (nonpar_SSR) are also listed in Table 2.Table 2
Sum of the squared residuals between the estimated parameters obtained using the three methods and the true values in Example 1.
Method
all_SSR (max)
all_SSR (min)
all_SSR (mean)
nonpar_SSR (max)
nonpar_SSR (min)
nonpar_SSR (mean)
LM
unSep
0.2930
0.0371
0.1715
0.1699
0.0113
0.0424
LM
VP + LS
0.2930
0.0371
0.1715
0.1699
0.0113
0.0424
LM
VP + CCVL
0.9767
0.1864
0.5134
0.1699
0.0113
0.0424We compare the fitting residual sum of squares and calculation process of the three methods in Table3.Table 3
Sum of the squared residuals between predicted values obtained using the three methods and the true observations in Example 1.
Method
y_SSR (max)
y_SSR (min)
y_SSR (mean)
MSE
Iteration times (mean)
Function count (mean)
LM
unSep
5.0763
3.0089
4.2680
0.0107
11.9
251.5
LM
VP + LS
5.0763
3.0089
4.2680
0.0107
9.8
144.9
LM
VP + CVCL
14.3523
7.3535
9.9158
0.0248
9.8
144.9Table2 indicates that, in terms of numerical value, the all_SSR values obtained by LMVP + LS and LMunSep methods are equal. The results of nonpar_SSR using LMVP + CCVL are equal to those of the other two methods; however, the results of all_SSR are higher, indicating that the estimation of linear parameters is worse than that of the other two methods.Table3 indicates that the sum of squares of residual errors between the predicted and true values obtained using LMVP + CCVL is the largest, and those of the other two methods are equal. The mean square errors of the LMunSep, LMVP + LS, and LMVP + CCVL methods are 0.0107, 0.0107, and 0.0248, respectively, which are relatively small, indicating that the results obtained by all three methods are reliable. Compared with the number of iterations and average number of calculations of the function, the method based on VP exhibits a considerable reduction compared with the method in which the parameters are not separated.The residual change of the objective function in the iterative process is shown in Figure4.Figure 4
Residual changes of the three methods in the iteration process for Example 1.
(a)
(b)
(c)According to Figure4, in the case of Gaussian function fitting, on the whole, the residual changes of the three methods are all reduced and the estimated results of the final parameters are very close to the truth value, suggesting that the three methods are convergent.In addition, according to the structural characteristics of the model, when the nonlinear parameters are fixed in theLMVP + LS and LMVP + CCVL methods, the nonlinear function matrix is full rank and the number of conditions is less than five for all ten calculations; thus, there is no ill-condition problem.Based on the experimental results, we draw two conclusions: (1) Because there is no rank deficiency or ill posing in the matrix of the nonlinear function, the LS method helps to attain more accurate linear parameters directly than the regularization method. (2) TheLMVP + LS method eliminates the linear parameters using the VP algorithm based on SVD, reduces the dimension of the parameters to be estimated, and improves the efficiency of the iteration process. The number of iterations required in the LMVP + LS method is less than that for the LMunSep method and produces the best effect under the same precision.


## 3.2. Example 2: Fractional Fitting Model

The fractional fitting model is also a commonly used curve-fitting method. The model is described by the following expression:(35)ya,b;t=a1+a2t+a3t2−a411+b1+0.5b2−t/b32+11+b1−0.5b2−t/b32−a511+b4+0.5b5−t/b62+11+b4−0.5b5−t/b62−a611+b7−t/b82+ξt,where a=a1,a2,a3,a4,a5,a6Τ is a linear parameter vector to be estimated and b=b1,b2,b3,b4,b5,b6,b7,b8Τ is the nonlinear parameter vector. The model in (35) is written in the form of a matrix as follows:(36)ya,b;t=Φ200×6b;ta6×1+ξ200×1,where Φu,θ;t=φ1,φ2,…,φ6 and φii=1,2,…,6 is the coefficient term corresponding to the linear parameter ai. The true values of the parameters in the model in (35) are a∗=5.0,0.5,0.3,0.4,0.5,0.7Τ and b∗=0.5,1,2,0.8,0.8,2.5,0.8,3Τ. Observation error ξ is a normal distribution with a standard deviation of 0.001. According to the above parameters, ten groups of observations are randomly generated, where each group has 200 data points. The initial values of the nonlinear and linear parameters are randomly selected according to the following uniform distribution:(37)a0=a∗+ξa,ξa∼N0,0.01,b0=b∗+ξb,ξb∼N0,0.01.The error equations listed are based on the model in (32) and 200 pairs of observed values ti,yi, and the parameters are solved using the LMunSep, LMVP + LS, and LMVP + CCVL methods. The parameter estimates obtained by taking the average of the ten calculations are listed in Table 4.Table 4
Parameters estimation and true values for Example 2.
Order
Parameter
True value
LM
unSep
LM
VP + LS
LM
VP + CCVL
1
a
^
1
5
4.9631
−14618.7461
4.9955
2
a
^
2
0.5
0.7870
−724.1217
0.5343
3
a
^
3
0.3
0.2143
672.7237
0.2999
4
a
^
4
0.4
0.4850
284.4777
0.3869
5
a
^
5
0.5
0.6116
13585.6419
0.5080
6
a
^
6
0.7
0.6835
−41575.4314
0.7272
7
b
^
1
0.5
0.5356
0.5784
0.5784
8
b
^
2
1
0.9094
0.9502
0.9502
9
b
^
3
2
1.2188
2.0133
2.0133
10
b
^
4
0.8
0.9804
0.8459
0.8459
11
b
^
5
0.8
0.6996
0.7919
0.7919
12
b
^
6
2.5
2.2681
2.5049
2.5049
13
b
^
7
0.8
0.8597
0.8080
0.8080
14
b
^
8
3
3.0277
3.0042
3.0042Because the estimated parameters obtained by the three methods do not differ significantly, and the fitting curves are basically coincident, only one fitting curve is drawn. The fitting curve is shown in Figure5.Figure 5
Observation points and fitting curves obtained using the three methods for Example 2.Table4 and Figure 5 indicate that the three methods produce similar true values for the nonlinear parameters. The calculation results in the LMVP + CCVL and LMVP + LS methods are the same because the separation of the objective function is equal, and the nonlinear parameters are all estimated by the LM algorithm. When the LMVP + LS method is used, significant deviation is observed between the true values of the linear parameters. Combined with Figure 5, the fitting curves of the three methods are similar, which shows that the LMVP + LS method cannot obtain the optimal parameters even though it satisfies the least squares principle in solving the linear parameters. This is because the condition number of the nonlinear function matrix is considerably greater than 100 in the process of solving the separation problem after VP. The number of conditions is above 107 for all ten calculations, and these were seriously ill-posed problems.Because the differences in the linear parameters solved by theLMVP + LS method are much higher than those of the nonlinear parameters, the original difference was reduced by 104 and then plotted. The differences between the parameter estimation solutions of the three methods and the true value are shown in Figure 6.Figure 6
Difference between the true value and the estimation of parameters obtained using the three methods in Example 2.Figure6 indicates that the difference between the parameter estimation and the true value obtained by the LMVP + CCVL method is the smallest among the three methods, that is, the result is closest to the true values. The difference obtained by the LMunSep method is next. The difference obtained by the LMVP + LS method is the largest, even if that is reduced by 104.To better understand the comparison of the parameter estimation results of the three methods, Table5 lists the maximum, minimum, and average values of the sum of squares of residuals between all the parameter estimation values and the true values (all_SSR) obtained from the ten calculations. Because the latter two methods eliminate linear parameters using the VP algorithm, the sums of squares of residuals between nonlinear estimated parameters and the true value (nonpar_SSR) are also listed in Table 5.Table 5
Sums of the squared residuals between the estimated parameters obtained using the three methods and the true values in Example 2.
Method
all_SSR (max)
all_SSR (min)
all_SSR (mean)
nonpar_SSR (max)
nonpar_SSR (min)
nonpar_SSR (mean)
LM
unSep
5.4669
0.9927
2.4807
0.0227
0.0159
0.01915
LM
VP + LS
1.44E + 11
6.59E + 08
2.58E + 10
0.0226
0.0155
0.01896
LM
VP + CCVL
0.307
0.0168
0.1232
0.0226
0.0155
0.01896Table5 indicates that the sum of the residuals of all the parameters using the LMVP + LS method is much larger than the other two methods. The main reason is that the matrix of the nonlinear function is seriously ill posed, and the results of direct solution using the least squares method are not the optimal parameters. Therefore, the LMVP + CCVL method using the iteration method by correcting characteristic values based on the L-curve obtained the optimal estimates that are similar to the true values. The sum of the residual squares of parameter estimates and true values is the smallest. The overall result is also superior to the LMunSep method when the parameters are not separated. For the sum of residual squares of nonlinear parameters, the LMVP + CCVL and LMVP + LS methods exhibit the same results, which are lower than the values obtained by the LMunSep method.A comparison of the three methods in terms of model fitting and calculation processes is presented in Table6.Table 6
Sums of the squared residuals between predicted value obtained by the three methods and the true observations in Example 2.
Method
y_SSR (max)
y_SSR (min)
y_SSR (mean)
MSE
Iteration times (mean)
Function count (mean)
LM
unSep
7.30E − 04
2.60E − 04
3.78E − 04
1.89E − 06
161.4
2516.2
LM
VP + LS
8.37E − 04
3.60E − 04
5.84E − 04
2.92E − 06
20.8
207.7
LM
VP + CCVL
4.12E − 04
3.17E − 05
1.90E − 04
9.52E − 07
20.8
207.7Table6 indicates that in terms of the sum of squared residues between the predicted value and the true value, irrespective of the maximum, minimum, or average value, the lowest values are obtained when the LMVP + CCVL method is used, followed by the LMunSep method; highest values are obtained using the LMVP + LS method. In terms of the mean square error, the LMVP + CCVL method provides the lowest values. The number of iterations and average number of functions considerably decrease for the VP algorithm; the average number of iterations was reduced from 161.4 to 20.8 and the average calculation number of functions was reduced from 2516.2 to 207.7.The residual change in the objective function of the iterative process is shown in Figure7.Figure 7
Residual changes of the three methods in the iteration process of Example 2.
(a)
(b)
(c)According to Figure7, in the case of Gaussian function fitting, on the whole, the residual changes of the three methods are all reduced and the estimated results of the final parameters are very close to the truth value, so it can be seen that the three methods are convergent.The experimental analysis indicates that the problem is seriously ill posed, and when considering separable nonlinear least squares, the results obtained by theLMVP + CCVL method afford the highest fitting precision. The computational efficiency is also much better than that for the LM method without separation.


## 3.3. Example 3: Decomposition of Full-Waveform LiDAR Data

This section of the experimental data originated from the full-waveform LiDAR data of a region collected in 2016. The full-waveform lidar system record the backscattered energy at different elevation points in a certain elevation range in the form of waveforms. The survey area mainly contains buildings, trees, and roads. The point cloud and waveform information are stored in the LAS file in the LAS1.3 standard format. The waveform data sampling interval is 1 ns and the number of samples is 287. The echo waveform is regarded as the superposition of several Gaussian functions. The model form is as follows:(38)ya,b;ti=∑j=14ajexp−ti−uj22θj2+ξti,i=1,2,…,287,where ξti is the Gaussian noise, a=a1,a2,a3,a4Τ is the linear parameter vector to be estimated, and b=μ1,θ1,μ2,θ2,μ3,θ3,μ4,θ4Τ is the nonlinear parameter vector to be estimated.The parameters were solved by theLMunSep, LMVP + LS, and LMVP + CCVL methods. The fitting curves are shown in Figures 8(a)–8(c).Figure 8
Fitting curves of theLMunSep, LMVP + LS, and LMVP + CCVL methods.
(a)
(b)
(c)It can be seen from Figure8 that, in the case of the LMunSep method, only the first decomposition waveform can be fitted to the observation value and the latter three decomposition waveforms are not well fitted. LMVP + LS and LMVP + CCVL provide an accurate fitting result.As there is no truth value in the actual measurement, the estimated result is no longer compared with the truth value. To better understand the comparison of the parameter estimation results of the three methods, the mean square error between the fitting results and the observed values, the maximum fit difference (Diff-max), the minimum fit difference (Diff-min), the number of iterations, and function count calculations of the three methods are presented in Table7.Table 7
Comparison of the solution results of the three methods.
Method
Diff-max
Diff-min
MSE
Iteration times
Function count
LM
unSep
48.0734
0.0191
285.8365
20
290
LM
VP + LS
1.4111
0.0076
0.2577
10
106
LM
VP + CCVL
1.8342
0.0004
0.2577
10
106It can be seen from Table7 that, like the simulation experiment results, owing to parameter separation, the reduced dimension of the parameter estimation improves the chance of convergence. LMVP + LS is the best of the three methods, with the least number of iterations and function calculations. LMVP + CCVL is second. However, the LMunSep method, in which parameters are not separated, does not produce the correct result.The residual changes of the objective function in the iterative process are shown in Figure9.Figure 9
Residual changes of theLMunSep, LMVP + LS, and LMVP + CCVL methods in the iteration process of Example 3.
(a)
(b)
(c)From Figure9, we can know that all three methods tend to converge, but LMunSep, whose parameters are not separated, does not converge to the optimal result, and the sum of squared residual convolutions is still large. The sum of the residuals of the remaining two methods tends to zero. The number of iterations is reduced from 20 to 10, and function count is also decreased from 290 to 106. The calculation is greatly simplified.In short, these experiments indicate that the estimation method of parameter separation shows a great improvement in terms of the number of iterations, number of function calculations, and calculation result.


## 4. Conclusion

In this study, linear and nonlinear parameters were separated by the VP algorithm based on SVD. Further, the separable least squares problem was transformed into a least squares problem with only nonlinear parameters, which reduced the dimension of the parameters, number of iterations, and function calculations and improved operation efficiency. In addition, to solve the ill-posed problem of the coefficient matrix comprising nonlinear functions when solving for linear parameters, an iteration method that corrects characteristic values based on the L-curve was adopted. This also helped to ensure the convergence of model parameter estimation and improved prediction accuracy. The parameter estimation method used in our study is suitable for cases with a large number of linear parameters in separable nonlinear least squares problems. Additionally, an iteration method that involves correcting characteristic values based on L-curve was used to solve the ill-conditioned coefficient matrix while solving linear parameters. However, one limitation is that the situation of rank deficit often occurred in the process of separable nonlinear least squares parameter estimation; this will be addressed in future research.



---

*Source: 9619427-2020-02-10.xml*