A new biased estimation method based on Neumann series for solving ill-posed problems

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.

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+\epsilon$

1

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_{\text{LSE}}={(A^{T}A)}^{-1}{A}^T\cdot y$

2

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

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

3

As stated before, it is known that the $x_{\text{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_{\text{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_{\text{RE}}={(B+k\cdot I_n)}^{-1}\cdot z$

4

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

$x_{\text{GRE}}={(B+\text{diag}(k_1,k_2,\cdots ,k_n))}^{-1}\cdot z$

5

in which $\text{diag}(k_1,k_2,\cdots ,k_n)$ is the diagonal matrix with the diagonal elements $k_i\ge 0$ , $i=1,2,\cdots ,n$ . Apparently, equations (3) and (4) are the special cases of equation (5) with the choices $k_1=k_2=\cdots =k_n=0$ and $k_1=k_2=\cdots =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$

6

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

$\text{cond}(B+K)<\text{cond}(B)$

7

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}+\cdots ]\cdot z$

8

The sufficient condition for convergence of equation (8) is

$\frac{\Vert K\Vert }{\Vert B+K\Vert }<1$

9

where $\Vert K\Vert$ and $\Vert B+K\Vert$ 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_{\text{BE}}^{\text{1}}={(B+K)}^{-1}\cdot z$

10

where $x_{\text{BE}}^{\text{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=\text{diag}(k_1,k_2,\cdots ,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_{\text{BE}}^2={(B+K)}^{-1}\cdot [I+Q]\cdot z$

11

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

12

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

13

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$

14

$r_0=\delta \times \text{mean}(\text{diag}(B))$

15

where $\text{mean}(\text{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 $\delta =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_{\text{BE}}^i$ and the residual vector $A{x}_{\text{BE}}^i-y$ . (3) Compute the performance index for each solution using the formula as

$p^i=\sqrt{{\left(\frac{\Vert x_{\text{BE}}^i\Vert }{\text{max}(\Vert x_{\text{BE}}^i\Vert )}\right)}^2+{\left(\frac{\Vert A{x}_{\text{BE}}^i-y\Vert }{\text{max}(\Vert A{x}_{\text{BE}}^i-y\Vert )}\right)}^2}$

16

where $\Vert x_{\text{BE}}^i\Vert$ and $\Vert A{x}_{\text{BE}}^i-y\Vert$ denote the 2-norms of the solution vector $x_{\text{BE}}^i$ and the residual vector $A{x}_{\text{BE}}^i-y$ , $\text{max}(\Vert x_{\text{BE}}^i\Vert )$ denotes the maximum value in $\Vert x_{\text{BE}}^i\Vert$ when i is taken from 1 to g, and $\text{max}(\Vert A{x}_{\text{BE}}^i-y\Vert )$ denotes the maximum value in $\Vert A{x}_{\text{BE}}^i-y\Vert$ 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_{\text{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=\left[\begin{array}{ccccc}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{array}\right],y=\left[\begin{array}{c}-10.5\\ 10.45\\ 1.4\\ 12\\ 14.1\\ -0.11\\ 21.2\\ 1.5\\ 9.01\\ 12\end{array}\right]$

17

In the example, the correlation matrix $B=A^{T}A$ exhibits a huge condition number $\text{cond}(B)=1.2892\times {10}^5$ . The true value of x can be obtained if noise is not considered as $x_{\text{TRUE}}={[\begin{array}{ccccc}1& 1& 1& 1& 1\end{array}]}^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_{\Delta x}=\frac{{\Vert x_{\text{estimate}}-x_{\text{TRUE}}\Vert }_2}{{\Vert x_{\text{TRUE}}\Vert }_2}$

18

Table 1.

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

$x_{\text{TRUE}}$	$x_{\text{LSE}}$	$x_{\text{RE}}$ = $x_{\text{BE}}^1$	$x_{\text{BE}}^{11}$ ( $g=20$ )	$x_{\text{BE}}^{13}$ ( $g=100$ )	$x_{\text{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_{\Delta x}=0$	$e_{\Delta x}=1.1618$	$e_{\Delta x}=0.7870$	$e_{\Delta x}=0.7289$	$e_{\Delta x}=0.7202$	$e_{\Delta x}=0.6752$

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

Table 2.

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

$x_{\text{TRUE}}$	$x_{\text{LSE}}$	$x_{\text{RE}}$ = $x_{\text{BE}}^1$	$x_{\text{BE}}^7$ ( $g=20$ )	$x_{\text{BE}}^{11}$ ( $g=100$ )	$x_{\text{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_{\Delta x}=0$	$e_{\Delta x}=1.1618$	$e_{\Delta x}=0.7781$	$e_{\Delta x}=0.7162$	$e_{\Delta x}=0.6904$	$e_{\Delta x}=0.6770$

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

Table 1 presents the results using $\delta =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_{\text{BE}}^{11}$ and the accuracy order is $x_{\text{BE}}^{11}>x_{\text{RE}}>x_{\text{LSE}}$ . When $g=100$ and $g=1000$ , one can see from Table 1 that the optimal biased estimates are $x_{\text{BE}}^{13}$ and $x_{\text{BE}}^{32}$ , respectively. Overall, the most accurate solution is $x_{\text{BE}}^{32}$ , followed by $x_{\text{BE}}^{13}$ and $x_{\text{BE}}^{11}$ .

Table 2 presents the results using $\delta =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_{\text{BE}}^7$ and the accuracy order is $x_{\text{BE}}^7>x_{\text{RE}}>x_{\text{LSE}}$ . When $g=100$ and $g=1000$ , one can see from Table 2 that the optimal biased estimates are $x_{\text{BE}}^{11}$ and $x_{\text{BE}}^{21}$ , respectively. Overall, the most accurate solution is $x_{\text{BE}}^{21}$ , followed by $x_{\text{BE}}^{11}$ and $x_{\text{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=\left[\begin{array}{ccccc}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{array}\right],y=\left[\begin{array}{c}-10.512\\ 10.443\\ 1.4485\\ 11.94\\ 14.085\\ -0.1535\\ 21.192\\ 1.6535\\ 8.9494\\ 11.865\end{array}\right]$

19

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_{\text{BE}}^3$ ) closest to the true value of x.

Table 3.

Results obtained by different methods.

$x_{\text{TRUE}}$	TLS	TLS with L-curve	LSE	RE with L-curve	VOM	$x_{\text{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_{\Delta x}$ = 0	$e_{\Delta x}$ = 6.7322	$e_{\Delta x}$ = 0.8290	$e_{\Delta x}$ = 1.3088	$e_{\Delta x}$ = 0.8547	$e_{\Delta x}$ = 0.8231	$e_{\Delta 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_{ij})}_{n\times n},\hspace{0.17em}\hspace{0.17em}{h}_{ij}=\frac{1}{i+j-1}$

20

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$

21

where the true value of unknown vector x is $x_{\text{TRUE}}={\left[1,\hspace{0.17em}1,\hspace{0.17em}\cdots ,\hspace{0.17em}1\right]}_{20\times 1}^T$ , and the observation vector y is assumed to be equal to $H_{20}\cdot x_{\text{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_{\text{TRUE}}$	LSE	$x_{\text{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_{\Delta x}$ = 0	$e_{\Delta x}$ = 87.9268	$e_{\Delta 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_{\text{BE}}^2$ is very close to the true value $x_{\text{TRUE}}$ .

Conclusions

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. In this study, a new BE method based on Neumann series is proposed to solve the ill-posed problems more effectively. The developed theory establishes the essential relationship between BE and unbiased estimation and can unify the existing unbiased and biased estimate formulas. Two numerical examples show that the proposed biased estimate has improved accuracy over the existing robust estimates to a certain extent. The proposed method may be very helpful to solve the ill-posed least squares problems in engineering practice.