A Subspace Preconditioned LSQR Gauss-Newton Method with a Constrained Line Search Path Applied to 3D Biomedical Microwave Imaging

Three contributions that can improve the performance of a Newton-type iterative quantitative microwave imaging algorithm in a biomedical context are proposed. (i) To speed up the iterative forward problem solution, we extrapolate the initial guess of the field from a few field solutions corresponding to previous source positions for the same complex permittivity (i.e., “marching on in source position”) as well as from a Born-type approximation that is computed from a field solution corresponding to one previous complex permittivity profile for the same source position. (ii)The regularized Gauss-Newton update system can be ill-conditioned; hence we propose to employ a two-level preconditioned iterative solution method. We apply the subspace preconditioned LSQR algorithm from Jacobsen et al. (2003) and we employ a 3D cosine basis. (iii) We propose a new constrained line search path in the Gauss-Newton optimization, which incorporates in a smooth manner lower and upper bounds on the object permittivity, such that these bounds never can be violated along the search path. Single-frequency reconstructions from bipolarized synthetic data are shown for various three-dimensional numerical biological phantoms, including a realistic breast phantom from the University of Wisconsin-Madison (UWCEM) online repository.


Introduction
During the last decade significant progress has been reported in fully 3D quantitative microwave imaging algorithms, for example, [1][2][3][4][5][6][7][8][9][10][11].The aim of such algorithms is to reconstruct a 3D spatially varying complex permittivity profile from measured scattered field data.Consequently they consist of a numerical inversion scheme which employs electromagnetic field equations that exactly account for the vectorial nature and scattering of the electromagnetic waves.Research efforts have focused in particular on biomedical applications [3,5,11], since it has been pointed out that the electromagnetic properties at microwave frequencies are sensitive to tissue types and various physiological and pathological parameters [12][13][14][15].The interest in these algorithms has been further stimulated by the possible benefits of microwave breast cancer screening as a supplemental diagnosis and/or therapy monitoring modality [4,[16][17][18][19][20]. Quantitative microwave imaging thus may open new perspectives as a moderate-cost, noninvasive, and nonionizing medical imaging technique.
It is well known that reconstructing the complex permittivity is a nonlinear and ill-posed inverse problem.The solution of this problem can be formulated as the minimizer of an appropriately regularized cost function that is composed of a data fit term and a regularization term or factor which incorporates a priori information about the object.Starting from an initial estimate, the permittivity profile is then updated progressively by means of an iterative optimization technique.Different approaches for choosing the cost function and the optimization technique have been proposed.Let us mention here the "modified gradient" type approach [3,4,6,21,22] and the "Newton-type" approach.In this paper we consider this last approach, where the cost function only depends on the complex permittivity and is optimized with a Newton-type technique.In this case a set of forward problems is solved in each iteration of the optimization.Whereas 2 International Journal of Antennas and Propagation the "modified gradient" approach avoids these computationally demanding forward problem solutions, the "Newtontype" approach needs far less iterations to converge.
In previous work we have reported reconstructions from experimental data with the Gauss-Newton method applied to various types of cost functions for a variety of 3D piecewisehomogeneous objects [8,23], including the 3D objects from the Fresnel database [24,25].In the present paper we deal with biological objects and we discuss three contributions that can improve the performance of a Newton-type algorithm with regard to the challenges posed by such objects.The reconstruction of a biological object is challenging when it is inhomogeneous with highly contrasted permittivities and when the number of reconstruction variables is large.Such conditions have an impact on the computational burden for the forward and update problems.In this paper we have chosen to implement the proposed contributions into the multiplicative smoothing (MS) regularized Gauss-Newton algorithm from De Zaeytijd et al. [8].Let us mention that we also have implemented them into the stepwise relaxed value picking regularized Gauss-Newton algorithm for the reconstructions of the 3D Fresnel targets in [24], but their details were not reported there.
Our first contribution focuses on the forward problem solution and addresses the choice of the initial guess of the total field inside the reconstruction domain for the iterative forward solver.Given the fact that, in each iteration of the reconstruction algorithm, for all permittivity profiles in the line search the forward problem repeatedly needs to be solved for all illuminations, a reduction in the number of iterations of the forward solution leads to important overall computational savings.In Tijhuis et al. [26,27] it was shown that a "marching on in anything" procedure, which extrapolates the initial guess from a few field solutions corresponding to previous values of a varying physical parameter, such as frequency and illumination angle, considerably enhances the speed of convergence.In [28] a "marching on in source position" extrapolation was applied in the first two iterations of a 2D distorted-wave iterative Born algorithm; the corresponding field solutions-for all sources and the two permittivity profiles-were stored and used to initialize a "marching on in profile" extrapolation for the remaining iterations of the reconstruction; a similar approach was employed in the line searches of a 2D quasi-Newton algorithm.In the 3D case a "marching on in profile" extrapolation is memory expensive; hence we applied the "marching on in source position" extrapolation for all iterations in [8].In this paper we extrapolate the initial guess from a few solutions corresponding to previous source positions for the same permittivity, as we did in [8], but we also include a Born-type approximation that is computed from a field solution corresponding to one previous permittivity profile for the same source position.We employ this technique with the stabilized biconjugate gradient fast Fourier transform (BiCGSTAB-FFT) volume integral equation (VIE) forward solver reported in [29], but it can be applied with other iterative forward solvers as well, for example, [30,31].
Our second contribution deals with the numerical solution of the update system in the Gauss-Newton optimization algorithm.In the early 2D microwave reconstruction algorithms, for example, the distorted Born iterative method (DBIM) [32] and the Newton-Kantorovich method [1], an update equation JΔ = Δe scat , with J the Jacobian matrix, is obtained from linearizing the scattered field e scat around the current permittivity iterate .The solution of this equation is formulated in the least squares sense, that is, the solution to the normal equations of the linearized (nonregularized) data error cost function-this corresponds to a Gauss-Newton update of the nonlinear data error cost function [33][34][35].The resulting "least squares" update system with a dense (approximate) Hessian matrix J  J is ill-conditioned and is regularized with a Tikhonov method-note that the regularization thus is applied to the permittivity update rather than to the permittivity.It is solved with a direct inversion method for the case of small objects in, for example, [1,34].More recently a Tikhonov regularized update system has been employed for breast imaging in 2D [36] and in 3D [17]; in the latter it is solved with a conjugate gradient technique.An equivalent approach consists in solving the nonregularized update system with the conjugate gradient for least squares (CGLS) algorithm, where regularization is achieved by terminating the CGLS algorithm after an adequate number of iterations.This approach has been tested in 2D with measured biomedical data [37,38] and in 3D with synthetic data for realistic Magnetic Resonance Imaging (MRI-) derived numerical breast phantoms [18] and with clinical breast data [19].
In the present paper we start from a regularized data error cost function, for example, [8, 11, 23-25, 28, 39-42], in which the regularization is defined on the permittivity (see, e.g., [8,43] for interpretations) and easily allows incorporating a priori information, for example, [42].Application of the Gauss-Newton technique then yields an update system where the (approximate) Hessian matrix J  J +  2 Σ is composed of contributions from the data error and from the regularizing function.We observed that this update system still can be ill-conditioned, even in the presence of the positive definite matrix  2 Σ.In [8,23] we solved the system with a BiCGSTAB routine.In [11] a CGLS iterative technique is combined with an implicit Jacobian calculation and a parallelized algorithm is tested with measured biomedical data.In this paper we employ a two-level preconditioned iterative method to obtain a computationally efficient and stable solution of the Gauss-Newton update.We apply the subspace preconditioned LSQR algorithm from Jacobsen et al. [44] with a 3D cosine basis.This algorithm extracts from the original update system a smaller system, which inherits its ill-conditioning but can be rapidly solved in a direct manner.The remainder of the system is better conditioned and can be solved iteratively with much less iterations than the original system.
Our third contribution is concerned with the use of permittivity constraints.Upper and lower bounds for the complex permittivity, which are based on a priori knowledge about the object, aim to improve the convergence of the algorithm.Often they are enforced, for example, by assigning after each update solution the constraint's values to those reconstruction variables that violate the constraints [34].In [45] such enforced update is reiterated a few times with the CGLS algorithm.In [40] constraints are incorporated in the Gauss-Newton framework by means of a nonlinear transformation that maps the constrained permittivity values on new, unconstrained optimization variables.In this paper we propose a technique that is inspired by the use of parameter transformations, but which is only applied in the line searches.We propose a new, constrained search path in the Gauss-Newton optimization, which incorporates in a smooth manner lower and upper bounds on the object permittivity, such that these bounds never can be violated along the search path.The technique can be implemented without modification of the original optimization scheme and line search algorithm.
The outline of the paper is as follows.In Section 2 the 3D regularized Gauss-Newton algorithm of De Zaeytijd et al. [8] is summarized, since the contributions of this paper are illustrated with this algorithm.In Section 3 the combined "marching on in source position-Born-type approximation" extrapolation procedure is detailed.The subspace preconditioned LSQR implementation is discussed in Section 4 and Section 5 treats the new constrained search path.In Section 6 singlefrequency reconstructions from synthetic data are shown for various numerical biological phantoms, including a realistic 3D MRI-derived numerical breast phantom from the University of Wisconsin-Madison (UWCEM) online repository.Section 7 presents the conclusions.

Formulation
The object with unknown relative complex permittivity is immersed in a homogeneous background medium with relative complex permittivity   , where   is the relative permittivity,  the conductivity,  0 the permittivity of vacuum, r = (, , ) the 3D position vector, and  = 2 the angular frequency.In the following, the   dependency of the time-harmonic fields is implicitly assumed and the -dependency of the complex permittivity is omitted, since a single frequency is adopted.In order to reconstruct the (complex) permittivity within a bounded investigation domain D (see Figure 1(a)), scattering data are collected by successively illuminating the object with   different known incident electric fields E inc  = [ inc , ,  inc , ,  inc , ], which are radiated by elementary electrical dipoles positioned in points r  exterior to D and with polarizations along unit vectors û ; for each illumination  the resulting scattered electric field components E scat  (r   ) ⋅ û  are measured in    positions r   along unit vectors û  [8].All these (complex scalar) scattered field components are stored in one   -dimensional measured field vector e meas , with   = ∑   =1    .Note that  (resp.,   ) refers to a unique combination of a position r  (resp., r   ) and a polarization û (resp., an orientation û  ).
For the numerical reconstruction, we approximate the unknown relative complex permittivity (r) with a piecewise constant function that assumes one value  ] in each cell ] of a uniform cuboid grid D  , which comprises   cubic cells with size Δ  (Figure 1(b)).These values are the reconstruction variables and they are collected in the   -dimensional permittivity vector .The solution of the inverse scattering problem is then defined as the minimizer of a regularized cost function F().In this paper we employ a MS regularized cost function where F LS = ‖e scat () − e meas ‖ 2 /N LS is the (nonlinear) least squares data error, with N LS = ‖e meas ‖ 2 a normalization International Journal of Antennas and Propagation constant and with e scat () the   -dimensional simulated scattered field vector, which contains the scattered field component values E scat  (r   )⋅ û  as computed with the forward solver [29] for a given permittivity vector ;  is a regularization parameter and F  is the smoothing regularizing function from [8]; see (A.1) in the Appendix.
For the minimization of F we apply a Gauss-Newton technique with approximate line searches.In [8] the permittivity vector is updated from iterate  to iterate  + 1 as where s  is the (modified) Gauss-Newton search direction computed from update system (4) and   is the descent step length, that is, a positive real number obtained by means of an inexact line search algorithm [46, pp. 34-38], such that F(  ) is approximately minimized along the direction s  .
In Section 5 of the present paper, we propose an alternative, constrained search path for (3).The (modified) Gauss-Newton search direction for cost function ( 2) is the solution of the system where (⋅)  stands for conjugate (⋅) * transpose (⋅)  , where  2  = N LS F LS  /(1+F   ) and where the subscript  indicates quantities evaluated in   .Furthermore, J  is the Jacobian matrix, with  ] =  scat  / ] ; note that we employ the (independent) variables  ] and  * ] (similar to [47]) and that we omit the subscript  wherever this is more convenient; the vector Ω  and the (positive definite) matrix Σ  contain first-order derivatives (Ω ] = F   / ] ) and second-order derivatives (Σ ] =  2 F   / ]  *  ), respectively, of the regularizing function F   ; see the Appendix for the explicit expressions.The expression between brackets in the RHS of (4) is the gradient with respect to the variables  * ] of cost function (2) and the factor J   J  + 2  Σ in the LHS is obtained after neglecting in the approximate Hessian matrix (i.e., the submatrix that remains after neglecting the elements  2 F  / ]   and  2 F  / * ]  *  ) a number of terms that are specific for the MS regularization (hence we call it a modified Gauss-Newton technique [8]).In [8] system (4) is solved iteratively with the BiCGSTAB routine [48].
Note that the search direction from ( 4) is also the least squares solution of the system or the minimizer of the regularized linear least squares problem where Δe meas  = e meas − e scat  and L  is the Cholesky factor of Σ  = L  L   .The computational effort in the reconstruction algorithm is mainly due to the multiple forward problem solutions and to the computation of the Gauss-Newton search directions.Indeed, for each iterate   the forward problem is solved for   illuminations, in order to compute the scattered field vector e scat  , which is needed to calculate the cost function F(  )-this is used in the stopping criterion and at the beginning of the line search-and to calculate the gradient in the RHS of (4), and in order to compute the total field on the grid, which is needed to calculate the Jacobian matrix (7) (this may require additional forward problem solutions); see further; the computation of the search direction s  involves the solution of a large and usually ill-conditioned linear system (4); furthermore an   multi-illumination forward problem F(  +   s  ) has to be solved for each trial value of   in the line search algorithm.Note that the elements of the Jacobian matrix are given by [8] where   is the background wave number, Φ ] is a 3D unity pulse function with support cell ], and E  is the total field on the computational grid resulting from dipole excitation  in the transmitting position r  and oriented along û (this quantity is available from the forward problem solution for   ); E   is the total field on the computational grid resulting from a dipole excitation in the receiver position r   and oriented along û  ; hence additional forward problem solutions are needed for those receiver positions and orientations that do not coincide with transmitting positions and orientations.
In the following sections we propose some ways to improve the computational efficiency.

