A simple algorithm to find the L-curve corner in the regularisation of ill-posed inverse problems
We propose a simple algorithm to locate the "corner" of an L-curve, a function often used to select the regularisation parameter for the solution of ill-posed inverse problems. The algorithm involves the Menger curvature of a circumcircle and the golden section search method. It efficiently finds the regularisation parameter value corresponding to the maximum positive curvature region of the L-curve. The algorithm is applied to some commonly available test problems and compared to the typical way of locating the l-curve corner by means of its analytical curvature. The application of the algorithm to the data processing of an electrical resistance tomography experiment on thin conductive films is also reported.
\section{Introduction}	
 
The solution $\bm{\hat x}$ of an ill-posed inverse problem is often formulated in terms of residuals
 
\begin{equation}
\begin{aligned}
\label{eqn:lsq}
\bm{\hat x} =\textup{arg }  \underset{\bm{x}}{\textup{min}} \left \{||\bm{A}\bm{x}-\bm{b}||^{2}\right\}, \quad \bm{A} \in \mathds{R}^{m \times n}, \quad m>n
\end{aligned}
\end{equation}
where the quantity $\bm{Ax}-\bm{b}$ is the vector of residuals between the experimental data vector $\bm{b}$ and the reconstructed data $\bm{A}\bm{x}$ for a given $\bm{x}$. Regularization techniques are applied to make the problem less sensitive to the noise of $\bm{b}$ and find a stable solution. 
One popular regularization method is to replace the system~\ref{eqn:lsq} with 
\begin{equation}
\begin{aligned}
\label{eqn:funct}
\bm{\hat x}_{\lambda}= \textup{arg }  \underset{\bm{x}}{\textup{min}} \left \{ ||\bm{A}\bm{x}-\bm{b}||^{2}+\lambda \bm{R}(\bm{x})\right \},\quad \lambda \in \mathds{R}, \quad \lambda \geq 0 
\end{aligned}
\end{equation}
where the regularization term $\bm{R}(\bm{x})$ represents a cost function, which may include some prior information about the solution. The scalar factor $\lambda$ is the  \emph{regularization parameter}, serving as weighing factor of $\bm{R}(\bm{x})$. The choice of $\lambda$ is crucial for the solution to be meaningful.
As an example, we consider the regularization method named after Tikhonov~\cite{tikhonov2013numerical}, in which $\bm{R}(\bm{x}) =  ||\bm{x}||^{2}$.
Several methods~(see \cite[p.~44]{almasy2009inverse} and references therein) have been developed in order to find an optimal tuning of $\lambda$ for a given problem. 
Of particular interest is the L-curve method~\cite{hansen1992analysis}. 
The L-curve is a two-dimensional curve, parametric in $\lambda$, defined by points with cartesian coordinates
\begin{equation}
\label{eqn:lcurve}
P(\lambda)=(\xi(\lambda), \eta(\lambda))\rightarrow\begin{cases}
& \xi(\lambda)= \log||\bm{Ax}-\bm{b}||^2\\
&\eta(\lambda)= \log||\bm{x}||^2
\end{cases}
\end{equation}
The point of maximum positive curvature $P(\lambda_\mathrm{opt})$, the ``corner'', can be associated to the optimal reconstruction parameter $\lambda_\mathrm{opt}$.  The underlying concept is that the ``corner'' represents a compromise between the fitting to the data and the regularization of the solution~\cite{hansen1993use}.
Numerical search algorithms have been proposed for the identification of $\lambda_\mathrm{opt}$; among them we mention the splines method~\cite{hansen1992analysis,hansen2007adaptive}, the triangle method~\cite{castellanos2002triangle} and the L-ribbon method~\cite{calvetti1999estimation}.
Here we propose an alternative method which is based on an estimation of the local curvature of the L-curve from three sampled points (which define a circle), and a sampling update rule based on the golden section search. The method attempts to reduce the computational effort by minimizing the number of sampled points of the L-curve explicitly computed. A detailed algorithm is described; an application to a reconstruction problem in electrical resistance tomography is reported.
\section{Algorithm}
\label{sec:alg}
The proposed algorithm~\ref{alg} is written in pseudo-code. The algorithm calls two functions. The function P = \texttt{l\_curve\_P($\lambda$)} is based on the the specific regularization problem being solved; it is assumed that at each call, given as input the regularization parameter $\lambda$ it solves the system (1) and provides as output the point $P(\lambda)$, that is the coordinates $\xi(\lambda)$ and $\eta(\lambda)$ of the L-curve. The function $C_k = \texttt{menger($P_j$, $P_k$, $P_\ell$)}$ is defined below in Sec.~\ref{sec:curv}.
The algorithm is iterative and identifies $\lambda_\mathrm{opt}$ with the golden section search, described in Sec.~\ref{sec:gss}.
\subsection{Curvature}
\label{sec:curv}
The function $C_k = \texttt{menger($P_j$, $P_k$, $P_\ell$)}$ is based on the definition of the curvature of a circle by three points given by Menger~\cite{menger1930untersuchungen} . In our case three values $\lambda_{j}<\lambda_{k}<\lambda_{\ell}$ of the regularization parameter identify three points $P(\lambda_{j})$, $P(\lambda_{k})$ and $P(\lambda_{\ell})$ on the L-curve.
Here we follow the notation of~\ref{eqn:lcurve} for the coordinates of a generic point $P(\lambda)$. For sake of simplicity of notation the substitution~\ref{eqn:subst} is made:
 %
