Reconstruction of 3-D Microwave Images based on a Block-BiCGStab Algorithm

In 3-D microwave imaging, gradient-based optimization algorithms usually make use of the so-called stabilized version of the biconjugate gradient iterative method (BiCGStab) in order to solve multiple linear systems. We propose to use a block version of BiCGStab to jointly solve the mutiple right-hand side linear systems. Illuminations are partitioned in subgroups, which makes the method more efficient. The reconstruction process is studied on realistic simulated data and illustrates the efficiency of the method compared to BiCGStab.


Introduction
The development of microwave imaging algorithms has received a lot of interest in the last few years with applications such as biomedical imaging and geoscience [1,2]. Solving the inverse scattering problem under realistic condtions presents several difficulties, a critical one being the reduction of the computational cost associated with three-dimensional problems. Nonlinear inversion methods usually rely on iterative algorithms to reconstruct the dielectric properties of the unknown object. A regularized data misfit cost function is minimized using a gradient-based optimization procedure, such as the Distorted Born iterative method (DBIM) [3] or the Gauss-Newton method [2,4]. At each iteration, the computation of a descent direction and of a step size requires the evaluation of the objective function and of its gradient. Such computations rely on the resolution of linear systems equal to the number of illuminations. In practice, exact inversion of such systems becomes computationally prohibitive for large 3-D problems and approximate solutions are obtained with iterative methods, such as the conjugate gradient method [5], the quasi minimal residual (QMR) method [6] or the stabilized version of the biconjugate gradient stabilized (BiCGStab) method [1,7].
Different methods have been proposed to solve efficiently such multiple systems, such as the marching-on-anything technique [4] or a massive parallelization scheme [2]. Since the linear systems share the same system matrix, block versions of iterative algorithms are appropriate, which simultaneously solve linear systems with multiple right-hand sides. In [8], a block-QMR algorithm was proposed for microwave imaging. In [9], a block version of BiCGStab was developed for random system inversion, which was shown to outperform the Block-QMR algorithm. In [10], we proposed a Partial-Block BiCGStab procedure to solve multiple forward problems in 3-D microwave imaging, which partitions the sources in several subgroups and uses Block-BiCGStab on each group. In this paper, the Partial-Block BiCGStab algorithm is applied for the computation of the objective function and its gradient at each iteration of the inversion process.
In Section 2, the 3-D formulation is presented. The inverse problem is introduced in Section 3 where the Block-BiCGStab algorithm and the Partial-Block version are detailed. Finally, the computational costs of the different methods are compared on simulated data in Section 4.