Initial Guess for the Forward Solver
In this section we further reduce the computational effort for the forward problem solution by extending the "marching on in source position" extrapolation procedure [27,28] with a Born-type approximation.We apply this technique to the BiCGSTAB-FFT iterative forward solver reported in [29,49], but it can be implemented with other iterative forward solvers as well.Let us first remind the reader of some main principles of this solver.The 3D domain integral equation for the electrical field is expressed in terms of the unknown electric flux density D  =  0 E  and of a mixed potential formulation for the contrast current densities.This expression is discretized with a Galerkin Method of Moments, by employing a number of   vectorial rooftop functions to test the equation and to expand the electric flux density.Note that the vector potential is not expanded directly as in [30] and that the 1/-singularity of the Green function is handled by singularity subtraction.
The rooftop functions are defined on a uniform cuboid grid D  , which is an integer subdivision of the inversion grid D  ; hence the forward problem cell size Δ   can be chosen as a fraction 1, 1/2, 1/3, . . . of the reconstruction cell size Δ  ; Δ   is usually around one-tenth of the wavelength in the considered medium.With using the normalized contrast function, defined as (r) = [(r) −   ]/(r), the resulting   ×   linear system is written as where the   -dimensional vectors e inc  , d  , and d ±  contain the tested incident field values, the unknown expansion coefficients  , of the electric flux density, and the products  ±   , , respectively, on the grid D  for illumination  and permittivity vector .The signs ± indicate one of both support cells of a rooftop function.The products W  d  and Z ± d ±  correspond to the tested total and scattered fields, respectively, on the grid.Note that the elements of the dense matrices Z ± do not depend on the permittivity; hence they need to be calculated only once for a series of scattering simulations with varying contrast.Equation ( 8) is solved for the coefficients d  with the BiCGSTAB method [48], starting from an initial guess d 0  .The FFT method [30,50,51] is used to speed up the matrix-vector multiplications Z ± d ±,  in each iteration  of BiCGSTAB while the multiplication W  d   is fast because W  is a sparse matrix.The resulting BiCGSTAB-FFT method has a memory requirement of O(  ) and a computational effort of O(  log   ), with  the number of iterations needed to reach the desired accuracy, expressed by a relative accuracy threshold   .
In our inversion scheme the forward problem has to be solved several times: for a varying permittivity vector  and with each of those for varying illuminations .In [27] it is shown that the number of iterations  to solve (8) can be significantly reduced by means of a "marching on in anything" technique provided that   is not much smaller than the relative error introduced by noise and the discretization.This technique proposes an adequate choice for the initial guess d 0  based on available solutions which correspond to slightly different illumination and/or object configurations.Suppose we have a few vectors x  ,  = 1, . . ., , that can be regarded as approximations for d  .The initial guess d 0  is then calculated as the linear combination which minimizes the error ‖L  d 0  − e inc  ‖ 2 between the LHS and RHS of (8).The coefficients   thus are a solution of the linear system which is a small system, since  = 3 usually is sufficient.
In this paper we propose to use as the approximations x  to d  on the one hand a few solutions d   , which were computed for the same permittivity vector but for nearby transmitter positions-that is, "marching on in source position"and on the other hand a Born-type approximation d   , which is calculated as with where the vectors d ,±  have elements In this expression, the vector e   (,   ) contains the tested values of the total field E   , which is generated by the applied current density and the Born-type contrast current density where (r) is the new permittivity function (gathered in the new permittivity vector ) and E   (r) is the total field corresponding to a previous permittivity function   (r) (gathered in the permittivity vector   ) for the same illumination .The coefficients   , are available from this previous solution.Since the matrix W  is sparse and since the multiplications can be done with FFTs, the calculation of d   is fast.Note further that in (10) only the multiplication of L  with d   needs to be done, since the other products L  d   = e inc   are the incident field vectors, which are available.The inclusion of the Born-type approximation in the "marching on" scheme thus allows for a simple extrapolation over the permittivity without having to store multiple solution vectors d   for a number of different permittivity profiles   as would be the case in a "marching on in permittivity" scheme [28].