\begin{equation}
\begin{aligned}
\label{eqn:subst}
\xi(\lambda_i)&\rightarrow \xi_i,\\
\eta(\lambda_i)&\rightarrow \eta_i,\\
P(\lambda_i)&\rightarrow P_i.
\end{aligned}
\end{equation}
So we define a signed curvature $C_{k}$ of the circumcircle as
\begin{equation}
\begin{aligned}
\label{eqn:menga}
C_k = \frac{2\cdot(\xi_{j}\eta_{k}+\xi_{k}\eta_{\ell}+\xi_{\ell}\eta_{j}-\xi_{j}\eta_{\ell}-\xi_{k}\eta_{j}-\xi_{\ell}\eta_{k})}
{\Big(\overline{P_{j}P_{k}}\cdot \overline{P_{k}P_{\ell}}\cdot \overline{P_{\ell}P_{j})}\Big)^{1/2}}, 
\end{aligned}
\end{equation}
where
\begin{equation}
\begin{aligned}
\label{eqn:defs}
\overline{P_{j}P_{k}}&=(\xb-\xa)^{2}+(\yb-\ya)^{2},\\
\overline{P_{k}P_{\ell}}&=(\xc-\xb)^{2}+(\yc-\yb)^{2},\\
\overline{P_{\ell}P_{j}}&=(\xa-\xc)^{2}+(\ya-\yc)^{2}, 
\end{aligned}
\end{equation}
are the euclidean distances between the sampled L-curve points. Note that we choose to index the curvature with the intermediate index ($k$) of the three points.
\subsection{Golden Section Search}
\label{sec:gss}
The algorithm is initialized by assigning the search interval $[\lambda_1,\lambda_4]$. Two other values $\lambda_{2}$ and $\lambda_{3}$ are calculated as
\begin{equation}
\begin{aligned}
\label{eqn:goldenratio}
\lambda_{2} &= (\lambda_{4}+\varphi\cdot\lambda_{1})/(1+\varphi), \\
\lambda_{3} &= \lambda_{1}+(\lambda_{4}-\lambda_{2}),\\
\end{aligned}
\end{equation}
where $\varphi = (1+\sqrt{5})/2$ is the golden section\cite{ kiefer1953sequential} .
 %
