# A New Method for TSVD Regularization Truncated Parameter Selection

Authors: Zemin Wu; Shaofeng Bian; Caibing Xiang; Yude Tong
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The truncated singular value decomposition (TSVD) regularization applied in ill-posed problem is studied. Through mathematical analysis, a new method for truncated parameter selection which is applied in TSVD regularization is proposed. In the new method, all the local optimal truncated parameters are selected first by taking into account the interval estimation of the observation noises; then the optimal truncated parameter is selected from the local optimal ones. While comparing the new method with the traditional generalized cross-validation (GCV) andL curve methods, a random ill-posed matrices simulation approach is developed in order to make the comparison as statistically meaningful as possible. Simulation experiments have shown that the solutions applied with the new method have the smallest mean square errors, and the computational cost of the new algorithm is the least.


---

## Body

## 1. Introduction

Ill-posed problem is widespread in the field of geophysical survey, such as GNSS rapid positioning, precise orbit solution of spacecraft, and downward continuation of airborne gravity [1–7]. The least square estimation is not stable in the ill-posed problem and other approaches have to be utilized. Tikhonov [8] proposed the Tikhonov regularization theory, Hanke [9] utilized iteration while Walker and Birch [10] and Hoerl and Kennard [11] applied ridge regression method. Hańćkowiak [12] and Lin et al. [13] also studied this problem. Among those approaches, methods based on singular value decomposition (SVD) have drawn extensive attention, thanks to the numerical stability of their solutions [7]. On that basis, Wiggins [14] and Xu [15] intensively studied the truncated singular value decomposition (TSVD) regularization method.The standard form of the observation equation is (nonstandard form can be transformed to a standard one) as follows:(1)L=AX+ε,
whereAis the m×n coefficient matrix, m>n,  rank(A)=n,X is the n×1 parameter vector, L is the m×1 observation vector, εis them×1residual vector,E(ε)=0, and E(εεT)=σ02I.The condition number of normal matrix (N=ATA) is generally used as the measurement of the ill-posed severity as follows:
(2)cond(N)=σ12σn2,
where σ1 and σn are the maximum and minimum singular values of A, respectively.According to application experiences, iflg(cond)>4, the system is considered as ill-posed [16].In ill-posed problem, some singular values of the coefficient matrix approximate to 0, the least square estimation will enormously amplify the observation noises and degrade the precision. In TSVD regularization, items containing these small singular values are discarded to maintain the stability of the solution. Setting a small positive numberη as the threshold value, items containing singular values smaller than η are discarded; thus, the amplification to the observation noises is restrained. If there are k singular values bigger than η, the TSVD solution is given as
(3)X^k=∑i=1kviuiTσiL.Normally, the mean square error (MSE) is used to evaluate the quality of ill-posed problem solution. The key point of the TSVD is how to select a proper truncated parameter to get the smallest mean square error of the solution. Many scholars have extensively studied this problem. Golub et al. [17] proposed GCV method, Hansen [18, 19] and Hansen and O’Leary [20] made use of the property of the L curve to select truncated parameter, and these methods will be introduced in Section 2. A new truncated parameter selecting method will be proposed in Section 3. In Section 4, a random ill-posed matrices simulation approach for simulation experiment will be developed, and the performances of GCV, L curve, and the new methods will be compared Section 2.


## 2. Traditional Truncated Parameter Selecting Method

### 2.1. GCV Method

