Methods for calculating the electrode position Jacobian for impedance imaging

Electrical impedance tomography (EIT) or electrical resistivity tomography (ERT) current and measure voltages at the boundary of a domain through electrodes. Significance: The movement or incorrect placement of electrodes may lead to modelling errors that result in significant reconstructed image artifacts. These errors may be accounted for by allowing for electrode position estimates in the model. Movement may be reconstructed through a first-order approximation, the electrode position Jacobian. A reconstruction that incorporates electrode position estimates and conductivity can significantly reduce image artifacts. Conversely, if electrode position is ignored it can be difficult to distinguish true conductivity changes from reconstruction artifacts which may increase the risk of a flawed interpretation. Objective: In this work, we aim to determine the fastest, most accurate approach for estimating the electrode position Jacobian. Approach: Four methods of calculating the electrode position Jacobian were evaluated on a homogeneous halfspace. Main results: Results show that Fréchet derivative and rank-one update methods are competitive in computational efficiency but achieve different solutions for certain values of contact impedance and mesh density.

In both geophysics and biomedical impedance imaging, electrode movement and boundary modelling errors can be the cause of significant image reconstruction artifacts , Zhou and Dhalin 2003, Boyle and Adler 2011. Reconstruction artifacts can sometimes be difficult to identify and can make a correct interpretation of a reconstructed conductivity distribution challenging. Data misfit can be informative in detecting such errors, though attributing data misfit to specific artifacts is generally not explored in any particular depth. Therefore, mismatches between the physical and model boundary, boundary errors, may increase the risk of an incorrect interpretation. Movement of the model boundary may be estimated from electrode movement by interpolating surface movement between the electrodes. Some artifacts may be addressed by adapting the model to allow for electrode position estimates. In two dimensions, movements that are conformal do not introduce artifacts (Boyle et al 2012a). In three dimensions and assuming full Robin-Neumann data at the boundary, for a convex domain the conductivity, shape, and contact impedance are unique up to a rigid transformation (a rotation or translation) because the other 4 dimensions of a conformal map, the three specifying a Möbius inversion and one for scaling, disappear from the null space (Kolehmainen et al 2007). Non-conformal shape changes may be explained by a specific anisotropic conductivity, isotropic conductivity and shape change, or anisotropic conductivity and shape change (Lionheart 1998). For an assumed isotropic conductivity, exactly those electrode displacements that are not conformal may be reconstructed by setting electrode positions as model parameters in the inverse problem. To reconstruct electrode position in an L2-norm Gauss-Newton framework, an electrode position Jacobian must be constructed to estimate the first derivatives of the measurements with respect to the electrode position.
In this work, we aim to determine the fastest, most accurate approach for estimating the electrode position Jacobian. The four Jacobian methods, detailed in the following sections, are: the naïve perturbation method, a minimal mesh perturbation method, a rank-one matrix update (Gómez-Laberge and Adler 2008), and the recently developed Fréchet derivative (Dardé et al 2012). We make use of a simplified four electrode homogeneous half-space model to investigate computational efficiency and deviations from anticipated solutions. Throughout this work, the conductivity is assumed to be isotropic. By evaluating these methods, we mean to explore the sensitivity of the perturbation methods, understand and implement the Fréchet derivative, and ask the Engineering question: Do these methods give the same answer for a simplified model?