Four values of $\lambda$ define four points on the L-curve and allow to calculate two curvatures, $C_{2}$ from $\{P(\lambda_{1})$, $P(\lambda_{2})$, $P(\lambda_{3})\}$ and $C_{3}$ from $\{P(\lambda_{2})$, $P(\lambda_{3})$, $P(\lambda_{4})\}$. 
The curvatures $C_{2}$ and $C_{3}$ are compared; consistent reassignment and recalculation are done in order to work at each iteration with four points $P(\lambda_1)\ldots P(\lambda_4)$.  
The algorithm terminates when the search interval 
$[\lambda_1,\lambda_4]$ is smaller than a specified threshold $\epsilon$ and returns $\lambda_\mathrm{opt}$.
 
Note that it may happen that the curvature $C_3$ associated to the right-hand circle is negative at the initial stage of the search, since $C_k$ is defined with sign in~\ref{eqn:menga}. Moreover by definition the corner must always correspond to a positive curvature and it lays on the left-side of the plot.  For these reasons the algorithm performs a check, and while $C_{3} < 0$ the search extreme $\lambda_{1}$ is kept fixed, $\lambda_{4}$ is shifted toward smaller values and   $\lambda_{2}$ and  $\lambda_{3}$ are recalculated. The condition on $C_3$ is strong enough that even in case of both negative curvatures it guarantees the convergence toward the corner.

Some considerations can be done: (a) according to the golden section search method, the algorithm needs to recalculate only one $P(\lambda)$ at each iteration (except for the first iteration), the other can be simply reassigned; this limit the calculation effort; 
(b) as $P(\lambda_{1})$ and $P(\lambda_{4})$ are distant at the first iterations, $C_2$ and $C_3$ are just rough approximations of the curvature of the L-curve in different regions, but become more accurate as the distance between the search extremes decreases.
\begin{algorithm}[H]
\caption{L-curve corner search}
\label{alg}
\begin{algorithmic}[1]

\STATE{Initialize  $\lambda_{1}$ and $\lambda_{4}$;}
	\COMMENT{search extremes}
\STATE{Assign $\epsilon$;}
	\COMMENT{termination threshold}
\STATE{$\varphi \leftarrow (1+\sqrt{5})/2$;}
	\COMMENT{golden section}
\STATE{$\lambda_{2} \leftarrow (\lambda_{4}+\varphi\cdot\lambda_{1})/(1+\varphi)$;}
\STATE{$\lambda_{3} \leftarrow \lambda_{1}+(\lambda_{4}-\lambda_{2})$;}
 \FOR{$i=1$ to $4$}
 	\STATE{$P_{i}\leftarrow$\texttt{l\_curve\_P}($\lambda_i$);}
		\COMMENT{func. \texttt{l\_curve\_P} returns~\ref{eqn:lcurve}}
 \ENDFOR
\REPEAT
\STATE{$C_2 \leftarrow$\texttt{menger}($P_1$,$P_2$,$P_3$);}
	\COMMENT{func. \texttt{menger} calls~\ref{eqn:menga}}
\STATE{$C_3 \leftarrow$\texttt{menger}($P_2$,$P_3$,$P_4$);}
	
	\REPEAT 
		\STATE{$\lambda_{4} \leftarrow \lambda_{3}$;\quad $P_{4} \leftarrow P_{3}$;}
		\STATE{$\lambda_{3} \leftarrow \lambda_{2}$;\quad $P_{3} \leftarrow P_{2}$;}
		\STATE{$\lambda_{2} \leftarrow (\lambda_{4}+\varphi\cdot\lambda_{1})/(1+\varphi)$;}
		\STATE{$P_{2}\leftarrow$\texttt{l\_curve\_P}($\lambda_2$);}
		\STATE{$C_3 \leftarrow$\texttt{menger}($P_2$,$P_3$,$P_4$);}
	\UNTIL{$C_3 >0$}