Subspace Preconditioned Gauss-Newton Update
The solution of system (4) for the Gauss-Newton search direction s by means of a direct inversion method would require O((  ) 3 ) operations and O(  (  ) 2 ) operations would be needed to compute the product J  J. Since the number of unknowns   in a 3D inverse problem is large, it is more efficient to solve (4) iteratively.This requires per iteration a multiplication of an   -dimensional vector with J, followed by a multiplication of an   -dimensional vector with J  .The computational complexity O(    ) is much less than a direct inversion provided that the number of iterations  can be kept small.However the condition number of J-and therefore of J  J-generally is very large.Even with the wellconditioned regularization term  2 Σ in (4) the conditioning of the system matrix remains problematic and worsens towards the end of the optimization in case of MS regularization, since  is proportional to the least squares data error.We therefore propose to solve (4) with the iterative subspace preconditioned LSQR algorithm (SPLSQR) of Jacobsen et al. [44], which is an implementation of a two-level iterative method [52] for the solution of large-scale regularized linear least squares problems such as (6).Since this algorithm is conceived for real system matrices, minimization problem (6) is reformulated as with x = [ R (s) where R and I stand for the real and imaginary parts, respectively, and where K is an   × 2  matrix with   = 2  + 2  , x is a vector of length 2  , and y is a vector of length   .As can be seen in Figure 2 for a generic (but relatively small) inverse problem, the singular value spectrum of the matrix K with  = 0, that is, without regularization, gradually decreases over a large number of singular values without showing a clear threshold where this spectrum could be truncated.However, when  ̸ = 0, the regularization introduces a platform at the lower end of the spectrum of K. Nonetheless the conditioning of K still is not very good, because of the rapidly decreasing singular values in the first part of the spectrum.
The principle idea of the SPLSQR algorithm [44] is a splitting of the solution space R 2  into two orthogonally complementary subspaces S V and S  of dimensions  V and   , respectively, with  V +   = 2  , spanned by the (orthogonal) columns of the matrices , respectively.This means we look for a solution x = Vk + Ww.It is furthermore desirable that S V has a small dimension  V and that the conditioning of KW is much better than the conditioning of K. Indeed, introducing the QR factorization where , and Because K has full rank, the same holds for R and the first term in the RHS of ( 18) can be made zero for every w.Therefore, ( 15) is equivalent to Since ( 20) is a small  V ×  V upper-triangular system, it is solved conveniently for k by back substitution.The solution of (19) for w is done with the LSQR algorithm [53], which although mathematically equivalent is numerically more stable than the CGLS method [44] and which can be formulated directly in terms of p = Ww, thereby avoiding the construction of and multiplications with the large matrix W. Calculating the QR factorization (17) and performing the multiplications with Z or Z  are relatively cheap if done with Householder transformations, because  V is small.With an appropriate choice of V (and W) the number of LSQR iterations to solve (19) can be kept small.If, for instance, V consists of the  V principal right singular vectors of K, then the spectrum of KV, denoted as spec(KV), coincides with spec(K) up to   V , the  V th singular value of K. Furthermore, since W is orthogonal to V, the subspace S  is spanned by the 2  −  V least significant right singular vectors of K; hence spec(KW) coincides with (  V +1 , . . .,  2  ).It can be seen from Figure 2 that if  V is large enough, spec(KW) practically completely lies within the plateau in spec(K), introduced by the regularization.This means that the conditioning of KW in this situation is very good.It also implies that KV inherits the ill-conditioning of K, but this does not pose a problem, since (20) is solved directly.However, it is computationally too expensive to construct this SVD subspace; hence a subspace that "resembles" it is chosen.It is well known that for most applications the right singular vectors of an ill-posed problem become more oscillatory as the corresponding singular values decrease.Therefore, we propose in this paper a truncated 3D discrete cosine basis (the DCT-2 basis from [54]), defined on the cubic grid D  , which is used for the real and imaginary part of the complex search direction vector s separately; that is, the matrix V now has the form where V  contains per column one of the  V DCT basis vectors with the lowest spatial frequencies in the -, -, and -directions.In this case the multiplication of K with V can be evaluated efficiently by performing 3D discrete cosine transforms [54] on the rows of K (actually on the first and the second half of these rows separately) and retaining only the  V =  V   V   V  components with the lowest spatial frequencies in the -, -, and -directions, respectively.This subspace can be considered as a coarse grid approximation to the actual solution space.Figure 2 shows spec(KV) for this subspace and it can be seen that it indeed coincides with spec(K) for large singular values.

A New Constrained Line Search Path
In order to improve the convergence of the optimization technique, but mainly to prevent the permittivity and conductivity values from becoming nonphysical or too high to handle with the chosen forward discretization cell size Δ   , including a priori knowledge concerning the expected upper and lower bounds on the complex permittivity in the biological object by means of constraints is recommended.Let us denote the upper and lower bounds by   max and   min , respectively, for the real part of the reconstruction variables and by   max and   min , respectively, for the imaginary part; that is, where R and I stand for the real and imaginary parts, respectively.Such constraints can be incorporated in the Gauss-Newton framework using a nonlinear transformation that maps the constrained permittivity values on new, unconstrained optimization variables [40].Here, we follow an approach which is inspired by the use of parameter transformations but which is only applied in the line search.Furthermore our transformation differs from the ones reported in [40].More specifically, we propose to replace the search path in line search (3) with a smooth, constrained path that entirely lies within the constraints, if the starting point   does so, and which starts along a descent direction.With this constrained path, it is sure that the cost function will be reduced if   is not a local minimizer and that the constraints will not be violated.This path is defined as where s  is the descent direction, solution of (4), and where we propose for the vector function f the expression Path (23) only considers the constraints that can be violated along path (3) with   > 0, hence the distinction between cases  ] ≥ 0 and  ] < 0 in (24).For small   (3) and ( 23) coincide, but in the vicinity of the constraints path ( 23) is bent away from (3) and it has a limit point on the constraints, as is illustrated in Figure 3(a) for a problem with only one complex optimization variable.Although theoretically the optimization variables can never reach their bounds with this procedure and the path always starts along a descent direction, it is possible that very little progress is made if one or more variables are very close to one or more of their bounds.Indeed, in such a situation the line search path deviates from its initial direction s already for very small -values and starts to run along the projection of s on the nearest boundaries of the constrained optimization domain as is illustrated in Figure 3(b).It is possible that although s, calculated with (4), is a descent direction, its projection on these boundaries is not and the cost function starts to increase again for very small values of .In this case the line search is terminated after only a negligible reduction of the cost function.If s is the steepest descent direction, then its projection s  on boundaries like the ones considered here (i.e., upper and lower bounds on the optimization variables), for example, s  = (k  s)k, with k a unit vector, remains a descent direction (or is zero), since the projection of s  on the gradient direction satisfies s   (−s) = −‖k  s‖ 2 ≤ 0 [46,49].Therefore, it is possible in our implementation to switch temporarily to the steepest descent direction when the abovementioned problems with the constraints are encountered.Note that the use of the new search path implements the constrained optimization in a smooth way, without requiring any adjustments to the update system and the line search algorithm.Indeed, the same inexact line search algorithm as in [8] is used, comprising bracketing and sectioning phases and polynomial interpolation; see [46, pp. 34-38].It also avoids sequences of time consuming minimizations of reduced problems, as is needed in active set methods.However, mapping (24) is highly nonlinear and a large step in  may correspond to a negligible step in the actual optimization domain when many optimization variables are close to their bounds.The line search algorithm then may require many steps, hence many forward problem simulations, before a local minimum along the search path is reached.Although the Born-type extrapolation from Section 3 results in a rapid solution of additional forward problems when the permittivity profile changes little, a large number of forward problem solutions may increase the total computation time significantly.Therefore we recommend avoiding too stringent constraints.For example, instead of letting the maximum bound on the imaginary part be at most   max = 0 as required by the passivity condition, we allow for a (sufficiently small) International Journal of Antennas and Propagation positive value.Note finally that although our implementation uses constant values throughout the grid for the constraints   min ,   max ,   min , and   max , these constraints could vary locally in the grid.