Applications of the electrode position jacobian
For impedance imaging, the electrode position Jacobian can be used to reconstruct electrode position. This is practical when the interior conductivity distribution is known and the forward model is otherwise accurate. In general, the interior conductivity is not known and a joint reconstruction of electrode position and conductivity can significantly reduce the occurrence of artifacts (Soleimani et al 2006, Gómez-Laberge and Adler 2008, Jehl et al 2015. One approach to achieving a joint reconstruction is an extension of the conductivity reconstruction over an L2-norm using the well known Gauss-Newton iterative method (Nocedal and Wright (2006), section 10.2). Over the course of a number of iterations, the conductivity and electrode position may be reconstructed based on the EIT measurements. An update δσ δx , ( ) is calculated (1) based on the combined Jacobian J, the inverse noise covariance of the measurements W, regularization R scaled by hyperparameter λ, the measurement misfit b, and the prior misfit c. A line search may be used to determine α the amount of the update to apply (2). (2) The measurement misfit b is calculated for each iteration using a forward model F from the conductivity σ n and movement x n estimates at iteration n. Likewise, the prior misfit is calculated as the difference between the current parameter estimates σ x , n n and the prior estimates σ x , . The combined Jacobian is simply the conductivity σ J and electrode position Jacobians J x concatenated. The regularization matrix R is a block diagonal matrix combining conductivity regularization σ R and movement regularization R x . The relative strength of the regularization is controlled by the ratio of the regularization hyperparameters β η / . The initial conductivity is generally taken to be a homogeneous estimate of the conductivity or the prior.
Nearly the same formulation may be used for a difference reconstruction where there is a baseline set of measurements and measurements at a later time. The change in conductivity and electrode position may then be computed as a single-step in the Gauss-Newton update (1) for a sufficiently linear change. A line search (2) is typically not necessary (α = 1). The prior is generally taken to be the background conductivity and initial electrode positions ( = c 0) giving where the Jacobian is calculated on a prior estimate of conductivity and electrode position. In the difference solution (3), mutual errors in the sets of measurements m m , 0 1 will generally cancel making the reconstruction problem slightly less challenging since the forward model does not need to match the real system with as much fidelity as the absolute reconstruction (1). Joint reconstructions, regularization, and balancing the two parameter types are not explored further in this work. Joint electrode position and conductivity reconstructions are an on-going topic of research in both the biomedical electrical impedance tomography (EIT) and geophysics electrical resistivity tomography (ERT) fields (Jehl et al 2015, Wilkinson et al 2016.

Electrode position jacobian
An 'electrode position Jacobian' is a first-order estimate of the effect of an electrode movement on the impedance imaging measurements. In impedance imaging, the electrode position Jacobian has a prominent role in some techniques for addressing the issue of boundary errors. An electrode position Jacobian matrix J x , captures how a small electrode movement δx affects each measurement δU (figure 1) where matrix row i, column j is the ith measurement's change in voltage given a small movement on electrode j. The electrode position Jacobian estimates the effect of an electrode's movement in the same way that a conductivity Jacobian reflects how a regional change in conductivity σ j affects each measurement U i . For the point electrode model (PEM), the Jacobian captures movement effects, but real electrodes modelled using the complete electrode model (CEM) have non-infinitesimal geometry and an associated contact impedance which affects current flow in the vicinity of the electrode (Somersalo et al 1992). In general, the size, shape, rotation, and contact impedance are assumed to be constants in calculating the Jacobian. Electrodes moved in such a way are treated as inflexible and non-rotating electrodes. This has some clear limitations in the context of curved surfaces where the electrode can only move a certain distance without some sort of rotation to maintain contact with the boundary. In practice, the components of a movement vector are typically applied as tangential and normal movements rather than translation in absolute coordinates. Contact impedance and full movement vectors could, in principle, be estimated as columns in the Jacobian. Contact impedance and the conductivity distribution have been simultaneously reconstructed for impedance images (Heikkinen et al 2002, Winkler and Rieder 2015. Sensitivity to normal movements are usually such that their effects are difficult to estimate accurately unless measurements are taken on driven electrodes Lionheart 2015, Crabb et al 2017).
For a minimizing the 2-norm of a residual using Gauss-Newton updates, the Jacobian J, regularization R, hyperparameter λ and data misfit r, determine the search direction δx. A line search then determines a distance α along that search direction. If the system is sufficiently linear, a line search is unnecessary and α = 1 can work well.
Errors of magnitude in the Jacobian may be circumvented by a reasonably accurate line search. On the other hand, any systematic errors in the electrode position Jacobian's direction estimate may entirely confound an attempt to address the electrode position issue. Therefore, the electrode position Jacobian plays a key role in estimating any electrode placement error in the model or accounting for electrode movements. ) with a measured potential difference V relative to a reference electrode (dark green ) is moved δx resulting in a change in the measured voltage

Movement artifacts
Surprisingly small amounts of electrode or surface movement can be a source of significant conductivity reconstruction artifacts when electrodes are incorrectly positioned. Electrode movements of as little as 5.7% of electrode spacing cause significant artifacts for elliptical 16 electrode (lung EIT) configurations . Electrode movements of 10% have been observed to give 20% resistivity artifacts for a collinear electrode array (2D geophysics ERT) configurations (Zhou and Dhalin 2003). To an experienced observer, an individual small electrode movement can be identified as a conductivity image artifact surrounding the erroneously placed electrode.
In figure 2, a two-dimensional half-space finite element method (FEM) model with background resistivity of 100 Ω · m, a conductive target (90 Ω · m), and a resistive target (110 Ω · m) were simulated to generate measurements on a 32 electrode linear array (5 m electrode spacing) with the Wenner stimulus pattern (Barker 1979, Loke andBarker 1996). Conductivity images were reconstructed after adding noise (40 dB SNR additive Gaussian white noise (AGWN)) and moving two electrodes (electrode #2 and #12) in opposite directions by up to 50% (2.5 m) of electrode spacing (5 m). Artifacts initially appear near the moved electrodes for small movements and, for larger movements, spread to contaminate the entire image by masking the true targets (figures 2(b)-(g)). When movements are reversed the artifacts near the electrodes are generally reversed: conductive artifacts become resistive and vise-versa (compare figures 2(g) and (h)). Larger movements or movements of multiple electrodes, where the effects interact, can lead to images that are very difficult to interpret. More extreme electrode movements may result in a reconstruction that fails to converge such that there is no useful image to interpret.
The electrode position issue is, in fact, more tractable than it first appears. Impedance imaging methods are quite good at detecting small electrode movements. The artifacts caused by these electrode movements can be hard to distinguish from each other and from other noise induced artifacts in a reconstruction, but by examining the data misfit, rather than the artifacts, correlations between electrodes can be uncovered. Data misfit that is strongly correlated with particular electrode movements will be identified in a joint Gauss-Newton reconstruction (1). Using the same forward model (figure 2(a)), the change in all measurements was plotted against a single electrode's movement over the range of −50% to +50% of electrode spacing (figure 3). The measurement changes were nearly linear for a single electrode's movement: generally within 1% of measurements. The Jacobian, therefore, was nearly linear for a relatively broad range of electrode movement. We note that the assumption of orthogonality in the Jacobian will be approximately correct for 'small' electrode movements or an individual large movement.

Perturbation Jacobian
Calculations to estimate the electrode position Jacobian by direct perturbation of an FEM mesh requires multiple forward solutions: one per electrode and position dimension. (In three dimensions, a 32 electrode array would require 96 additional forward solutions per Gauss-Newton iteration.) The method is a direct implementation of equation (4): the partial derivatives are replaced with a small perturbation δx j of electrode j's location where δU i is the difference between two forward model F solutions and the only change between forward models is the location of electrode j. The method inherently captures all relevant movement effects in the model because the model is completely rebuilt to incorporate the updated electrode position. Figure 4(a) illustrates a naïve perturbation implementation for movement of a single CEM electrode's nodes in an FEM mesh. Perturbations such as the one illustrated here must select a small enough perturbation to approximate the partial derivative while using a large enough perturbation to remain numerically stable. The perturbation method is fundamentally a calculation of the difference ∆u between two inverted FEM system matrices − A 1 given the same stimulus B Figure 2. Electrode movement artifacts; simulated reconstructions based on forward model (a), each with two electrodes (electrode #2 and #12, circled, of 32 electrodes numbered left-to-right at 5 m intervals) having electrode displacements of (b) 0% (c) 1%, (d) 5%, (e) 10%, (f) 25% and (g) 50% of electrode spacing on a 2-dimensional halfspace reconstruction (40 dB SNR, λ = 0.01, Laplace regularization, Wenner stimulus pattern). Note that when electrode movement is reversed (h) −50%, the conductivity artifacts near the moved electrodes are, for the most part, reversed. Single or well separated electrode location errors introduce characteristic 'ringing' artifacts that can overwhelm conductivity-based information. (i) Resistivity (Ω · m).
where the operator T selects the difference between electrodes from all FEM nodal voltages X to form a vector of all measurements. Numerical stability issues due to a perturbation are typically observed when the difference in floating point errors between two system matrices overwhelm the change in measured voltage at the electrodes. The electrode position Jacobian is a function of the system matrix which is dependent on the conductivity distribution, boundary, electrode placement, electrode size and electrode contact impedance. The selection of a stable perturbation magnitude is then dependent on variables that are typically unknown or uncertain to a significant degree. We observe that a difference solution will cancel common errors in two inverted system matrices. Movement of all of an electrode's nodes affects the shape of all elements directly connected to the electrode. A piece-wise linear forward model of potential within the domain usually requires smaller elements near the electrode boundary to accurately approximate the electrical discontinuity at the edge of an electrode. Therefore, a small perturbation to the electrode location will affect a large number of mesh elements exactly where the mesh is most sensitive to numerical noise.

Improved perturbations
To improve compute time for the perturbation method, we introduced an alternate 'minimal' perturbation strategy (figure 4(b)). Rather than moving all nodes associated with an electrode, only nodes at the boundary of the electrode were shifted. The two solutions, naïve and minimal perturbation, are analytically equivalent because the boundary conditions are consistent. None the less, numerical differences exist due to the order of operations: at which point the problem is converted from a continuous domain to the discrete one. The choice of when and how to discretise implicitly selects a subset of the possible solutions. The minimal perturbation method should suppress instabilities in the naïve perturbation method because fewer elements of the difference in system matrices are subject to non-cancelling floating point errors.
The minimal perturbation method may be extended to handle large electrode movements (figure 5). Assuming sufficient element refinement over the surface, the method may be applied by selecting nodes on the mesh for a new electrode location. The new electrode's boundary nodes are then perturbed to match the exact electrode location and diameter. First, a new electrode location x is identified. The local node spacing nearby is measured h and all nodes within the electrode diameter plus a margin + d h 2 e / are assigned to replace the former CEM nodes. The new electrode's boundary is identified and nodes on that boundary are perturbed radially with respect to the electrode centre so that the electrode boundaries exactly match the prescribed shape and location. For iterative solutions, perturbations start with a common mesh so that mesh quality is not successively degraded by perturbations.
Finally, by identifying the elements affected by nodal perturbations and CEM connectivity changes, the system matrix A may be updated by only recalculating modified elements in S e ( ) and reusing unmodified elements. The system matrix may then be reassembled at little relative cost. Typically, we observed a speed-up of between two and six times the naïve perturbation approach.

Rank-one update perturbation
One technique for reducing the quantity of calculations required for a mesh perturbation is to develop a rank-one matrix update to the FEM system matrix. The rank-one matrix update for are those elements connected to the perturbed nodes which must be updated in the FEM system matrix along with CEM connectivity and contact impedance distribution. electrode position (Gómez-Laberge and Adler 2008) has generally been used for single-step Gauss-Newton solutions where the electrode movement is small; typically less than 1% of electrode spacing.
The rank-one update is an application of the Sherman-Morrison formula (Sherman and Morrison 1950) which computes the inverse sum of an invertible matrix A (our FEM system matrix) and an outer product uv T : the movement perturbation. The particular node to perturb is selected by setting u v , as zero column vectors and non-zero at the row u i and column v j of the matrix entry (i,j) that is to be perturbed. The determinant of the perturbed matrix A is given by the well known matrix determinant lemma (Ding 2007) for a perturbation uv T . The FEM matrices are calculated for a system matrix A decomposed into a connectivity matrix C which associates local and global node numbering in the element mesh, the element shape functions S, and the conductivity per element D. The conductivity Jacobian σ J for element n is calculated as a difference in potentials T , given potential distributions X (Yorkey et al 1987, Adler and. We note that T is an operator which selects difference measurements from all calculated forward solutions and returns them as a vector. For a single stimulus, T can be constructed as a linear matrix followed by a vectorization, but this does not generalize unless all stimulus patterns use the same sequence of measurements. When the shape functions are first-order piecewise linear functions of potential with piece-wise constant conductivity, they can be calculated in two dimensions (n D = 2) as where S is a block diagonal matrix with each element e defined by its node locations S e ( ) , E 1 \ is the matrix E with the first row removed, and | | E det gives the area of the 2D element. The variables in the matrix E are illustrated by an example of a triangular element (blue) with nodes located at = p p p , ) . Applied to electrode position in impedance imaging, the process of Gómez-Laberge and Adler (2008) is to use a decomposition of the FEM system matrix A to determine the impact of perturbing a node. The electrode position Jacobian requires the application of the product of derivatives rule where x n refers to a global node numbered n and affects all element shape functions S e ( ) connected to that node so that some number of local nodes p involved in the matrices E are perturbed. The partial derivatives of the first-order interpolatory shape function may then be approximated using the rank-one perturbations vectors u and v to select the row and column to manipulate by a small perturbation (Gómez-Laberge and Adler 2008).
The rank-one update avoids the need to reassemble an entire FEM system matrix and recalculate the matrix inverse in the forward solution for each electrode position. The previously described 'minimal perturbation' technique might be applied to reduce the number of nodes moved in the rank-one update. In practice, the benefit of a minimal perturbation is small because the rank-one update is already computationally inexpensive and numerically stable for perturbations involving a small number of nodes.

Fréchet derivative for tangential movement
The Fréchet derivative for tangential electrode movement (Dardé et al 2012) accounts for contact impedance z c effects when using the CEM for calculating the electrode position Jacobian in a manner similar to the adjoint method for calculating the conductivity Jacobian.
The adjoint method for conductivity takes the dot product of the stimulus and measurement vector fields to calculate the conductivity Jacobian (Wang et al 2001, Borsic et al 2012. The stimulus field E m stim ( ) is the current vector for each element N j of the FEM mesh when a current is applied across stimulus electrodes used for measurement m according to a stimulus and measurement protocol. The measurement field E m meas ( ) is found by swapping measurement and stimulus electrodes to create a new stimulus driven with unit current and re-solving the forward problem for each measurement.
The adjoint method is nearly always a better choice when available; it is more efficient to calculate and more numerically stable than perturbation methods. Avoiding perturbations to the forward system matrix A supports the reuse of the computationally expensive decomposition (QR, LU, SVD).
The Fréchet derivative for electrode position takes a similar form as (21); an integral over the dot product of two fields constructed from the stimulus field E m stim ( ) and the field due to measurement electrodes used as stimulus E m meas ( ) . Examining the region near each electrode we observe that each measurement m is calculated as the difference between potentials at two electrodes. For the CEM, each electrode has a measured voltage U and a varying potential u along the contact surface with the domain (figure 6(a)). Between these two potentials, lies a distributed contact impedance layer which allows current to flow along the electrode. In calculating the electrode position Jacobian, we are interested in the change in the measured electrode voltage δU assuming the reference electrode's measured voltage is unaffected by the primary electrode's movement.
For a small electrode movement h, the outward tangential component υ ∂E of the movement Jacobian J x (the component normal to the edge of the electrode along its circumference) may be calculated according to the main result derived in (Dardé et al 2012), equation (11) where the measured voltage V and varying potential along the contact surface v are calculated when used as a stimulus electrode ( figure 6(b)). The 2D electrode is of length and contact impedance z c , with contact impedance units of [Ω.m]. We note that our statement of this equation differs from that presented in Dardé et al (2012) by adding the 1/ denominator term, where the units of contact impedance were not explicitly stated in the prior work. We find that this restatement gives consistent results for large contact impedances between the various methods while, at the same time, conforming to our definition of the contact impedance units used elsewhere in this work. Integrating over the surface of the electrode ∂E gives the electrode position Jacobian. The implementation of the Fréchet derivative for electrode position in two dimensions is particularly straightforward. For two-dimensional models, point evaluation at the edges of the electrode gives where electrode potential when used as a measurement electrode U and stimulus electrode V are combined with the potential at each edge of the electrode. The potential under the electrode when used as a measurement electrode, at both left u 1 and right u 2 edges of the electrode are combined with those as if used as a stimulus electrode, at both left v 1 and right v 2 edges of the electrode. The potentials v v V , , were calculated by reversing the stimulus and measurement electrodes in the stimulus and measurement pattern and solving the forward problem again. In three dimensions, the Jacobian calculation involves the line integral along the electrode edge. These forward solutions would be available without additional calculations for a 'complete' stimulus sequence where all combinations of stimulus electrodes were already calculated.
We implement the Fréchet derivative for tangential movement in two dimensions and compare it to the three perturbation methods: naïve perturbation, minimal perturbation, and the matrix rank-one update.

Methods
A rectangular two-dimensional model with four CEM electrodes on its upper surface was constructed in electrical impedance and diffuse optics reconstruction software (EIDORS) using NetGen (figure 7). Many models were generated that conformed to the same geometry but with variations in the electrode diameter, mesh density and contact impedance. Mesh density, electrode diameter and contact impedance were varied over orders of magnitude. Simulations were run by generating a model for each parameter in the range: mesh density, contact impedance, and electrode diameter. A single parameter was varied over its range while the other two parameters were held fixed. On each model, the electrode position Jacobian was calculated for each electrode using our implementation of the naïve perturbation, minimal perturbation, rank-one update, and Fréchet derivative. The results were plotted to demonstrate the points of agreement and illustrate numerical instabilities between methods.
The model's conductivity was set to unity (1 Ω · m). The size of the model was set large enough so that expanding the sides or bottom did not significantly affect the measurements. An alternate approach would have been to use a mixed boundary condition that approximates the homogeneous conductivity in the ±x and −z directions to infinity. The expanded domain is substantially easier to implement and achieves the same result. The four electrodes were spaced equidistant at 1 m intervals centre-to-centre. The electrode diameter d e was varied from 0.9 m (nearly touching) down to 2 mm. Contact impedance z c was varied between 10 −14 Ω · m to 10 +14 Ω · m. Mesh density, as measured by the maximum element height h max , was varied between 0.10 m and 1 m. Node density at the electrodes was scaled in proportion to the overall node density ( figure 8).
The geometry and electrode parameters were selected so that the results should be general to any equally spaced linear electrode array. For example, contact impedance effects are a function of the ratio of background resistivity to contact impedance.
In figure 7 , current flow density is indicated by stream lines (red) and voltage potential by the background colour map (blue-white-yellow). Inset in figure 7, are a close up of the positive stimulus electrode s+ (upper inset), and measurement electrode m+ (lower inset). The diversion of current flow through the measurement electrode can be observed for low electrode contact impedances in the lower inset.

Simulations
Plots illustrating variations in the Jacobian with respect to contact impedance (figure 9(a)), mesh density (figure 9(b)) and electrode diameter (figure 9(c)) follow. Plots show the naïve perturbation, minimal perturbation, matrix rank-one update, and Fréchet derivative for the positive and negative measurement electrodes (m+, m−). The dashed vertical line (red) indicates the configuration when a variable was fixed for other plots in figure 9.
Similar plots for all electrodes were observed. We observed that the measurement electrode m− plot's Jacobian were flipped vertically ( , we can see that this is the expected behaviour. Given a difference measurement between two electrodes (m+, m−), in a fixed and smoothly varying field caused by stimulus electrodes (s+, s−), we can consider movement of a single electrode. Moving a measurement electrode m+ in the +x direction reduces the distance between the measurement electrodes, leading to a reduced difference measurement. Conversely, moving the other measurement electrode m− in the +x direction increases the distance between measurement electrodes, increasing the difference in measured potential. A similar thought experiment for the stimulus electrodes leads to the conclusion that these plots show the generally expected trends for electrode movement. The measurement electrodes in this model are separated by a third of the distance of the stimulus electrodes, so that the scaling between calculated stimulus and measurement Jacobians is approximately correct as well. Our two independently developed implementations of the Fréchet derivative, based on Dardé et al (2012), were found to give the same solution.
For large contact impedances, all methods performed well and gave similar results ( figure  9(a)). For small contact impedances, the naïve perturbation method became unstable for our selected perturbation step (∆ = − 10 6 m). Changing the perturbation step size modified the threshold at which the naïve perturbation became unstable but never completely removed the instability (figure 10). This unavoidable instability was expected because there are practical upper and lower limits on the perturbation step size. The mesh density and electrode diameter were varied but these did not noticeably affect the plot of contact impedance versus perturbation size. For large step sizes, the 'minimal perturbation for large steps' method gave 'true' results to the limit where adjacent electrodes collided because the technique limits mesh element deformations to half an element's area. The rank-one update and naïve methods failed when nodal perturbations moved nodes beyond the enclosing elements leading to incorrect FEM meshes and an erroneous Jacobian. The rank-one perturbation is accurate when contact impedance was large relative to the conductivity of the adjacent domain because, in the limit ∞ z lim c ( → ), the contact impedance does not influence current flow through the electrode. On the other hand, the rank-one perturbation lacks derivatives of the contact impedance with respect to position, as the contact impedance and electrode shape are key components of the CEM, and these may become significant at some threshold. The missing derivatives from the CEM may be one source of small discrepancies in the naïve and minimal perturbations when compared to the rank-one perturbation results for small contact impedances. Finite precision numeric accuracy will set a lower limit on contact impedance at which point a shunt electrode model should be adopted. Given that the electrode shape was being held constant for the implementation of the rank-one perturbation used here, we believe these effects should be small.
For small contact impedances, we observed that the Fréchet derivative tends to zero. In the context of (23), this behaviour makes sense. The method uses the difference in potential across the electrode to calculate the effect of electrode movement which is subject to the limits of finite precision calculations in the FEM model of the voltage distribution. A contact impedance that tends towards zero will have a nearly zero potential difference across the electrode as it approaches a shunt electrode model. It turns out, in these simulations, that the Fréchet derivative will then tend to zero as the potential across the electrode drops more rapidly than the contact impedance. A similar reduction in potential difference across the electrode would occur for very small electrodes or weak stimulus. Indeed, for PEM electrodes the movement Jacobian breaks down under the Fréchet derivative , pp 103-5). (In the continuum model Fréchet with a Dirac delta input current, the sensitivity to normal perturbations are unbounded on driven electrodes and zero elsewhere. This indicates that there would likely be numerical problems if one tried to implement a PEM normal electrode position Jacobian.) We note that this behaviour is at odds with the perturbation method results.
The Fréchet derivative approaching zero for small contact impedances is inaccurate: the other methods presented here can correctly estimate an electrode position Jacobian under these conditions while the Fréchet derivative fails. This is somewhat disappointing, as the original proposition for the Fréchet derivative, beyond its computational efficiency, is that it accurately captures the CEM effects which are much more apparent for low contact impedances. The other methods succeed under low contact impedance conditions by sampling the rate of change in the electric fields beyond the electrode's boundary. This suggests a relatively simple idea for 'fixing' the Fréchet derivative: to sample just beyond the electrode's boundary. Clearly, this would have deep implications on the mathematical derivation of the Fréchet derivative: modifications are not considered further in this work.
All methods exhibit the same stable and uniform behaviour for the majority of variations in electrode diameter ( figure 9(b)). For electrode diameters where the electrodes nearly touch, the absolute magnitude of the Jacobian grows. The case where electrodes nearly cover the boundary was not explored further, as all methods seemed to be in strong agreement.
For very coarse meshes, the Fréchet derivative was found to underestimate the Jacobian by some form of systematic but noisy offset. This noise was possibly due to the sensitivity of the Fréchet derivative to errors in the FEM estimated voltages under the electrodes. The other perturbation-based methods rely solely on the change in measured electrode voltage, rather than a difference between measured voltage and voltage under the electrode. The naïve perturbation, minimal perturbation and rank-one matrix update gave consistent results for all tested mesh densities. We note that the sign was correct in all cases: a line search could correct for any error introduced by the Fréchet derivative method.
Calculation of the Jacobian can take a considerable fraction of the calculation total time for a single-step Gauss-Newton conductivity solution (Boyle et al 2012b). Fast calculation of the electrode position Jacobian can improve productivity in algorithm development and throughput for Gauss-Newton iterations where the Jacobian is updated at each iteration. For a particular forward model ( , the average calculation times for each electrode position Jacobian implementation were measured in Matlab 5 (Intel Core i5-2500K, 3.30 GHz, 32 GB mem) (figure 11). Results are shown with and without EIDORS caching enabled. The EIDORS caching feature was not disabled because functions were called multiple times for certain operations, knowing that the cache will filter out most of the Figure 11. Electrode position Jacobian compute time for tangential movement; (a) total run times compared between the different methods with and without function memoization (caching), (b) break down of timing where, for the rank-one update and Fréchet derivative, the clear boxes indicate a 'sunk cost' that is already available if caching is enabled since the system matrix has not been perturbed. computational cost of these calls. EIDORS uses a function memoization implementation where a cached function result is returned when function inputs match a cache entry (Michie 1968). To control the cache hit-rate, the ground node in the forward model was perturbed slightly so that the top-level function memoization would initially be defeated. Assembly and a single forward solution were not included in the timings, as these would already have been calculated in a Gauss-Newton iteration, and were considered a sunk cost. Subsequent calls to memoized results within a method benefited from the efficiency of a cached result if the forward model was not modified. When the cache was enabled or disabled, results showed that the Fréchet derivative method was fastest, while the rank-one update method was a close second.
Observing the profiling results, it was apparent that the vast majority of the computational time for the perturbation Jacobian (up to 96%) was consumed in recomputing the element shape functions E and then reassembling the system matrix A. With EIDORS' memoized caching mechanism, the entire system matrix was recomputed when a nodal perturbation was detected. For the specific case of computing the electrode position Jacobian, a finer grained caching of the inverted element matrices E, prior to being combined into the system matrix, considerably sped up the nodal perturbation methods because a handful of element matrices were recomputed. These improvements were implemented for the 'min. perturb' method.
We note that the rank-one update, essentially, performs this optimization by only modifying elements of the system matrix that are directly connected to the perturbed node. The Fréchet derivative method avoids the system matrix computation costs altogether by not modifying the system matrices. It is possible, for a 'complete' stimulus set where all possible combinations of stimulus electrodes are used, that no further forward solutions are required than have already been computed to that point in the Gauss-Newton algorithm. Such a computational saving is hard to quantify in a limited test bench such as the one presented here, but could reduce the incremental cost of calculating a electrode position Jacobian to a few multiplyand-accumulate operations if the forward solution nodal voltages are already available from the conductivity Jacobian.

Discussion
Boundary movement or imprecise electrode placement results in modelling errors that can cause significant reconstructed conductivity artifacts. Modelling boundary and electrode position estimates to address modelling errors and dynamic movement can correct some artifacts by incorporating electrode position as inverse problem parameters to be resolved. The Jacobian calculation contributes to determining what parameter changes may be beneficial in minimizing a cost function.
In this work, four methods for calculating the electrode position Jacobian were compared through simulation on a simplified two-dimensional homogeneous half-space. These methods were: the perturbation method, the rank-one matrix update (Gómez-Laberge and Adler 2008), and the Fréchet derivative (Dardé et al 2012). Simulations demonstrated that: • Variations in electrode diameter gave stable results for all methods up to diameters where adjacent electrodes nearly touched. • The perturbation method became unstable for small contact impedances depending on perturbation magnitude. • The Fréchet derivative method was sensitive to coarse meshes and was not appropriate for small contact impedances. • The Fréchet derivative and rank-one update methods were very fast, especially when accounting for previously solved forward problems in an iterative solution.
Based on these findings we recommend the Fréchet derivative for calculating the electrode position Jacobian when contact impedances are greater than the surrounding resistivity distribution and the rank-one matrix update otherwise. When using the Fréchet derivative we recommend evaluation at two or more mesh densities to confirm that the Jacobian solution has converged. We observed that from a computational point of view, the Fréchet derivative is most efficient, particularly when forward solutions are already available from calculating the conductivity Jacobian via the adjoint method or due to forward solutions previously calculated in an iterative reconstruction. The rank-one update is more efficient than the perturbation methods, but slower than the Fréchet derivative.
For conductivity changes it is known that the formal Fréchet derivative and discretisation operations commute as demonstrated in , p 63). Thus for conductivity changes one can use either formulation to estimate the Jacobian matrix with the same accuracy, with the preferred choice, for a given measurement protocol and discretisation, being the one with minimal computational cost. For boundary shape and electrode position changes we have no formal proof that the derivative of discretisation (rank-one perturbation method) and discretisation of derivative (Fréchet derivative method) agree for a given discretisation. Indeed the results demonstrated in this paper suggest that the two methods diverge as contact impedance approaches zero z 0 c ( → ). In the absence of a formal understanding, we recommend the derivative of discretisation (rank-one perturbation method) for synthetic (simulation) data, and the discretisation of the Fréchet derivative (Fréchet derivative method) for measured data. If one has access to voltage data from a simulated boundary shape perturbation (synthetic data), this implies working at a given discretisation level, and it seems plausible that the rank-one perturbation method would be most accurate. If one has access to voltage data from a continuum boundary shape deformation (as with real data) it seems appropriate to use the Fréchet derivative method to accurately model the continuum deformation.
If a large change in electrode displacement is required to explain the data, a perturbationbased update will fail: overlapping element shape functions may be calculated. Overlapping element shape functions may lead to invalid FEM solutions. A superficial 'fix' is to apply greater regularization to restrict the step size but this does not help achieve an accurate answer if larger movements are required to explain the measurements. In short, Gauss-Newton iterative solutions may make use of the rank-one update in calculating the Jacobian but will need an alternate strategy to safely update the FEM system matrix when movements exceed a fraction of the minimum element height of the mesh. Gauss-Newton updates and line search test points require an updated forward model where the electrode locations are updated. For tangential movements, we recommend a method where the new electrode location would be used to identify candidate nodes on the FEM mesh. Nodes from this set that are on the boundary of the electrode would then be perturbed to align them exactly with the new electrode location. This method is preferable to a mechanism where the whole mesh is distorted to push an electrode to its new location because it maintains the quality of the mesh elements, as measured by in-circle to out-circle ratio or other mesh quality criteria (Shewchuk 1996, Knupp 2003, Pébay and Baker 2003, which is important to the final FEM solution error (Berzins 1999). Potentially, the system matrix may be updated more efficiently than a complete recalcul ation by only updating shape sub-matrices S e ( ) containing nodes that were perturbed and any electrode-to-node connectivity.
Removing the mesh node movements from the perturbation in calculating an electrode position Jacobian is only possible when the electrode movement precisely agree with existing node locations. This is not generally true, but for a fine enough mesh provides a means of testing the theory that perturbations are the cause of numerical instability for the naïve perturbation method in contrast to the improved stability of the minimal perturbation method. Experiments applying a 'no node perturbation' change precisely agree with the minimal perturbation method over a range of mesh densities when the step change agrees with node placements. These experiments show that the perturbation of the electrode edge nodes is not significantly affecting the numerical accuracy of the solutions: a better perturbation strategy is not likely to be possible.
Unlike the conductivity Jacobian, the electrode position Jacobian may be nearly sparse for many common electrode configurations. Even for electrodes with very low contact impedance, effectively shunting current across their surface, the gaps between inactive electrodes must be small to have a measurable effect on the electrode position Jacobian. The sparsity of the electrode position Jacobian can be exploited to: reduce the number of calculations to build the Jacobian, reduce storage (for example, Compressed Sparse Row or Compressed Sparse Column storage), and expedite calculations using the Jacobian (for example, sparse matrix algebra).
Movement estimates tangential to the surface were faithfully calculated for many scenarios in this work. Using the same framework, movement normal to the surface was simulated with a naïve perturbation implementation, but appears to suffer from severe numerical instabilities that render the resulting Jacobian unusable in its current form. A Fréchet derivative for the normal component has also been developed (Dardé et al 2013) and offers the possibility of a stable Jacobian for normal boundary movement. The methods in this work can be extended to normal and tangential movement in three dimensions; preliminary results for three-dimensional tangential perturbation methods show similar outcomes.