Two Multigrid Algorithms for Inverse Problem in Electrical Impedance Tomography

    A direct partial differential problem involves an interior differential equation and a set of initial and/or boundary conditions which stably determines a unique solution. An inverse problem is one in which the differential equation and/or initial and/or boundary conditions are not fully given and instead the results of a set of solution observations (measurements) are known. They may contain errors, and even without errors the problem is usually ill-posed: the known data may be approximated by widely different solutions, but the ill-posedness of an inverse problem does not necessarily imply more expensive solution process. On the contrary: once the nature of the ill-posedness has been generally understood, to solve an inverse problem may even be much less expensive than to solve a corresponding direct problem. I am going to demonstrate this methodological point on an example of inverse problem for Electrical Impedance Tomography (EIT).

    An EIT device for medical use consists of a set of electrodes attached to the chest of a patient. A small known current is passed between two  driver electrodes. In each measurement a current is passed through a different electrode pair, while the voltage drops at all the electrodes are recorded. The collected data are used in order to compute the conductivity distribution in a part of the patient's chest and then to display it on a screen in order to detect anomalies, such as tumors. The electrical potential satisfies the equation Ñ×(sÑu)=0, where s is an electrical conductivity. The set of measurements gives ideally (in the limit of many small electrodes and as many measurements) the Neumann to Dirichlet mapping: the Dirichlet u boundary condition resulting from any Neumann ¶u/¶n condition. The inverse problem is to evaluate s from this mapping. The conductivity depends on the EIT data in a very weak way. Therefore the inverse problem of EIT is ill-posed.

    There exist some works on numerical methods for the relevant problems, but their number is rather sparse and even those papers do not consider the question of numerical efficiency, despite its importance for applications. In this specific inverse EIT problem, employing local Fourier decompositions, we have shown that all conductivity-function components of wavelenght l are ill-posed at distances r>>l from the boundary. Hence there is no need to use at such distances fine solution grids, since all we can know about the solution can be calculated with grids whose meshsizes increase proportionality to r.

    Two multigrid solvers for this problem are presented:

    The first multigrid solver uses the Tikhonov regularization technique. The inverse problem is reformulated as a variational minimization problem. The resulting Euler equations form a PDE system (for u, s and  Lagrange multipliers), which makes the problem suitable in principle for an effective numerical solution by multigrid methods. The multigrid solver starts with large and then procceds to progressively smaller regularization parameters. Many features of classical multigrid were adapted to this particular problem (such as intergrid communications, boundary condition treatment and coarse grid solution).

    In the case of a large regularization parameter, numerical results demonstrate a good convergence of the developed solver, but the obtained solution is too smeared and doesn't approximate the real conductivity too well. At small regularization the final approximation is much better, especially near the boundary. However, in this case the system is no longer elliptic and for efficiency the multigrid solver must use more sophisticated relaxation scheme, which effectively decomposes the system into its scalar factors. Although the multigrid cycles assymptotically slow down, the final approximation to the conductivity is obtained by just one multigrid cycle per grid refinement even for the smallest regularization for which solution still exists.

    In the second multigrid solver a regularization in the classical sense is replaced by a careful choice of grids. In order to avoid introducing ill-posed components into the current approximation, changes to the solution far from the measurements are calculated only on coarse enough grids (with meshsizes at least comparable to the distance from the measurements), which only define large-scale averages of the solution. On the other hand, near the boundary of measurements we introduce both high and low oscillatory components to the solution. For this we use semi-coarsening in the direction parallel to the boundary of measurements and full-coarsening in the perpendicular direction.

    The algorithm  is based on three main principles: 1) due to the ill-posedness, the computational work on the grid of meshsize h should only be done near the boundary, at distances O(h) from it; 2) similar changes introduced into s at the different O(h) distances from the boundary of measurements cause very different changes in u on this boundary Þ for smooth Fourier components in y a full coarsening (both in x and y) does not resolve these differences, hence we must use a semi coarsening in y direction (i.e. without coarsening the meshsize in the x direction); 3) it is necessary to start always from the better-posed changes Þ in each cycle (FMG, V-cycle, W-cycle, semi-cycle) we introduce changes to the current approximation first on the coarsest (or semi-coarsest) grid, where those changes are well-posed, and only after that changes from the finer grids.

    Numerical results show that the second algorithm approximates s better than the regularized method with its many artificial parameters (the regularization coefficients, which should change over the domain), for less work (no Lagrange multipliers). The accuracy of both algorithm, in particular near the boundary of measurements, is nearly the same. But the second algorithm does not use any regularization in a classical sense, and as a result significant space and time saving is achieved. In this algorithm there is no need to solve a large system of discretized differential equations and boundary conditions; also, no Lagrange-multiplier equations need to be treated.