A Multiscale Optimization Framework for Reconstructing Binary Images using Multilevel PCA-based Control Space Reduction
AA Multiscale Optimization Framework for BiomedicalApplications through Multilevel PCA-based ControlSpace Reduction
Priscilla M. Koolman
College of Engineering & Science, Florida Institute of Technology, Melbourne, FL 32901, USA
Vladislav Bukshtynov ∗ Department of Mathematical Sciences, Florida Institute of Technology, Melbourne, FL 32901, USA
Abstract
An efficient computational approach for optimal reconstructing parameters of binary-type physical properties for models in biomedical applications is developed and vali-dated. The methodology includes gradient-based multiscale optimization with mul-tilevel control space reduction by using principal component analysis (PCA) coupledwith dynamical control space upscaling. The reduced dimensional controls are usedinterchangeably at fine and coarse scales to accumulate the optimization progress andmitigate side effects at both scales. Flexibility is achieved through the proposed proce-dure for calibrating certain parameters to enhance the performance of the optimizationalgorithm. Reduced size of control spaces supplied with adjoint-based gradients ob-tained at both scales facilitate the application of this algorithm to models of highercomplexity and also to a broad range of problems in biomedical sciences. This tech-nique is shown to outperform regular gradient-based methods applied to fine scaleonly in terms of both qualities of binary images and computing time. Performance ofthe complete computational framework is tested in applications to 2D inverse prob-lems of cancer detection by the electrical impedance tomography (EIT). The resultsdemonstrate the efficient performance of the new method and its high potential forminimizing possibilities for false positive screening and improving the overall qualityof the EIT-based procedures.
Keywords:
PDE-constrained optimal control ◦ gradient-based method ◦ controlspace reduction ◦ multiscale parameter estimation ◦ principal component analysis ◦ electrical impedance tomography ◦ cancer detection problem ∗ Corresponding author: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] J u l Introduction
In this work, we propose and validate a computational approach for optimal reconstructingthe physical properties of the media based on any available, possibly incomplete and noisy,measurements. In particular, this approach is useful in various applications in biomedical sci-ences to operate with physical models characterized by near binary distributions observed forsome of their physical properties, e.g. heat or electrical conductivity. The proposed compu-tational framework is using gradient-based multiscale optimization techniques supplied withmultilevel control space reduction over both fine and coarse scales used interchangeably.Proper “communication” established between scales in terms of projecting the solution fromone scale onto another benefits in both quality and computational efficiency of the obtainedresults.As seen in many practical applications, fine scale optimization performed on fine meshesis able to provide high resolution images for the searched distributions. Fine meshes alsocontribute enormously towards increased sensitivity by enforcing accuracy in computing ad-joint states and constructing adjoint-based gradients if in use. The size of the control spacedefined over fine scales may be significantly decreased by applying any types of parame-terization, for instance by using linear transformations based on available sample solutions(realizations) when applying principal component analysis (PCA). However, fine scale op-timization may still suffer from over-parameterization if the problem is under-determined,i.e. the number of controls overweighs the size of available data (measurements). On theother hand, optimization performed on coarse meshes could arrive at a solution much fasterdue to the size of the control space. Usually, the solutions obtained at coarse meshes are ofa low quality and less accurate due to sensitivity naturally “coarsened”. In addition to thiscoarse scale optimization may suffer from being over-determined if the available data andthe size of the control space are not properly balanced.Various techniques of multiscale modeling have been used for decades with proven successfor numerous applications in computational mathematics, engineering and computer science.See [22,29,36] for a comprehensive description of this general area and an extensive literaturereview. There is also a recently growing interest in using multiscale techniques in applicationsto biological and biomedical sciences [5, 14, 31, 35] with limited involvement of methodsbroadly used in optimization and control theory.The proposed multiscale optimization framework utilizes all advantages mentioned abovewhile using fine and coarse scales. Moreover, using them both in one process helps usmitigate their side effects. For example, fine scale solution images may not provide clearboundaries between regions identified by different physical properties in space. As a result,a smooth transition cannot provide an accurate recognition of shapes, e.g. of cancer-affectedregions while solving an inverse problem of cancer detection (IPCD). In our computations,fine scale optimization is used to approximate the location of regions with high and lowvalues of a physical parameter, namely electrical conductivity. Projecting solutions onto thecoarse scale provides a dynamical (sharp-edge) filtering to the fine scale images optimizedto better match the available data. The filtered images then projected back onto the finescale preserving some information on recent changes obtained at the coarse scale. In fact,2e developed a computationally efficient procedure for automated scale shifting in order toaccumulate optimally progress obtained at both scales. At the extent of how the boundariesof the cancerous spots are recovered by projections between scales with assigned controls atthe coarse scale, this approach may be related to a group of level–set methods which utilizemultiscale techniques and adaptive grids [11,17,23,25,27,30,32]. We also use some notations,main ideas and governing principles of multiscale parameter estimation (MPE), refer to[15, 18, 19, 23] for some details. In the current paper, we keep the main focus on applyingour new computational approach to IPCD by the Electrical Impedance Tomography (EIT)technique, however, this methodology could be easily extended to a broad range of problemsin biomedical sciences, also in physics, geology, chemistry, etc.EIT is a rapidly developing non-invasive imaging technique gaining popularity by enablingvarious medical applications to perform screening for cancer detection [1, 3, 8, 21, 24, 33, 37].A well-known fact is that the electrical properties, e.g. electrical conductivity or permit-tivity, of different tissues are changing if the status of a tissue is changing from healthyto cancer-affected. This physical phenomenon allows EIT to produce images of biologicaltissues by interpreting their response to applied voltages or injected currents [8, 12, 21]. Amathematical model for solving different types of EIT problems and performing both an-alytical and computational analysis of such solutions was suggested in [13]. The inverseEIT problem deals with reconstructing the electrical conductivity by measuring voltages orcurrents at electrodes placed on the surface of a test volume. This so-called Calderon typeinverse problem [10] is highly ill-posed, refer to topical review paper [7]. Since 1980s variouscomputational techniques have been suggested to solve this inverse problem computation-ally. Recent paper [4] reviews the current state of the art and the existing open problemsassociated with EIT and its applications.This paper proceeds as follows. In Section 2 we present the mathematical description ofthe inverse EIT as an optimal control problem to be solved at both fine and coarse scales byapplying control space reduction using PCA (fine scale) and upscaling via dynamical parti-tioning (coarse scale). Procedures for performing multiscale optimization with interchangingfine and coarse phases are discussed in Section 3. Model descriptions and detailed compu-tational results including optimization parameter calibration are presented in Section 4.Concluding remarks are provided in Section 5.
In the recent paper [2] the inverse EIT problem is formulated as a PDE-constrained optimalcontrol problem in the Besov spaces framework for which the Fr´echet gradient and optimalityconditions are derived. The authors also developed the projective gradient method in Besovspaces and provided extensive numerical analysis for 2D models by implementing PCA-based solution space re-parameterization and Tikhonov-type regularization. In our currentdiscussion of the inverse EIT model we use the same notations as established in [2].3et Ω ⊂ R n , n = 2 ,
3, be an open and bounded set representing body and we assume thatfunction σ ( x ) : Ω → R + represents electrical conductivity at point x ∈ Ω. Electrodes ( E (cid:96) ) m(cid:96) =1 with contact impedances ( Z (cid:96) ) m(cid:96) =1 ∈ R m + are attached to the periphery of the body ∂ Ω. If theso-called “current–to–voltage” model is used, electrical currents ( I (cid:96) ) m(cid:96) =1 ∈ R m are appliedto the electrodes and induce constant voltages (electrical potentials) U = ( U (cid:96) ) m(cid:96) =1 ∈ R m onthe same electrodes. In this paper, however, we use the “voltage–to–current” model wherevoltages U (cid:96) are applied to electrodes E (cid:96) to initiate electrical currents I (cid:96) . In either model, itis assumed that both electrical currents and voltages satisfy the conservation of charge andground (zero potential) conditions, respectively m (cid:88) (cid:96) =1 I (cid:96) = 0 , m (cid:88) (cid:96) =1 U (cid:96) = 0 . (1)We formulate the inverse EIT (conductivity) problem [10] as a PDE-constrained optimalcontrol problem [2] by considering minimization of the following cost functional J ( σ ) = m (cid:88) (cid:96) =1 ( I (cid:96) − I ∗ (cid:96) ) , (2)where ( I ∗ (cid:96) ) m(cid:96) =1 ∈ R m are measurements made for electrical currents I (cid:96) . The latter may becomputed as I (cid:96) = (cid:90) E (cid:96) σ ( x ) ∂u ( x ) ∂n ds, (cid:96) = 1 , . . . , m (3)based on conductivity field σ ( x ) set here as a control variable. A distribution of electricalpotential u ( x ) : Ω → R is obtained as a solution of the elliptic problem ∇ · [ σ ( x ) ∇ u ( x )] = 0 , x ∈ Ω (4a) ∂u ( x ) ∂n = 0 , x ∈ ∂ Ω − m (cid:91) (cid:96) =1 E (cid:96) , (cid:96) = 1 , . . . , m (4b) u ( x ) + Z (cid:96) σ ( x ) ∂u ( x ) ∂n = U (cid:96) , x ∈ E (cid:96) , (cid:96) = 1 , . . . , m (4c)in which n is an external unit normal vector on ∂ Ω. Here we have to mention a well-knownfact that the inverse EIT problem to identify electrical conductivity σ ( x ) in discretizeddomain Ω with available input data ( I ∗ (cid:96) ) m(cid:96) =1 of size m is highly ill-posed. Therefore, weformulate an optimal control problem which is adapted to the situation when the size ofinput data can be increased through additional measurements while keeping the size of theunknown parameters, i.e. elements in the discretized description for σ ( x ), fixed. Followingthe discussion in [2] related to the “rotation scheme” we set U = U, I = I and consider m − U j = ( U j , . . . , U m , U , . . . , U j − ) , j = 2 , . . . , m (5)4pplied to electrodes E , E , . . . , E m respectively. Using the “voltage–to–current” modelallows us to measure associated currents I j ∗ = ( I j ∗ , . . . , I j ∗ m ). In addition, the total numberof available measurements could be further increased from m up to Km by applying (5)to K different permutations of potentials within set U . Having a new set of Km input data( I j ∗ ) Kmj =1 and in light of the Robin condition (4c) used together with (3), we now consider theoptimal control problem on minimization of the updated cost functional J ( σ ) = Km (cid:88) j =1 m (cid:88) (cid:96) =1 β j(cid:96) (cid:34)(cid:90) E (cid:96) U j(cid:96) − u j ( x ; σ ) Z (cid:96) ds − I j ∗ (cid:96) (cid:35) , (6)where each function u j ( · ; σ ) , j = 1 , . . . , Km , solves elliptic PDE problem (4a)–(4c). Addedweights β j(cid:96) in (6) in general allow setting the importance of measurement I j ∗ (cid:96) (when β j(cid:96) > β j(cid:96) = 0) from cost functional J computations. We alsocould note that the forward EIT problem (4a)–(4c) together with (3) may be used to generatevarious model examples (synthetic data) for inverse EIT problems to adequately mimiccancer related diagnoses seen in reality.Finally, the solution of the optimal control problemˆ σ ( x ) = argmin σ J ( σ ) (7)to minimize cost functional (6) subject to PDE constraint (4) could be obtained by thefirst-order optimality condition which requires the directional differential of cost functional δ J ( σ ; δσ ) to vanish for all perturbations δσ . By invoking the Riesz representation theorem [6]in L functional space δ J ( σ ; δσ ) = (cid:104) ∇ σ J , δσ (cid:105) L = (cid:90) Ω ∇ σ J δσ d Ω , (8)an iterative algorithm is proposed in [2] to solve problem (7) by means of cost functionaladjoint gradients (with respect to control σ ) ∇ σ J = − Km (cid:88) j =1 ∇ ψ j ( x ) · ∇ u j ( x ) (9)computed based on solutions ψ j ( · ; σ ) : Ω → R , j = 1 , . . . , Km , of the adjoint PDE problem ∇ · [ σ ( x ) ∇ ψ ( x )] = 0 , x ∈ Ω (10a) ∂ψ ( x ) ∂n = 0 , x ∈ ∂ Ω − m (cid:91) (cid:96) =1 E (cid:96) (10b) ψ ( x ) + Z (cid:96) ∂ψ ( x ) ∂n = 2 β (cid:96) (cid:20)(cid:90) E (cid:96) u ( x ) − U (cid:96) Z (cid:96) ds + I ∗ (cid:96) (cid:21) , x ∈ E (cid:96) , (cid:96) = 1 , . . . , m (10c)5 .2 Fine Scale: PCA-based Control Space Reduction A well-known problem in numerical optimization is that a spatially discretized form of theoptimal control problem discussed in Section 2.1 is over-parameterized even for small size2D models. To overcome ill-posedness due to over-parameterization of discretized σ ( x ) alongwith the fact that the obtained solutions should also honor any available prior information(such as available images, etc.), we implement re-parameterization of the control space basedon principal component analysis, also known as Proper Orthogonal Decomposition (POD)or Karhunen–Lo`eve (KL) expansion.More specifically, PCA enables us to represent control σ ( x ) in terms of uncorrelatedvariables (components of vector ξ ) mapping σ ( x ) and ξ by σ = Φ ξ + ¯ σ, (11a) ξ = ˆΦ − ( σ − ¯ σ ) , (11b)where Φ is the basis (linear transformation) matrix, ˆΦ − is the pseudo-inverse of Φ, and ¯ σ isthe prior mean. In the PCA used in our numerical experiments, the truncated singular valuedecomposition (TSVD) of a (centered) matrix, containing N r sample solutions (realizations)( σ ∗ n ) N r n =1 , as its columns, is used to construct the basis matrix Φ. The prior mean is given by¯ σ = (1 /N r ) (cid:80) N r n =1 σ ∗ n , see [2, 9] for details on constructing a complete PCA representation.The optimal control problem initially defined in Section 2.1 is now restated in terms of newmodel parameters ξ ∈ R N ξ used in place of control σ ( x ) as followsˆ ξ = argmin ξ J ( ξ ) (12)subject to discretized PDE model (4) and using control mapping (11). New gradients ∇ ξ J of cost functional J ( σ ) with respect to new control ξ can be expressed as ∇ ξ J = Φ T ∇ σ J (13)to define projection of gradients ∇ σ J obtained by (9) from initial (physical) σ -space ontothe reduced-dimensional ξ -space. An optimum in employing PCA-based control space reduction on a fine mesh discussed inSection 2.2 will be achieved by finding the minimal (optimal) size of new control ξ to enablehonoring prior information from available sample solutions [2, 9]. As often seen in practicalapplications, this optimal size cannot prevent the optimal control problem (12) from beingstill over-parameterized. Therefore, one will be interested in further re-parameterization byfinding a new control space defined by a reasonably small number of parameters, and thus,having fewer local minima.As a motivation for our new multiscale approach we used an idea of gradient-basedmultiscale (parameter) estimation (GBME) raised from the general MPE principles [15].6igure 1 illustrates the general concept of GBME. It employs various approaches for gradient-based refinement of the control space for dynamical space upscaling, i.e. control grouping,by analysis of changes in the gradient structure. For example, “noncompetitive” controlsidentified by relatively small components in the gradient, shown in red in Figures 1(a,b),are grouped into a new control. The associated (cumulative) component of the upscaledgradient, added as a red bar in Figure 1(c), makes the new control competitive and “visible”by other controls shown in blue and green. (a) old (fine scale) controls g r a d i e n t c o m p o n e n t s (b) new (coarse scale) controls up s c a l e d g r a d i e n t c o m p o n e n t s (c) Figure 1:
Schematic illustrating the general concept of gradient-based multiscale estima-tion (GBME). Plots in (a,b) show “noncompetitive” vs. “competitive” controls identified byassociated small (red bars) and big (blue and green bars) components of the gradient. Anew (cumulative) gradient component obtained via control grouping is depicted by the bigred bar in (c) with dimensions compounded by summing the respective dimensions of smallred bars in (b). For simplicity, all gradients shown in (a,b,c) have only positive components.Keeping aside for a while a principle by which spatial controls are grouped (partitioned) inour new approach, let us focus on obtaining upscaled gradients for new controls ( ζ j ) N ζ j =1 ∈ R N ζ + .We assume that control σ ( x ) in problem (7) is discretized over the fine mesh containing N elements and represented by a finite set of controls ( σ i ) Ni =1 ∈ R N + . This set is then partitionedinto N ζ subsets C j , j = 1 , . . . , N ζ , by selecting (without repetition) N j controls for j -thsubset and defining a map, i.e. fine–to-coarse partition, M : ( σ i ) Ni =1 → N ζ (cid:91) j =1 C j , C j = { σ i : P i,j = 1 , i = 1 , . . . , N } , N ζ (cid:88) j =1 | C j | = N ζ (cid:88) j =1 N j = N, (14)where the partition indicator function is defined as P i,j = (cid:40) , σ i ∈ C j , , σ i / ∈ C j . (15)7o proceed with gradients, we now consider discretized directional differential δ ˜ J , ob-tained from (8) by the first-order scheme, which is consistent with the discretized form ofdomain Ω decomposed into N spatial elements ( δ Ω i ) Ni =1 each of area (or volume in 3D) ∆ i δ J ( σ ; δσ ) ≈ δ ˜ J = N (cid:88) i =1 ∂ J ∂σ i ∆ i δσ i , (16)where δσ i perturbs controls σ i . Whenever control grouping is in place, we assume that allcontrols σ i within the same j -th group, σ i ∈ C j , are perturbed equally, i.e. δσ i = δζ j if P i,j = 1 for all i = 1 , . . . N and j = 1 , . . . N ζ . Then one could easily show that the spatialgrouping is fully consistent with the Riesz theorem δ ˜ J = N ζ (cid:88) j =1 N (cid:88) i =1 P i,j ∂ J ∂σ i ∆ i δσ i = N ζ (cid:88) j =1 ∂ J ∂ζ j δζ j = (cid:104) ∇ ζ J , δζ (cid:105) (17)and upscaled gradients ∇ ζ J are computed by summing up those components of discretizedgradients ∇ σ J related to controls σ i ∈ C j , i.e. ∂ J ∂ζ j = N (cid:88) i =1 P i,j ∂ J ∂σ i ∆ i . (18)In general, a gradient-based framework to perform optimization on multiple, fine and coarse,meshes will benefit from the following. • Forward simulations on fine N-element meshes allow constructing highly accurateadjoint-based gradients ∇ σ J . • Reasonably small number of controls N ζ (cid:28) N defined at coarse scales tends to lessenthe number of local minima. • Upscaling gradients at coarse scales by (18) allows dynamical control space relaxationwithout interrupting iterative optimization to follow changes in upscaling map (14).Figure 2 illustrates the general concept of such dynamical control space relaxation and ob-taining upscaled gradients for two new controls at a coarse scale. As shown in Figures 2(a,b),gradient components for all fine mesh controls may have the same order of magnitude. Thus,these controls, shown in red and blue, could be grouped following the idea which is differentfrom the analysis of the gradient structure used in GBME. We discuss this in detail as a newgrouping (partitioning) approach in Section 3.3.
The proposed approach for optimization utilizing multilevel control space reduction overmultiple scales, both fine and coarse, is motivated by a range of possible applications in8 upscaled control (a) fine scale controls g r a d i e n t c o m p o n e n t s (b) coarse scale controls up s c a l e d g r a d i e n t c o m p o n e n t s (c) Figure 2:
Schematic illustrating the general concept of dynamical control space relaxationand obtaining upscaled gradients at coarse scale by making 2-group assignments for fine scalecontrols, N ζ = 2. Plots in (a,b) show all controls on a fine mesh. Different colors are used toidentify controls to be included into (blue) upscaled control σ ( x ) in EIT. Fig-ure 3(left) shows an example of the histogram typical for such distributions.A common approach to solve optimal control problem (7) is to solve problem (12) insteadby applying PCA-based space reduction by mapping N -element (fine-mesh) σ -space andreduced-dimensional ξ -space. An optimal solution ˆ σ ( x ) obtained after applying map (11a)to ˆ ξ is of a Gaussian type. In case one of two modes is relatively small, a histogram for thesolution image is hardly recognized as being bimodal, e.g. as shown in Figure 3(middle), andpossible conversions to binary images may be very inaccurate.To provide a remedy and obtain the optimal solution ˆ σ ( x ) of a required binary typethe proposed approach employs multiscale optimization at both fine and coarse scales eachwith their own sets of controls by using them interchangeably. While various schemes areavailable for switching between scales, here we consider a simple one: scales are changedevery n s optimization iterations. We also define the coarse scale indicator function χ c ( k ) = (cid:40) , (2 k s − n s < k ≤ (2 k s − n s , (fine scale)1 , (2 k s − n s < k ≤ k s n s , (coarse scale) (19)where k s = 1 , , . . . and k = 0 , , , . . . denote respectively the counts for switching cyclesand optimization iterations. We choose to terminate the optimization run once the following9riterion is satisfied (cid:12)(cid:12)(cid:12)(cid:12) J ( σ k ) − J ( σ k − ) J ( σ k ) (cid:12)(cid:12)(cid:12)(cid:12) < (1 − χ c ) (cid:15) f + χ c (cid:15) c , k (cid:54) = k s n s + 1 (20)subject to chosen tolerances (cid:15) f , (cid:15) c ∈ R + . σ f r e qu e n c y “true” model: binary σ l σ h σ f r e qu e n c y fine scale optimization σ l σ h f r e qu e n c y coarse scale control σ σ l σ h coarse scale optimization coarse–to–fine fine–to–coarse projection projection Figure 3:
Schematic illustrating the general concept of the multiscale optimization frame-work. (left) A typical histogram representing binary distribution of true electrical conduc-tivity σ ( x ) used in EIT. In all three plots σ l and σ h values represent two modes associatedrespectively with healthy and cancer-affected regions within domain Ω. (middle) An exam-ple of the Gaussian-type histogram typical for solution σ k ( x ) obtained after k iterations ata fine scale. (right) A binary histogram for solution σ k ( x ) obtained at k -th iteration at acoarse scale. Positions of blue and red bars are associated with current values of σ klow and σ khigh controls, while their heights are computed based on the fine scale representation σ ( ξ k )cut off by the current value of the coarse scale threshold control σ kth . See Sections 3.2–3.3 fordetails. Coarse–to–fine and fine–to–coarse projections are defined respectively by (21)–(23)and (25)–(29).Our multiscale optimization framework is shown schematically in Figures 3(middle,right).Here, we would like to emphasize that all solutions obtained at both fine and coarse scalesare represented on (projected onto) the same N -element fine mesh. The word “multiscale”,in fact, refers to updates provided to discretized control σ ( x ). At a fine scale all elements σ i , i = 1 , . . . , N , are updated by applying PCA-based transformation, while at a coarse scalethese elements are sorted into two groups and then updated within each group by followingthe same rule. We discuss optimization phases at both scales as well as switching proceduresin the following two sections. 10 .2 Fine Scale Phase We denote the solution for control σ k = σ k ( x ) obtained at a fine scale at k -th iteration as σ ( ξ k ). During the fine scale optimization phase, χ c ( k ) = 0, control σ k is updated by solvingoptimization problem (12) in the reduced-dimensional ξ -space and by using map (11a) asdescribed in Section 2.2, i.e. σ k = σ ( ξ k ). Alternatively, during the coarse scale optimizationphase ( χ c ( k ) = 1) σ ( ξ k ), updated last time at the end of the fine scale phase, is used incoarse scale control grouping discussed in Section 3.3.To start the first switching cycle, k s = 1, and optimization itself, k = 0, the initial guess σ may be taken, for instance, corresponding to any approximate theoretical prediction σ ( x ).It is obvious that every time when coarse scale switches to the fine one, k = 2 k s n s , k s ≥ ξ k should be updated to ensure it receives as much as possible informationrelated to recent changes in σ k during the coarse scale phase. On the other hand, this updateshould not worsen the results σ ( ξ k ) previously obtained at the fine scale. Here, we proposethe following scheme to project solution σ k obtained at the end of the coarse scale phaseonto ξ -space by using a convex combination of σ ( ξ k ) and σ k ∀ k = 2 k s n s , k s ≥ σ ( ξ k ) = α c → f σ ( ξ k ) + (1 − α c → f ) σ k , α c → f ∈ [0 , . (21)Control ξ k then could be re-initialized from ¯ σ ( ξ k ) by using map (11b). As ¯ σ ( ξ k ) and σ k havedifferent distributions, namely of Gaussian and binary types, the coarse–to–fine projectionscheme (21) may benefit from projecting σ k first to its PCA equivalent σ kP CA = Φ ˆΦ − ( σ k − ¯ σ ) + ¯ σ, (22)see [9, 34] for details. σ kP CA then could be used in (21) in place of σ k . Also, an optimal choiceof relaxation parameter α c → f could be made by solving the following optimization problemin 1D α c → f = ˆ α = argmin α ≤ α ≤ J (¯ σ ( ξ k )) ≤ J ( σ ( ξ k )) (23)which appears to be highly nonlinear due to the inequality constraint to control the qualityof fine scale solutions σ ( ξ k ) in transition between subsequent switching cycles. To run optimization at a coarse scale we define a new 3-component control vector ζ = ( ζ j ) j =1 in which the first two entries are the low and high values of (binary) electrical conductivity σ ( x ) associated with healthy and cancer-affected regions in domain Ω, i.e. ζ = σ low , ζ = σ high . (24)They are shown schematically as respectively blue and red bars in Figure 3(right). Thethird component, ζ = σ th , takes responsibility for the shape of those regions (healthy and11ancer-affected) and is set as a separation threshold to define a boundary between low andhigh conductivity regions as shown in green in Figure 3(right). Such a simple structure ofcontrol ζ allows us to create a simplified representation of the coarse scale solution ζ k forcontrol σ k at k -th iteration based on the current fine scale representation σ ( ξ k ) = ( σ i ( ξ k )) Ni =1 σ ki = (cid:40) σ klow , σ i ( ξ k ) < σ kth ,σ khigh , σ i ( ξ k ) ≥ σ kth , i = 1 , . . . N, (25)where 0 < σ klow < σ khigh , min i σ i ( ξ k ) < σ kth < max i σ i ( ξ k ) , i = 1 , . . . N. (26)Simply, (25) provides a rule for creating fine–to–coarse partition M in (14) when N ζ = 2based on the current state of control ζ . During the coarse scale optimization phase, χ c ( k ) = 1,control σ k is updated by solving a 3D optimization problem in the ζ -spaceˆ ζ = argmin ζ J ( ζ ) (27)subject to constraints (bounds) provided in (26), and then σ k = σ ( ζ k ). When solvingproblem (27) during the first switching cycle, k = n s , ζ k could be initially approximated bysome constants, for example σ kth = 12 (cid:104) max i σ i ( ξ k ) + min i σ i ( ξ k ) (cid:105) ,σ klow = mean i (cid:8) σ i ( ξ k ) : σ i ( ξ k ) < σ kth (cid:9) ,σ khigh = mean i (cid:8) σ i ( ξ k ) : σ i ( ξ k ) ≥ σ kth (cid:9) , i = 1 , . . . N. (28)Switching from fine scale to coarse one when k = (2 k s − n s , k s >
1, could be even morestraightforward by utilizing the components of control ζ obtained at the previous coarse scalephase, i.e. ζ k = ζ k − n s . (29)In order to solve (27) by any approach which requires computing a gradient, its first twocomponents ∂ J ( ζ ) ∂ζ = ∂ J ∂σ low , ∂ J ( ζ ) ∂ζ = ∂ J ∂σ high could be easily obtained by using gradient summation formula (18) after completing parti-tioning map M (14)–(15) by employing (25). On the other hand, the third component maybe approximated by a finite difference scheme, e.g. of the first order, ∂ J ( ζ ) ∂ζ = ∂ J ∂σ th = J (cid:0) σ k ( ζ , ζ , ζ + δ ζ ) (cid:1) − J (cid:0) σ k ( ζ , ζ , ζ ) (cid:1) δ ζ + O ( δ ζ ) . (30)Parameter δ ζ in (30) is to be set experimentally pursuing trade-off between being reasonablysmall to ensure accuracy and large enough to protect numerator from being zero. In fact,formulas (25)–(29) provide a complete description of the fine–to–coarse projection for control σ ( x ) used in our approach. A summary of the complete computational scheme to performour new PCA-based multilevel optimization over multiple scales is provided in Algorithm 1.12 lgorithm 1 Computational workflow for PCA-based multiscale optimization k ← χ c ← σ ← initial guess σ ( x ) ξ ← σ by (11b) repeat u k ← σ k by solving forward problem (4) ψ k ← u k , σ k by solving adjoint problem (10) ∇ σ J ( σ k ) ← u k , ψ k by (9) if χ c = 1 then σ ( ξ k ) ← ξ k by (11a) ∇ ζ J ( ζ k ) ← ζ k , σ ( ξ k ) , ∇ σ J ( σ k ) by (18), also using (14)–(15) and (25), and (30) else ∇ ξ J ( ξ k ) ← ∇ σ J ( σ k ) by (13) end if update ζ k +1 and ξ k +1 by using, depending on χ c ( k ), descent directions D ζ ( ∇ ζ J ) or D ξ ( ∇ ξ J ) obtained respectively from ∇ ζ J ( ζ k ) or ∇ ξ J ( ξ k ) ζ k +1 = ζ k − χ c ( k ) τ k D ζ (cid:0) ∇ ζ J ( ζ k ) (cid:1) , (31a) ξ k +1 = ξ k − (1 − χ c ( k )) τ k D ξ (cid:0) ∇ ξ J ( ξ k ) (cid:1) (31b) if χ c = 1 then σ k +1 ← ζ k +1 , σ ( ξ k +1 ) by (25) else σ k +1 ← ξ k +1 by (11a) end if k ← k + 1 χ c ← k by (19) if χ c ( k ) (cid:54) = χ c ( k − thenif χ c = 1 then σ k ← ζ k , σ ( ξ k ) by (25) else ξ k ← σ k , σ ( ξ k ) by (21)–(23) σ k ← ξ k by (11a) end ifend ifuntil termination criterion (20) is satisfied to given tolerances (cid:15) f and (cid:15) c Main Results
Our optimization framework integrates computational facilities for solving forward PDEproblem (4), adjoint PDE problem (10), and evaluation of the gradients according to (9),(13), and (18). These facilities are incorporated mainly by using
FreeFem++ , see [20] fordetails, an open–source, high–level integrated development environment for obtaining nu-merical solutions of PDEs based on the Finite Element Method (FEM). For solving numer-ically forward PDE problem (4), spatial discretization is carried out by implementing FEMtriangular finite elements: P2 piecewise quadratic (continuous) representation for electricalpotential u ( x ) and P0 piecewise constant representation for conductivity field σ ( x ). Systemsof algebraic equations obtained after such discretization are solved with UMFPACK , a solver fornonsymmetric sparse linear systems [16]. The same technique is used for numerical solutionsof adjoint problem (10). All computations are performed using 2D domainΩ = (cid:8) x ∈ R : x + x < r (cid:9) (32)which is a disc of radius r Ω = 0 . m = 16 equidistant electrodes E (cid:96) with half-width w = 0 .
12 rad covering approximately 61% of boundary ∂ Ω as shown in Figure 4(a). Electricalpotentials U (cid:96) , see Figure 4(b), are applied to electrodes E (cid:96) as seen in (5) following the“rotation scheme” discussed in Section 2.1. We also consider adding up to three additionalpermutations within the set of potentials U , and, by choosing K ∈ { , , , } , we increasethe total number of measurements from m = 256 ( K = 1) to Km = 1024 ( K = K max =4). The potentials are chosen to be consistent with the ground potential condition (1).Determining the Robin part of the boundary conditions in (4c) we equally set the electrodecontact impedance Z (cid:96) = 0 .
1. Figure 4(c) also shows an example of the distribution of flux σ ( x ) ∇ u ( x ) of electrical potential u ( x ) in the interior of domain Ω during EIT.Physical domain Ω is discretized using mesh created by specifying 176 vertices overboundary ∂ Ω and totaling N = 7726 triangular FEM elements inside Ω. This fine meshis then used to construct gradients ∇ σ J , ∇ ξ J , and ∇ ζ J to perform the optimizationprocedure as described in Algorithm 1. To solve problems (27) and (12) iteratively as seenin (31), our framework is utilizing respectively the Steepest Descent (SD) and ConjugateGradient (CG) approaches [26] to obtain descent directions D ζ and D ξ . Stepsize parameters τ k in (31) are obtained by applying line minimization search [28].The actual (true) electrical conductivity σ true ( x ) we seek to reconstruct will be givenanalytically for each model by σ true ( x ) = (cid:40) σ c , x ∈ Ω c ,σ h , x ∈ Ω h , Ω = Ω c ∪ Ω h , Ω c ∩ Ω h = ∅ (33)and setting σ c = 0 . c (up to 4 spots of different size depending onthe model’s complexity) and σ h = 0 . h . In terms of the initial guessfor control σ ( x ) we take a constant approximation to σ true given by σ = ( σ h + σ c ) = 0 . E E E E E E E E E E E E E E E (a) (b) I I I I I I I I I I I I I I I I (c) Figure 4: (a) Equispaced geometry of electrodes E (cid:96) , (cid:96) = 1 , , . . . ,
16, placed over bound-ary ∂ Ω. (b) Electrical potentials U (cid:96) . (c) Electrical currents I (cid:96) (positive in red, negative inblue) measured at electrodes E (cid:96) . Black arrows show the distribution of flux σ ( x ) ∇ u ( x ) ofelectrical potential u ( x ) in the interior of domain Ω.In order to avoid early termination at a coarse scale, termination tolerances in (20) are setas (cid:15) c = 0 and (cid:15) f = 10 − .For all computations mentioned in the rest of Section 4 we use a PCA-based map (11) be-tween N -dimensional discretized control σ ( x ) and reduced-dimensional ξ -space as describedin Section 2.2. A set of N r = 1000 realizations ( σ ∗ n ) n =1 is created using a generator ofuniformly distributed random numbers. Each realization σ ∗ n “contains” from one to seven“cancer-affected” areas with σ c = 0 .
4. Each area is located randomly within domain Ω andrepresented by a circle of randomly chosen radius 0 < r ≤ . r Ω . To perform TSVD wechoose the number of principal components N ξ by retaining 662, 900, and 965 basis vec-tors in the PCA description. These values correspond to the preservation of respectively r ξ = 99%, 99 . .
99% of the “energy” in the full set of basis vectors, see [2, 9] fordetails.
We created our (benchmark) model I ∗ (cid:96) are recordedby running (4) and (3) with σ ( x ) = σ true represented by P2 FEM elements. Then, as men-15 .10.150.20.250.30.350.40.450.5 (a) (b) (c) Figure 5:
Model σ true ( x ). (b) A binary-type histogramfor model σ P CA ( x ) for model σ true ( x ) to its PCA equivalent with r ξ = 99% by using (22).tioned in Section 4.1, the actual reconstruction ˆ σ ( x ) is obtained in P0 finite element space.Taking also into account that the “measurement” data is not projected into its PCA equiv-alent, we compared cost functional J ( σ ) evaluated at σ true with added noise and at σ P CA ,as shown in Figure 5(c), obtained after applying PCA projection. The equivalent noise isestimated to be at level 0 . − . S = { K, r ξ , j max , n s } of 4 majorparameters to be calibrated: • K ∈ { , , , } to define the total number of measurements, Km , respectively as 256,512, 768, 1024, • r ξ ∈ { , . , . } to set the number of principal components, i.e. number ofcontrols N ξ at a fine scale, respectively to 662, 900, 965, • j max ∈ { , Km } to consider two cases for optimization at the coarse scale: β j(cid:96) = 1 , j =1 , . . . , Km , vs. β j(cid:96) = 0 except j = 1 to consider respectively full data vs. m = 16measurements to avoid over-parameterization, and • n s ∈ { , } frequency of switching between fine and coarse scales.Figure 6 shows the results obtained after performing the calibration procedure, i.e. afterrunning optimization for model S . We have chosentwo critical factors to evaluate the performance in each case: the L -norm difference betweenoptimal solution ˆ σ ( x ) and true conductivity field σ true ( x ) ( x -axis), and the absolute errorin recovered ˆ σ high by comparing it with the known value σ h = 0 . y -axis). We note thatfor all experiments ˆ σ low is recovered very accurately. Five outliers inside the dotted boxin Figure 6(a) are excluded from the entire calibration statistics. Then open circles inFigure 6(a) provide the results averaged over four cases, K = 1 , , ,
4, with K = 4 (red circle)appeared as the best one. Figures 6(b,c,d) show the same 43 outcomes (less 5 outliers) then16 K = 1, 256 measK = 2, 512 measK = 3, 768 measK = 4, 1024 meas (a) r = 99%, N = 662r = 99.9%, N = 900r = 99.99%, N = 965 (b) j max = Km, Km measj max = 1, m = 16 meas (c) n s = 5 itn s = 10 it (d) Figure 6:
Model K = 1 , , ,
4, (b) r ξ = 99% , . , . j max = 1 , Km , and (d) n s =5 ,
10. Open circles in (a) provide the results averaged over four cases K = 1 , , ,
4. Opencircles and triangles in (b,c,d) represent averages obtained respectively for K = 1 and K = 4.Five outliers shown inside the dotted box in (a) are removed from plots in (b,c,d). The resultobtained with parameter schedule S ∗ in (34) is shown by hexagons in all four plots.colored according to values of the rest parameters, namely r ξ , j max , and n s . Open circles andtriangles there represent averages obtained respectively for K = 1 and K = 4. We concludethat simple models, like our model S ∗ = { K = 4 , r ξ = 99% , j max = 1 , n s = 5 } . (34) (a) -6 -4 -2 (b) = low2 = high3 = th (c) (d) (e) -6 -4 -2 (f) Figure 7:
Model σ true ( x ) in Figure 5(a). Graphs in (b,f) presentnormalized cost functional J ( σ k ) / J ( σ ) as a function of iteration count k evaluated at fine(red dots) and coarse (blue dots) scales. Pink dots in (b) show values α c → f as solutions forproblem (23). Changes in the coarse scale controls ζ k = [ σ klow σ khigh σ kth ] are shown in (c) with σ c = 0 . σ h = 0 . S ∗ . Figure 7(a) shows the outcome of this reconstruction with a dottedcircle representing the location of the cancer-affected region taken from known σ true ( x ),see Figure 5(a). The location of this region is captured accurately. In addition, as seenin Figure 7(c), reconstructed values of both coarse scale controls ˆ ζ = ˆ σ low = 0 . ζ = ˆ σ high = 0 . ( σ k ) / J ( σ ) as a function of iteration count k evaluated at both fine (red dots) and coarse(blue dots) scales. Pink dots show values α c → f as solutions for problem (23). Furtheranalysis of changes in α c → f ( k ) aligned with the tracked behavior of J ( σ k ) may suggest moredevelopment for termination conditions to provide even better performance.Figures 7(d,e,f) show also the results of performing optimization for the same model σ h = 0 . (cid:15) f = 10 − .Finally, we conclude here that our multiscale computational framework, when properlycalibrated, is able to provide binary images consistent with the obtained measurements withsignificant reduction in computational time. We now present results obtained using our new multiscale optimization framework appliedto models with an increased level of complexity. The added complications are the numberof cancer-affected regions (more than one) and the variations in the size of those regions.Although additional calibration for optimization parameters might be seen as useful to obtainbetter results, we use the same parameter schedule S ∗ in (34) obtained using our benchmarkmodel σ true ( x ) for our model (cid:15) f = 10 − . As confirmed by the solution image and the associated histogram seenin Figures 8(b,c), we arrived at the same conclusion as for model α c → f (in pink) different from 1 identifies weightedprojections of the coarse scale solutions onto the fine scale performed with weight 1 − α c → f .Updated fine scale solutions then are used for constructing fine–to–coarse partitions M in(14) and (25) by changing values of the coarse scale threshold control ζ = σ th (in black).We also acknowledge the error in recovering high conductivity part ˆ ζ = ˆ σ high = 0 . σ true ( x ) is shown in Figure 9(a). As in the previous case, we runa fine scale only optimization terminated with (cid:15) f = 10 − which gives a solution image and19 .10.150.20.250.30.350.40.450.5 (a) (b) (c) (d) -5 -4 -3 -2 -1 (e) = low2 = high3 = th (f) Figure 8:
Model σ true ( x ). Outcomes for optimization:(b,c) when performed only at a fine scale, and (d,e,f) after applying multiscale framework byAlgorithm 1. Plots in (b,d) show the obtained images with added dotted circles to representthe location of four cancer-affected regions taken from known σ true ( x ). The fine scale solutionhistogram is presented in (c). Graph in (e) shows normalized cost functional J ( σ k ) / J ( σ ) asa function of iteration count k evaluated at fine (red dots) and coarse (blue dots) scales. Pinkdots show values α c → f as solutions for problem (23). Changes in the coarse scale controls ζ k = [ σ klow σ khigh σ kth ] are shown in (f) with σ c = 0 . σ h = 0 . .10.150.20.250.30.350.40.450.5 (a) (b) (c) (d) -5 -4 -3 -2 -1 (e) = low2 = high3 = th (f) Figure 9:
Model σ true ( x ) with added structure ofboundary electrodes (in black). Outcomes for optimization: (b,c) when performed only ata fine scale, and (d,e,f) after applying multiscale framework by Algorithm 1. Plots in (b,d)show the obtained images with added dotted circles to represent the location of three cancer-affected regions taken from known σ true ( x ). The fine scale solution histogram is presentedin (c). Graph in (e) shows normalized cost functional J ( σ k ) / J ( σ ) as a function of iterationcount k evaluated at fine (red dots) and coarse (blue dots) scales. Pink dots show values α c → f as solutions for problem (23). Changes in the coarse scale controls ζ k = [ σ klow σ khigh σ kth ]are shown in (f) with σ c = 0 . σ h = 0 . ζ , ˆ ζ = ˆ σ high = 0 . σ th as a single control for allthree regions may significantly contribute to amplify these errors.Our last model σ true ( x ) isshown in Figure 10(a). This model contains only one circular-shaped cancer-affected regionof the same size as the smallest region in model σ ( x ) = σ h , ∀ x ∈ Ω) is very close to the order of noise that appeared naturally in21 .10.150.20.250.30.350.40.450.5 (a) (b) (c) (d) -6 -4 -2 (e) = low2 = high3 = th (f) Figure 10:
Model σ true ( x ) with added structure ofboundary electrodes (in black). Outcomes for optimization: (b,c) when performed only ata fine scale, and (d,e,f) after applying multiscale framework by Algorithm 1. Plots in (b,d)show the obtained images with added dotted circles to represent the location of one smallcancer-affected region taken from known σ true ( x ). Note that the color bar in (d) is re-scaledfor better shape capturing. The fine scale solution histogram is presented in (c). Graph in (e)shows normalized cost functional J ( σ k ) / J ( σ ) as a function of iteration count k evaluatedat fine (red dots) and coarse (blue dots) scales. Pink dots show values α c → f as solutions forproblem (23). Changes in the coarse scale controls ζ k = [ σ klow σ khigh σ kth ] are shown in (f) with σ c = 0 . σ h = 0 . (cid:15) f = 10 − .After comparing fine scale and multiscale optimization images, it is true to say that thelatter could provide more assistance in concluding on possible abnormal changes in tissuesand navigating the surgeons. On the other hand, the fine scale image may be misleading assome (yellowish) regions may be deceptively interpreted as being cancerous. The multiscaleoptimization, however, allows keeping its images devoid of this problem and, as such, showshigh potential in minimizing possibilities for false positive screening and improving the overallquality of the EIT-based procedures. Similarly to our previous models, the same conclusion22rrives after observing continuous information exchange between fine and coarse scales seen inFigures 10(e,f) which results in creating an image with clear binary resolution ready to locatethe cancer-affected region. We also note that its size appears larger than in the true modeland the error in recovering the coarse scale control ˆ ζ = ˆ σ high = 0 . In this work, we presented an efficient computational approach for optimal reconstructingthe physical properties of the media characterized by distributions close to binary. In partic-ular, this approach could be useful in various applications in biomedical sciences to operatewith physical models supplied with some, possibly incomplete and noisy, measurements. Theproposed computational framework includes an array of gradient-based multiscale optimiza-tion techniques supplied with multilevel control space reduction over both fine and coarsescales used interchangeably. Quality and computational efficiency of the obtained results areensured by developing a methodology for establishing an effective “communication” betweenscales by projecting the solutions from one scale onto another and accumulating optimallyprogress obtained at both scales. Such multiscale optimization paired with multilevel controlspace reduction allows using computational advantages seen at both scales and to mitigatetheir negative impacts.We investigated the performance of our complete computational framework in applica-tions to the 2D inverse problem of cancer detection by the Electrical Impedance Tomographytechnique. Our first benchmark model mimics a simple biological tissue case with confirmedpresence of one circular-shaped area affected by cancer. The proposed procedure for cali-brating certain parameters was applied to this model to ensure the enhanced performance ofour optimization framework. We also presented results obtained by applying the calibratedframework to multiple models of an increased level of complexity, namely with three and forcancer-affected regions of various sizes. For every model, we obtained clear images with anice binary resolution enabled to locate all cancer-affected regions. In addition, our multi-scale optimization framework proved its high efficiency by completing computations 3 timesfaster than in cases when only fine scale was in use.We also check the applicability of our framework in applications to procedures for cancerrecognition at the early stages by a model containing one tiny cancerous spot with a diametercomparable with the size of the boundary electrodes. Despite the errors in recovering thetrue shape and the values of electrical conductivity, we conclude that obtained images ofthat quality will provide valuable assistance in recognizing possible abnormal changes intissues and further navigating medical professionals with their decisions. We conclude thatthe properly calibrated multiscale optimization framework is able to provide binary imagesconsistent with the provided measurements using significantly reduced computational time.23n general, we see a high potential of the proposed computational framework in minimizingpossibilities for false positive screening and improving the overall quality of the EIT-basedprocedures.There are many ways in which our multiscale optimization algorithm can be tested andextended. We provided an example of a simple calibration procedure, but we expect theperformance may be further improved by extending the list of calibrated parameters andapplying more advanced minimization techniques to perform local and global searches whileperforming optimization at both fine and coarse scales. Given that we used data providedby a specific electrode configuration, it will be of interest to apply a further analysis ofthe measurement structure, for example considering a 32-electrode scheme and improvingsensitivity by optimizing the structure of available data. We also plan to investigate theuse of flexible schemes for switching between scales including new approaches for projectingsolutions. The impact of the noise present in measurements should be also systematicallyanalyzed. Also of interest is the extension of our multiscale optimization approach to includepossibility of using various PCA sample structures, multiple coarse scale controls associatedwith different spatial regions, and applicability to bimodal distributions.Finally, it is important to test our new approach in various applications to real data anddifferent types of cancerous tissues, as this would certainly suggest areas in which furtherdevelopments may be required. Despite the fact that this approach was initially tested withsynthetic EIT-related problems, we believe that this methodology could be easily applied toa broad range of problems in biomedical sciences, also in physics, geology, chemistry, etc.
References [1] Abascal, J.F.P.J., Lionheart, W.R.B., Arridge, S.R., Schweiger, M., Atkinson, D.,Holder, D.S.: Electrical impedance tomography in anisotropic media with known eigen-vectors. Inverse Problems (6), 1–17 (2011)[2] Abdulla, U.G., Bukshtynov, V., Seif, S.: Breast cancer detection through electricalimpedance tomography and optimal control theory: Theoretical and computationalanalysis. arXiv:1809.05936 [3] Adler, A., Arnold, J., Bayford, R., Borsic, A., Brown, B., Dixon, P., Faes, T.J., Frerichs,I., Gagnon, H., G¨arber, Y., Grychtol, B., Hahn, G., Lionheart, W., Malik, A., Stocks, J.,Tizzard, A., Weiler, N., Wolf, G.: GREIT: towards a consensus EIT algorithm for lungimages. In: 9th EIT conference 2008, 16-18 June 2008, Dartmouth, New Hampshire.Citeseer (2008)[4] Adler, A., Gaburro, R., Lionheart, W.: Handbook of Mathematical Methods in Imaging,chap. Electrical Impedance Tomography, pp. 701–762. Springer New York, New York,NY (2015)[5] Alber, M., Tepole, A.B., Cannon, W.R., De, S., Dura-Bernal, S., Garikipati, K., Kar-niadakis, G., Lytton, W.W., Perdikaris, P., Petzold, L., Kuhl, E.: Integrating machine24earning and multiscale modeling – perspectives, challenges, and opportunities in thebiological, biomedical, and behavioral sciences. Digital Medicine (115) (2019)[6] Berger, M.S.: Nonlinearity and Functional Analysis. Acad. Press, New York (1977)[7] Borcea, L.: Electrical impedance tomography. Inverse Problems , 99–136 (2002)[8] Brown, B.: Electrical impedance tomography (EIT): A review. Journal of MedicalEngineering and Technology (3), 97–108 (2003)[9] Bukshtynov, V., Volkov, O., Durlofsky, L., Aziz, K.: Comprehensive framework forgradient-based optimization in closed-loop reservoir management. Computational Geo-sciences (4), 877–897 (2015)[10] Calderon, A.P.: On an inverse boundary value problem. In: Seminar on Numeri-cal Analysis and Its Applications to Continuum Physics, pp. 65–73. Soc. Brasileira deMathematica, Rio de Janeiro (1980)[11] Chen, W., Cheng, J., Lin, J., Wang, L.: A level set method to reconstruct the discon-tinuity of the conductivity in EIT. Science in China Series A: Mathematics , 29–44(2009)[12] Cheney, M., Isaacson, D., Newell, J.: Electrical impedance tomography. SIAM Review (1), 85–101 (1999)[13] Cheng, K.S., Isaacson, D., Newell, J., Gisser, D.G.: Electrode models for electric currentcomputed tomography. IEEE Transactions on Biomedical Engineering (9), 918–924(1989)[14] Clancy, C.E., An, G., Cannon, W.R., Liu, Y., May, E.E., Ortoleva, P., Popel, A.S.,Sluka, J.P., Su, J., Vicini, P., Zhou, X., Eckmann, D.M.: Multiscale modeling in theclinic: Drug design and development. Annals of Biomedical Engineering (9), 2591–2610 (2016)[15] Cominelli, A., Ferdinandi, F., De Montleau, P., Rossi, R.: Using gradients to refineparameterization in field-case history-matching projects. SPE Reservoir Evaluationand Engineering (3), 233–240 (2007)[16] Davis, T.A.: Algorithm 832: UMFPACK V4.3 – an unsymmetric-pattern multifrontalmethod. ACM Transactions on Mathematical Software (TOMS) (2), 196–199 (2004)[17] Gibou, F., Fedkiw, R., Osher, S.: A review of level-set methods and some recent appli-cations. Journal of Computational Physics , 82–109 (2018)[18] Grimstad, A.A., Mannseth, T.: Nonlinearity, scale, and sensitivity for parameter esti-mation problems. SIAM Journal on Scientific Computing (6), 2096–2113 (2000)2519] Grimstad, A.A., Mannseth, T., Nævdal, G., Urkedal, H.: Adaptive multiscale perme-ability estimation. Computational Geosciences , 1–25 (2003)[20] Hecht, F.: New development in FreeFem++. Journal of Numerical Mathematics (3-4), 251–265 (2012)[21] Holder, D.S.: Electrical Impedance Tomography. Methods, History and Applications.Taylor & Francis (2004)[22] Horstemeyer, M.F.: Practical Aspects of Computational Chemistry, chap. MultiscaleModeling: A Review. Springer (2009)[23] Lien, M., Berre, I., Mannseth, T.: Combined adaptive multiscale and level-set parameterestimation. Multiscale Modeling & Simulation (4), 1349–1372 (2005)[24] Lionheart, W.: EIT reconstruction algorithms: Pitfalls, challenges and recent develop-ments. Physiological Measurement (1), 125–142 (2004)[25] Liu, D., Khambampati, A.K., Du, J.: A parametric level set method for electricalimpedance tomography. IEEE Transactions on Medical Imaging (2), 451–460 (2018)[26] Nocedal, J., Wright, S.J.: Numerical Optimization, 2nd edn. Springer (2006)[27] Osher, S., Sethian, J.: Fronts propagating with curvature-dependent speed: Algorithmsbased on Hamilton-Jacobi formulations. Journal of Computational Physics (1), 12–49(1988)[28] Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.: Numerical Recipes:The Art of Scientific Computing, 3rd edn. Cambridge University Press (2007)[29] Steinhauser, M.: Computational Multiscale Modeling of Fluids and Solids: Theory andApplications, 2nd edn. Springer (2017)[30] Tai, X.C., Chan, T.: A survey on multiple level set methods with applications foridentifying piecewise constant functions. International Journal of Numerical Analysisand Modeling (1), 25–47 (2004)[31] Tawhai, M., Bischoff, J., Einstein, D., Erdemir, A., Guess, T., Reinbolt, J.: Multiscalemodeling in computational biomechanics. IEEE Engineering in Medicine and BiologyMagazine (3), 41–49 (2009)[32] Tsai, R., Osher, S.: Level set methods and their applications in image science. Com-munications in Mathematical Sciences (4), 1–20 (2003)[33] Uhlmann, G.: Electrical impedance tomography and Calder´on’s problem. Inverse Prob-lems (12), 123,011 (2009) 2634] Volkov, O., Bukshtynov, V., Durlofsky, L., Aziz, K.: Gradient-based Pareto optimalhistory matching for noisy data of multiple types. Computational Geosciences (6),1465–1485 (2018)[35] Walpole, J., Papin, J.A., Peirce, S.M.: Multiscale computational models of complexbiological systems. Annual Review of Biomedical Engineering , 137–154 (2013)[36] Weinan, E.: Principles of Multiscale Modeling. Cambridge University Press (2011)[37] Zou, Y., Guo, Z.: A review of electrical impedance techniques for breast cancer detec-tion. Medical Engineering and Physics25