Numerical Examples
In this section we discuss reconstructions from synthetic data of a number of simplified 3D biological models as well as of a realistic 3D MRI-derived breast phantom from the UWCEM online repository.All computations are performed on a 64bit computer with 2 GHz Dual Core AMD Opteron processor and 8 GB RAM.FFTs are calculated using the FFTW library [55].The LSQR solver is a self-written nonoptimized code.
In some of the examples we added white Gaussian noise to the simulated scattered fields.We define the signal-to-noise ratio (SNR) as SNR = 10 log 10     e meas where  2 is the variance of the noise, with  0 the exact permittivity profile.We also define the noise level   as the least squares data error obtained for  0 :   = F LS ( 0 ).
In this paper we did not focus on optimizing the transmitter and receiver configurations for each example, by trying to use the smallest possible number of transmitters and receivers in the most appropriate locations for an adequately chosen number of permittivity unknowns.In [56] estimates  are given for the number of degrees of freedom in each component of a single-view 3D scattered field; the corresponding information is accessible by uniformly positioning  receivers over a measurement sphere around the scatterer.However it is stressed in [56] that their results do not apply to the case where transmitters and/or receivers are located in the scatterer's reactive zone.In all our following examples transmitters and receivers are in the scatterer's near field; they are positioned on aspect-limited surfaces and the numbers of receiver positions surpass estimates  provided in [56]; two polarizations are employed for each transmitter and receiver.If the number of permittivity unknowns is larger than the information content in the data-which also can occur with large data sets-then it is the aim of the regularization to reduce the dimension of the solution space.
Let us mention that our fast forward solver, which employs a cuboid grid located completely inside the measurement surface, has a drawback with respect to the spherical and circular measurement surfaces and targets employed in the examples below, in the fact that the distance between these measurement surfaces and targets cannot be reduced beyond a limit determined by the "radius" of the cuboid; in our examples we even added some extra spacing; the circularly and elliptically cylindrical reconstruction domains employed with FEM and FDTD solvers in [17,33] are more convenient in this respect.

