An extended L-curve method for choosing a regularization parameter in electrical resistance tomography

To deal with the ill-posedness of the inverse problem in ERT, Tikhonov regularization has been widely used [6–8]. The method replaces the original ill-posed problem by a minimization problem whose general form is

$$\begin{eqnarray}&&\min \,F(x)=\min \left(\|Ax-b\|_{2}^{2}+\lambda \|L\left(x-{{x}_{0}}\right)\|_{2}^{2}\right),\end{eqnarray} \tag{ 7 }$$

where $\lambda$ is a positive parameter called the regularization parameter, $L$ is chosen to incorporate desirable properties on the solution such as smoothness, ${{x}_{0}}$ is the estimated solution according to some prior information. Here and below $\|\centerdot \|$ denotes the Euclidean norm.

In general, choosing ${{x}_{0}}$ as zero and $L$ as an identity matrix gives a standard Tikhonov regularization form:

$$\begin{eqnarray}&&\min \,F(x)=\min \left(\|Ax-b\|_{2}^{2}+\lambda \|x\|_{2}^{2}\right).\end{eqnarray} \tag{ 8 }$$

The solution of (8) under the regularization parameter $\lambda$ can be expressed as:

$$\begin{eqnarray}&&{{x}_{\lambda}}={{\left({{A}^{T}}A+\lambda I\right)}^{-1}}{{A}^{T}}b.\end{eqnarray} \tag{ 9 }$$

The effect of $\lambda I$ makes the inverse of ${{A}^{T}}A+\lambda I$ exists [4].

We decompose the matrix $A$ by the means of singular value decomposition (SVD) to give a further survey to the solution.

$$\begin{eqnarray}&&A=U \Sigma {{V}^{T}}=\underset{i=1}{\overset{n}{\sum}}\,{{u}_{i}}{{s}_{i}}v_{i}^{T},\end{eqnarray} \tag{ 10 }$$

where $\Sigma \in {{R}^{n\times n}}$ is a diagonal matrix with the singular values, satisfying: $\Sigma =\text{diag}\left({{s}_{1}}\cdots {{s}_{n}}\right),{{s}_{1}}\geqslant {{s}_{2}}\geqslant \cdots \geqslant {{s}_{n}}\geqslant 0$ and the matrices $U\in {{R}^{m\times n}}$ and $V\in {{R}^{\text{n}\times \text{n}}}$ consist of left and right singular vectors $U=\left({{u}_{1}}\cdots {{u}_{n}}\right),\,V=\left({{v}_{1}}\cdots {{v}_{n}}\right)$.

The solution of equation (9) can now be written as

$$\begin{eqnarray}&&{{x}_{\lambda}}=\underset{i=1}{\overset{n}{\sum}}\,{{f}_{i}}\frac{u_{i}^{T}b}{{{s}_{i}}}{{v}_{i}},\end{eqnarray} \tag{ 11 }$$

where ${{f}_{i}}=\frac{s_{i}^{2}}{s_{i}^{2}+\lambda},i=1\cdots n$ is called filter factor.

No matter which regularization method is used to solve the inverse problem, regularization parameter should be chosen properly so that the method can work sucessfully. The regularization parameter controls the degree of smoothing or regularization applied to the problem. A well-known method to choose the regularization parameter is L-curve method [12–16].

The L-curve method chooses the parameter corresponding to a corner of the plot (always in log–log scale) which is constructed by the norm of regularized solution $\|{{x}_{\lambda}}\|_{2}^{{}}$ and the norm of corresponding residual $\|A{{x}_{\lambda}}-b\|_{2}^{{}}$. The corner for determining the regularization parameter traditionally is called the global corner. The norm of regularization solution and corresponding residual can be calculated as

