Sage Journals: Discover world-class research

Abstract

The ill-posed least squares problems often arise in many engineering applications such as machine learning, intelligent navigation algorithms, surveying and mapping adjustment model, and linear regression model. A new biased estimation (BE) method based on Neumann series is proposed in this article to solve the ill-posed problems more effectively. Using Neumann series expansion, the unbiased estimate can be expressed as the sum of infinite items. When all the high-order items are omitted, the proposed method degenerates into the ridge estimation or generalized ridge estimation method, whereas a series of new biased estimates can be acquired by including some high-order items. Using the comparative analysis, the optimal biased estimate can be found out with less computation. The developed theory establishes the essential relationship between BE and unbiased estimation and can unify the existing unbiased and biased estimate formulas. Moreover, the proposed algorithm suits for not only ill-conditioned equations but also rank-defect equations. Numerical results show that the proposed BE method has improved accuracy over the existing robust estimation methods to a certain extent.

Keywords

Ill-posed problem biased estimate Neumann series least squares estimate ridge estimate

Introduction

Many engineering problems need to solve linear equations. Least squares estimation (LSE) is the most commonly used method to solve linear equations. It is also called unbiased estimation since it satisfies the optimal linear unbiasedness. But when the coefficient matrix of the equation system is ill-conditioned, the calculation results obtained by the LSE often have large errors or even complete distortion. This phenomenon is called the ill-posed least squares problem. As is well known, the ill-posed least squares problems often arise in many engineering applications such as machine learning, intelligent navigation algorithms, surveying and mapping adjustment model, and linear regression model.^1

–5 Without loss of generality, consider a multiple linear regression model

y = A \cdot x + ε

where y is the $m \times 1$ vector of observation data, A is the $m \times n$ coefficient matrix of full column rank, x is the $n \times 1$ vector of unknown parameters, and ε is the $m \times 1$ error vector. The LSE solution of x is

x_{LSE} = {(A^{T} A)}^{- 1} A^{T} \cdot y

Letting $B = A^{T} A$ and $z = A^{T} y$ , equation (2) reduces to

x_{LSE} = B^{- 1} \cdot z

As stated before, it is known that the $x_{LSE}$ is not the best choice if the correlation matrix $A^{T} A$ is ill-conditioned. In this case, little variation in y may lead to very large change in $x_{LSE}$ . To solve the ill-posed problem, ridge estimate (RE)⁶ and generalized ridge estimate (GRE)⁷ are proposed in the last decades to obtain more accurate estimator of x. The RE of x is

x_{RE} = {(B + k \cdot I_{n})}^{- 1} \cdot z

where I_n is the $n \times n$ identity matrix and k is the ridge parameter. The GRE of x is

x_{GRE} = {(B + diag (k_{1}, k_{2}, \dots, k_{n}))}^{- 1} \cdot z

in which $diag (k_{1}, k_{2}, \dots, k_{n})$ is the diagonal matrix with the diagonal elements $k_{i} \geq 0$ , $i = 1, 2, \dots, n$ . Apparently, equations (3) and (4) are the special cases of equation (5) with the choices $k_{1} = k_{2} = \dots = k_{n} = 0$ and $k_{1} = k_{2} = \dots = k_{n} = k$ , respectively. Note that RE and GRE destroy the equivalence relationship of the linear equations, and the solutions are biased estimations (BEs). The main problem existing in RE or GRE is how to choose a suitable value of the ridge parameter. Many methods, such as L-curve method and ridge trace method, have been proposed for choosing the ridge parameters in the last decades.^{8

–24} But so far, these processes for choosing the ridge parameters all require very complex calculations, and the ridge parameters obtained are often not the optimal solutions. It is very necessary to develop a new BE method for solving the ill-posed problems more effectively.

BE method based on Neumann series

In this section, a new BE based on Neumann series is proposed to solve the ill-posed least squares problems. Using Neumann series expansion, the unbiased estimate can be expressed as the sum of infinite items. When all the high-order items are omitted, the proposed method degenerates into the RE or GRE method, whereas a series of new biased estimates can be acquired by including some high-order items. Using the comparative analysis, the optimal biased estimate can be found out with less computation. Moreover, the proposed algorithm suits for not only ill-conditioned equations but also rank-defect equations. The main formulas of the proposed method are derived as follows.

