CN105975759A - L curve-based spectrum correction iteration method - Google Patents

Abstract

The invention discloses an L curve-based spectrum correction iteration method and belongs to the field of biased estimation of ill-conditioned problems. An iteration frequency m corresponding to a maximum curvature of an L curve is obtained by utilizing the L curve generated in a spectrum correction iteration process, and a calculation result when the iteration frequency is m serves as a parameter estimation result. According to the method, the advantage that a spectrum correction model is simple is inherited; an estimation result superior to a least square method can be obtained for the ill-conditioned problem; and the method has relatively high calculation efficiency.

Description

A kind of spectrum correcting iteration method based on L-curve

Technical field

The invention belongs to the biased estimation field of ill-conditioning problem, be specifically related to a kind of spectrum correcting iteration method based on L-curve.

Background technology

During geodetic surveying data processes, ill-conditioning problem is one of problem demanding prompt solution.Embody at gravitational inversion, Engineering Control net adjusted data, the aspect such as Rapid precision locating of GPS network.Coefficient matrix morbid state, can cause parametric solution unstable, and estimates of parameters differs bigger problem with true value.For ill-conditioning problem, biased estimation method has obtained studying widely, and biased estimation method the earliest is that the Stein of homogeneous compaction estimates, subsequently, the biased estimation method such as Principal Component Estimation, ridge estimaion, generalized ridge regres sion, Spectrum correction iteration proposes in succession.

At present, ridge estimaion is to apply biased estimation method widely, has higher computational accuracy, but there is ridge parameter and ask for the problem of difficulty, and conventional ridge parameter acquiring method has GCV method, L-curve method, Hoerl-Baldwin method etc..Generalized ridge regres sion is the extension to ridge estimaion, more sophisticated on theoretical model, but also increase generalized ridge parameter ask for difficulty.Thus, in the biased estimation field of ill-conditioning problem, not yet there is the method for estimation of optimum to disclosure satisfy that all practical application request.

Spectrum correction iteration is a kind of iterative estimate method, its abundant iteration result is consistent with method of least square result, belong to unbiased esti-mator result, and valuation belongs to biased estimator in the middle of its iteration, thus, if the optimal estimation in iterative process can rationally be determined, it is possible to give full play to the simple advantage of its model, for improving the biased estimation computational accuracy of ill-conditioning problem and efficiency further, have very important significance.

Summary of the invention

For above-mentioned technical problem present in prior art, the present invention proposes a kind of spectrum correcting iteration method based on L-curve, reasonable in design, and convenience of calculation improves efficiency, has good promotional value.

To achieve these goals, the present invention adopts the following technical scheme that

A kind of spectrum correcting iteration method based on L-curve, is carried out in accordance with the following steps:

Step 1: utilize formula (1) to parameter estimationCarry out spectrum and revise iterative computation, it is judged that whether parameter estimation variable quantity is less than threshold epsilon, i.e.Whether less than threshold epsilon；

If: judged result is that parameter estimation variable quantity is less than threshold epsilon, i.e.Then terminate iterative computation, perform step 2；

Or judged result is that parameter estimation variable quantity is not less than threshold epsilon, i.e.Then continue iterative computation, until parameter estimation variable quantity is less than threshold epsilon, then perform step 2；

${\hat{X}}^{\left(i\right)}={(B^{T}PB+E)}^{-1}(B^{T}Pl+{\hat{X}}^{(i-1)})---\left(1\right);$

In formula, B is the factor arrays of equation, and P is power battle array, and l is constant term, and E is unit battle array, and i is iterative computation number of times；

Step 2: ask for discrete point set during iterations i in spectrum correction iterative process And calculate adjacent discrete point (x_i,y_iInterpolation polynomial f (x between)_i,y_i), obtain spectrum and revise iterative process L-curve；

Step 3: ask for spectrum and revise the point of maximum curvature, the discrete point (x that this point of detection range is nearest on iterative process L-curve_m,y_mIterations m corresponding to)；

Step 4: by the parameter estimation corresponding to iterations mAs final estimated result.

Preferably, in step 3, specifically include