Numerical Arm Phantoms.
We consider two simplified arm phantoms: an adult's arm and a child's arm.The adult's arm is from [3] and consists of muscle with permittivity  muscle = 49.6 − 40.4 and bone with permittivity  bone = 8.0 − 3.2 and is immersed in water with permittivity   = 77.3− 21.2.For the child's arm we employ the same permittivities.The working frequency is 1 GHz, which yields a background wavelength of   = 3.38 cm.The investigation domain D is a cube with side 10 cm (≈3  ).The initial guess for the reconstructions is the background medium.

Adult's Arm
Phantom.This phantom is an 8.7 cm wide finite cylinder of muscle containing a 3.3 cm wide finite cylinder of bone.As with the example in [3], which was inspired by the simulated antenna configurations in [57], we employ three parallel rings containing transmitter and receiver antennas and surrounding domain D; see Figure 4. We also use 18 transmitter and 90 receiver positions evenly distributed over the three rings, but in each transmitter position we apply two polarizations, along û and the azimuthal direction û , to an electrical elementary dipole and for each illumination, the scattered field is computed in all 90 receiver positions along these two directions.The data vector thus contains   = 6480 field components or complex numbers; note that due to reciprocity 153 × 4 = 612 of these data values are redundant, which is the number of transmitter position pairs,  2  18 , times the number of data values per transmitter position.The synthetic data are generated with a discretization D  of 30 × 30 × 30 cells with size 3.33 mm (≈  /10) and Gaussian noise yielding a SNR of 30 dB is added (the noise level is   = 10 −3 ).The exact permittivity profile is depicted in Figures 5 and 6, where the permittivity cell is twice as large as that employed for the discretization D  of the fields.
For the inversion we use a grid with 15 × 15 × 15 cells with size 6.67 mm (≈  /5); hence the number of unknowns is   = 3375.Constraints on the complex permittivity are necessary for the convergence of the reconstruction and are implemented in the line search as described in Section 5.The stopping criterion is set to F LS < 4 × 10 −3 , which accounts for the noise and the discretization error.From the reconstructions in Figures 5 and 6 we observe that the shape, dimensions, and position of the muscle and bone are correct.The permittivity values are rather close to the exact values, in particular for the real part in Figure 5.Some smoothing is visible due to the smoothing regularization.These results were obtained after 8 iterations (i.e., 8 solutions of update system (4)), including 21 multiview forward problem solutions, where the "marching on in source position" technique combined with a Born-type approximation from Section 3 was used for the initial guess of the BiCGSTAB-FFT iterative forward solver.The total execution time was 2 hours and 30 minutes.