Equation (3) can be rewritten as

x = B^{- 1} \cdot z = {(B + K - K)}^{- 1} \cdot z = {((B + K) - K)}^{- 1} \cdot z

where the $n \times n$ squares matrix K is called the regularization matrix, which is designed to reduce the condition number for matrix inversion. It must be satisfied

cond (B + K) < cond (B)

The selection of the regularization matrix K is the first important issue in the proposed method, which will be discussed in the next section. Using Neumann series,^25

–28 equation (6) can be expanded as

x = {(B + K)}^{- 1} \cdot [I + K {(B + K)}^{- 1} + K {(B + K)}^{- 1} K {(B + K)}^{- 1} + \dots] \cdot z

The sufficient condition for convergence of equation (8) is

\frac{‖ K ‖}{‖ B + K ‖} < 1

where $‖ K ‖$ and $‖ B + K ‖$ denote the norms of matrices K and $B + K$ .

Equation (8) is of great significance. It contains the existing LSE, ridge estimation, and generalized ridge estimation formulas. In other words, it establishes the essential relationship between BE and unbiased estimation and then can unify the unbiased and biased estimate formulas. Based on equation (8), a series of new BE forms can be derived. The details are as follows:

If all items are considered in equation (8), the solution by equation (8) is the unbiased estimate of x, which is equivalent to the least squares estimate.

If we ignore all the high-order items in equation (8), we have

x_{BE}^{1} = {(B + K)}^{- 1} \cdot z

where $x_{BE}^{1}$ denotes the first-order biased estimate derived from equation (8). Obviously, equations (4) and (5) are the special cases of equation (10) with the choices $K = k \cdot I_{n}$ and $K = diag (k_{1}, k_{2}, \dots, k_{n})$ , respectively. In other words, the RE and GRE are the special cases of equation (8) when the high-order items are ignored.

From equation (8), one can see that the biased estimate is essentially a part of the unbiased estimate by ignoring some high-order items. Letting $Q = K {(B + K)}^{- 1}$ , a series of new biased estimates can be derived from equation (8) as

x_{BE}^{2} = {(B + K)}^{- 1} \cdot [I + Q] \cdot z

x_{BE}^{3} = {(B + K)}^{- 1} \cdot [I + Q + Q^{2}] \cdot z

x_{BE}^{i} = {(B + K)}^{- 1} \cdot [I + Q + Q^{2} + \dots + Q^{i - 1}] \cdot z

There is an optimal biased estimate in the above series of biased estimates as shown in equations (10) to (13). Searching the optimal biased estimate is the second important issue in the proposed method, which will also be discussed in the next section.

Discussion on the special issues

As stated before, there are two important issues in the proposed method that need to be addressed to get the optimal solution more quickly. The first issue is the selection of the regularization matrix K. Technically, the regularization matrix K can be chosen arbitrarily on the premise that equation (7) is satisfied. But it is worth noting that the computation cost in searching the optimal biased estimate will depend on the selected regularization matrix. In this section, a simple formula to obtain the regularization matrix is given as

K = r_{0} \cdot I_{n}

r_{0} = δ \times mean (diag (B))

where $mean (diag (B))$ denotes the mean value of diagonal elements in matrix B and δ is a coefficient that reflects the noise level in data measurement. For example, if the data noise level is 1%, the coefficient δ can be chosen as $δ = 0.01$ .

The second issue of searching the optimal biased estimate will be addressed by a comparative analysis. The main steps are as follows. (1) Compute the series solutions of x by equation (13) when i is taken from 1 to g. The value of g can be determined by experience; for example, $g = 50$ , $g = 100$ , $g = 200$ , and so on. Generally, it can get a better solution by taking a larger value of g. But this will also lead to an increase in computational effort. An empirical way to determine the suitable value of g is to take the multiple of the dimension of matrix B. As will be shown in numerical examples, more accurate solutions can be obtained when the value of g increases. (2) For each solution, calculate the 2-norms of the solution vector $x_{BE}^{i}$ and the residual vector $A x_{BE}^{i} - y$ . (3) Compute the performance index for each solution using the formula as

p^{i} = \sqrt{{(\frac{‖ x_{BE}^{i} ‖}{max (‖ x_{BE}^{i} ‖)})}^{2} + {(\frac{‖ A x_{BE}^{i} - y ‖}{max (‖ A x_{BE}^{i} - y ‖)})}^{2}}