Step 3.1: according to adjacent discrete point (x in step 2_i,y_iInterpolation polynomial f (x between)_i,y_i), utilize formula k_i=| f_i″|/(1+f_i ^{′ 2})^3/2Calculate in discrete point curvature k a little_i, wherein f ' is f (x_i,y_i) first derivative, f is " for f (x_i,y_i) second dervative；

Step 3.2: the some k of maximum curvature on search L-curve_max, determine apart from the discrete point (x that this point is nearest_m,y_mIterations m corresponding to).

The Advantageous Effects that the present invention is brought:

The present invention proposes a kind of spectrum correcting iteration method based on L-curve, compared with prior art, the present invention inherits the spectrum simple advantage of correction model, improves computational accuracy and computational efficiency, and be easily achieved and extend, be conducive to popularization and the application of spectrum correction model.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of a kind of spectrum correcting iteration method based on L-curve of the present invention.

Fig. 2 is the curve chart that spectrum revises iterative process.

Fig. 3 is that spectrum revises the L-curve figure produced in iterative process.

Fig. 4 be respectively adopted spectrum revised law, L-curve ridge estimaion method, L-curve spectrum revised law add up the calculating time diagram of 10 times continuously.

Detailed description of the invention

Below in conjunction with the accompanying drawings and the present invention is described in further detail by detailed description of the invention:

According to least square indirect adjustment principle, the valuation making parameter X isThen error equation is:

$V=B\hat{X}-l---\left(1\right)$

Wherein, B is factor arrays, and l is constant term.The least square solution of error equation is:

$N\hat{X}=W---\left(2\right)$

N=B in formula^TPB, W=B^TPl.Wherein, P is power battle array.

Add at formula (2) two ends simultaneously:

$(N+E)\hat{X}=W+\hat{X}---\left(3\right)$

In formula, E is unit battle array.

Can be composed correction iterative formula by formula (3) is:

${\hat{X}}^{\left(i\right)}={(N+E)}^{-1}(W+{\hat{X}}^{(i-1)})---\left(4\right)$

In formula,I & lt for X composes correction iterative estimate value.

A kind of spectrum correcting iteration method (as shown in Figure 1) based on L-curve, is carried out the most in accordance with the following steps:

Step 1: utilize formulaIterative computation goes out parameter estimationI is iterative computation number of times, it is judged that whether parameter estimation variable quantity is less than threshold epsilon, i.e.Whether less than threshold epsilon；If parameter estimation variable quantity is less than threshold epsilon, then terminate iterative computation, perform step 2；If parameter estimation variable quantity is not less than threshold epsilon, then continue iterative computation, until judging that parameter estimation variable quantity, less than threshold epsilon, then performs step 2；

Step 2: according to the parameter estimation solved in step 1Asking for iterations successively is discrete point set during iCalculate adjacent discrete point (x_i,y_iInterpolation polynomial f (x between)_i,y_i), it is thus achieved that about the L-curve of iterations i；

Step 3: according to the interpolation polynomial f (x between adjacent discrete point_i,y_i), calculate first derivative f ', second dervative f ", according to formula k_i=| f_i″|/(1+f_i ^{′ 2})^3/2Calculate the curvature of discrete point, the point of maximum curvature on search L-curve, and determine apart from the iterations m corresponding to the discrete point that this point is nearest；

Step 4: by the parameter estimation corresponding to iterations mAs final estimated result.

The present invention is tested, simulates an ill-conditioning problem example.

Unknown parameter has 5, X=[x₁ x₂ x₃ x₄ x₅]^T, parameter true value isInitial parameter values is X₀=[0.5 0.5 0.5 0.5 0.5]^T, add the design matrix after random error and observation vector be shown in formula (5), conditional number cond (A^TA)=20838, in Very Ill-conditioned.

The main contents of this experiment include:

(1) utilize common method of least square (LS method) and Spectrum correction iteration to carry out ill-condition equation to solve, calculate parameter estimationTable 1 lists each method valuationWith true valueThe norm of differenceSize.

Table 1 method of least square revises solution results contrast with spectrum

(2) utilizing Spectrum correction iteration, ridge estimaion method based on L-curve and Spectrum correction iteration based on L-curve to carry out ill-condition equation to solve, table 2 lists each method valuationWith true valueThe norm of differenceSize.

Fig. 3 is that spectrum revises the L-curve figure produced in iterative process.