\IF{$C_{2} > C_{3}$}

	\STATE{$\lambda_\mathrm{opt}\leftarrow\lambda_{2}$;}
		\COMMENT{store $\lambda_\mathrm{opt}$}
	\STATE{$\lambda_{4} \leftarrow \lambda_{3}$; \quad $P_{4} \leftarrow P_{3}$;}
	\STATE{$\lambda_{3} \leftarrow \lambda_{2}$; \quad $P_{3} \leftarrow P_{2}$;}
	\STATE{$\lambda_{2} \leftarrow (\lambda_{4}+\varphi\cdot\lambda_{1})/(1+\varphi)$;}
	\STATE{$P_{2}\leftarrow$\texttt{l\_curve\_P}($\lambda_2$);}	
		\COMMENT{only $P_2$ is recalculated}
\ELSE
	\STATE{$\lambda_\mathrm{opt}\leftarrow\lambda_{3}$}
	\STATE{$\lambda_{1} \leftarrow \lambda_{2}$; \quad $P_{1} \leftarrow P_{2}$;}
	\STATE{$\lambda_{2} \leftarrow \lambda_{3}$; \quad $P_{2} \leftarrow P_{3}$;}
	\STATE{$\lambda_{3} \leftarrow \lambda_{1}+(\lambda_{4}-\lambda_{2})$;}
	\STATE{$P_{3}\leftarrow$\texttt{l\_curve\_P}($\lambda_3$);}
		\COMMENT{only $P_3$ is recalculated}	
\ENDIF
\UNTIL{ $(\lambda_{4}-\lambda_{1})/\lambda_{4}<\epsilon$}
\RETURN $\lambda_\mathrm{opt}$

\end{algorithmic}
\end{algorithm}

\begin{figure}[h]
\centering 
\subfloat[Iteration 1]{\label{fig:lcurves1}\includegraphics[height=9cm]{L1.eps}} 
\subfloat[Iteration 6]{\label{fig:lcurves2}\includegraphics[height=9cm]{L2.eps}} 
\subfloat[Iteration 28]{\label{fig:lcurves3}\includegraphics[height=9cm]{L3.eps}} 
\caption{The proposed algorithm at various iterations. The reference L-curve is reported as a solid line. Solid circles represent the points $P(\lambda)$ evaluated at the present iteration. Empty circles represent the points evaluated at past iterations. In~\ref{fig:lcurves1} arrows indicate the initial search extremes. In~\ref{fig:lcurves3} the four points are no more resolved, see Fig.~\ref{fig:zoom} for clarity.}
\label{fig:lcurves}
\end{figure}
\section{Example of application}
\label{sec:experiments}

In the following we show an example of application of the algorithm to electrical resistance tomography~\cite{holder2004electrical, seo2013electrical}. This example aims to image the sheet conductance spatial distribution of tin-oxide films on glass substrates having circular geometry and with electrical contacts on the sample boundary as shown in Fig.~\ref{fig:reca}. 
Four-terminal resistance measurements among the contacts are performed with a scanning setup; the values are the elements of the data vector $\bm{b}$. A detailed description of the experiment is given in~\cite{cultrera2016Electrical}.  

The reconstructed image (shown in Fig.~\ref{fig:recb}) is obtained by solving a discretized Laplace equation. Tikhonov regularization is applied. 
Two software packages of common usage in the solution of ERT problems (EIDORS~\cite{adler2005eidors}, Regularisation Tools~\cite{hansen1994regularization} ) are here employed running under Matlab environment. EIDORS routines are used to generate a two-dimensional circular mesh (2304 elements) with 16 contact points (corresponding to a $\bm{b}$ of size 208) on the mesh boundary, to discretize the Laplace equation and thus obtain the matrix $\bm{A}$. 
The library Regularisation Tools is employed to implement the function \texttt{L\_curve\_P($\lambda$)}. Other EIDORS routines are used to graphically render the solution $\bm{x}$ as in Fig.~\ref{fig:recb}.

The proposed algorithm is implemented in Matlab as well. Figures~\ref{fig:lcurves} and~\ref{fig:zoom} show the evolution of the search for $\lambda_\mathrm{opt}$ performed by the algorithm. These figures display also a full L-curve obtained by dense sampling of \texttt{L-\_curve\_P($\lambda$)} as a reference.