where $‖ x_{BE}^{i} ‖$ and $‖ A x_{BE}^{i} - y ‖$ denote the 2-norms of the solution vector $x_{BE}^{i}$ and the residual vector $A x_{BE}^{i} - y$ , $max (‖ x_{BE}^{i} ‖)$ denotes the maximum value in $‖ x_{BE}^{i} ‖$ when i is taken from 1 to g, and $max (‖ A x_{BE}^{i} - y ‖)$ denotes the maximum value in $‖ A x_{BE}^{i} - y ‖$ when i is taken from 1 to g. Note that many optimization problems can be transformed into convex optimization models. Thus, equation (16) is one of the reasonable models to find the optimal biased estimate since it is also a convex model. The following numerical examples also show that this formula is reasonable to some extent. It must be pointed out that there may be a better optimization model than equation (16), which will be further studied in the future work. (4) Find the minimum value in all the pⁱ and the corresponding $x_{BE}^{i}$ is the desired optimal biased estimate of x. The above comparative analysis will be further illustrated in the following numerical examples.

Numerical examples

An ill-posed surveying adjustment model

Consider the ill-posed surveying adjustment model²⁹ $y = A \cdot x$ with

A = [\begin{matrix} 2 & - 5 & 1 & 1 & - 9.5 \\ - 2 & 4 & 1 & - 1.05 & 8.5 \\ - 2 & 1 & 1 & - 1 & 2.4 \\ - 1 & 2.5 & 4 & - 0.5 & 7 \\ - 1 & 3.2 & 4 & - 0.5 & 8.4 \\ 1 & 1 & - 3 & 0.4 & 0.49 \\ 3 & 7 & - 3 & 1.5 & 12.7 \\ 5 & - 1 & - 2 & 2.5 & - 3 \\ 4 & 2 & - 2 & 2.01 & 3 \\ 4 & 3 & - 2 & 2 & 5 \end{matrix}], y = [\begin{matrix} - 10.5 \\ 10.45 \\ 1.4 \\ 12 \\ 14.1 \\ - 0.11 \\ 21.2 \\ 1.5 \\ 9.01 \\ 12 \end{matrix}]

In the example, the correlation matrix $B = A^{T} A$ exhibits a huge condition number $cond (B) = 1.2892 \times 10^{5}$ . The true value of x can be obtained if noise is not considered as $x_{TRUE} = {[\begin{matrix} 1 & 1 & 1 & 1 & 1 \end{matrix}]}^{T}$ .

Now we test the proposed algorithm by using the contaminated observation vector. The contaminated observation vector y^c is generated by adding a random number to each element of the exact vector y in equation (17). Without loss of generality, we assume the contaminated observation vector is

y^{c} = {[- 10.1781, 10.6404, 1.3724, 12.05, 13.6077, - 0.1122, 20.9422, 1.554, 9.3286, 12.1123]}^{T} .

Results obtained by LSE (equation (3)), RE (equation (4)), and the proposed BE (equation (13)) are given in Tables 1 and 2 for different values of δ and g. In these tables, the relative error between the estimate and true value is described by

e_{Δ x} = \frac{{‖ x_{estimate} - x_{TRUE} ‖}_{2}}{{‖ x_{TRUE} ‖}_{2}}

Table 1.

Results of LSE, RE, and BE when $δ = 0.01$ with different values of g.

$x_{TRUE}$	$x_{LSE}$	$x_{RE}$ = $x_{BE}^{1}$	$x_{BE}^{11}$ ( $g = 20$ )	$x_{BE}^{13}$ ( $g = 100$ )	$x_{BE}^{32}$ ( $g = 1000$ )
1	1.4153	1.1690	1.2414	1.2462	1.2884
1	3.1332	0.4406	0.5298	0.5493	0.7270
1	1.5281	0.7992	0.8841	0.8890	0.9334
1	0.2334	0.5883	0.5634	0.5540	0.4716
1	−0.0768	1.2610	1.2190	1.2092	1.1207
$e_{Δ x} = 0$	$e_{Δ x} = 1.1618$	$e_{Δ x} = 0.7870$	$e_{Δ x} = 0.7289$	$e_{Δ x} = 0.7202$	$e_{Δ x} = 0.6752$