Child's Arm
Phantom.This phantom is narrower than the previous phantom, it is positioned obliquely (Figure 7), and we employ larger sets of data and unknowns.It consists of a 6 cm diameter finite cylinder of muscle and a bone structure, which is modeled as a 2 cm diameter finite cylinder and a 4 cm diameter sphere.We only consider the portion of the arm that falls within the investigation domain D. The successive illuminations come from 120 electrical dipoles along two orientations, û and the azimuthal direction û , in 60 positions that are evenly distributed over 5 horizontal circles with radius 10.14 cm (=3  ) and with vertical spacing 2.36 cm (=0.7  ).For each illumination, the scattered field is calculated in all 60 positions along the same directions, which leads to   = 14400 complex numbers, nearly half of which are redundant due to reciprocity.The synthetic data are generated with a discretization D  of 30 × 30 × 30 cells with size 3.33 mm (≈  /10); no noise is added this time.
For the inversion we use a slightly coarser grid-in order to avoid committing an inverse crime-with 25 × 25 × 25 cells with size 4 mm (≈0.12  ), which yields   = 15625 complex permittivity unknowns.From   <   -this inequality is even stronger if the redundant data due to reciprocity are removed-it follows that the rank of the Jacobian matrix J is smaller than   ; hence update system (4) would be singular if no regularization were applied.When the iterations proceed, parameter   in (4) becomes smaller and the system evolves toward a singular system, which motivates the use of the subspace preconditioning of Section 4. Note that in practice a perfect data fit is not possible due to measurement and discretization errors or-in our simulation study-due to the misfit of the simulation grids for the data generation and the reconstruction; hence system (4) does not become singular.The SPLSQR algorithm is used with  V = 512; that is, only the discrete cosine basis vectors with the 8 lowest spatial frequencies in the -, -, and -directions are retained.
The result of Figures 8 and 9 is obtained in 7 iterations (i.e., 7 solutions of update system (4)) after a total execution time of about 5 hours and 45 minutes if the relative accuracy of the BiCGSTAB-FFT iterative forward solver is set to   = 10 −3 .The data error at this point is F LS = 2.1 × 10 −4 and is not significantly reduced by further proceeding with the iterations.The constraints that were imposed on the permittivity unknowns are 1.0 < R( ] ) < 85.0 and −50.0 < I( ] ) < 1.0.The fact that we allow for slightly positive imaginary parts is motivated by the observation that too severe constraints can stall or even stop the convergence into a local minimum.The regularization parameter in (2) was set to  = 10 −6 -note that with MS regularization the choice of this parameter is not very critical [8].From Figures 8 and 9 it can be concluded that a nice reconstruction is obtained, where the structures in the arm are clearly visible in the correct permittivity ranges.It deviates from the actual profile by its smoother appearance, which is due to the type of regularization we employed.
The "marching on" technique of Section 3 has been used with  = 4; that is, the initial guess for a forward problem solution is obtained as a linear combination of Born-type approximation (11) and the solutions for three previous transmitter positions on the same horizontal circle in Figure 7.There is no extrapolation over different dipole circles.Without the "marching on" scheme the total computation time for the present example increases to about 8 hours, an increase of about 40%.

Simple Breast Phantoms.
We consider two numerical malignant breast phantoms: one without and one with a skin layer.The breast is a hemisphere with a radius of 5 cm that is filled with a homogeneous, averaged breast tissue with permittivity  breast = 9.99 − 2.82.This value is obtained with a single-pole Debye dispersion model by adopting the Debye parameters from a numerical breast phantom in [58].For the phantom with skin, the breast is covered with a 2.5 mm skin layer with permittivity  skin = 40.94− 16.17(dry skin).We chose a lossless background medium with permittivity   = 10, which matches the real part of the fatty breast tissue  breast .Some guidelines for selecting a matching medium are given in [22].The operating frequency is 1 GHz, yielding a background wavelength   = 9.48 cm.The tumor is a 2 cm diameter sphere with permittivity  tumor = 49.93−14.43[58] at the position (−2 cm, 2 cm, 0 cm).The antenna configuration is depicted in Figure 10note that an experimental hemispherical radar antenna array is reported in [59].There are 72 antenna positions on a hemisphere with radius 9 cm surrounding the front side of the breast, with 12 positions evenly spaced on each of the 6 meridians.Only the 48 central positions indicated in bold are used for the transmitting dipole, which is subsequently oriented along a meridional û and an azimuthal û polarization.The investigation domain D, with dimensions 5 cm × 10 cm × 10 cm, is also indicated.The scattered field is computed in all 72 antenna positions for both orientations and Gaussian noise yielding a SNR of 30 dB is added (the noise level is   = 10 −3 ).The resulting dimension of the data vector is   = 13824; note that due to reciprocity 1128 × 4 = 4512 of these data values are redundant, which is the number of transmitter position pairs,  2  48 , times the number of data values per transmitter position.The cell size in D  used for the data generation is 5 mm (=  /19, ≈ tumor /8) for the phantom without skin (see Figure 11 (top)) and 2.5 mm (=  /38, ≈ tumor /17), that is, the thickness of the skin, for the phantom with skin, in order to facilitate modeling of the skin layer (see Figure 12 (top)).We also employ these respective cell sizes for the forward problem solutions during the reconstructions.

Simple Breast Phantom without Skin.
For the reconstruction of this phantom the cell size of the grid D  is 5 mm; hence the number of complex permittivity unknowns is   = 4000.The relative accuracy threshold for the BiCGSTAB-FFT iterative forward solver is set to   = 10 −3 , since there is no point in much lowering it given the noisy data.The initial guess is the background medium everywhere in D. The regularization parameter is  = 10 −4 and the permittivity constraints are 1.0 < R( ] ) < 80.0 and −50.0 < I( ] ) < 1.0.The preconditioned LSQR algorithm with  V = 4×6×6 = 144 is used to solve for the Gauss-Newton update to an accuracy of 10 −5 .
After 6 iterations the least squares data error reaches the noise level, that is, F LS ≈ 10 −3 , and the optimization is stopped.With further iterations, the reconstruction errorthat is, the error between the reconstructed and exact permittivity profiles-increases again, while the cost function is only marginally reduced.The resulting images are shown in Figure 11 (bottom), where the tumor is clearly visible at the right location and with a higher relative permittivity and conductivity than the surrounding breast tissue.Considering the limited aperture data and the small dimensions of the tumor with respect to the background medium wavelength (<   /4), this result is quite satisfactory.The total execution time for this reconstruction was 31 minutes.When the acceleration techniques proposed in this paper were not applied, that is, without the "marching on" technique and by using a nonpreconditioned BiCGSTAB solver for the Gauss-Newton update, the optimization still was not finished after 160 minutes and in that time only 3 iterations were completed.The main cause of this is the large number of iterations needed by the nonpreconditioned BiCGSTAB solver to converge to the desired accuracy, especially after the third iteration when the data error is already close to the noise level.

