The Electrical Impedance Tomography (EIT) problem involves collecting electrical measurements on the smooth boundary ∂ to determine the spatially varying electrical conductivity distribution γ within the bounded region ⊆ R d (d = 2 , 3) . We assume γ has a strictly positive, isotropic and bounded conductivity distribution with no current sources inside . We represent the EIT forward problem with the following elliptic differential equation, div (γ ∇u ) = 0 in , (1) where u ∈ H 1 ( ) is the electric potential and γ is known.
An EIT experiment involves applying an electrical current (Neumann data) I on ∂ and measuring the resulting electrical potential (Dirichlet data) U on ∂ . The data collected then provides information on the Neumann-to-Dirichlet (NtD) map. Partial knowledge of the NtD map is used to approximate γ from a set of EIT experiments using different input currents [1][2][3] . EIT can be applied to nondestructive testing [4] , the monitoring of oil and gas mixtures in oil pipelines [5] , noninvasive medical imaging [6,7] , etc. For example, EIT has been used in civil engineering to monitor water infiltration into soil [8][9][10][11] . Stacey [12] also studied the EIT's feasibility by monitoring moisture movement in Berea sandstone. Hou and Lynch [13] applied EIT to loaded cementitious composites that were fiber reinforced.
EIT is well known to suffer from a high-degree of nonlinearity and severe ill-posedness [14,15] . Therefore, regularization is required to produce reasonable electrical impedance images. Most reconstruction methods are deterministic, such as the factorization method [17] , d-bar method [18] , and variational type methods for least-squares fitting. The variational type methods minimize a certain discrepancy functional using Tikhonov-based regularization or an iterative type method of a linearized model or fully nonlinear model such as sparisty cosnstraints [14,15] , iteratively regularized Gauss Newton method [16] just to mention a few examples. These analytical methods can be effective in determining specific conductivity distributions, but statistical inversion methods [19] can offer an alternative approach. In [20] , Kaipio et al. optimizes the current patterns based on criteria regarding functionals of the posterior covariance matrix. The Bayesian approach has also been used to study the errors from model reduction and partially unknown geometry [21,22] . In [23] , Bardsley quantified uncertainty of MCMC-based image reconstruction. Statistically, the sparsity regularization enforces the p prior on the expansion coefficients for a certain basis like the deterministic approaches for EIT [14,15,24] . In this paper, we apply the Pilot Adaptive Metropolis algorithm proposed in [25][26][27] to investigate the algorithm's efficacy with synthetic data. We use a TV prior and a combination of a TV prior and a 1 -type prior to solve a Bayesian inverse problem for which the dimension of the problem is neither too high neither too low (approximately 500 unknowns). In [40] , Lassas and Siltanen describe how the solution to a Bayesian inverse problem diverges when the dimension of the problem converges to infinity; however, for the dimension (approximately 500 unknowns) considered in our problem, we found experimentally that our method converges to a solution. Similarly, in [41] , Lucka found that, for 1 -type priors, the efficiency of MCMC methods dramatically decreased when the sparsity or dimension of the problem increased, which we confirmed in numerical experiments were only 1 is used as standalone priors. In this study, we find that the our MCMC method converged efficiently when we use a combination of 1 -type and TV priors.
The paper is organized as follows. In Section 2 , we describe the complete electrode model (CEM) for the EIT problem. In Section 3 , we formulate the inverse problem. We then describe the IRGN method in Section 4 . In Section 5 , we formulate the statistical inverse problem, and we discuss choices for regularizing the prior density. In Section 6 , we discuss the Markov Chain Monte Carlo (MCMC) method and Pilot Adaptive Metropolis algorithm. We compare the inversion results from the statistical approach with the traditional IRGN method in Section 7 to demonstrate the efficacy of our coupled model. Lastly, we conclude with Section 8.
Suppose L electrodes have been fixed around the surface of an object. Current is applied to a subset of these electrodes, and the resulting voltage is measured at all other electrodes. This EIT experiment is repeated several times with different electrode configurations to efficiently characterize the imaged object. The inverse problem then is to reconstruct the conductivity distribution inside the object from a finite set of surface point measurements.
The Eq. ( 1) is used to solve the electric potential at the electrodes and inside for some applied current I . We obtained (1) by a scaling analysis of Maxwell's equation for electromagnetic fields inside of an object [1,28] . We denote the class of admissible conductivities as
Currents are applied to electrodes on the surface ∂ of the body. Let j denote the inward unit normal component for the current density that is produced by current on the surface. Then, the Neumann boundary condition is
(1) and ( 2) are known as the continuum model . In practice, the current density j at the electrodes is unknown, but E l γ ∂u ∂n dS = I l is known, where n is the unit outward normal to , E l is the surface area of the l th electrode and I l is the current injected into E l . Thus the Neumann condition (2) can be rewritten as
Furthermore, we know j = 0 for the current density on the boundary between the electrodes, such that
The potential on the electrodes is considered to be constant: U l = constant . This property is known as the shunting effect , which is represented by
There is also an electro-chemical effect due to the formation of a thin and highly resistive layer between the electrodes and the body. Electrical impedance from this layer, z l , is called the effective contact impedance or surface impedance at E l . This effect changes (5) to
The complete electrode model (CEM) consists of (1), ( 3), ( 4) and ( 6) , together with the following conditions.
Although the CEM has a unique solution, the accuracy of it is determined by the predicted experimental measurements [29] .
For the FEM solution to the forward CEM model, we reformulate the variational form of CEM (see [28] for details) as
where f (v , V ) = L l=1 I l V l is a functional that maps H → R . Existence and uniqueness of a solution ( u , U ) ∈ H can be shown using the Lax-Milgram lemma. In [29] , it is shown that ( u , U ) satisfies the CEM conditions if and only if it satisfies (9) .
Following the typical FEM approach, let T = { T 1 , . . . , T | T | } be the triangulation of , which has N mesh points for the finite dimensional subspace
. We must determine the coefficients α i and β k in this formulation. Choosing v = φ i and V = ν k when the set of test functions is of the form (φ 1 , 0) , . . . , (φ N , 0) , (ν 1 , 0) , . . . , (ν L -1 , 0) results in the following system of equations in matrix form: ˆ A θ = f, (10) where θ = (α, β ) T ∈ R N+ L -1 , the matrix ˆ A is the standard FEM system matrix (see [28] for details) and the right hand side f is given by,
For a fixed current vector I and positive contact impedances
to be the forward operator that maps the conductivity γ to the solution of the forward problem obtained by FEM described above. [30] showed that this operator is Frèchet differentiable. Recall that Frèchet differentiability of F in γ means that lim
The forward problem uses the known conductivity distribution γ to find the boundary data associated with a given source.
Our goal is to estimate the conductivity distribution γ from all pairs of current vectors I ∈ ˜ R L and resulting voltage vectors U ∈ ˜ R L . We minimize the cost functional to recover γ since the inverse problem is severely ill-posed from the finite set of measurements.
where u δ approximates the exact data u with the accuracy δ, i.e.,
However, regularization is needed to improve the ill-posed problem and instead, we minimize,
where λ is the regularization parameter, R ( • ) is the regularization term and γ * is the known background. There are several choices for R ( • ). The p regularization R p (γγ * ) is defined as
where 0 < p ≤ 2 is a constant, [19] . The p regularization enforces sparsity for 0 < p ≤ 1 and smoothness when p ≥ 2. Total variation regularization is used for most practical applications to obtain smooth images. Total variation is defined as [42,43] ,
where | • |( ) is a finite Radon measure. In particular, for γ ∈ W 1,1 ( ), it reduces to the standard notation for TV regularization,
The regularization function is represented by a norm for most analytical methods. In this paper we used one of the most successful methods for solving the ill-conditioned problem, R 2 known as Tikhonov regularization. The cost functional from Tikhonov regularization is
where W is weight function and upon discretization becomes a weight matrix. There are several iterative approaches to minimize (13) . In this paper, we used a modified iteratively regularized Gauss-Newton (IRGN) method for the minimization, which is described in the next section.
Suppose λ k is some sequence of regularizing parameters satisfying the conditions
Let the unique global minimum of (13) be denoted by ˜ γ . Assume ˜ γ satisfies the invertibility conditions and the conditions for F as described in [31] . Then, the unique global minimum of ( 13) is explicitly given by ˜
where γ k is the k th approximation to γ and F (γ k ) is the Jacobian matrix at the k th iteration, and W 2 = W T W . The above algorithm is generalized further using a line search procedure from Smirnova et al. [31] . A variable step size,
The modified IRGN algorithm is then
The line search parameter s k is chosen to minimize the scalar objective function
where p k is the search direction, which solves
This step is accomplished through a backtracking strategy until either one of the strong Wolfe conditions,
is satisfied [33] , or the maximum number of backtracking steps has been reached. We use the theoretically derived values of c 1 = 0 . 0 0 01 and c 2 = 0 . 9 , derived in [33] . Due to the inexact nature of u δ , we adopt a stopping rule from Bakushinsky and Smirnova [32] to terminate the iterations at the first index K = K(δ) , such that
Computing the Jacobian of the forward operator is required for most iterative Newton methods. If the regularity assumptions on the domain and the coefficients are satisfied, the forward operator is differentiable. The operator F is Frèchet differentiable. It maps γ ∈ int( Q ) to the solution ( u , U ) ∈ H of the forward problem with current vector
Inverse problems are typically written in terms of a minimization problem. In this section, we write the inverse problem in terms of the posterior density of the conductivity γ given the measurements U on ∂ to obtain the expected value of the conductivity given the surface measurements. This estimate is a reasonable point estimate to the ill-posed inverse problem. In the following, we formulate the finite dimensional posterior density for the conductivity distribution given the finite dimensional measurement on ∂ . Let δ L represent the L dimensional measurement noise corresponding to the measurements U L . Denote γ * M , U * L and δ * L as the random variables for γ M , U L and δ L , respectively.
For the posterior density π γ * M (γ M | U L , δ L ) , we first assume the EIT forward operator ( γ M , δ L ) is well posed, such that
with the property that R m (x )
from ( 24) and π X (x
π Y (y ) . Here, the x and y are values of the random variables X and Y , respectively. By integrating the measurement noise δ L out of Eq. ( 25) , we can then obtain the joint density of γ * M and
We assume that the measurement noise
, and that it is additive, that is U L = (γ M ) + δ L . Eq. ( 27) can now be simplified to
The posterior density of γ * M given the measurements
from the conditional probability formula. Note that π
The posterior density can be defined by choosing a density for the measurement noise δ * L and a prior density function for γ * M . We assume data with random Gaussian noise,
where C is a positive definite covariance matrix. The density of γ * M is the prior density because it contains prior knowledge about γ * M . We also assume
where R ( • ) is a regularization function, α > 0 is a constant, and χ A ( γ M ) a indicator function with A = [0 , ∞ ) n .
In this section, we discuss several choices for the regularization function R ( • ). We can apply most regularization functions in both the statistical and the analytical settings. However, the statistical setting typically allows to chose from more general regularizations then its analytical counterpart.
The p regularization R p (γ M ) is defined as
where c i represent the weights, γ b M is a typical background conductivity and 0 < p ≤ 2 is a constant [19] . R p (γ M ) is a norm if p ≥ 1 and only defines a metric when 0 < p < 1. The regularization function is represented by a norm for most analytical methods, but the statistical reconstruction can also handle cases when 0 < p < 1. In theory, a good choice for the weights c i are large values at the boundary and exponentially decreasing values towards the center of because the variance is smaller on the boundary than in the center [34] . The p regularization enforces sparsity when 0 < p ≤ 1 and enforces smoothness when p ≥ 2.
Total variation regularization is used to obtain smooth images, which is needed for most practical applications. The total variation regularization for an infinite dimensional space is defined as R T V c (γ ) as previously defined with γ as the infinite dimensional version of γ M . The discrete analogue for a two-dimensional body of the total variation regularization R T V c [19] is
where l i is edge length of the i th adjacent triangle and * i = (0 , 0 , . . . , 0 , 1 a (1 ,i ) , 0 , . . . , 0 , -1 a (2 ,i ) , 0 , . . . , 0) , (38) where a = { a ( j,i ) } h i =1 , j∈{ 1 , 2 } is the set of the numbers for all adjacent triangle tuples ( a (1, i ) , a (2, i ) ).
In general, we can choose π γ * M (γ M ) to be any kind of prior density on instead of using the classical analytical choices. We recommend choosing a continuous regularization function for γ M that is integrable over for proper reconstructions.
The p and total variation prior can provide better inversion results when combined together [26] ,
where α 1 , α 2 > 0 are constants. For proper choice of α 1 , α 2 > 0, this regularization function can provide smooth images that give the correct background conductivity.
We have discussed the posterior density with several meaningful prior densities or regularization functions. Now, we can find the Bayesian estimate for γ * M based on the measurements
There is no direct method for finding the Bayesian estimate
have a closed form. However, we can use the Markov Chain Monte Carlo Method (MCMC) to generate a large set of random result from Gelman et al. [37] . For a normal target and proposal distribution, the optimal acceptance ratio is approximately .45 in the one dimensional case and .234 when the number of dimension converges to infinity [37] . Suppose we perform t 0 adaptions, one every t r iterations, where 1 < t r t 0 < B * < N . Let c i denote a variable recording if the i th iteration of this algorithm has been accepted, c i := 1 , i th iteration has been accepted, 0 , else.
The estimator for the acceptance ratio for the k th proposal distribution is
Informally speaking, the algorithm modifies the covariance matrix in the pilot time t r t 0 to converge closely to the optimal acceptance ratio. Then, the standard Metropolis-Hastings algorithm begins with the latest state and proposal distribution of the pilot time.
In Algorithm 1 , the pilot adaptive Metropolis algorithm is recapitulated with an arbitrary starting state x (0) ∈ E and a starting guess for the positive definite covariance matrix C 0 . The chain only satisfies the Markov property after the last adaption at time t r t 0 in Algorithm 1 . The chain usually still moves towards the high probability areas of the target distribution during the pilot time [25][26][27]35] .
We generate synthetic data for our simulations with a known conductivity distribution γ . We proceed to solve the forward problem using γ in the finite element Galerkin method. For each case, we had two different discrete conductivity distributions inside a circular domain . We used | T | = 4128 triangles and N = 2129 linear basis functions for u and γ on . There are L = 16 equally spaced electrodes placed along the boundary with each electrode covering a surface area of 5mm. Next, we added random Gaussian noise with a standard deviation that is 1% and 3% of the maximum measurement, respectively. In order to avoid the inverse crime, we reconstructed γ on with a smaller mesh size consisting of 1032 triangles with 549 mesh points.
Our simulations specifications are listed in the following table, which includes the known γ with a residual δ in (22) and a stopping rule for choosing ρ in (12) .
The sequence of step lengths s k is chosen through a backtracking strategy, s 1 , s 1 / 2 , . . . until either the strong Wolf condition from (20) or ( 21) is satisfied. The maximum number of backtracking steps was set to 16. We also imposed
)10 -3 to prevent singularity in the numerical computation.
No backtracking is possible if the sequence of regularization parameter λ k decreases too quickly. We define λ k = λ 1 c c+ k -1
with λ 1 = 1 , c = 4 , to provide us with ˆ d = 1 . 25 .
The true conductivity consists of a homogeneous background and one circular inclusion of radius 0.01 mm centered at (0.0242 mm, 0.015 mm). The conductivities of the background and the inclusions are 7 • 10 -4 Ohm and 10 -8 Ohm , respectively, as shown in Fig. 1 . The reconstructions with 1% and 3% noise levels are shown in Fig. 2 (a) and 2 (c) for the IRGN method and Fig. 2 (b) and 2(d) for the statistical inversion method.
Inversion results from the IRGN method are smooth from Tikhonov regularization. The inclusion location is effectively captured, but its support is slightly larger than the true inclusion. In particular, its significantly extended towards the center of when the noise level is higher at 3%. The magnitude of the inclusion conductivity is overestimated by 10 -4 Ohm compared to the true value of 10 -8 Ohm . In contrast, the reconstruction by the statistical inversion method is more localized at the true location with a reasonably homogeneous background. However, with a higher noise level, the background conductivity starts to get distorted. We computed the 1 and 2 reconstruction errors as
, where γ T is the true conductivity distribution.
The true conductivity consists of a homogeneous background and double circular inclusions of radius 0.01mm centered at ( ± 0.036mm, 0mm). The conductivities of the background and the inclusions are 7 • 10 -4 Ohm and 10 -8 Ohm , respectively, as shown in Fig. 3 . The reconstructions with 1% and 3% noise levels are shown in Fig. 4 (a) and 4 (c) for the IRGN method and Fig. 4 (b) and 4 (d) for the statistical inversion method. We observe that both methods are able to retrieve the inclusions. In IRGN method, the supports of the inclusions are extended towards the center of , and the magnitude of the background conductivity is slightly overestimated. However, the statistical approach provides a more localized solution with a sharper background. We list the reconstruction errors for γ in Table 3 , which shows that the errors e 1 and e 2 are smaller for the statistical inversion method.
The true conductivity consists of a homogeneous background and double circular inclusions of radius 0.01 mm centered at ( ± 0.036 mm, 0 mm) and (0 mm, ± 0.036 mm). The conductivities of the background and the inclusions are 7 • 10 -4 Ohm and 10 -8 Ohm , respectively, as shown in Fig. 5 . Multiple inclusions are challenging for most numerical algorithms. However, both the approaches we used in this paper produce reasonable reconstructions of γ for different noise levels. As expected, with higher noise level and more inclusions, their support becomes larger in IRGN method, eventually making them barely observable. On the other hand, the statistical inversion does produce significantly sharper images with more localized inclusions even for a higher noise level.
We constructed the EIT forward model using the CEM as the Galerkin finite element approximation and then discussed the setup for the linearized inverse problem. We discussed several regularization functions for the statistical inverse problem. We compared these general regularization functions for the statistical inverse problem to the regularization function used in the deterministic setting. We obtained a Bayes' estimate for the electrical conductivity from the Pilot Adaptive Metropolis algorithm that implements a smoothing criteria and enforces sparsity in the prior distribution. We used both the p and total variation priors for regularization in the Pilot Adaptive Metropolis algorithm. We found that coupling sparsity regularization and smoothness regularization improved performance more than the case when only sparsity constraints or smoothness constraints were used.
The smoothness regularization provided regularity in the solution while the sparsity approach reduced the ill-posed nature of the inverse problem. Both criteria resulted in electrical impedance images with higher resolution. We determined that sparsity and smoothness regularization are needed during EIT inversion for improved image reconstruction. For this study, the statistical algorithm provided better reconstructions compared to IRGN in terms of both 2 and 1 errors. One drawback to the statistical approach is that it can be computationally expensive to run. For low dimensional cases, the running times were about 40-45 min when we used the deterministic method as an initial guess.
Additional disadvantages of the statistical algorithm and IRGN method are their dependencies on the (1) location of the inclusions and (2) data noise. Reconstruction errors are higher as data noise increases and the location of inclusions move further inside the object. This result is a direct consequence of the ill-posed inverse problem. We avoided this effect by creating inclusions closer to the surface for higher data resolution and using a relatively low noise level. Most experimental noise for EIT does not exceed 10%, but we do not know the noise distribution. Our EIT data is synthetic and produced from a controlled environment. We will continue our work to determine if we can reconstruct high resolution images when the noise distribution is unknown. Our IRGN method and statistical algorithm will be tested further with experimental data. Overall, we found that the statistical approach is a promising and computationally efficient algorithm when it uses the inverted results from the IRGN as an initial starting point.