The results contrast of 2 three kinds of different solutions of table

(3) in order to compare the computational efficiency of algorithms of different, be respectively adopted spectrum revised law, L-curve ridge estimaion method, L-curve spectrum revised law calculate.Statistics each method continuously repeats the calculating time of 10 times, and Fig. 4 is shown in statistical result.Experiment computer processor used is Intel (R) Core (TM) i3M 380 2.53GHz, and internal memory 4GB, operating system is Windows 7, and algorithm implementation tool is MATLAB 7.8.0 (R2009a).

By result above it can be seen that

(1) combine shown in table (1), Fig. 2, the result of calculation of method of least square and the Spectrum correction iteration result after carrying out abundant iteration 400 times is consistent, demonstrating the Spectrum correction iteration result of calculation after abundant iteration is unbiased, and the parameter estimation that tradition Spectrum correction iteration result of calculation is when being iteration 138 times, result is better than method of least square, thus traditional Spectrum correction iteration is biased estimation, and when there is ill-conditioning problem when, biased estimation better than unbiased esti-mator.

(2) combine shown in table (2), Spectrum correction iteration, ridge estimaion method based on L-curve and Spectrum correction iteration result of calculation based on L-curveIt is respectively 1.2923,0.8547 and 0.8171, it is seen that ridge estimaion method result of calculation based on L-curve is better than Spectrum correction iteration, and Spectrum correction iteration computational solution precision based on L-curve is the highest.

(3) combine shown in Fig. 4, the relatively computational efficiency of three kinds of methods, wherein, it is 0.3640s that Spectrum correction iteration calculates time average, it is 0.9398s that ridge estimaion method based on L-curve calculates time average, it is 0.5077s that Spectrum correction iteration based on L-curve calculates time average, and therefore Spectrum correction iteration based on L-curve inherits the spectrum simple advantage of revised law model, and computational efficiency is better than ridge estimaion method based on L-curve.

Certainly, described above is not limitation of the present invention, and the present invention is also not limited to the example above, change that those skilled in the art are made in the essential scope of the present invention, retrofits, adds or replaces, and also should belong to protection scope of the present invention.

Claims (2)

1. a spectrum correcting iteration method based on L-curve, it is characterised in that: carry out in accordance with the following steps:

Step 1: utilize formula (1) to parameter estimationCarry out spectrum and revise iterative computation, it is judged that whether parameter estimation variable quantity Less than threshold epsilon, i.e.Whether less than threshold epsilon；

If: judged result is that parameter estimation variable quantity is less than threshold epsilon, i.e.Then terminate iteration meter Calculate, perform step 2；

Or judged result is that parameter estimation variable quantity is not less than threshold epsilon, i.e.Then continue iteration meter Calculate, until parameter estimation variable quantity is less than threshold epsilon, then perform step 2；

${\hat{X}}^{\left(i\right)}={(B^{T}PB+E)}^{-1}(B^{T}Pl+{\hat{X}}^{(i-1)})---\left(1\right);$

In formula, B is the factor arrays of equation, and P is power battle array, and l is constant term, and E is unit battle array, and i is iterative computation number of times；

Step 2: ask for discrete point set during iterations i in spectrum correction iterative processAnd calculate adjacent discrete point (x_i,y_iInterpolation polynomial f (x between)_i,y_i), obtain spectrum and repair Positive iterative process L-curve；

Step 3: ask for spectrum and revise the point of maximum curvature, the discrete point (x that this point of detection range is nearest on iterative process L-curve_m,y_m) Corresponding iterations m；

Step 4: by the parameter estimation corresponding to iterations mAs final estimated result.

Spectrum correcting iteration method based on L-curve the most according to claim 1, it is characterised in that: in step 3, tool Body includes

Step 3.1: according to adjacent discrete point (x in step 2_i,y_iInterpolation polynomial f (x between)_i,y_i), utilize formula k_i=| f_i″|/(1+f_i′²)^3/2Calculate in discrete point curvature k a little_i, wherein f ' is f (x_i,y_i) first derivative, f " is f(x_i,y_i) second dervative；

Step 3.2: the some k of maximum curvature on search L-curve_max, determine apart from the discrete point (x that this point is nearest_m,y_m) corresponding Iterations m.