Simple
Breast Phantom with Skin.The breast phantom with skin is more realistic but it renders the inverse problem more challenging: there is a nonnegligible scattering from the skin interface, which can obscure the signal from the tumor, and the small thickness of the skin layer results in a finer discretization Δ   = 2.5 mm.For this phantom we demonstrate possibilities and limitations of the reconstruction algorithm with regard to a partial reconstruction based on a priori information.We assume the skin layer to be known and we perform the optimization only for the permittivity cells inside the inner skin contour.Furthermore, we reduce the number of permittivity unknowns by optimizing for aggregates of permittivity cells that are cubes with side 5 mm, that is, 8 cells of the forward grid D  , or portions of cubes next to the skin.This yields   = 1965 complex permittivity unknowns.Note that the preconditioned LSQR algorithm with the 3D discrete cosine basis defined on the cuboid grid cannot be used for this partial optimization example.It may be possible to develop another subspace of breast specific basis functions for this problem.The nonpreconditioned BiCGSTAB solver thus is used instead, but with the relatively small number of reconstruction variables the iteration count remains acceptable.The accuracy for this solver is set to 10 −5 and the regularization parameter is  = 10 −6 .The resulting images after 5 iterations and about 2.5 hours are shown in Figure 12 (bottom).The tumor is again clearly visible at the right location and with a higher relative permittivity and conductivity than the surrounding breast tissue, but these values again are lower than the exact values.Some higher values also appear in a number of spots near the (high permittivity) skin interface, which might be a compensating behavior for reconstruction errors elsewhere in the image.The longer computation time results from the finer discretization in the forward problem.

Reconstruction of Effective Permittivities.
Inspired by a homogenization strategy reported in [58] for 2D time domain ultrawideband detection of breast tumors, we use our inversion algorithm to estimate a homogenized effective breast permittivity from the scattering data of the malignant breast phantom without skin.This value could be used, for example, as an approximation for the healthy tissue permittivity in other-qualitative-inversion methods, such as the linear sampling algorithm [22,60,61], where anomalies are detected against a known background by solving only a linear problem.It also could be used as a guideline for adapting the matching medium   or as an initial estimate for a more detailed quantitative reconstruction.We use only one (aggregated) reconstruction variable for this reconstruction-we keep the same discretization D  as before for the forward problem solutions-and its support coincides with the breast For the malignant breast phantom with skin we estimate two effective permittivities, one for the skin layer and one for the breast tissues, in a similar way as presented for the first phantom.Starting from the background permittivity   , the result after 3 iterations and 52 minutes is  skin,eff = 39.03 − 14.65 and  breast,eff = 10.6 − 3.12, values that are within 6% of the actual complex permittivities  skin = 40.94− 16.17 and  breast = 9.99 − 2.82 in the exact phantom.

Realistic 3D MRI-Derived Numerical
Breast Phantom.The UWCEM Numerical Breast Phantom Repository contains a number of anatomically realistic MRI-derived numerical breast phantoms for breast cancer detection and treatment applications [62].These phantoms capture the structural heterogeneity of normal breast tissue and incorporate the realistic dispersive dielectric properties of normal breast tissue from 0.5 to 20 GHz reported by [14,15].We selected for this paper Phantom 1 from ACR class 1 (ID 071904), which is a mostly fatty breast phantom with some glandular and fibroconnective inhomogeneities.The complex permittivity in a sagittal slice through this phantom at a frequency of 2 GHz is depicted in Figure 13 [63].For the background medium we chose a material with permittivity   = 10.0 − 2.0 (the corresponding conductivity is   = 0.223 S/m), which differs from the background medium employed for the simple breast phantoms by the addition of some losses.The background wavelength is   = 4.72 cm.
Since the cell size of this MRI-derived phantom is very small, only 0.5 mm, and since such a high resolution is not needed at the considered frequency nor desirable due to the high memory and computational costs, we derived the -, -, and -directions is employed, that is, a coarse grid approximation with roughly one-third of the resolution of the full permittivity grid D  .The regularization parameter is  = 10 −5 .To test the abilities of the reconstruction method we perform a blind reconstruction; that is, we do not perform reconstruction solely within a known breast contour as in [17,18], and the initial estimate is the background medium everywhere in D. The constraints on the permittivity are 1.0 < R( ] ) < 55.0 and −50.0 < I( ] ) < 1.0.The "marching on" scheme is used with  = 4; that is, the initial guess for the BiCGSTAB-FFT iterative forward solver is obtained as a linear combination of Born-type approximation (11) and the solutions for three previous transmitter positions.Figure 16 shows the result after 13 iterations (or 18 hours and 25 minutes); the changes in the permittivity profile then are less than 1 percent.The shape of the breast and the overall structure of the inhomogeneities are clearly visible in the reconstruction.The tumor is located correctly as well.Even the small lump of fibroglandular tissue near the nipple can be resolved.We observe some smoothing that might be due to the type of regularization.However, the reconstructed relative permittivity and conductivity values of the inhomogeneities are lower than the exact values: in the center of the tumor the reconstructed permittivity is 21 − 6.5 instead of  tumor = 50.0− 10.0.In [17,18] quantitative reconstructions of numerical breast phantoms from the UWCEM repository are presented for a smaller 1 cm diameter inclusion.The reconstructions are from multiple frequency data values and employ a finer discretization grid than we did.From a visual inspection of Figure 8 in [17] at 1.5 GHz, the inclusion hardly can be distinguished in the permittivity curves, but it clearly appears in the difference image of reconstructions from data with and without inclusion (Figure 9 in [17]), the quantitative values still being too low though.Similar observations hold for [18] at 2.5 GHz, where a contrast agent was applied to the tumor.In [19] a 1.2 cm size tumor is reconstructed in a mostly fatty breast from single-frequency clinical data at 1.3 GHz and the identification of this tumor is facilitated by comparison with a reconstruction of the other-healthy-breast where no tumor is present.Finally we employ this example to illustrate the effect of the subspace preconditioning.We compare the efficiency of the SPLSQR algorithm to the efficiency of a conventional iterative solver without preconditioning in solving (4).The conventional solver is the BiCGSTAB routine [48].For the first ( = 0) and the final ( = 12) Gauss-Newton iterations we let the SPLSQR algorithm with the same parameters as mentioned earlier solve (19) to a relative accuracy of 10 −4 and we calculate the resulting accuracy on original system (4).We then let BiCGSTAB solve (4) to that accuracy and compare the number of iterations and the time needed by both methods.The comparison is summarized in Table 1.In the beginning of the optimization the SPLSQR algorithm is more than 4 times faster than BiCGSTAB and towards the end of the optimization, when the preconditioning becomes more crucial due to the decrease of parameter  in (4), the speedup factor approaches 20.The reduction in solution time is less than the reduction in the number of iterations in the case  = 0, since the computation of QR factorization ( 17) is also included in the former.Note that the time needed to solve (4) a single time without preconditioning at  = 12 is 90% of the total reconstruction time with preconditioning.Note furthermore that we use the BiCGSTAB routine from the PIM library [48], which is an optimized Fortran library, International Journal of Antennas and Propagation  while we implemented the SPLSQR ourselves in C; thus the speedup might even be more significant with an optimized implementation.