$$\begin{eqnarray}&&\|{{x}_{\lambda}}\|_{2}^{2}=\underset{i=1}{\overset{n}{\sum}}\,{{\left(\,{{f}_{i}}\frac{u_{i}^{T}b}{{{s}_{i}}}\right)}^{2}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\|A{{x}_{\lambda}}-b\|_{2}^{2}=\underset{i=1}{\overset{n}{\sum}}\,{{\left(\left(1-{{f}_{i}}\right)s_{i}^{T}b\right)}^{2}}.\end{eqnarray} \tag{ 13 }$$

A standard L-curve with its corner is demonstrated in figure 1. However the L-curve method sometimes cannot provide a proper parameter to solve the inverse problem even using the parameter associated with the global corner correctly [8, 22–24].

Zoom In Zoom Out Reset image size

3.3.1. Simulation example of ERT.

A simulation example where the L-curve method fails to select the proper regularization parameter for ERT is given here for analysis. The example is built in COMSOL (a multi-physical software) and MATLAB. The reconstructed field is a circle region with 16 electrodes around equally (see in figure 2). The parameters of the field, electrodes and object are listed in table 1.

Table 1. Parameters of reconstructed field and electrodes.

Field		Electrodes			Object
Radius	Conductivity	Length	Thickness	Material	Radius	Location	Conductivity
5 cm	1 S m−1	10 mm	3 mm	Titanium	15 mm	(0,0.035)	2 S m−1

Zoom In Zoom Out Reset image size

The true distribution and the reconstructed image obtained from L-curve method based on Tikhonov regularization are displayed in figure 3. The result indicates that it fails to get the reasonable distribution information using a typical L-curve method selecting regularization parameter.

Zoom In Zoom Out Reset image size

3.3.2. Analysis of the example.

The convexity property of the L-curve [25] is introduced in order to check when and why the typical L-curve method breaks down for inverse problem of ERT.

According to the convexity properties of L-curve, define ${{\theta}_{1}}(\lambda )=\frac{\|A{{x}_{\lambda}}-b{{\|}_{2}}}{\|{{x}_{\lambda}}{{\|}_{2}}}$, then if and only if all $\lambda$ in a vicinity hold the condition:

$$\begin{eqnarray}&&\theta _{1}^{\prime}(\lambda )\leqslant \frac{{{\theta}_{1}}(\lambda )}{\lambda},\end{eqnarray} \tag{ 14 }$$

the L-curve is convex in the vicinity. Let $\chi (\lambda )={{\theta}_{1}}(\lambda )/\lambda$ and $C(\lambda )=\theta _{1}^{\prime}(\lambda )-\chi (\lambda )$, then the L-curve is concave or convex in the vicinity depending on whether the sign of $C(\lambda )$ is either positive or negative.

The L-curve and associated $z=\chi (\lambda )$, $z=\theta _{1}^{\prime}(\lambda )$ and $C(\lambda )$ for the example mentioned above, are depicted in figure 4. It can be seen that the sign of $C(\lambda )$ is first positive, then negative, again positive, again negative, and finally positive along the whole definition domain, which means the convexity of the L-curve is changed. Every time the curve turns to convex from concave, a convex corner will appear. Thus two convex corners appear on the curve of the example–one is the global corner traditionally used for the L-curve method to determine the regularization parameter and the other is a new corner. The situation where there exist several convex local corners may cause the failure of the L-curve method in determining a proper regularization parameter for ERT.

Zoom In Zoom Out Reset image size

The quantitative analysis on criteria used for evaluating the method is discussed. These evaluation criteria are called relative error and correlation coefficient, which can be computed as [4]

$$\begin{eqnarray}&&\text{relative} \_\text{error}=\frac{\|\hat{x}-x{{\|}_{2}}}{\|x{{\|}_{2}}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\text{corr} \_\text{coefficient}=\frac{\underset{i=1}{\overset{n}{\sum}}\,\left({{{\hat{x}}}_{i}}-\overline{{\hat{x}}}\right)\left({{x}_{i}}-\bar{x}\right)}{\sqrt{\underset{i=1}{\overset{n}{\sum}}\,{{\left({{{\hat{x}}}_{i}}-\overline{{\hat{x}}}\right)}^{2}}\underset{i=1}{\overset{n}{\sum}}\,{{\left({{x}_{i}}-\bar{x}\right)}^{2}}}},\end{eqnarray} \tag{ 16 }$$

where $x\in {{R}^{n\times 1}}$ is the true conductivity distribution of the model, $\hat{x}\in {{R}^{n\times 1}}$ is the reconstructed distribution, ${{x}_{i}}$ and ${{\hat{x}}_{i}}$ are $i\text{th}$ elements of $x$ and $\hat{x}$ respectively, $\bar{x}$ and $\overline{{\hat{x}}}$ are the mean values of $x$ and $\hat{x}$, respectively. The smaller the relative error and the larger the correlation coefficient are, the better reconstructed image can be achieved.

For the example mentioned above, the relationship between regularization parameters and corresponding relative errors in semi-log coordinate is displayed in figure 5 (left). It has been found that there is a reasonable range of parameters, which is marked in the red box, with acceptable small relative errors. The right plot of the figure 5 presents the L-curve in which the reasonable range has also been marked in the red box corresponding to the one in the left. And then an inspiring conclusion can be obtained from figure 5 that the new corner point locates in the reasonable range of parameters. This provides a new guidance for modifying the L-curve method to choose an effective parameter. Thus an extended L-curve method is proposed by determining the regularization parameter associated with either global corner or the new corner.

Zoom In Zoom Out Reset image size

3.3.3. Realization of the extended method.

Simulation results show when the typical L-curve method fails the regularization parameter corresponding to the global corner is smaller than the regularization parameter corresponding to the new corner and when the typical L-curve method works well the regularization parameter corresponding to the global corner is larger than the regularization parameter corresponding to the new corner. So in the global algorithm the larger one between the regularization parameters corresponding to the global corner and the new corner is chosen as the final parameter.

The global algorithm for the extended method of the choosing regularization parameter in ERT can be summarized as follows:

Start:

Require:$A$,$b$;

Initialize:${{x}^{0}}=0$;

Determine the regularization parameter:

Compute a series of solution norm and residual norm under different value of regularization parameters

Compute the parameters corresponding to both global corner and new corner,

Select the larger one between the two as the final parameter;

Calculate the solution using Tikhonov regularization with the parameter determined;

End.

Two strategies are put forward to determine the new corner point and further the extended method can be implemented.

3.3.3.1. Strategy using the second-order differential of L-curve.

Since there is a sudden change of the second-order differential at the corner of L-curve, a strategy based on the second-order differential of L-curve is proposed to determine the new corner.

As the L-curve is often displayed in log–log coordinate, the second-order differential of the L-curve can be computed as

$$\begin{eqnarray}&&\hat{k}=\frac{{{\text{d}}^{2}}\hat{\xi}}{\text{d}{{{\hat{\rho}}}^{2}}}=\frac{{{{\hat{\rho}}}^{\prime}}{{{\hat{\xi}}}^{\prime\prime}}-{{{\hat{\rho}}}^{\prime\prime}}{{{\hat{\xi}}}^{\prime}}}{{{\left({{{\hat{\rho}}}^{\prime}}\right)}^{3}}},\end{eqnarray} \tag{ 17 }$$

where $\hat{\rho}=\log \left(||A{{x}_{\lambda}}-b||_{2}^{2}\right)$ and $\hat{\xi}=\log \left(||{{x}_{\lambda}}||_{2}^{2}\right)$. Then the new corner point can be found according to the second largest peak of the second-order differential while the largest peak of the second-order differential refers to the global corner, as seen in figure 6 (left).

Zoom In Zoom Out Reset image size

3.3.3.2. Strategy using the curvature of L-curve.

This strategy uses the curvature of L-curve and selects the new corner at the point with the second largest peak of curvature. The curvature of the L-curve can be computed as

$$\begin{eqnarray}&&\hat{c}=\frac{{{{\hat{\rho}}}^{\prime}}{{{\hat{\xi}}}^{\prime\prime}}-{{{\hat{\rho}}}^{\prime\prime}}{{{\hat{\xi}}}^{\prime}}}{{{\left({{\left({{{\hat{\rho}}}^{\prime}}\right)}^{2}}+{{\left({{{\hat{\xi}}}^{\prime}}\right)}^{2}}\right)}^{3/2}}},\end{eqnarray} \tag{ 18 }$$

Figure 6 (right) shows the relationship between regularization parameter and curvature of L-curve.

Solving the inverse problem of ERT by the extended L-curve method is a very time-consuming process. To reduce the running time, a projection method based on Krylov subspace is added to the extended L-curve method.

The Krylov subspace is defined as:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{K}_{k}}\left({{A}^{T}}A,{{A}^{T}}b\right)\equiv \text{span}\left\{{{A}^{T}}b,\left({{A}^{T}}A\right){{A}^{T}}b,{{\left({{A}^{T}}A\right)}^{2}}{{A}^{T}}b,\cdots,{{\left({{A}^{T}}A\right)}^{k-1}}{{A}^{T}}b\right\}, \end{array}\end{eqnarray} \tag{ 19 }$$

where $k$ is the dimension of the subspace. The subspace contains important information on the inverse problem, and omits the useless information like noise. Lanczos bidiagonalization algorithm is utilized to realize the orthonormal decomposition of ${{K}_{k}}$ and achieve the orthonormal subspace ${{W}_{k}}$ [27]. The algorithm constructs a bidiagonal matrix ${{B}_{k}}\in {{R}^{(k+1)\times k}}$ and two orthogonal matrices ${{\bar{V}}_{k}}\in {{R}^{n\times k}}$ and ${{\bar{U}}_{k+1}}\in {{R}^{m\times (k+1)}}$. ${{\bar{V}}_{k}}$ and ${{\bar{U}}_{k+1}}$ form the bases of the Krylov subspace ${{K}_{k}}\left({{A}^{T}}A,{{A}^{T}}b\right)$ and ${{K}_{k+1}}\left(A{{A}^{T}},b\right)$, respectively. They hold the relationship like:

$$\begin{eqnarray}&&A{{\bar{V}}_{k}}={{\bar{U}}_{k+1}}{{B}_{k}}.\end{eqnarray} \tag{ 20 }$$

Then the objective function of Tikhonv regularization equation (8) can be reformed as:

$$\begin{eqnarray}&&x={{W}_{k}}\theta,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\min \,{{F}_{K}}(x)=\min \left(\|A{{W}_{k}}\theta -b\|_{2}^{2}+\lambda \|{{W}_{k}}\theta \|_{2}^{2}\right).\end{eqnarray} \tag{ 22 }$$

Let ${{W}_{k}}={{\bar{V}}_{k}}$, and the equation (22) can be further reformulated as

$$\begin{eqnarray}\begin{array}{*{35}{rl}} \|A{{W}_{k}}\theta -b\|_{2}^{2}+\lambda \|{{W}_{k}}\theta \|_{2}^{2}= & \|{{{\bar{U}}}_{k+1}}{{B}_{k}}\theta -b\|_{2}^{2}+\lambda \|{{V}_{k}}\theta \|_{2}^{2} \\ = & \|\bar{U}_{k+1}^{T}\left({{{\bar{U}}}_{k+1}}{{B}_{k}}\theta -b\right)\|_{2}^{2}+\lambda \|{{V}_{k}}\theta \|_{2}^{2} \\ = & \|{{B}_{k}}\theta -\bar{U}_{k+1}^{T}b\|_{2}^{2}+\lambda \|{{V}_{k}}\theta \|_{2}^{2}. \end{array}\end{eqnarray} \tag{ 23 }$$

The alternative form of the function is:

$$\begin{eqnarray}&&\min \left(\|{{B}_{k}}\theta -\bar{U}_{k+1}^{T}b\|_{2}^{2}+\lambda \|{{V}_{k}}\theta \|_{2}^{2}\right).\end{eqnarray} \tag{ 24 }$$

The extended L-curve method with Krylov subspace takes three steps to solve the inverse problem of ERT: firstly, project the inverse problem into a subspace; secondly, solve the projected problem by the extended L-method; thirdly, perform the inverse projection transformation to the solution in order to get the real one.

The L-curve with its global and new corners for the example mentioned above is displayed in figure 7, and two strategies of the second-order differential and the curvature are used. The reconstructed results are shown in figure 8, while relative error and correlation coefficient are illuminated in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be seen from figures 7 and 8 that typical the L-curve method delivers a regularization parameter associated with the sharper L-corner of the maximum curvature and yields an unsatisfactory solution corresponding to a small regularization parameter for the example of ERT. From figure 7, the parameter computed by the strategy of second-order differential is similar to the one computed by the strategy of curvature, and both of them locate just near the new corner. The results in figures 8 and 9 show that the extended L-curve method selecting regularization parameter associated with the new corner obtains an effective solution for the inverse problem and achieves a satisfactory image for the reconstructed field. The most important is that it can fix the problem in which the traditional L-curve method fails and can extend the application of the L-curve.

To further verify the effectiveness of the proposed method, more simulation models are constructed. The conductivity of the objects is 2 S m⁻¹. Figures 10 and 11 exhibit the true distributions, reconstructed images for simple models and complex models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For these models, the traditional L-curve method may fail to reconstruct image, but the extended method with any of these two strategies can reconstruct the distributions effectively. Especially, in dealing with the complex models in figure 11, the extended method achieves a very attractive performance.

Choose Model 1 and Model 8 to analyse the convexity and the corners of the L-curves. Figure 12 displays the curves $z=\chi (\lambda )$ and $z=\theta _{1}^{\prime}(\lambda )$ (left) and the L-curves (right) of the two models. The values of $z=\chi (\lambda )$ and $z=\theta _{1}^{\prime}(\lambda )$ demonstrate that the L-curve changes its convexity in the whole definition domain and this will lead to two corners appearing on L-curves (which be seen in the right). This means that in the situations where the traditional L-curve method fail to reconstruct the distribution, two corners can be found in the L-curve and the parameter corresponding to the new corner can work out a satisfactory result.

Zoom In Zoom Out Reset image size

Comparing the two strategies for computing the new corner, even though both of them can realize the image reconstruction of ERT with high quality, the strategy of curvature outperforms the strategy of differential sometimes.

These models are also used to verify the extended L-curve method with Kryolv subspace. The imaging results and the corresponding evaluation criteria are shown in figures 13 and 14, respectively. Figure 13 shows the extended L-curve method with Krylov subspace can achieve almost the same good reconstructed results as the extended L-curve method without Krylov subspace. And the conclusion is further confirmed by the relative errors demonstrated in figure 14. Comparison between the running time for these two methods (shown in figure 14) fully prove the extended L-curve method with Kryolv subspace improves the speed of the extended method and reduces the computational time. In conclusion, the extended L-curve method with Kryolv subspace can effectively reconstruct the distribution for ERT in a much shorter time without losing the accuracy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the practical use of the extended L-curve method, some static experiments are carried out with an ERT system previously developed by Tianjin University [28]. The schematic diagram of the system is shown in figure 15. The parameters of the system are listed in table 2. Four experimental models are constructed and the parameters of the models are listed in table 3.

Table 2. Parameters of the system.

Electrodes	Material	Titanium
	Shape	Square
	Size	$20\times 10\times 1$ mm
Pipeline	Material	Perspex
	Inner diameter	186 mm
	Outer diameter	200 mm
Objects	Material	Perspex
	Shape	Rod

Zoom In Zoom Out Reset image size

Experimental results are given in figure 16. The results demonstrate that the extended L-curve method can achieve effective reconstructed images. The evaluation criteria for experimental data are given in figure 17. The relative errors are small while the correlation coefficients are large, what's more, the value of both the two criteria are similar for the extended method with or without Krylov subspace. From the running time, it can be seen that the method with Krylov subspace improves the speed of the computation. In general, the extended method performs well in dealing with experimental data from practical system. And, combined with Krylov subspace, the method can not only achieve effective reconstruction, but also save the running time at the same time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size