Based on statistical point of view, Golub et al. [17] proposed GCV method for truncated parameter selection. In linear system AX=L, assuming that A# is the matrix mapping observation vector L to regularized solution X^k, A#L=X^k. GCV method is equal to selecting a truncated parameter k to minimize the GCV function, which is defined as
(4)GCV(k)=∥(I-AA#)L∥22[trace(I-AA#)]2,
where “trace(·)” denotes the trace of a matrix.The disadvantage of GCV method is that the GCV function converges very slowly in some cases. It may lead to the GCV function minimization whilek is equal to n, in which situation, the truncated parameter cannot be selected, and the GCV method fails.


### 2.2.L Curve Method

Hansen [18, 19] extensively studied the application of L curve method in regularization. With the variation of k, the discreet points (ln(∥X^k∥2) and ln(∥L-AX^k∥2))(k=1,2…,n) can be fitted with a cubic spline curve ω=Ω(ρ). He defined the maximum curvature point of this curve by (5) [18] as follows:
(5)maxκ=ρ′ω′′-ρ′′ω′[(ρ′)2+(ω′)2]3/2,
where ρ and ω are the independent and dependent variables of the cubic spline function, respectively.The original discreet point which lies closest to the maximum curvature point from the left side is selected. Its truncated serial numberk is the truncated parameter asked for.Further descriptions of theL curve method can be seen in Hansen [18, 19]. Hansen [19] indicated that the truncated parameter selected by the L curve method is not the optimum parameter but the asymptotically optimum one.


## 2.1. GCV Method

Based on statistical point of view, Golub et al. [17] proposed GCV method for truncated parameter selection. In linear system AX=L, assuming that A# is the matrix mapping observation vector L to regularized solution X^k, A#L=X^k. GCV method is equal to selecting a truncated parameter k to minimize the GCV function, which is defined as
(4)GCV(k)=∥(I-AA#)L∥22[trace(I-AA#)]2,
where “trace(·)” denotes the trace of a matrix.The disadvantage of GCV method is that the GCV function converges very slowly in some cases. It may lead to the GCV function minimization whilek is equal to n, in which situation, the truncated parameter cannot be selected, and the GCV method fails.


## 2.2.L Curve Method

Hansen [18, 19] extensively studied the application of L curve method in regularization. With the variation of k, the discreet points (ln(∥X^k∥2) and ln(∥L-AX^k∥2))(k=1,2…,n) can be fitted with a cubic spline curve ω=Ω(ρ). He defined the maximum curvature point of this curve by (5) [18] as follows:
(5)maxκ=ρ′ω′′-ρ′′ω′[(ρ′)2+(ω′)2]3/2,
where ρ and ω are the independent and dependent variables of the cubic spline function, respectively.The original discreet point which lies closest to the maximum curvature point from the left side is selected. Its truncated serial numberk is the truncated parameter asked for.Further descriptions of theL curve method can be seen in Hansen [18, 19]. Hansen [19] indicated that the truncated parameter selected by the L curve method is not the optimum parameter but the asymptotically optimum one.


## 3. A New Truncated Parameter Selecting Method

In this section, we develop a new approach to selecting the truncated parameter. We focus on the MSE of the solution, making some proper estimation, and finally derive the conditions that the optimal truncated parameter should meet. In this section and the subsequent sections,k is called the optimal value if MSE(X^k) is the smallest in all the MSE(X^i)(i=1,2…,n); k is called the local optimal value if MSE(X^k) is the smallest just in MSE(X^k-1), MSE(X^k), and MSE(X^k+1). Before the optimal value selection, we first select all the local optimal values as the candidates.


### 3.1. Local Optimal Truncated Parameter Selection

The MSE of TSVD solution is(6)MSE(X^k)=E[(X^k-X)T(X^k-X)].Taking the singular values decomposition into consideration, the least square solution of the observation equation can be written as(7)X^LS=∑i=1nviuiTσiL=X+∑i=1nviuiTσiε.According to (1) and (3), we rewrite the TSVD solution as
(8)X^k=∑i=1kviviTX+∑i=1kviuiTσiε.Inserting (8) into (6) yields
(9)MSE(X^k)=∑i=k+1nviTXXTvi+∑i=1kσ02σi2.As mentioned above, ifk is the local optimal value, the following inequalities are satisfied:
(10)MSE(X^k)-MSE(X^k-1)=σ02σk2-vkTXXTvk<0,MSE(X^k+1)-MSE(X^k)=σ02σk+12-vk+1TXXTvk+1>0.DenotingY=VTX=[y1,y2…,yn]T, inequalities (10) can be rewritten as
(11)σ02<σk2yk2,σ02>σk+12yk+12.Asσ02 and Y are unknown, we should estimate them first. Premultiplying matrix VT to both sides of (7), we obtain a new vector denoted asY^LS as follows:
(12)Y^LS=[u1TLσ1,u2TLσ2,…,unTLσn]T,(13)σi2y^iLS2=(uiTL)2,
where y^iLS is the ith element of Y^LS.Inserting (1) into (13) yields
(14)E[(uiTL)2]=uiTE(εεT)ui+σi2(viTX)2=σ02+σi2yi2.Equation (14) is an estimator of yi; taking (14) into account, inequalities (11) can be approximately written as
(15)2σ02<(ukTL)2,2σ02>(uk+1TL)2.Then, we estimateσ02; an unbiased estimator of σ02 can be obtained from least-square estimation as
(16)σ^0LS2=(AX^LS-L)T(AX^LS-L)(m-n).We have to treatσ^0LS2 cautiously, because it is just an approximate value of σ02. Replacing σ02 with σ^0LS2 may cause some misjudgments. For example, in the situation shown in Figure 1, k is the local optimal value since f(k+1)<σ02<f(k) (denoting f(i)=(uiTL)2-σ02). But, if we replace σ02 with σ^0LS2, the opposite judgment will be made since σ^0LS2<f(k+1)<f(k).Figure 1
Scheme ofσ02 estimation.In order to overcome this shortcoming, the interval estimation ofσ02 is used instead of σ^0LS2 in this algorithm. Assuming that [t1,t2] is the confidence interval of σ02 under the confidence level 1-α, we replace σ02 with t1 and t2 in the first and second inequality in inequalities (15), respectively. By this means, a fault-tolerance interval between the judgment thresholds is given. As seen in Figure 1, by properly selecting t1 and t2, there are t1<f(k+1)<f(k)<t2, and so, k is judged as the local optimal value.To get a proper interval estimation ofσ02, analysis goes as follows.The estimation of the observation noises is(17)ε^=[ε^1,ε^2…,ε^m]T=AX^LS-L,
where ε^i~N(0,σ02) (note that among the m elements, onlym-nof them are independent).Thus, according to the definition ofχ2 distribution, ∑i=1mε^i2/σ02~χ2(m-n). Since [t1,t2] is the confidence interval ofσ02under the confidence level 1-α, we yield
(18)P{t1<σ02<t2}=1-α.Considering theχ2 distribution, (18) can be transformed as follows:
(19)P{ε^Tε^t2<χ2<ε^Tε^t1}=1-α.SelectingP{χ2≤ε^Tε^/t2}=P{χ2≥ε^Tε^/t1}=α/2 yields
(20)t1=ε^Tε^χ1-α/22(m-n),t2=ε^Tε^χα/22(m-n),Equation (20) is the estimator of t1 and t2, the final local optimal truncated parameter selecting inequalities reads as follows:
(21)2ε^Tε^χ1-α/22(m-n)<(ukTL)2,2ε^Tε^χα/22(m-n)>(uk+1TL)2.There is an undetermined parameter in inequalities (21): the confidence level is 1-α. This parameter adjusts the size of the fault-tolerance interval. Enlarging α will decrease the fault-tolerance interval and the risk of leaking local optimal value increases. While decreasing α will increase the fault-tolerance interval, it will cause the judgment efficiency to decline. In order to balance the fault-tolerance and the judgment efficiency of inequalities (21), we first set α=0.5 in this algorithm. If no truncated parameter matches inequalities (21), the value of α will be decreased until one or more truncated parameters can match it.


### 3.2. Optimal Truncated Parameter Selection

If there is only one local optimal value, we regard it as the optimal value. If not, assuming that there arer local optimal values k1,k2…,kr, It is essential to select the optimal value from them. Analysis goes as follows.The mean square error function can be transformed as(22)MSE(X^k)=tr(cov(X^k))+∥E(X^k)-X∥22,
where “cov(·)” denotes the covariance matrix. Inserting (1) and (3) into (22), the first and second items of which can be simplified as follows:
(23)tr(cov(X^k))=E(X^kTX^k)+E(X^kT)E(X^k)=σ02∑i=1k1σi2,∥E(X^k)-X∥22=∥∑i=1kviviTX-X∥22=∑i=k+1n(XTvi)2.Assuming thatkp and kq are both local optimal values, kp<kq, MSE(X^kp) minus MSE(X^kq) is denoted as Δ(kp,kq), we derive
(24)Δ(kp,kq)=∑i=kp+1kq(XTvi)2-σ02∑i=kp+1kq1σi2.According to inequalities (23), to determine the sign of Δ(kp,kq), we take proper estimations of ∑i=kp+1kq(XTvi)2 and σ02∑i=kp+1kq1/σi2 simultaneously.If we regardX as an unknown random vector, the elements of which are independent and have the same statistics weight. That is, E(XXT)=σx2I. Then, it yields
(25)∑i=kp+1kq(XTvi)2≈σx2(kq-kp).Askp and kq are both local optimal values, the regularization solutions X^kp and X^kq are much more stable than the least square solution. We estimate σx2 by X^kp and σ02 by its unbiased estimator σ^0LS2, which yields
(26)Δ(kp,kq)≈X^kpTX^kp(kq-kp)n-σ^0LS2∑i=kp+1kq1/σi2.According to (26), if X^kpTX^kp(kq-kp)/n-σ^0LS2∑i=kp+1kq1/σi2<0, we regard kp as the better truncated parameter and vice versa. After comparing all the local optimal values k1,k2,…,kr utilizing (26), we can get the minimized MSE of all the solutions, and the optimal truncated parameter is determined.


### 3.3. Algorithm Flow

In Sections3.1 and 3.2, we have discussed the conditions that the optimal truncated parameter should meet. And then, the equations to select the local optimal and the optimal truncated parameter are deduced in detail. Although the derivation process is somehow complicated and trivial, the execution of this algorithm is quite simple. In this subsection, we summarize the steps of the algorithm flow to make it clear and more readable.BeginStep1. Set the confidence level of interval estimation 1-α=0.5.Step2. Calculate the two thresholds t1 and t2 using  (20).Step  3. For each truncated parameters k from 1 to n-1, calculate inequalities (21). If the inequalities are satisfied, select the corresponding truncated parameter as a local optimal truncated parameter.Step4. If no truncated parameter matches inequalities (21), set α=α-0.1, return to Step2;else, skip this step.Step5. Estimate the variance of observation noise σ02 by (16).Step  6. If there is only one local optimal parameter, it is the optimal truncated parameter simultaneously;else, make pairwise comparisons of the local optimal truncated parameters according to (26), until getting the one minimizing the MSE as the optimal truncated parameter.End


## 3.1. Local Optimal Truncated Parameter Selection

The MSE of TSVD solution is(6)MSE(X^k)=E[(X^k-X)T(X^k-X)].Taking the singular values decomposition into consideration, the least square solution of the observation equation can be written as(7)X^LS=∑i=1nviuiTσiL=X+∑i=1nviuiTσiε.According to (1) and (3), we rewrite the TSVD solution as
(8)X^k=∑i=1kviviTX+∑i=1kviuiTσiε.Inserting (8) into (6) yields
(9)MSE(X^k)=∑i=k+1nviTXXTvi+∑i=1kσ02σi2.As mentioned above, ifk is the local optimal value, the following inequalities are satisfied:
(10)MSE(X^k)-MSE(X^k-1)=σ02σk2-vkTXXTvk<0,MSE(X^k+1)-MSE(X^k)=σ02σk+12-vk+1TXXTvk+1>0.DenotingY=VTX=[y1,y2…,yn]T, inequalities (10) can be rewritten as
(11)σ02<σk2yk2,σ02>σk+12yk+12.Asσ02 and Y are unknown, we should estimate them first. Premultiplying matrix VT to both sides of (7), we obtain a new vector denoted asY^LS as follows:
(12)Y^LS=[u1TLσ1,u2TLσ2,…,unTLσn]T,(13)σi2y^iLS2=(uiTL)2,
where y^iLS is the ith element of Y^LS.Inserting (1) into (13) yields
(14)E[(uiTL)2]=uiTE(εεT)ui+σi2(viTX)2=σ02+σi2yi2.Equation (14) is an estimator of yi; taking (14) into account, inequalities (11) can be approximately written as
(15)2σ02<(ukTL)2,2σ02>(uk+1TL)2.Then, we estimateσ02; an unbiased estimator of σ02 can be obtained from least-square estimation as
(16)σ^0LS2=(AX^LS-L)T(AX^LS-L)(m-n).We have to treatσ^0LS2 cautiously, because it is just an approximate value of σ02. Replacing σ02 with σ^0LS2 may cause some misjudgments. For example, in the situation shown in Figure 1, k is the local optimal value since f(k+1)<σ02<f(k) (denoting f(i)=(uiTL)2-σ02). But, if we replace σ02 with σ^0LS2, the opposite judgment will be made since σ^0LS2<f(k+1)<f(k).Figure 1
Scheme ofσ02 estimation.In order to overcome this shortcoming, the interval estimation ofσ02 is used instead of σ^0LS2 in this algorithm. Assuming that [t1,t2] is the confidence interval of σ02 under the confidence level 1-α, we replace σ02 with t1 and t2 in the first and second inequality in inequalities (15), respectively. By this means, a fault-tolerance interval between the judgment thresholds is given. As seen in Figure 1, by properly selecting t1 and t2, there are t1<f(k+1)<f(k)<t2, and so, k is judged as the local optimal value.To get a proper interval estimation ofσ02, analysis goes as follows.The estimation of the observation noises is(17)ε^=[ε^1,ε^2…,ε^m]T=AX^LS-L,
where ε^i~N(0,σ02) (note that among the m elements, onlym-nof them are independent).Thus, according to the definition ofχ2 distribution, ∑i=1mε^i2/σ02~χ2(m-n). Since [t1,t2] is the confidence interval ofσ02under the confidence level 1-α, we yield
(18)P{t1<σ02<t2}=1-α.Considering theχ2 distribution, (18) can be transformed as follows:
(19)P{ε^Tε^t2<χ2<ε^Tε^t1}=1-α.SelectingP{χ2≤ε^Tε^/t2}=P{χ2≥ε^Tε^/t1}=α/2 yields
(20)t1=ε^Tε^χ1-α/22(m-n),t2=ε^Tε^χα/22(m-n),Equation (20) is the estimator of t1 and t2, the final local optimal truncated parameter selecting inequalities reads as follows:
(21)2ε^Tε^χ1-α/22(m-n)<(ukTL)2,2ε^Tε^χα/22(m-n)>(uk+1TL)2.There is an undetermined parameter in inequalities (21): the confidence level is 1-α. This parameter adjusts the size of the fault-tolerance interval. Enlarging α will decrease the fault-tolerance interval and the risk of leaking local optimal value increases. While decreasing α will increase the fault-tolerance interval, it will cause the judgment efficiency to decline. In order to balance the fault-tolerance and the judgment efficiency of inequalities (21), we first set α=0.5 in this algorithm. If no truncated parameter matches inequalities (21), the value of α will be decreased until one or more truncated parameters can match it.


## 3.2. Optimal Truncated Parameter Selection

If there is only one local optimal value, we regard it as the optimal value. If not, assuming that there arer local optimal values k1,k2…,kr, It is essential to select the optimal value from them. Analysis goes as follows.The mean square error function can be transformed as(22)MSE(X^k)=tr(cov(X^k))+∥E(X^k)-X∥22,
where “cov(·)” denotes the covariance matrix. Inserting (1) and (3) into (22), the first and second items of which can be simplified as follows:
(23)tr(cov(X^k))=E(X^kTX^k)+E(X^kT)E(X^k)=σ02∑i=1k1σi2,∥E(X^k)-X∥22=∥∑i=1kviviTX-X∥22=∑i=k+1n(XTvi)2.Assuming thatkp and kq are both local optimal values, kp<kq, MSE(X^kp) minus MSE(X^kq) is denoted as Δ(kp,kq), we derive
(24)Δ(kp,kq)=∑i=kp+1kq(XTvi)2-σ02∑i=kp+1kq1σi2.According to inequalities (23), to determine the sign of Δ(kp,kq), we take proper estimations of ∑i=kp+1kq(XTvi)2 and σ02∑i=kp+1kq1/σi2 simultaneously.If we regardX as an unknown random vector, the elements of which are independent and have the same statistics weight. That is, E(XXT)=σx2I. Then, it yields
(25)∑i=kp+1kq(XTvi)2≈σx2(kq-kp).Askp and kq are both local optimal values, the regularization solutions X^kp and X^kq are much more stable than the least square solution. We estimate σx2 by X^kp and σ02 by its unbiased estimator σ^0LS2, which yields
(26)Δ(kp,kq)≈X^kpTX^kp(kq-kp)n-σ^0LS2∑i=kp+1kq1/σi2.According to (26), if X^kpTX^kp(kq-kp)/n-σ^0LS2∑i=kp+1kq1/σi2<0, we regard kp as the better truncated parameter and vice versa. After comparing all the local optimal values k1,k2,…,kr utilizing (26), we can get the minimized MSE of all the solutions, and the optimal truncated parameter is determined.


## 3.3. Algorithm Flow

In Sections3.1 and 3.2, we have discussed the conditions that the optimal truncated parameter should meet. And then, the equations to select the local optimal and the optimal truncated parameter are deduced in detail. Although the derivation process is somehow complicated and trivial, the execution of this algorithm is quite simple. In this subsection, we summarize the steps of the algorithm flow to make it clear and more readable.BeginStep1. Set the confidence level of interval estimation 1-α=0.5.Step2. Calculate the two thresholds t1 and t2 using  (20).Step  3. For each truncated parameters k from 1 to n-1, calculate inequalities (21). If the inequalities are satisfied, select the corresponding truncated parameter as a local optimal truncated parameter.Step4. If no truncated parameter matches inequalities (21), set α=α-0.1, return to Step2;else, skip this step.Step5. Estimate the variance of observation noise σ02 by (16).Step  6. If there is only one local optimal parameter, it is the optimal truncated parameter simultaneously;else, make pairwise comparisons of the local optimal truncated parameters according to (26), until getting the one minimizing the MSE as the optimal truncated parameter.End


## 4. Random Simulation Experiment

### 4.1. The Approach of Random Simulation

Although we can use a particular coefficient matrixA to compare different truncated parameter selecting methods, such comparisons are of limited value [21]. Because the doubt to the conclusions drawn from them will be threefold: (i) whether the conclusions are still valid if the dimensions of the coefficient matrix have changed; (ii) whether the conclusions are still valid if the condition number and singular values distribution of the coefficient matrix have changed; (iii) whether the conclusions are valid if they are applied to another coefficient matrix which has the same dimensions and condition number with A.In order to make the comparison of the three methods as persuasive as possible, we have to simulate the coefficient matrices by taking into account the three issues mentioned above. According to (2), A=UΛVT; thus, the simulation problem has turned into how to design the matrices U, Λ, and V and their dimensions.Theoretically, the dimensionsm and n can be set any arbitrary positive integer which satisfy: m>n>1. Considering the computing burden, we set 5≤n≤25 and n<m≤40, each of them is randomly generated in their interval. As it is mentioned in Section 1, a coefficient matrix is regarded as ill-posed if lg(cond)>4. But the condition number cannot be infinite; otherwise, the ill-posed matrix will be turned into rank defect matrix. In the simulation, the mathematical expectation of lg(cond) is uniformly distributed over[4,7]. Furthermore, about the singular values, we assume that the mathematical expectation of the sum of the logarithms of all the singular values is 0. Finally, the equations set can be written as follows:
(27)E(lgσ12-lgσn2)=ψ,lgσi2-lgσi+12=d+ξ,E(∑i=1nlgσi2)=0,
where the random variables ψ and ξ conform to uniform distribution, ψ~R[4,7]; d=ψ/(n-1) is also a random variable, ξ~R[-d,d].Utilizing the Gram-Schmidt method, unitary matrixU (or V) can be derived from an m-dimension (or n-dimension) random full rank matrix. For instance, assume that B is an m-dimension random full rank matrix and b1,b2,…,bm are column vectors of B; thus, the column vectors of U are as follows:
(28)u1=b1|b1|,uk=bk-∑i=1k-1(bk,ui)ui|bk-∑i=1k-1(bk,ui)ui|,(k=2,…,m).To make the simulation meaningful, we sample 100 random generated coefficient matrices for our simulation experiment. The dimensions and condition numbers of these matrices are shown in Figures2(a) and 2(b), respectively.Dimensions and condition numbers of the 100 coefficient matrices randomly generated, respectively.
(a)
(b)


### 4.2. Comparisons of Truncated Parameter Selecting Methods

In this section, we will compare the performances of TSVD solutions applied with the GCV,L curve, and the new methods proposed in this paper. The criterion is the mean square errors of the solutions. In the simulation, we set X as n-dimension vector, all elements of which are 1. The elements of residual vector ε are Gaussian white noises, generated by randomizer, ε~N(0,σ02I). To enrich our comparison, two observation noise levels are used in our simulation, the variances of which are LTL/(100m) and LTL/(1000m). In consideration of the randomness of ε, objective comparisons cannot be obtained if we only do the experiment for once. In our simulation, 1000 times independent repeated experiments have been done for each coefficient matrix and each observation level. The mean square errors which are used in our comparisons are the mean value of those in the 1000 experiments.(i) Results Comparison between GCV and the New Methods. The comparisons between GCV and the new methods under the observation noise levels σ02=LTL/(1000m) and σ02=LTL/(100m) are shown in Figures 3(a) and 3(b), respectively. In order to gain further details of the comparisons, we draw the D-value curves of their mean square errors. The D-value curves of Figures 3(a) and 3(b) are shown in Figures 3(c) and 3(d), respectively; the negative value means that the new method has smaller mean square errors. Since the GCV method may fail in some cases, the solutions shown in Figures 3(a) to 3(d) do not contain those failure cases.Results comparison between GCV and the new methods in two observation noise levels.
(a)
(b)
(c)
(d)Seen in Figures3(a) and 3(b), the solid lines are always beneath the dotted lines. In Figures 3(c) and 3(d), most values are negative. That means that the solutions of TSVD regularization applied with the new method have smaller mean square errors.(ii) Results Comparison between
L
Curve and the New Methods. Similar to previous comparisons, the comparisons between L curve and the new methods under two observation noise levels are shown in Figures 4(a) and 4(b), respectively. The D-value curves of Figures 4(a) and 4(b) are shown in Figures 4(c) and 4(d), respectively; the negative value means that the new method has smaller mean square errors.Results comparison betweenL curve and the new methods in two observation noise levels.
(a)
(b)
(c)
(d)As seen in Figures4(a) and 4(b), the solid lines are also always beneath the dotted lines; and in Figures 4(c) and 4(d), most values are negative. So we can verdict that the new method has a better performance in truncated parameter selection than L curve method.(iii) Time Consuming Comparisons. Besides the quality of results, computational cost is another important aspect to evaluate the effectiveness of an algorithm. In the simulations, we record the time consuming of the algorithms for each coefficient matrix. The average time consuming of these three algorithms are shown in Figure 5(a). (Note that these time consuming values are the average values of 1000 times independent repeated experiments.) Figure 5(b) has shown the time consuming of the new algorithm under two different noise levels. The observation noise levels σ02=LTL/(1000m) and σ02=LTL/(100m) are, respectively, labeled “Noise level 1” and “Noise level 2” in Figure 5(b).(a) time consuming comparisons of the three methods; (b) time consuming of the new algorithm in two noise levels.
(a)
(b)From Figures5(a) and 5(b), we can see that the computational costs of these algorithms are all small. Nevertheless, the cost of the new method is obviously smaller than that of the other two by dozens of times. With the increasing of the matrices’ dimension or the count of times, this advantage will be more and more obvious. Another conclusion we can get from Figure 5(b) is that the observation noises only have a little influences with the computational costs.


## 4.1. The Approach of Random Simulation

Although we can use a particular coefficient matrixA to compare different truncated parameter selecting methods, such comparisons are of limited value [21]. Because the doubt to the conclusions drawn from them will be threefold: (i) whether the conclusions are still valid if the dimensions of the coefficient matrix have changed; (ii) whether the conclusions are still valid if the condition number and singular values distribution of the coefficient matrix have changed; (iii) whether the conclusions are valid if they are applied to another coefficient matrix which has the same dimensions and condition number with A.In order to make the comparison of the three methods as persuasive as possible, we have to simulate the coefficient matrices by taking into account the three issues mentioned above. According to (2), A=UΛVT; thus, the simulation problem has turned into how to design the matrices U, Λ, and V and their dimensions.Theoretically, the dimensionsm and n can be set any arbitrary positive integer which satisfy: m>n>1. Considering the computing burden, we set 5≤n≤25 and n<m≤40, each of them is randomly generated in their interval. As it is mentioned in Section 1, a coefficient matrix is regarded as ill-posed if lg(cond)>4. But the condition number cannot be infinite; otherwise, the ill-posed matrix will be turned into rank defect matrix. In the simulation, the mathematical expectation of lg(cond) is uniformly distributed over[4,7]. Furthermore, about the singular values, we assume that the mathematical expectation of the sum of the logarithms of all the singular values is 0. Finally, the equations set can be written as follows:
(27)E(lgσ12-lgσn2)=ψ,lgσi2-lgσi+12=d+ξ,E(∑i=1nlgσi2)=0,
where the random variables ψ and ξ conform to uniform distribution, ψ~R[4,7]; d=ψ/(n-1) is also a random variable, ξ~R[-d,d].Utilizing the Gram-Schmidt method, unitary matrixU (or V) can be derived from an m-dimension (or n-dimension) random full rank matrix. For instance, assume that B is an m-dimension random full rank matrix and b1,b2,…,bm are column vectors of B; thus, the column vectors of U are as follows:
(28)u1=b1|b1|,uk=bk-∑i=1k-1(bk,ui)ui|bk-∑i=1k-1(bk,ui)ui|,(k=2,…,m).To make the simulation meaningful, we sample 100 random generated coefficient matrices for our simulation experiment. The dimensions and condition numbers of these matrices are shown in Figures2(a) and 2(b), respectively.Dimensions and condition numbers of the 100 coefficient matrices randomly generated, respectively.
(a)
(b)


## 4.2. Comparisons of Truncated Parameter Selecting Methods

In this section, we will compare the performances of TSVD solutions applied with the GCV,L curve, and the new methods proposed in this paper. The criterion is the mean square errors of the solutions. In the simulation, we set X as n-dimension vector, all elements of which are 1. The elements of residual vector ε are Gaussian white noises, generated by randomizer, ε~N(0,σ02I). To enrich our comparison, two observation noise levels are used in our simulation, the variances of which are LTL/(100m) and LTL/(1000m). In consideration of the randomness of ε, objective comparisons cannot be obtained if we only do the experiment for once. In our simulation, 1000 times independent repeated experiments have been done for each coefficient matrix and each observation level. The mean square errors which are used in our comparisons are the mean value of those in the 1000 experiments.(i) Results Comparison between GCV and the New Methods. The comparisons between GCV and the new methods under the observation noise levels σ02=LTL/(1000m) and σ02=LTL/(100m) are shown in Figures 3(a) and 3(b), respectively. In order to gain further details of the comparisons, we draw the D-value curves of their mean square errors. The D-value curves of Figures 3(a) and 3(b) are shown in Figures 3(c) and 3(d), respectively; the negative value means that the new method has smaller mean square errors. Since the GCV method may fail in some cases, the solutions shown in Figures 3(a) to 3(d) do not contain those failure cases.Results comparison between GCV and the new methods in two observation noise levels.
(a)
(b)
(c)
(d)Seen in Figures3(a) and 3(b), the solid lines are always beneath the dotted lines. In Figures 3(c) and 3(d), most values are negative. That means that the solutions of TSVD regularization applied with the new method have smaller mean square errors.(ii) Results Comparison between
L
Curve and the New Methods. Similar to previous comparisons, the comparisons between L curve and the new methods under two observation noise levels are shown in Figures 4(a) and 4(b), respectively. The D-value curves of Figures 4(a) and 4(b) are shown in Figures 4(c) and 4(d), respectively; the negative value means that the new method has smaller mean square errors.Results comparison betweenL curve and the new methods in two observation noise levels.
(a)
(b)
(c)
(d)As seen in Figures4(a) and 4(b), the solid lines are also always beneath the dotted lines; and in Figures 4(c) and 4(d), most values are negative. So we can verdict that the new method has a better performance in truncated parameter selection than L curve method.(iii) Time Consuming Comparisons. Besides the quality of results, computational cost is another important aspect to evaluate the effectiveness of an algorithm. In the simulations, we record the time consuming of the algorithms for each coefficient matrix. The average time consuming of these three algorithms are shown in Figure 5(a). (Note that these time consuming values are the average values of 1000 times independent repeated experiments.) Figure 5(b) has shown the time consuming of the new algorithm under two different noise levels. The observation noise levels σ02=LTL/(1000m) and σ02=LTL/(100m) are, respectively, labeled “Noise level 1” and “Noise level 2” in Figure 5(b).(a) time consuming comparisons of the three methods; (b) time consuming of the new algorithm in two noise levels.
(a)
(b)From Figures5(a) and 5(b), we can see that the computational costs of these algorithms are all small. Nevertheless, the cost of the new method is obviously smaller than that of the other two by dozens of times. With the increasing of the matrices’ dimension or the count of times, this advantage will be more and more obvious. Another conclusion we can get from Figure 5(b) is that the observation noises only have a little influences with the computational costs.


## 5. Conclusions

A new truncated parameter selecting method is proposed in this paper. Focusing on how to decrease the MSE of the solution, we divide the truncated parameter selection into two parts: local optimal and optimal truncated parameter selection. Detailed analyses are made and finally the equations to select the local optimal and the optimal truncated parameter are obtained. Although the derivation of the new algorithm is somehow trivial, its execution is quite simple and efficient.The key point of local optimal truncated parameter selection is the estimation of the observation noises. To overcome the inaccuracy of point estimation, we apply the interval estimation with the confidence level1-α. The parameter α balances the fault-tolerance and the judgment efficiency of local truncated parameter selection. In this algorithm, α is firstly set as 0.5, it can adjust automatically according to the situation of local truncated parameter selection.The selection of optimal truncated parameter is based on the local optimal truncated parameters. We regard the parameter vector as a random variant vector, making use of the stability of its regularization solutions, and finally derive the equations of optimal truncated parameter selection.A random ill-posed matrices simulation approach is developed and a great deal of experiments is made in Section4. The results have shown that the solutions applied with the new method have the smallest mean square errors. Meanwhile, the computational cost of the new method is the least of the three.



---

*Source: 161834-2013-11-17.xml*