Formulation of the 3-D Problem
An unknown inhomogeneous object of complex permittivity (r) is included in a volume V , where r denotes the spatial position in V . We assume a homogeneous background medium with complex permittivity b and magnetic permeability µ 0 , that is the vacuum permeability. We suppose that all objects and media are nonmagnetic. Successive known harmonic incident electric fields E inc i (r), i = 1, . . . , N S illuminate the object with angular frequency ω. The scattered fields E scat i are measured at locations r R j , j = 1, . . . , N R . In the following, the time dependence term e −jωt is omitted. The observation equation links the data with the total electric field E i in the domain V . It reads [11]: where k b and χ respectively denote the wavenumber k b = ω 2 µ 0 b and the contrast function The ∇∇· term represents the gradient of the divergence. The homogeneous Green's function g(r, r ) is defined as g(r, r ) = exp jk b r − r /4π r − r . Note that the integral in (1) is a convolution with the kernel g(r, r ). The total field E i (r) in the domain V is ruled by the electric field integral equation (EFIE) or domain equation: In order to obtain a discrete model, we assume that the volume V is a cuboid divided into N cubic voxels of dimensions N x , N y and N z . Let e inc i denote the vector of length 3N containing the components of the incident field discretized on each voxel n with center r n , n = 1, . . . , N . Similarly, let e i denote the discretized total electric field. The contrast function χ(r) is assumed to be constant on each voxel. The vector x of length N contains the unknown contrast on all voxels. The domain equation (2) is discretized using the method of moments [12], which yields: where X is a 3N × 3N diagonal matrix whose diagonal contains three replicas of the vector x and G is a 3N ×3N convolution matrix corresponding to the discretization of (k 2 b +∇∇·)g(r, r ). The discretization procedure is detailed in [13]. Eq. (3) can be written as a linear system: where I 3N is the identity matrix of size 3N . Determining the total field e i in (4) corresponds to solving the forward problem. The observation equation (1) is discretized similarly. The three components of the scattered fields E scat i measured at the N R receivers are stored in the vector y i of size 3N R . The discretized observation equation reads The N R × N matrix G o is constructed similarly to G. The vector b is a perturbation term accounting for noise and model errors. Combining the domain equation (4) and the observation equation (5), the measurements y i can be expressed as functions of the contrast x according to: 3. Three-dimensional reconstruction 3.1. Optimization of a regularized objective function Solving the inverse problem consists in determining the contrast function of an unknown object embedded in the volume V from measurements of the scattered field for N S incident fields. Using the discretized model (6), the problem amounts to estimating the discretized contrast x from data y 1 , . . . , y N S . This problem is ill-posed, therefore we define the solution as the minimizer of a regularized cost function: where is the quadratic data misfit and R(x) is a regularization term operating on the differences between adjacent voxels: where {X k, ,m } k, ,m represents the 3-D arrangement of x. We use the " 2 1 " function ϕ(u) = δ 2 + |u| 2 , which tends to preserve edges in the reconstructed 3D-image. Moreover, for any δ = 0, the cost function is continuously differentiable on R N and gradient-based optimization methods can be used to solve (7). The parameter λ controls the trade-off between the data misfit and the regularization function. In this paper, parameters λ and δ are set empirically. We consider the local minimization of criterion (7) with the L-BFGS algorithm [14]. Each iteration of such algorithm involves one computation of the gradient and possibly several computations of the objective function for the line search step. The gradient reads: where the symbol † denotes the Hermitian adjoint, the overline denotes the complex conjugation and diag{u} stands for the diagonal matrix with diagonal u.

Block algorithms for linear system resolutions
The evaluation of the cost function F requires the computation of N S forward problems L −1 x e inc i , i = 1, . . . , N S . The computation of the gradient also requires the resolution of N S systems with matrix L x . For realistic 3-D problems, exact inversion of such linear systems is computationally prohibitive: L x is a full matrix with 3N × 3N elements. Consequently, the systems are generally solved by iterative methods, such as the BiCGStab algorithm [1,2], which must be called N S times for the computation of F and N S times for the resolution of systems with normal matrix L x in (9). The most computationally demanding steps in BiCGStab are the matrix-vector multiplications with the matrix L x , which essentially correspond to convolution operations, that can be efficiently implemented by Fast Fourier Transform algorithms.
In this paper, we propose a method to solve all the linear systems jointly. Indeed, the operator L x does not depend on the illuminations and the linear systems can be aggregated. The forward problems involved in the computation of both the cost function and its gradient read: The computation of the gradient also requires a second set of system resolutions with matrix L x , that can also be aggregated as: with Y = [y 1 , . . . , y N R ]. We propose to use the Block-BiCGStab algorithm, which is a generalization of the BiCGstab algorithm to multiple right-hand side systems [9]. It is described in Algorithm 1 for the joint resolution of the forward problems in (10). The notation T, S F denotes the Frobenius product between matrices T and S, that is, the trace of matrix T † S. If N S = 1, Block-BiCGStab reduces to BiCGStab. In the following, the implementation of the computation of the total fields in (10) is detailed. The resolution of (11) is performed similarly.

Algorithm 1 Block-BiCGStab algorithm 1: For an initial guess
k ← k + 1 15: until ∀i, r Run the Block-BiCGStab algorithm on E inc (p) and E (p) instead of E inc and E.