LSE: least squares estimation; RE: ridge estimate; BE: biased estimation.

Table 2.

Results of LSE, RE, and BE when $δ = 0.005$ with different values of g.

$x_{TRUE}$	$x_{LSE}$	$x_{RE}$ = $x_{BE}^{1}$	$x_{BE}^{7}$ ( $g = 20$ )	$x_{BE}^{11}$ ( $g = 100$ )	$x_{BE}^{21}$ ( $g = 1000$ )
1	1.4153	1.1946	1.2485	1.2668	1.3079
1	3.1332	0.4455	0.5588	0.6348	0.8143
1	1.5281	0.8308	0.8914	0.9104	0.9552
1	0.2334	0.5951	0.5496	0.5138	0.4334
1	−0.0768	1.2598	1.2045	1.1667	1.0772
$e_{Δ x} = 0$	$e_{Δ x} = 1.1618$	$e_{Δ x} = 0.7781$	$e_{Δ x} = 0.7162$	$e_{Δ x} = 0.6904$	$e_{Δ x} = 0.6770$

LSE: least squares estimation; RE: ridge estimate; BE: biased estimation.

Table 1 presents the results using $δ = 0.01$ with different values of series number g. When $g = 20$ , one can see from Table 1 that the optimal biased estimate is $x_{BE}^{11}$ and the accuracy order is $x_{BE}^{11} > x_{RE} > x_{LSE}$ . When $g = 100$ and $g = 1000$ , one can see from Table 1 that the optimal biased estimates are $x_{BE}^{13}$ and $x_{BE}^{32}$ , respectively. Overall, the most accurate solution is $x_{BE}^{32}$ , followed by $x_{BE}^{13}$ and $x_{BE}^{11}$ .

Table 2 presents the results using $δ = 0.005$ with different values of series number g. When $g = 20$ , one can see from Table 2 that the optimal biased estimate is $x_{BE}^{7}$ and the accuracy order is $x_{BE}^{7} > x_{RE} > x_{LSE}$ . When $g = 100$ and $g = 1000$ , one can see from Table 2 that the optimal biased estimates are $x_{BE}^{11}$ and $x_{BE}^{21}$ , respectively. Overall, the most accurate solution is $x_{BE}^{21}$ , followed by $x_{BE}^{11}$ and $x_{BE}^{7}$ .

From the above results, one can conclude that (1) the proposed BE method has improved accuracy over RE and LSE to a certain extent, (2) more precise optimal solutions can be obtained with the series number g increasing, and (3) there is little difference between the optimal biased estimates if g is large enough. It is important to note that the computation cost of the proposed method will not increase significantly even if g is very large since only simple multiplication and sums of matrices are needed in the computation of series. It is apparent that the computation cost of the proposed method is far less than that of the L-curve method since the latter needs multiple inverse operations of matrices in choosing a suitable value of the ridge parameter.

To further illustrate the superiority of the proposed method, a more serious contaminated model³⁰ derived from equation (17) is used in the next discussion, whose coefficient matrix and observation vector are both contaminated as

A = [\begin{matrix} 1.8812 & - 5.1186 & 1.0129 & 1.0806 & - 9.5331 \\ - 2.2202 & 3.8944 & 1.0656 & - 1.0268 & 8.4156 \\ - 1.9014 & 1.1472 & 0.8832 & - 1.099 & 2.4498 \\ - 1.0519 & 2.5056 & 3.9539 & - 0.366 & 7.1488 \\ - 0.9673 & 3.0783 & 3.9738 & - 0.471 & 8.3454 \\ 1.0234 & 0.9959 & - 3.1213 & 0.5479 & 0.4053 \\ 3.0021 & 6.8872 & - 3.1319 & 1.6138 & 12.675 \\ 4.8996 & - 1.1349 & - 1.9069 & 2.4316 & - 2.9337 \\ 3.9053 & 1.9739 & - 1.9989 & 1.8808 & 2.9146 \\ 3.9626 & 3.0953 & - 2.0645 & 1.9927 & 4.8799 \end{matrix}], y = [\begin{matrix} - 10.512 \\ 10.443 \\ 1.4485 \\ 11.94 \\ 14.085 \\ - 0.1535 \\ 21.192 \\ 1.6535 \\ 8.9494 \\ 11.865 \end{matrix}]