The algorithm is run by choosing as initial search extremes $\lambda_{1} = 1\cdot10^{-10}$ and $\lambda_{4} = 1\cdot10^{-3}$.
Fig.~\ref{fig:lcurves} reports three snapshots of the running algorithm. Empty circles represent points visited at previous iterations, while filled circles represent the four points $P_1 \ldots P_4$ of the given iteration.
Fig.~\ref{fig:zoom} shows a detail of the last iteration. The reader can notice in Fig.~\ref{fig:zoom1} the L-curve the ``corner''. Fig.~\ref{fig:zoom2} shows the $P(\lambda_i)$ points of the last iteration, where the magnification factor does not allow anymore to appreciate the curvature.  
\begin{figure}[h]
\centering 
\subfloat[]{\label{fig:zoom1}\includegraphics[height=7cm]{L3Z1.eps}} 
\subfloat[]{\label{fig:zoom2}\includegraphics[height=7cm]{L3Z2.eps}}  
\caption{Detail of the last iteration. In~\ref{fig:zoom1} it is visible the ``corner'' of the present L-curve. In~\ref{fig:zoom2} the four $P(\lambda)$ are resolved. The label indicates $\lambda_{opt}$.}
\label{fig:zoom}
\end{figure}
Fig.~\ref{fig:conv} shows the evolution of the algorithm towards convergence. For sake of generality, the search extremes were chosen very far from each other; when knowledge about the problem allows to restrict the $\lambda$ search interval, the number of iterations can be reduced.
\begin{figure}[tbhp]
\centering 
\includegraphics[width=8cm]{convergence.eps} 
\caption{Behaviour of the algorithm versus iteration number for the present data with a threshold $\epsilon = 1\%$. The steep slope corresponds to lines 9-18 of the algorithm~\ref{alg}. The last point (iteration 28) corresponds to $\lambda_\mathrm{opt}$.}
\label{fig:conv}
\end{figure}
The $\lambda_\mathrm{opt}$ returned by the algorithm with a $\epsilon = 1\%$ is $1.84\cdot10^{-6}$. This result is used to perform the final image reconstruction shown in Fig.~\ref{fig:recb}.
\begin{figure}[tbhp]
\centering 
\subfloat[Scheme of the sample. Colors represent nominal values of sheet conductance.]{\label{fig:reca}\includegraphics[width=6cm]{reca.eps}} 
\qquad
\subfloat[Sheet conductance map reconstructed from the data with $\lambda=1.84\cdot10^{-6}$.]{\label{fig:recb}\includegraphics[width=6cm]{recb.eps}} 
\caption{Sheet conductance map obtained with the optimal regularization parameter $\lambda_\mathrm{opt}$ returned by our algorithm.}
\label{fig:rec}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}
The proposed algorithm allows, for a given inverse problem having the form \ref{eqn:lsq},  the determination of the regularization parameter $\lambda_\mathrm{opt}$ corresponding to the maximum positive curvature of the L-curve. The algorithm is designed for maximum simplicity of implementation on already existing solvers. In a test on real electrical resistance tomography problem, convergence is achieved in a couple of dozen iterations, each iteration corresponding to a redetermination of \ref{eqn:funct}.
\bibliographystyle{plain}

2016-08-17T02:08:15.000
1608
1608.04571
en
https://arxiv.org/abs/1608.04571
We propose a simple algorithm to locate the "corner" of an L-curve, a function often used to select the regularisation parameter for the solution of ill-posed inverse problems. The algorithm involves the Menger curvature of a circumcircle and the golden section search method. It efficiently finds the regularisation parameter value corresponding to the maximum positive curvature region of the L-curve. The algorithm is applied to some commonly available test problems and compared to the typical way of locating the l-curve corner by means of its analytical curvature. The application of the algorithm to the data processing of an electrical resistance tomography experiment on thin conductive films is also reported.
Numerical Analysis (math.NA)
A simple algorithm to find the L-curve corner in the regularisation of ill-posed inverse problems
Qwen/Qwen-72B
1. YES
2. YES