4: end for
The Block-BiCGStab method is an iterative algorithm, where each iteration requires the resolution of two linear systems of size N S with matrix R † 0 V at steps 6 and 12 of Algorithm 1. These systems are relatively small, therefore the corresponding matrices can be computed and system resolutions can be performed exactly. However, during the iterative process, the matrix R † 0 V may become ill-conditioned, which leads to instabilities. This may happen in particular when two sources, say with indices i and j, are very close. In that case, the incident fields e inc i and e inc j are highly correlated and this property is propagated into the quantities R (k) and V. To overcome the instabilities, we propose a procedure called Partial-Block BiCGStab. To avoid ill-conditioned matrices, the N S sources are divided into P subgroups such that all sources in one subgroup p are as far as possible from each other. Then, matrix E inc is split into P submatrices E inc (p) in which the columns are weakly correlated. Similarly, matrix E is split into submatrices E (p) . Then, the Block-BiCGStab is applied independently on the P subproblems. This procedure is summarized in Algorithm 2. For P = 1, Partial-Block BiCGStab reduces to Block-BiCGStab. For P = N S , it reduces to the BiCGStab algorithm. The choice of P is important, more information on its influence are given in [10]. The resolution of the multiple forward problems is faster for Partial-Block BiCGStab than for BiCGStab and is even better for high contrasted objects [10]. Here, the Partial-Block BiCGStab procedure is used to compute the cost function and the gradient as well in order to speed up the reconstruction algorithm using L-BFGS.

Simulation Results
In this section, we evaluate the computation time and the convergence of our method compared to BiCGStab and Block-BiCGStab. First, the resolution of the multiple forward problems are studied. Then, results on the reconstruction process are presented. We consider a simulated experiment presented in Fig. 1 (left) that can be found in the literature, with N S = 96 sources at frequency 2.45 GHz embedded in the air of permittivity b = 1. The sources are placed on three coaxial circles at heights z = −λ b , 0 and λ b , and radius 1.5λ b , where λ b is the wavelength of the incident fields. On each circle, N S /3 = 32 sources are equally distributed. A lossy cube of size 0.5λ b with contrast χ 1 = 6 + 1.2j is embedded in a cubic volume V of size λ b . Inside the lossy cube, a cuboid is inserted with square cross-section of 0.25λ b , height 0.5λ b and contrast χ 2 = 8 + 4j. The volume V is discretized into 8 × 8 × 8 voxels. The 96 receivers are placed on three coaxial circles at the same heights as the sources and with radius 1.3λ b . Synthetic data are simulated on a grid of 16 × 16 × 16 voxels and zero-mean white Gaussian noise is added with signal-to-noise ratio equal to 20 dB. The reconstruction procedure is implemented with Matlab on a personal computer with 8 Go RAM and eight cores clocked at 3.40 GHz (four hyperthreaded cores). The reconstruction process is parallelized by simply executing loops with Matlab parfor instruction. For BiCGStab, the 2N S system resolutions are sent in parallel to the different processors. The Block-BiCGStab algorithm is not easily parallelizable and therefore is not parallelized. For the Partial-Block BiCGStab algorithm, the resolutions of the different subproblems are also sent in parallel to the different processors. For Partial-Block BiCGStab, we divide the set of sources into four subgroups (see Fig. 1 left). For Block-BiCGStab and Partial-Block BiCGStab, R 0 is set to L x R (0) . Fig. 2 plots the evolution of the relative residual errors for the resolution of the 96 forward problems using the three algorithms: BiCGStab, Block-BiCGStab and Partial-Block BiCGStab. The tolerance is set to 10 −6 , which is a strict value. Two initializations are tested: E (0) = E inc , and E (0) close to the exact solution of (10) (which can be computed directly for such reasonable size problem). We note that regardless of the initialization, BiCGStab and Partial-Block BiCGStab converge in nearly the same number of iterations, but Partial-Block BiCGStab is the most efficient. On the contrary, the Block-BiCGStab algorithm converges much faster if it is initialized close to the solution. The computation times for the three resolutions with In the reconstruction process, the three algorithms were implemented to solve (10) and (11) involved in the cost function and its gradient. The hyperparameter λ was set to 10 −5 and δ to 10 −2 . At iteration k of the L-BFGS algorithm, the resolutions of (10) and (11) are initialized with the solutions obtained at iteration k − 1. Implementations based on the BiCGStab, Block-BiCGStab and Partial-Block BiCGStab algorithms converge respectively in 855 s, 1004 s and 573 s. The tolerance set to 10 −6 on the residual errors for the three algorithms leads to good accuracy in the resolutions of the linear systems.
A perspective is to solve inverse problems of larger size. We will also study the influence of the tolerance: a weak tolerance generates strong approximations in the system resolutions and could introduce instabilities in the convergence of the reconstruction algorithm.