In reference,³⁰ the authors presented the computation results of this contaminated model obtained by some existing methods such as LSE, RE with L-curve, total least squares estimate (TLS), TLS with L-curve, and virtual observation method. For ease of comparison, Table 3 gives the result obtained by the proposed BE method and the results in reference.³⁰ Compared with the existing methods, one can see from Table 3 that the proposed method can obtain the solution ( $x_{BE}^{3}$ ) closest to the true value of x.

Table 3.

Results obtained by different methods.

$x_{TRUE}$	TLS	TLS with L-curve	LSE	RE with L-curve	VOM	$x_{BE}^{3}$ (δ = 0.1, g = 20)
1	3.3041	1.1899	1.3944	1.2157	1.1567	1.1706
1	−2.8033	0.3903	0.1223	0.3728	0.4028	0.4040
1	0.0600	0.8122	0.7791	0.8280	0.7771	0.7903
1	3.5874	0.6109	0.2628	0.5983	0.6028	0.6153
1	2.9027	1.3060	1.4414	1.3157	1.2980	1.3016
$e_{Δ x}$ = 0	$e_{Δ x}$ = 6.7322	$e_{Δ x}$ = 0.8290	$e_{Δ x}$ = 1.3088	$e_{Δ x}$ = 0.8547	$e_{Δ x}$ = 0.8231	$e_{Δ x}$ = 0.8169

TLS: total least squares estimate; LSE: least squares estimation; VOM: virtual observation method; RE: ridge estimate; BE: biased estimation.

Hilbert ill-conditioned matrix

The Hilbert ill-conditioned matrix is often used to test the performance of the various robust estimation methods. The typical Hilbert matrix is defined as²²

H_{n} = {(h_{i j})}_{n \times n}, h_{i j} = \frac{1}{i + j - 1}

With the increase of its order n, the Hilbert matrix becomes more seriously ill-conditioned. In this example, the 20-order Hilbert matrix $H_{20}$ is used as the coefficient matrix of a linear equation as

H_{20} \cdot x = y

where the true value of unknown vector x is $x_{TRUE} = {[1, 1, \dots, 1]}_{20 \times 1}^{T}$ , and the observation vector y is assumed to be equal to $H_{20} \cdot x_{TRUE}$ . That is to say, the observation vector is assumed to be error-free. Using LSE and the proposed method, Table 4 gives the calculation results of this Hilbert ill-conditioned equation.

Table 4.

Solutions of the Hilbert ill-conditioned equation.

$x_{TRUE}$	LSE	$x_{BE}^{2}$ (δ = 0.001, g = 10)
1	1.0000	1.0001
1	1.0000	0.9982
1	1.0007	1.0035
1	1.0205	1.0023
1	0.8828	0.9990
1	−0.8125	0.9966
1	11.0000	0.9959
1	−7.7500	0.9967
1	11.5000	0.9984
1	−23.7500	1.0004
1	−4.2500	1.0023
1	5.5000	1.0038
1	22.0000	1.0047
1	12.0000	1.0049
1	−8.0000	1.0044
1	−40.0000	1.0031
1	−20.0000	1.0010
1	−60.0000	0.9982
1	18.0000	0.9948
1	3.3750	0.9908
$e_{Δ x}$ = 0	$e_{Δ x}$ = 87.9268	$e_{Δ x}$ = 0.0166

BE: biased estimation; LSE: least squares estimation.

From Table 4, one can see that the LSE solution is completely distorted even with the use of error-free data. This is due to the equation having a very serious morbidity problem since the condition number of $H_{20}$ is 1.9084 ×10¹⁸. But the proposed BE method is still very effective and reliable. The optimal solution $x_{BE}^{2}$ is very close to the true value $x_{TRUE}$ .

Conclusions

Footnotes

Acknowledgement

The author thanks the reviewers for a thorough and careful reading of the original article. Their comments are greatly appreciated and have helped to improve the quality of the article. In addition, the author would like to thank Dr Wang W, Master Bai ZC, and Master Li Na for their help in programming and text checking.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was funded by the National Natural Science Foundation of China (11202138).

ORCID iD

QW Yang

References

Eldén

. Algorithms for the regularization of ill-conditioned least squares problems. BIT 1977; 17(2): 134–145.

Chen

Billings