Conclusions
We have employed a three-dimensional quantitative reconstruction algorithm in the context of biomedical microwave imaging.The algorithm is based on a Gauss-Newton optimization scheme with line searches that is applied to a regularized nonlinear cost function.The forward solver has been improved by the combination of a "marching on in source position" and a "Born-type approximation" extrapolation procedure for choosing the initial guess of the total field on the computational grid.A two-level subspace preconditioned iterative method has been employed for the solution of the complex permittivity update system.In particular, we have applied the subspace preconditioned LSQR algorithm from Jacobsen et al. ( 2003) with a 3D cosine basis.A significant improvement in the computational efficiency with respect   to using the BiCGSTAB solver is observed here, but it would be worthwhile to compare its performance also with a CGLS iterative scheme.A modification of the line search path based on a parameter transformation has been proposed to allow for a smooth incorporation of constraints in the optimization method.Numerical experiments have been conducted at one single frequency, thereby avoiding difficulties with the dispersive nature of body tissues, and with bipolarized antenna configurations.They have shown that the method yields promising results for biomedical imaging, in particular for the challenging problem of breast tumor detection.Multifrequency schemes, dedicated antenna positioning configurations, and the incorporation of more a priori information in the reconstruction strategy may further improve the resolution of the images and the estimation of the permittivity.Possible benefits of using bipolarized versus monopolarized transmitters and/or receivers could be explored further, although a bipolarized configuration may substantially affect the complexity of a measurement apparatus.Future developments also should include the use of a problem-specific regularization strategy to cope with the smearing effect shown by the smoothing regularization.with Ω ] = F  / ] and where Σ is a real and constant matrix with Σ ] =    ,,ℎ in these expressions represents the set of neighboring cells of cell (, , ℎ), that is, the six cells that share a face with cell (, , ℎ).These also include the virtual neighboring cells just outside D  when (, , ℎ) is on the border of D  .

Figure 1 :
Figure 1: The 3D scattering configuration (a) and the permittivity grid D  with the piecewise constant approximation for the permittivity function (b).

Figure 2 :
Figure 2: The singular value spectra of the matrix K without regularization ( = 0) and with regularization ( ̸ = 0) for a generic inverse scattering problem with 2000 permittivity unknowns.The spectrum of the matrix KV where V contains the truncated 3D discrete cosine basis corresponding to  V = 432.

Figure 3 :
Figure 3: Illustration of the constrained line search path: (a) the starting point   lies far from the boundaries of the optimization domain, such that the new search path from (23) coincides with the original path from (3) up to large   ; (b) the starting point   lies close to a boundary of the optimization domain, such that the new search path from (23) starts running along the projection of (3) on this boundary already for small   .

Figure 4 :
Figure 4: A view on the investigation domain D and the dipole configuration for the numerical adult's arm phantom.Dipole positions are indicated with dots and orientations with arrows.All dipoles act as receivers; a subset acts as transmitters.

Figure 5 :
Figure 5: Real part of the exact permittivity used to simulate the data at 1 GHz with a SNR of 30 dB (left) and of the reconstructed permittivity (right) for the adult's arm phantom, in the -plane (top) and in three perpendicular planes (bottom).Distances are in m.

Figure 6 :
Figure 6: −Imaginary part of the exact permittivity used to simulate the data at 1 GHz with a SNR of 30 dB (left) and of the reconstructed permittivity (right) for the adult's arm phantom, in the -plane (top) and in three perpendicular planes (bottom).Distances are in m.

Figure 7 :
Figure 7: A view on the numerical child's arm phantom, the investigation domain D, and the dipole configuration.Dipole positions are indicated with dots and orientations with arrows.All dipoles act as transmitters and receivers.Distances are expressed in background wavelengths.

Figure 8 :
Figure 8: Real part of the exact permittivity used to simulate the data at 1 GHz without added noise (left) and of the reconstructed permittivity on the coarser grid (right) for the child's arm phantom, in the -plane (top) and -plane (bottom).Distances are in m.

Figure 10 :
Figure 10: Views on the investigation domain D and the dipole configuration for the simplified breast phantoms: front view (a), side view (b).The arrows visualize the dipole orientations and the large black dots indicate positions with both transmitters and receivers.

Figure 11 :
Figure 11: Exact permittivity (top) and reconstructed permittivity (bottom) for the numerical simple breast phantom with a 2 cm diameter tumor and without skin.Real part (left) and −imaginary part (right).The SNR of the data is 30 dB.Distances are in m.

Figure 12 :
Figure 12: Exact permittivity (top) and partially reconstructed permittivity inside the skin layer (bottom) for the numerical simple breast phantom with a 2 cm diameter tumor and with skin.Real part (left) and −imaginary part (right).The SNR of the data is 30 dB.Distances are in m.

Figure 13 :Figure 14 :
Figure 13: The permittivity in a sagittal slice through full-resolution MRI-derived numerical breast phantom 1 from ACR class 1 at 2 GHz.The cell size is 0.5 mm.Real part (a) and imaginary part (b).

Figure 15 :
Figure 15: A view of the dipole configuration for the realistic breast phantom.Dipole positions are indicated with dots and orientations with arrows.All dipoles act as receiver and transmitting dipoles are indicated with larger dots.The gray shape represents the outline of the breast phantom.Distances are in m.

Figure 16 :
Figure 16: The reconstructed permittivity in a sagittal slice at  = 0. Real part (a) and imaginary part (b).Distances are in m.

Table 1 :
A comparison between the SPLSQR algorithm and the BiCGSTAB iterative method for the solution of update system (4) in case of the UWCEM breast phantom.The number of iterations, the solution time, and the resulting accuracy are given for the first ( = 0) and the last ( = 12) iterations in the reconstruction algorithm.