Luo

. Orthogonal least squares methods and their application to non-linear system identification. Int J Control 1989; 50(5): 1873–1896.

Neumaier

. Solving ill-conditioned and singular linear systems: a tutorial on regularization. SIAM Rev 1998; 40(3): 636–666.

Mercorelli

. Biorthogonal wavelet trees in the classification of embedded signal classes for intelligent sensors using machine learning applications. J Frankl Inst 2007; 344(6): 813–829.

Mercorelli

. Denoising and harmonic detection using nonorthogonal wavelet packets in industrial applications. J Syst Sci Complex 2007; 20(3): 325–343.

Hoerl

Kennard

. Ridge regression: biased estimation for nonorthogonal problems. Technometrics 1970; 12(1): 55–67.

Hemmerle

. An explicit solution for generalized ridge regression. Technometrics 1975; 17(3): 309–314.

Mead

Renaut

. A Newton root-finding algorithm for estimating the regularization parameter for solving ill-conditioned least squares problems. Inverse Probl 2008; 25(2): 025002.

El Ghaoui

Lebret

. Robust solutions to least-squares problems with uncertain data. SIAM J Matrix Anal A 1997; 18(4): 1035–1064.

10.

Hansen

. Regularization tools: a MATLAB package for analysis and solution of discrete ill-posed problems. Numer Algorithms 1994; 6(1): 1–35.

11.

Hansen

. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992; 34(4): 561–580.

12.

Marquaridt

. Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation. Technometrics 1970; 12(3): 591–612.

13.

Brezinski

Redivo-Zaglia

Rodriguez

, et al. Multi-parameter regularization techniques for ill-conditioned linear systems. Numer Math 2003; 94(2): 203–228.

14.

Ding

Shu

Yeo

, et al. Development of least-squares-based two-dimensional finite-difference schemes and their application to simulate natural convection in a cavity. Comput Fluids 2004; 33(1): 137–154.

15.

Carson

Higham

. A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems. SIAM J Sci Comput 2017; 39(6): A2834–A2856.

16.

Silva

Ribeiro

Periçaro

. A new accelerated algorithm for ill-conditioned ridge regression problems. Comput Appl Math 2017; 37: 1–18.

17.

Gholami

Gheymasi

. Regularization of geophysical ill-posed problems by iteratively re-weighted and refined least squares. Computat Geosci 2016; 20(1): 19–33.

18.

Liu

Wang

. On the weighting method for mixed least squares -total least squares problems. Numerical Linear Algebra with Applications 2017; 24(5): e2094.

19.

Prikopa

Gansterer

Wimmer

. Parallel iterative refinement linear least squares solvers based on all-reduce operations. Parallel Comput 2016; 57: 167–184.

20.

Lee

Cheon

. Estimation for the multi-way error components model with ill-conditioned panel data. J Korean Stat Soc 2017; 46(1): 28–44.

21.

Zhang

Liew

. An improved moving least-squares Ritz method for two-dimensional elasticity problems[J]. Appl Math Comput 2014; 246: 268–282.

22.

Deng

Yin

Peng

, et al. An iterative algorithm for solving ill-conditioned linear least squares problems. Geodesy and Geodynamics 2015; 6(6): 453–459.

23.

Zouzias

Freris

. Randomized extended Kaczmarz for solving least squares[J]. SIAM J Matrix Anal A 2013; 34(2): 773–793.

24.

Jun

Park

. Generation of optimal correlations by simulated annealing for ill-conditioned least-squares solution. J Nucl Sci Technol 2015; 52(5): 670–674.

25.

Suzuki

. On the convergence of Neumann series in Banach space. Math Ann 1976; 220(2): 143–146.

26.

Benaroya

Rehak

. The Neumann series/Born approximation applied to parametrically excited stochastic systems. Probabilist Eng Mech 1987; 2(2): 74–81.

27.

Cho

Kim

. A refined semi-analytic design sensitivity based on mode decomposition and Neumann series. Int J Numer Meth Eng 2005; 62(1): 19–49.

28.

Yang

. Model reduction by Neumann series expansion. Appl Math Model 2009; 33(12): 4431–4434.

29.

Ning

. Total least squares adjustment theory and its applications[M]. Beijing: China Science and Technology Press, 2011, pp. 89–102. (in Chinese).

30.