Three-Dimensional Microwave Imaging: Fast and Accurate Computations with Block Resolution Algorithms

Abstract

This paper considers the microwave imaging reconstruction problem, based on additive penalization and gradient-based optimization. Each evaluation of the cost function and of its gradient requires the resolution of as many high-dimensional linear systems as the number of incident fields, which represents a large amount of computations. Since all such systems involve the same matrix, we propose a block inversion strategy, based on the block-biconjugate gradient stabilized (BiCGStab) algorithm, with efficient implementations specific to the microwave imaging context. Numerical experiments performed on synthetic data and on real measurements show that savings in computing time can reach a factor of two compared to the standard, sequential, BiCGStab implementation. Improvements brought by the block approach are even more important for the most difficult reconstruction problems, that is, with high-frequency illuminations and/or highly contrasted objects. The proposed reconstruction strategy is shown to achieve satisfactory estimates for objects of the Fresnel database, even on the most contrasted ones.

1. Introduction

Microwave imaging aims at estimating the dielectric properties (permittivity and conductivity) of objects illuminated by incident electromagnetic fields [1]. This technique has been used in a wide range of applications such as biomedical imaging [2,3], geophysics [4,5] and nondestructive testing [6,7]. Quantitative inverse scattering methods estimate the complex permittivity map inside the object from a set of measured scattered electromagnetic fields around it. When the object size is comparable to the wavelength, diffraction occurs and the forward scattering model, which defines the scattered field as a function of the dielectric properties at any point of the object, becomes non-linear. The resulting reconstruction problem is ill-posed [8,9], and its resolution becomes particularly difficult with high-frequency illuminations and/or for highly contrasted objects—that is, if the permittivity inside the object strongly differs from that of the background—where the Born approximation is not valid anymore.

Imaging methods usually rely on solving an optimization problem, whose complexity is strongly impacted by the non-linearity of scattering. In order to circumvent this difficulty, several works introduced additional variables, corresponding to the total fields inside the object, which are jointly estimated together with the contrast function (i.e., the relative complex permittivity). Within this formulation, the domain equation, which links the total fields to the contrast function, is not exactly solved, but the resulting bilinear model allows the low-cost computation of iterates in the optimization procedure. The Modified Gradient Method [10] and Contrast Source Inversion methods [11,12] fall within this category, which achieve satisfactory compromise between computation time and reconstruction quality in cases where diffraction remains limited [13]. Although some works addressed 3-D reconstruction based on this formulation [7,14,15,16], most of them were limited to 2-D problems, where the number of unknowns remained relatively small [13,17,18,19,20,21,22,23,24].

A second class of methods, which dominates in the recent 3-D microwave imaging literature [25,26,27,28,29,30,31,32,33,34,35]. directly incorporates the domain equation into the problem formulation. The price to pay for using such an exact form of the domain equation is a much more complex evaluation of the cost function and of its gradient: both require a numerical solver for computing the forward problem, based on the resolution of very large linear systems, especially in the 3-D case. The related computational effort then represents the major cost of the reconstruction procedure. Therefore, developing efficient strategies for such system resolutions remains one of the current challenges for improving the efficiency of microwave imaging techniques [33,36,37]. In order to reduce the size of the problem, some works focused on the reconstruction of strongly constrained, parameterized, objects. Efficient numerical methods were developed with global convergence properties, e.g., based on evolutionary algorithms [38] or Markov Chain Monte-Carlo methods, coupled with surrogate models [39]. In the particular case of point-like scatterers, the two-step procedure in [40] efficiently accounts for the sparsity of the object.

In this paper, we focus on methods that aim to estimate any permittivity distribution (that is, without any strong prior assumption on the inspected object), based on the exact problem formulation. Previous works based on multiplicative [33] or additive [13], edge-preserving or sparsity-enhancing [41], regularization fall within this category. All these methods rely on repeated computations of forward scattering problems. One generally resorts to iterative system solving algorithms, which are stopped as soon as an acceptable residual error is obtained. Several such algorithms have been used in microwave imaging, e.g., the conjugate gradient (CG) algorithm [42,43], the biconjugate gradient (BiCG) algorithm [44], the quasi-minimal residual (QMR) algorithm [45,46] or the biconjugate gradient stabilized (BiCGStab) algorithm [33,36,37,47]. BiCGStab is now one of the most popular algorithms in this field and it was shown to converge faster than CG, BiCG and QMR [47].

In microwave imaging, as many forward scattering models as the number of illuminating sources must be computed. Such computations are usually performed sequentially or via massively parallel procedures [33,48]. Note that all such linear systems involve the same system matrix, which only depends on the object dielectric properties and on the discretization scheme. A classical solution relies on the preliminary LU-decomposition of the system matrix [4]. However, this possibility becomes prohibitive in 3-D imaging [49]. Block iterative algorithms specifically address the resolution of such linear systems with multiple right-hand sides (RHSs), by jointly solving all systems. Block versions of several standard iterative methods have been proposed in the applied mathematics literature, including the block-conjugate (and bi-conjugate) gradient algorithms [50], the block-QMR algorithm [51] and the block-BiCGStab algorithm [52]. Such block algorithms were shown to achieve better performance than their sequential counterparts in generic contexts. Block-QMR has been used in microwave imaging [49]. The block-BiCGStab algorithm proposed in [52] was shown to outperform block-QMR. Its standard implementation, however, is rather unsatisfactory for microwave imaging forward problems, where the system matrices are much larger, and where there are many more RHS terms which are highly correlated [53].

Here, we propose a reconstruction procedure based on additive penalization, where the evaluation of both the cost function and its gradient are performed with a dedicated block-BiCGStab implementation. Specific tunings of the block-BiCGStab algorithm are discussed, and the overall efficiency of the inversion method is evaluated as a function of the problem difficulty (illuminating frequency and contrast). We show in particular that improvements are even greater as the problem complexity increases. We provide reconstruction results obtained on the Fresnel database [54], where the proposed implementation allows the accurate reconstruction of the most difficult objects.

In Section 2, the 3-D forward scattering model is built using an integral formulation. A discretization procedure based on the method of moments is used, and our inversion formulation based on additive regularization is detailed. In Section 3, we propose efficient adaptations of block-BiCGStab for forward scattering computations, and for their implementation within the reconstruction procedure. In Section 4, algorithmic performance is evaluated with simulated data, as a function of the object complexity, with different illuminating frequencies and contrast levels. In Section 5, our method is applied to real measurements extracted from the Fresnel database [54], where reconstruction quality and computation times are discussed. Concluding comments are given in Section 6.

2. Formulation of 3-D Forward and Inverse Scattering Problems

2.1. Forward Scattering Problem

The forward scattering model defines the total electric field in a given volume V containing the scattering object, as a function of the spatial distribution of the object permittivity. We adopt a standard approach based on the integral formulation of the domain equation [1], followed by a discretization procedure based on the method of moments [55]—see for example [14,36,37,43,44] for similar choices. We consider a background medium with constant complex permittivity ${ϵ}_b$, containing the object under study with complex permittivity $ϵ\left(\mathit{r}\right)$ where $\mathit{r}$ denotes the vector of spatial coordinates. The object is illuminated by an incident wave with angular frequency $\omega$. In the sequel, the time dependency $exp(-j\omega t)$ of all quantities is omitted. The electric field integral equation gives an implicit relation between the unknown total electric field at any point inside the domain under study and the object permittivity distribution [1]:

$${\overrightarrow{E}}^{\mathrm{tot}}\left(\mathit{r}\right)={\overrightarrow{E}}^{\mathrm{inc}}\left(\mathit{r}\right)+\left(k_b^2+\nabla \nabla ·\right){\int }_{{\mathit{r}}^{\prime }\in V}g(\mathit{r}-{\mathit{r}}^{\prime })\chi \left({\mathit{r}}^{\prime }\right){\overrightarrow{E}}^{\mathrm{tot}}\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime },\forall \mathit{r}\in V,$$

(1)

where ${\overrightarrow{E}}^{\mathrm{tot}}\left(\mathit{r}\right)$ and ${\overrightarrow{E}}^{\mathrm{inc}}\left(\mathit{r}\right)$ respectively denote the total and the known incident 3-D electric field vectors at point $\mathit{r}\in$V, ${\mu }_0$ is the vacuum permeability, $k_b=\sqrt{{\omega }^2{\mu }_0{ϵ}_b}$ is the wavenumber in the background medium, $g\left(\mathit{r}\right)=exp(j{k}_b\parallel \mathit{r}\parallel )/4\pi \parallel \mathit{r}\parallel$ is the Green function and $\chi \left(\mathit{r}\right)=(ϵ\left(\mathit{r}\right)-{ϵ}_b)/{ϵ}_b$ is the contrast function. The term $\nabla \nabla ·$ represents the gradient of the divergence with respect to $\mathit{r}$.

Using the discretization procedure based on the method of moments (see [14] for the 3-D case), the volume V is divided into N cubic voxels and the integral Equation (1) turns into:

$${\mathit{e}}^{\mathrm{tot}}={\mathit{e}}^{\mathrm{inc}}+{\mathbf{G}}_D\mathbf{X}{\mathit{e}}^{\mathrm{tot}}.$$

(2)

The size-$3N$ vectors ${\mathit{e}}^{\mathrm{inc}}$ and ${\mathit{e}}^{\mathrm{tot}}$ respectively contain the three spatial components of the 3-D incident and total electric field vectors, discretized on the N voxels that are sorted according to the lexicographic order. The size-N vector $\mathit{x}$ corresponds to the contrast function discretized on the N voxels, put in vector form, and $\mathbf{X}$ is the $3N\times 3N$ diagonal matrix whose diagonal contains three replicas of $\mathit{x}$. The $3N\times 3N$ matrix ${\mathbf{G}}_D$ corresponds to the discretization of the $3\times 3$ Green tensor $\left(k_b^2{\mathcal{I}}_3+\nabla \nabla \right)g(\mathit{r}-{\mathit{r}}^{\prime })$ on the N voxels, where ${\mathcal{I}}_3$ denotes the $3\times 3$ identity tensor (see for example [53] for details). From (2), the discretized total field is the solution to the following linear system:

$${\mathbf{L}}_x{\mathit{e}}^{\mathrm{tot}}={\mathit{e}}^{\mathrm{inc}},\mathrm{with}{\mathbf{L}}_x={\mathbf{I}}_{3N}-{\mathbf{G}}_D\mathbf{X},$$

(3)

where ${\mathbf{I}}_{3N}$ denotes the $3N\times 3N$ identity matrix.

In microwave imaging, the object is successively illuminated with a set of $N_S$ sources at different spatial locations, creating incident fields ${\mathit{e}}_i^{\mathrm{inc}},i=1,\dots ,N_S$. Then, computing the corresponding total fields ${\mathit{e}}_i^{\mathrm{tot}}$ amounts to solving $N_S$ linear systems ${\mathbf{L}}_x{\mathit{e}}_i^{\mathrm{tot}}={\mathit{e}}_i^{\mathrm{inc}}$, that is:

$${\mathbf{L}}_x{\mathbf{E}}^{\mathrm{tot}}={\mathbf{E}}^{\mathrm{inc}},\mathrm{with}{\mathbf{E}}^{\mathrm{tot}}=\left[{\mathit{e}}_1^{\mathrm{tot}},\dots ,{\mathit{e}}_{N_S}^{\mathrm{tot}}\right]\mathrm{and}{\mathbf{E}}^{\mathrm{inc}}=[{\mathit{e}}_1^{\mathrm{inc}},\dots ,{\mathit{e}}_{N_S}^{\mathrm{inc}}],$$

(4)

where all linear systems involve the same system matrix ${\mathbf{L}}_x$. Such a property is the key point for our block resolution strategy.

2.2. Observation Model

The scattered field ${\overrightarrow{E}}^{\mathrm{scat}}$ at any point $\mathit{r}$ outside V is defined by the volume integral observation equation:

$${\overrightarrow{E}}^{\mathrm{scat}}\left(\mathit{r}\right)=\left(k_b^2+\nabla \nabla ·\right){\int }_{{\mathit{r}}^{\prime }\in V}g(\mathit{r}-{\mathit{r}}^{\prime })\chi \left({\mathit{r}}^{\prime }\right){\overrightarrow{E}}^{\mathrm{tot}}\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }.$$

(5)

Using the same discretization procedure as in Section 2.1, considering a set of $N_M$ measurement points, it becomes:

$${\mathit{e}}^{\mathrm{scat}}={\mathbf{G}}_O\mathbf{X}{\mathit{e}}^{\mathrm{tot}},$$

(6)

where the size-$N_M$ vector ${\mathit{e}}^{\mathrm{scat}}$ collects the measurements of the scattered fields at the receiver locations, and the $N_M\times 3N$ matrix ${\mathbf{G}}_O$ corresponds to the discretization of the Green tensor $\left(k_b^2{\mathcal{I}}_3+\nabla \nabla \right)g(\mathit{r}-{\mathit{r}}^{\prime })$, for ${\mathit{r}}^{\prime }$ at the N voxels in V and $\mathit{r}$ at the measurement points [56]. Combining both domain (3) and observation (6) equations finally yields the observation model:

$${\mathit{e}}_i^{\mathrm{scat}}={\mathbf{G}}_O\mathbf{X}{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}={\mathbf{G}}_O\mathbf{X}{({\mathbf{I}}_{3N}-{\mathbf{G}}_D\mathbf{X})}^{-1}{\mathit{e}}_i^{\mathrm{inc}},\mathrm{for}\mathrm{any}\mathrm{incident}\mathrm{field}{\mathit{e}}_i^{\mathrm{inc}}.$$

(7)

2.3. Inverse Problem Formulation

The microwave imaging inverse problem consists of estimating the contrast function from the measured scattered fields at the receivers. This problem is ill-posed [8,9], notably because the relation between the contrast function and the data, obtained from the two coupled Equations (1) and (5) in the continuous formulation (Equations (2) and (6) in the discrete case), is highly nonlinear.

Here, we adopt a classical approach to inverse problems [13,57], which consists of minimizing the misfit between the data ${\tilde{\mathit{e}}}_i^{\mathrm{scat}}$ and the simulated scattered fields ${\mathit{e}}_i^{\mathrm{scat}}$, penalized by an additive regularization term. More precisely, we define the reconstructed contrast function $\widehat{\mathit{x}}$ by:

$$\widehat{\mathit{x}}=\underset{\mathit{x}}{\arg \min }\mathcal{F}\left(\mathit{x}\right),\mathrm{with}\mathcal{F}\left(\mathit{x}\right)=\frac{1}{2}\sum _{i=1}^{N_S}{\left|\left|{\tilde{\mathit{e}}}_i^{\mathrm{scat}}-{\mathbf{G}}_O\mathbf{X}{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right|\right|}^2+\mathcal{R}\left(\mathit{x}\right),$$

(8)

where the regularization term operates on differences of the contrast function at neighboring voxels:

$$\mathcal{R}\left(\mathbf{x}\right)=\gamma \sum _{i,j,i\sim j}\phi (x_i-x_j)\mathrm{with}\gamma >0,$$

(9)

where notation $i\sim j$ selects indices of neighbor voxels and $\phi$ is the “${\ell }_2{\ell }_1$” function $\phi \left(u\right)=\sqrt{{\delta }^2+u^2}$, where $\delta$ is a small positive constant ensuring the differentiability of the cost function, which was set to ${10}^{-2}$ in all our experiments. Such a penalization favors edge-preserving solutions and leads to a smooth cost function $\mathcal{F}$. Therefore, optimization can be addressed with gradient-based iterative algorithms. The penalization weight $\gamma$ trades off between fidelity to data and regularization, which is necessary since the problem is ill-posed [8,9]. Its value may depend on the noise level and on the problem size.

Note that alternate regularization strategies, such as multiplicative regularization [33,58], have also been proposed. Here, we chose additive regularization because it was previously shown to be particularly efficient in the case of highly contrasted objects [13]. However, the techniques presented in this paper may also be applied to reconstruction methods based on multiplicative regularization, since both the corresponding objective function and its gradient present similar algebraic characteristics.

2.4. Optimization and Computational Issues

We propose to solve the optimization problem (9) with the L-BFGS algorithm, which is an acknowledged method for solving high-dimensional convex optimization problems  ([59], see for example [13] in application to microwave imaging). Each iteration of L-BFGS consists of two steps: computing a descent direction, based on the computation of the gradient at the current point, and finding a convenient step size, for which several evaluations of the cost function may be necessary. The gradient of $\mathcal{F}$ reads:

$$\nabla \mathcal{F}\left(\mathit{x}\right)=-\mathbf{Q}\sum _{i=1}^{N_S}{\mathbf{J}}_i^{†}\left({\tilde{\mathit{e}}}_i^{\mathrm{scat}}-{\mathbf{G}}_O\mathbf{X}{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right)+\nabla \mathcal{R}\left(\mathit{x}\right)$$

(10)

where $\mathbf{Q}=[{\mathbf{I}}_N,{\mathbf{I}}_N,{\mathbf{I}}_N]$, the symbol ${}^{†}$ denotes the conjugate transpose and ${\mathbf{J}}_i$ is the Jacobian matrix associated to the i-th data misfit term in (8). It can easily be shown that:

$${\mathbf{J}}_i={\mathbf{G}}_O{\left({\mathbf{L}}_x^T\right)}^{-1}\mathrm{diag}\left\{{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right\},$$

(11)

where $\mathrm{diag}\left\{\mathit{u}\right\}$ stands for the diagonal matrix with diagonal $\mathit{u}$, such that the gradient reads:

$$\nabla \mathcal{F}\left(\mathit{x}\right)=-\mathbf{Q}\sum _{i=1}^{N_S}\mathrm{diag}{\left\{{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right\}}^{†}{\overline{{\mathbf{L}}_x}}^{-1}{\mathbf{G}}_O^{†}\left({\tilde{\mathit{e}}}_i^{\mathrm{scat}}-{\mathbf{G}}_O\mathbf{X}{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right)+\nabla \mathcal{R}\left(\mathit{x}\right),$$

(12)

where $\overline{{\mathbf{L}}_x}$ denotes the conjugate of ${\mathbf{L}}_x$.

Inspection of (8) and (12) reveals that the bulk of the computational effort required for evaluating $\mathcal{F}$ and $\nabla \mathcal{F}$ lies in the evaluation of matrix-vector products involving matrices $\mathbf{X}$ and ${\mathbf{G}}_O$, and in the resolution of linear systems with matrices ${\mathbf{L}}_x$ and $\overline{{\mathbf{L}}_x}$. Note that ${\mathbf{G}}_O$ is an $N_M\times 3N$ matrix and this manageable size makes it possible to explicitly store it for direct computations; similarly, $\mathbf{X}$ is a $3N\times 3N$ diagonal matrix, which results in small memory and computational requirements. Indeed, the major difficulty lies in the resolution of the $N_S$ systems of size $3N\times 3N$ involving matrix ${\mathbf{L}}_x$ (computation of ${\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}$), and in the resolution of another $N_S$ similar systems involving matrix $\overline{{\mathbf{L}}_x}$ (computation of ${\overline{{\mathbf{L}}_x}}^{-1}{\mathit{v}}_i$, with ${\mathit{v}}_i={\mathbf{G}}_O^{†}\left({\tilde{\mathit{e}}}_i^{\mathrm{scat}}-{\mathbf{G}}_O\mathbf{X}{\mathbf{L}}_x^{-1}{\mathit{e}}_i^{\mathrm{inc}}\right)$). Note that, in the frequent case where the number of sources exceeds the number of receivers, by virtue of the Lorentz reciprocity theorem (see for example [60]), one can equivalently switch the roles of sources and receivers, then reducing the number of systems to be solved.

In Section 3, we propose a dedicated block resolution strategy in order to jointly solve such multiple systems. This iterative method requires the computation of high-dimensional matrix-vector products involving matrix ${\mathbf{L}}_x$ (resp. $\overline{{\mathbf{L}}_x}$), that is, involving the $3N\times 3N$ matrix ${\mathbf{G}}_D$ (resp. $\overline{{\mathbf{G}}_D}$). Note that these products are convolution products, which can be efficiently computed with 3-D Fast Fourier Transform (FFT) algorithms [36,47].

3. Simultaneous Resolution of Multiple Forward Scattering Problems

We now address the joint resolution of the multiple forward scattering problems (4) by the block-BiCGStab algorithm proposed in [52]. We present the algorithm and we focus on implementation issues for solving one set of multiple scattering problems.

3.1. Description of the Block-BiCGStab Algorithm

The pseudo-code implementing the block-BiCGStab algorithm for the joint resolution of (4) is given in Algorithm 1. The notation ${\langle \mathbf{T},\mathbf{S}\rangle }_F$ represents the Frobenius product of two matrices. For comparison purposes, the sequential implementation with the BiCGStab algorithm, commonly used in microwave imaging, is recalled in Algorithm 2.

Algorithm 1 Joint resolution of $N_S$ linear systems ${\mathbf{L}}_x{\mathbf{E}}^{\mathrm{tot}}={\mathbf{E}}^{\mathrm{inc}}$ with block-BiCGStab.

1: For an initial matrix guess ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$, set ${\mathbf{R}}^{\left(0\right)}={\mathbf{E}}^{\mathrm{inc}}-{\mathbf{L}}_x{{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$, ${\mathbf{P}}^{\left(0\right)}={\mathbf{R}}^{\left(0\right)}$   2: Choose an arbitrary $3N\times N_S$ matrix ${\tilde{\mathbf{R}}}_0$   3: $k=0$   4: repeat 5:  $\mathbf{V}={\mathbf{L}}_x{\mathbf{P}}^{\left(k\right)}$   6:  solve $\left({\tilde{\mathbf{R}}}_0^{†}\mathbf{V}\right)\mathbf{A}={\tilde{\mathbf{R}}}_0^{†}{\mathbf{R}}^{\left(k\right)}$   7:  $\mathbf{S}={\mathbf{R}}^{\left(k\right)}-\mathbf{V}\mathbf{A}$   8:  $\mathbf{T}={\mathbf{L}}_x\mathbf{S}$   9:  $\omega ={\langle \mathbf{T},\mathbf{S}\rangle }_F/{\langle \mathbf{T},\mathbf{T}\rangle }_F$   10:  ${{\mathbf{E}}^{\mathrm{tot}}}^{(k+1)}={{\mathbf{E}}^{\mathrm{tot}}}^{\left(k\right)}+{\mathbf{P}}^{\left(k\right)}\mathbf{A}+\omega \mathbf{S}$   11:  ${\mathbf{R}}^{(k+1)}=\mathbf{S}-\omega \mathbf{T}$   12:  solve $\left({\tilde{\mathbf{R}}}_0^{†}\mathbf{V}\right)\mathbf{B}=-{\tilde{\mathbf{R}}}_0^{†}\mathbf{T}$   13:  ${\mathbf{P}}^{(k+1)}={\mathbf{R}}^{(k+1)}+({\mathbf{P}}^{\left(k\right)}-\omega \mathbf{V})\mathbf{B}$   14:  $k\leftarrow k+1$   15: until$\forall i,\parallel {\mathit{r}}_i^{\left(k\right)}\parallel /\parallel {\mathit{e}}_i^{\mathrm{inc}}\parallel <\mathrm{tolerance}$

Algorithm 2 Sequential resolution of $N_S$ linear systems ${\mathbf{L}}_x{\mathit{e}}_i^{\mathrm{tot}}={\mathit{e}}_i^{\mathrm{inc}}$, $i=1,\dots ,N_S$, with BiCGStab.

0: for$i=1,\dots ,N_S$do 1:  For an initial guess ${{\mathit{e}}_i^{\mathrm{tot}}}^{\left(0\right)}$,  set ${\mathit{r}}^{\left(0\right)}={\mathit{e}}_i^{\mathrm{inc}}-{\mathbf{L}}_x{{\mathit{e}}_i^{\mathrm{tot}}}^{\left(0\right)}$, ${\mathbf{p}}^{\left(0\right)}={\mathit{r}}^{\left(0\right)}$   2:  Choose an arbitrary vector ${\tilde{\mathit{r}}}_0$   3:  $k=0$   4:  repeat 5:   $\mathbf{v}={\mathbf{L}}_x{\mathbf{p}}^{\left(k\right)}$   6:   $\alpha =\left({\tilde{\mathit{r}}}_0^{†}{\mathit{r}}^{\left(k\right)}\right)/\left({\tilde{\mathit{r}}}_0^{†}\mathit{v}\right)$   7:   $\mathit{s}={\mathit{r}}^{\left(k\right)}-\alpha \mathit{v}$   8:   $\mathit{t}={\mathbf{L}}_x\mathit{s}$   9:   $\omega ={\mathit{t}}^{†}\mathit{s}/{\mathit{t}}^{†}\mathit{t}$   10:   ${{\mathit{e}}_i^{\mathrm{tot}}}^{(k+1)}={{\mathit{e}}_i^{\mathrm{tot}}}^{\left(k\right)}+\alpha {\mathit{p}}^{\left(k\right)}+\omega \mathit{s}$   11:   ${\mathit{r}}^{(k+1)}=\mathit{s}-\omega \mathit{t}$   12:   $\beta =\left(\alpha {\tilde{\mathit{r}}}_0^{†}{\mathit{r}}^{(k+1)}\right)/\left(\omega {\tilde{\mathit{r}}}_0^{†}{\mathit{r}}^{\left(k\right)}\right)$   13:   ${\mathit{p}}^{(k+1)}={\mathit{r}}^{(k+1)}+\beta ({\mathit{p}}^{\left(k\right)}-\omega \mathit{v})$   14:   $k\leftarrow k+1$   15:  until $\parallel {\mathit{r}}^{\left(k\right)}\parallel /\parallel {\mathit{e}}_i^{\mathrm{inc}}\parallel <\mathrm{tolerance},$   16: end for

In both algorithms, the most time-consuming operations concern the computation of matrix-vector products of the form ${\mathbf{L}}_x\mathit{p}=\mathit{p}-{\mathbf{G}}_D\mathbf{X}\mathit{p}$ in steps 5 and 8. Note that the same number ($2N_S$) of such products are performed by one iteration of block-BiCGStab and $N_S$ iterations of BiCGStab. In the remaining steps, a set of $N_S$ dot products and scalar-vector products (one per system inversion) in BiCGStab are replaced by one $N_S$-dimensional matrix-vector product in block-BiCGStab. Similarly, a set of $N_S$ scalar divisions are replaced by one $N_S\times N_S$ system inversion (steps 6 and 12). Therefore, such block operations generate a slight extra cost in block-BiCGStab compared to BiCGStab. Because of the small size of the involved systems, this remains negligible compared to the computation time required by multiplications with ${\mathbf{L}}_x$. Most of all, block operations introduce some coupling in the auxiliary variables and, subsequently, in the directions that are computed in order to update the solution at step 10. The rationale behind block versions is that coupling the descent directions yields more efficient update steps, and therefore results in a decrease in the total number of iterations and computing time.

3.2. Efficient Implementation for Solving Multiple Forward Problems

In its original description [52], the block-BiCGStab algorithm was evaluated on problems involving only a small number (5 to 10) of RHS vectors which were drawn independently, in moderate dimensions (up to $2500\times 2500$). In microwave imaging, the number of RHS vectors—the number of measurements—is much higher, and vectors are highly correlated, since they correspond to the different incident illuminating fields, which are spatially close to each other. Moreover, matrices involved in 3-D problems are usually much bigger than those used in [52].

Numerical experiments, whose salient results are reported in Section 4.2, revealed that the conditioning of matrix ${\tilde{\mathbf{R}}}_0^{†}\mathbf{V}$, which enters steps 6 and 12 of Algorithm 1, has a critical impact on the behavior and convergence of the block-BiCGStab algorithm. Inspection of steps 1 and 5 clearly indicates that the conditioning of $\mathbf{V}$ strongly depends on that of ${\mathbf{R}}^{\left(0\right)}$. With the standard initialization ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}={\mathbf{E}}^{\mathrm{inc}}$, when two sources i and j are close, the corresponding incident fields ${\mathit{e}}_i^{\mathrm{inc}}$ and ${\mathit{e}}_j^{\mathrm{inc}}$ which make up the corresponding columns of ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$, become highly correlated; hence, ${\mathbf{R}}^{\left(0\right)}$ becomes poorly conditioned and ${\tilde{\mathbf{R}}}_0^{†}\mathbf{V}$ is near-singular, which may lead to numerical errors propagating in the computations of matrices $\mathbf{A}$ and $\mathbf{B}$ (steps 6 and 12). This pathological behavior does not occur in the simulated experiments reported in [52], since only statistically independent RHS vectors are considered.

Here, we propose to add a small random quantity to ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$ in the initialization step 1. This reduces the correlation between columns of ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$, and thus improves the conditioning of matrix ${\tilde{\mathbf{R}}}_0^{†}\mathbf{V}$. Empirically, we have found that an initial perturbation with a signal-to-noise ratio ranging from 40 to 100 dB yielded roughly the same favorable effect, while a higher level of perturbation would slow down the convergence, and a smaller one would amount to no perturbation. In all experiments in this paper including such noisy initialization, the signal-to-noise ratio was set to 50 dB. The results reported in Section 4.2 show that the algorithm behaves in a satisfactory manner with such an initialization. Regarding the choice of matrix ${\tilde{\mathbf{R}}}_0$, various options were studied; the best results were obtained by using the choice proposed in [52], that is, ${\tilde{\mathbf{R}}}_0={\mathbf{R}}^{\left(0\right)}={\mathbf{E}}^{\mathrm{inc}}-{\mathbf{L}}_x{{\mathbf{E}}^{\mathrm{tot}}}^{\left(0\right)}$.

3.3. Implementation within the Reconstruction Procedure

We now consider the use of the block-BiCGStab algorithm for solving the two multiple linear systems involved in the evaluation of the cost function (8) and of its gradient (12) at each iteration of the optimization procedure. The first system reads ${\mathcal{S}}^{\left(t\right)}:{{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}={\mathbf{L}}_{x^{\left(t\right)}}^{-1}{\mathbf{E}}^{\mathrm{inc}}$, and the second one is of the form ${\overline{\mathcal{S}}}^{\left(t\right)}:{\mathbf{V}}^{\left(t\right)}={\overline{\mathbf{L}}}_{x^{\left(t\right)}}^{-1}{\mathbf{U}}^{\left(t\right)}$, where notation ${}^{\left(t\right)}$ indexes the iteration number. According to the expression of the gradient (12), the RHS term ${\mathbf{U}}^{\left(t\right)}$ in ${\overline{\mathcal{S}}}^{\left(t\right)}$ depends on the solution of ${\mathcal{S}}^{\left(t\right)}$, so the two systems can only be solved sequentially and not jointly. In the following, details are given considering ${\mathcal{S}}^{\left(t\right)}$, but similar arguments hold for ${\overline{\mathcal{S}}}^{\left(t\right)}$.

As discussed in Section 3.2, initialization of block-BiCGStab plays an important role in the algorithm efficiency. Our implementation is based on the fact that two consecutive iterates of the optimization algorithm are close to each other, especially when the reconstruction algorithm is near to convergence. Therefore, the solution ${{\mathbf{E}}^{\mathrm{tot}}}^{(t+1)}$ should be close to that obtained at the former iteration ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}$, and ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}$ could be used as an initialization of block-BiCGStab for solving ${\mathcal{S}}^{(t+1)}$. In particular, ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}$ should be much closer to ${{\mathbf{E}}^{\mathrm{tot}}}^{(t+1)}$ than the matrix formed by the incident fields ${\mathbf{E}}^{\mathrm{inc}}$, that was used for initializing block-BiCGStab for the resolution of a single set of forward problems in Section 3.2. We noticed, however, that implementing this strategy led to poor convergence rates, because the total fields ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}$ are highly correlated with each other, which generated conditioning issues similar to that explained in Section 3.2. Therefore, we propose to initialize iteration $t+1$ using a random perturbation added to ${{\mathbf{E}}^{\mathrm{tot}}}^{\left(t\right)}$.

4. Performance of Block-BiCGStab on Synthetic Data

The results presented in this section were obtained on synthetic data representing a typical experimental setup. First, we illustrate the impact of the implementation choices discussed in Section 3.2. Then, we investigate the performance of the proposed block-BiCGStab implementation for computing multiple forward scattering problems as a function of the problem difficulty, by varying the object contrast and the illuminating frequency. Finally, the computational efficiency of block-BiCGStab is evaluated on the full reconstruction process.

4.1. Description of the Numerical Example and Implementation Details

The object used in our experiments was composed of two nested cubes embedded in air. This example was taken from [33] and the contrast was increased by a factor of two in order to make the reconstruction problem more challenging. The object domain V was a cube with a 30-cm edge size, which was discretized onto a $30\times 30\times 30$ regular Cartesian grid ($1{\mathrm{cm}}^3$ voxels). It is usually acknowledged that the discretization grid must have at least four voxel sides per wavelength for accurate estimation of the scattered field. Here, the highest illuminating frequency is 3 GHz (see Section 4.3.1), that is, $\lambda =10$ cm, so that the voxel side is $\lambda$/10. The external part of the object was a cube with a 20-cm edge size, the contrast of which was set to ${\chi }_1=0.6+j0.8$. The internal part was a smaller cube with a 10-cm edge size and contrast ${\chi }_2=1.2+j0.4$. The object was illuminated by $N_S=160$ sources (vertical electric dipoles) that were regularly spaced around the object on five circles with a 60 cm diameter, at heights $z=-20$ cm, $-10$ cm, 0, 10 cm and 20 cm. Measurements were simulated by considering 160 receivers at the same locations as the sources. The object and the acquisition setup are represented in Figure 1.

Data were generated according to model (7) using a grid twice as fine as the one used in our numerical experiments ($60\times 60\times 60$ voxels) for more accurate computations, and circular Gaussian white noise was added in order to obtain a 20-dB SNR. Only the z component of the scattered field was used in our tests. Then, the reconstruction problem had ${30}^3$ = 27,000 complex unknowns and ${160}^2$ = 25,600 complex-valued data points.

All computations were performed with an Intel Core i7-5960X (eight cores) clocked at 3 GHz using Matlab, with multithreaded operations. For BiCGStab, a marching-on-in-source procedure [37] was used when it improved convergence. This procedure allows for better initialization of the solutions by using an interpolation based on the total fields previously computed for the nearest sources. The tolerance on the relative residuals imposed for convergence of both BiCGStab and Block-BiCGStab (line 15 in Algorithms 1 and 2) was set to ${10}^{-6}$. Therefore, the solutions obtained by the two algorithms are the same up to such numerical precision.

4.2. Impact of Initialization

The solution to the 3D forward problem corresponding to the experimental setup presented in Section 4.1 was computed with both BiCGStab and block-BiCGStab algorithms. First, the solutions were initialized to the incident fields (${{\mathit{e}}_i^{\mathrm{tot}}}^{\left(0\right)}={\mathit{e}}_i^{\mathrm{inc}}$). Evolution of the iterated relative residuals $\parallel {\mathit{e}}_i^{\mathrm{inc}}-{\mathbf{L}}_x{{\mathit{e}}_i^{\mathrm{tot}}}^{\left(k\right)}\parallel /\parallel {\mathit{e}}_i^{\mathrm{inc}}\parallel$ is shown in Figure 2a,b. BiCGStab converges in between 18 and 21 iterations for all problems, whereas block-BiCGStab residuals first slowly decrease, and then increase, leading to divergence.

Then, as described in Section 3.2, noise with 50 dB SNR was added to the initial solution of block-BiCGStab. The corresponding iterates are shown in Figure 2c—BiCGStab iterates, which are not represented, behave similarly to those represented in Figure 2a. On the contrary, block-BiCGStab iterates now converge much faster: in this example, block-BiCGStab converges in 40% less iterations than BiCGStab, corresponding to a similar saving in terms of computation time.

In all remaining experiments in the paper, the block-BiCGStab algorithm was initialized by adding such a noise to the incident fields.

4.3. Performance for Different Forward Scattering Configurations

We now compare BiCGStab and block-BiCGStab for several configurations of the 3-D forward problem based on the previous setup. The tests consisted of solving the $N_S$ systems $({\mathbf{I}}_{3N}-{\mathbf{G}}_D\mathbf{X}){\mathit{e}}_i^{\mathrm{tot}}={\mathit{e}}_i^{\mathrm{inc}}$, and we focused on the influence of the illuminating frequency and of the contrast values in $\mathbf{X}$, since both factors impact the problem difficulty—the Born approximation, for example, is only valid for low-contrast objects with small electric size.

4.3.1. Impact of the Frequency

We first study the algorithmic performance as a function of the illuminating frequency. The object defined in Section 4.1 was used, with frequencies equal to 1, 2 and 3 GHz. The corresponding wavelengths are respectively $\lambda =30$ cm, $\lambda =15$ cm and $\lambda =10$ cm, and the size of the reconstructed volume respectively corresponds to $\lambda$, $2\lambda$ and $3\lambda$: the higher the frequency, the bigger the electric size of the object.

Table 1. CPU time and number of iterations required for computing the 160 forward simulations with BiCGStab and block-BiCGStab for three different frequencies. Results are averaged over the 160 resolutions for BiCGStab and over 20 realizations of the random initialization for block-BiCGStab.

Frequency	CPU Time		Number of Iterations		CPU Time Ratio Block/BiCGStab
			[min, Average, max]
	BiCGStab	Block	BiCGStab	Block
1 GHz	116 s	110 s	[7 8.0 10]	[7 7.9 9]	0.95
2 GHz	187 s	138 s	[13 13.4 15]	[10 10 10]	0.74
3 GHz	245 s	164 s	[17 17.9 21]	[12 12 12]	0.67

gives the corresponding central processing unit (CPU) time and the number of iterations required for BiCGStab and block-BiCGStab to converge. First, in accordance with the discussion in Section 3.1, the CPU time is almost proportional to the number of iterations performed by both algorithms. Second, block-BiCGStab always yields a shorter computation time than BiCGStab, and the improvement becomes more significant as the frequency increases: at 3 GHz, block-BiCGStab is 33% faster than BiCGStab for solving the 160 systems.

4.3.2. Impact of the Contrast

We now evaluate the computation time as a function of the contrast. A scaling factor ranging from $0.25$ to $2.5$ was applied to the synthetic object defined in Section 4.1. The illuminating frequency was set to $3$GHz.

Table 2. CPU time and number of iterations required for computing the 160 forward simulations with BiCGStab and block-BiCGStab for different contrast levels. Results are averaged over the 160 resolutions for BiCGStab and over 20 realizations of the random initialization for block-BiCGStab.

Contrast	CPU Time		Number of Iterations		CPU Time Ratio Block/BiCGStab
			[min, Average, max]
	BiCGStab	Block	BiCGStab	Block
0.25	105 s	86 s	[7 7.01 8]	[6 6 6]	0.82
0.5	152 s	120 s	[10 10.3 12]	[8 8.5 9]	0.79
1	245 s	164 s	[17 17.9 21]	[12 12 12]	0.67
1.5	361 s	217 s	[25 26.9 29]	[15 16.1 18]	0.60
2	482 s	271 s	[34 36.1 39]	[19 20.2 21]	0.56
2.5	610 s	331 s	[44 45.9 50]	[23 24.8 26]	0.54

gives the corresponding CPU time and the number of iterations required for BiCGStab and block-BiCGStab. Here again, the CPU time is almost proportional to the number of iterations. In all cases, block-BiCGStab is faster than BiCGStab and, more importantly, the relative saving in CPU time increases with the contrast of the object: for a contrast factor of 2.5, block-BiCGStab runs approximately twice as fast as BiCGStab.

The results reported in Section 4.3.1 and Section 4.3.2 therefore suggest that block-BiCGStab is best suited for the resolution of difficult problems, i.e., involving large and/or highly contrasted objects.

4.4. Inversion Procedure and Object Reconstruction

We now consider the reconstruction problem with the simulation setup described in Section 4.1. Since the object was highly contrasted, illuminations at frequencies 1, 2 and 3 GHz were used and the following frequency hopping technique was implemented (see [28,61,62] for example):

• the initial solution was set to zero, and $N_{\mathrm{it}}$ iterations of the inversion algorithm were performed with the 1 GHz data;
• another $N_{\mathrm{it}}$ iterations were then performed with the 2 GHz data, with an initial solution set to the output of the previous step;
• the reconstruction algorithm was finally applied to 3 GHz data, with an initial solution set to the output of the previous step, and iterations are run until the ${\ell }_{\infty }$-norm of the gradient becomes lower than ${10}^{-6}$.

The number $N_{\mathrm{it}}$ of iterations for each intermediate frequency is a trade-off between reconstruction quality and computation time. In this experiment, we found that $N_{\mathrm{it}}=20$ iterations were sufficient to achieve satisfactory results (that is, more iterations at intermediate frequencies would not improve the final reconstruction and would unnecessarily increase the computation time). The regularization weight in (9) was set to $\gamma ={10}^{-6}$ for the three frequencies.

Reconstruction was performed using BiCGStab and block-BiCGStab, following the procedure proposed in Section 3.3. Since both algorithms used the same convergence thresholds, evaluations of the objective function and of its gradient at any point by the two algorithms were nearly identical; therefore, the minimization process followed roughly the same sequence, and the reconstructed objects obtained with both approaches are comparable. Reconstruction results are shown in Figure 3, together with the true object. The shapes of the two cubes are well retrieved and the contrast of the outer cube (${\chi }_1=0.6+j0.8$) is very close to that of the true object. The boundaries of the reconstructed cubes slightly exceed the true ones, especially for the inner cube. The contrast value in the inner cube is also slightly misestimated (underestimation of the real part and overestimation of the imaginary part). Note that the inner cube is highly contrasted (${\chi }_2=1.2+j0.4$), which makes it very difficult to reconstruct.

For this problem, the reconstruction lasted about 21.3 h when computations were performed with BiCGStab. The reconstruction time dropped to 16.8 hours when block-BiCGStab was used, which represents a saving of more than 20%.

5. Reconstruction of the Fresnel Database Objects

We finally consider the reconstruction of the five objects of the Fresnel database [54]. Such objects have already been reconstructed in many papers [27,28,29,61,63,64]—see also a summary of reconstruction results in [65]. Due to space limitations, we restrict the presentation to two objects, which can be considered as the two most difficult ones: the TwoSpheres and the Cylinder objects.

The 3-D experimental setup for the Fresnel data set is detailed in [54]. Using the reciprocity theorem and the two polarizations of the experimental data [29], we consider that the setup was equivalently composed of 36 sources placed on a circle at $z=0$ with a 1.796-m radius. Data were equivalently acquired with 81 receivers placed on a sphere with a 1.796-m radius [54]. The incident fields were plane waves polarized along the z-axis and acquisitions were performed at frequencies ranging from 3 to 8 GHz. The two polarizations were used (co- and cross-polarization). Since the measurements were performed in the far-field zone of the scattering object, the electric field at the receivers is transverse (that is, its longitudinal component is zero). When projected into the Cartesian coordinate system, the three spatial components of the measured field are generally non-zero. Thus, each data set is made up of $36\times 81\times 3=8748$ complex-valued measurement points. For all objects, the volume of interest V was defined as a cube with a 10 cm edge size divided into $30\times 30\times 30$ voxels, so that the object is fully enclosed in V in each case. The maximum frequency is 8 GHz, that is, $\lambda =37.4$ mm, so that the coarser discretization scheme corresponds to 11.3 voxel sides per wavelength. We then have 27,000 complex-valued unknowns. In all experiments, the regularization parameter $\gamma$ in (9) was set to $\gamma ={10}^{-6}$.

Here also, as discussed in Section 4.4, a frequency hopping procedure was used to carry out the reconstructions as follows:

• Perform $N_{\mathrm{it}}$ iterations of the reconstruction algorithm with the 3 GHz data and initial contrast function set to zero.
• Perform $N_{\mathrm{it}}$ iterations of the reconstruction algorithm with the 4 GHz data and initial contrast function set to the final iterate of the previous step.
• Repeat the previous step with the 5 GHz, 6 GHz and 7 GHz data.
• Perform reconstruction of the object with the 8 GHz data, initial contrast function set to the final iterate of the previous step and a stopping rule on the ${\ell }_{\infty }$-norm of the gradient set to ${10}^{-6}$.

For all objects except Cylinder, computing $N_{\mathrm{it}}=10$ iterations at each frequency empirically led to a satisfactory compromise between reconstruction quality and time. For Cylinder, which was bigger and more contrasted, $N_{\mathrm{it}}=20$ intermediate iterations were necessary in order to achieve the best reconstruction.

The TwoSpheres object was composed of two identical spheres in contact embedded in air, 50 mm in diameter and with uniform, real-valued, relative permittivity equal to $2.6$, that is, contrast $\chi =1.6$. Figure 4 (top) shows 30 slices of the true object and of our corresponding reconstruction, and Figure 4 (bottom) represents the isosurface of the reconstruction at the contrast value $\chi =1$. Note that we only present the real part of the estimated contrast function because the estimated imaginary part was very small. The latter is presented in the supplementary data files. It can be observed that the shape of the object was well reconstructed, even at the contact point. The contrast value was slightly overestimated at the contact point (the estimated contrast was about 1.8), and slightly underestimated around it (the estimated contrast was about 1.3).

The Cylinder object was composed of a cylinder of 80 mm in length and 80 mm in diameter, with uniform relative permittivity equal to $3.05$ (contrast $\chi =2.05$). It was the most difficult object of the database because of its large size and high contrast, and several methods failed to reconstruct this object [28,33,63,66], while many others produced poor quality results [27,29,61,62,64]. Figure 5 shows our reconstruction results. The shape and the contrast were adequately reconstructed: the contrast function was homogeneous inside the object and the estimated value at the center voxels was close to the true one, since the estimated permittivity is about 1.9. To the best of our knowledge, this reconstruction of Cylinder outperformed all other results previously reported in the literature.

On a broader perspective,

Table 3. Subjective comparison of the obtained reconstruction results with previously published works on Fresnel database objects: “=”, “+” and “++” respectively mean that our reconstruction is of comparable, better, and much better quality. N/A means that reconstruction of the object is not presented in the corresponding paper.

Paper Ref.	TwoCubes	IsocaSphere	CubeSpheres	TwoSpheres	Cylinder
[63]	++	=	=	++	N/A
[29]	++	=	+	++	++
[28]	+	+	=	+	N/A
[64]	++	+	+	++	++
[61]	+	=	+	=	+
[27]	+	=	=	=	++
[33]	+	=	N/A	N/A	N/A
[66]	+	N/A	N/A	=	N/A
[62]	N/A	N/A	+	=	++

provides a heuristic assessment of the reconstruction quality for all objects in the Fresnel database. Each cell in the table contains a subjective comparison of our result with previously obtained ones for each object, based on the figures in the corresponding papers. We consider that the proposed method provided results at least similar to previous ones, and often better for the most difficult objects. The reader is referred to supplementary files associated with the paper, that depict graphical representations of all reconstructed objects and provide all reconstruction results as Matlab files.

Table 4. Reconstruction time for the five objects of Fresnel database using BiCGStab and block-BiCGStab for linear system inversions.

Object from	Time	Time	Time Ratio
Database	BiCGStab	Block-BiCGStab	Block/BiCGStab
TwoCubes	2.33 h	2.07 h	0.89
IsocaSphere	4.38 h	3.58 h	0.82
CubeSpheres	3.21 h	2.60 h	0.81
TwoSpheres	8.60 h	6.41 h	0.75
Cylinder	60.0 h	28.3 h	0.47

finally gives the reconstruction time for the five objects of the Fresnel database using either BiCGStab or block-BiCGStab. The time saved by block-BiCGStab becomes more important as the difficulty of the reconstruction problem increases (i.e., as objects become larger and/or more contrasted). For the most difficult Cylinder object, the reconstruction time drops from 60 h to 28.3 h, which means that block-BiCGStab is more than twice faster than BiCGStab.

6. Discussion

We proposed a reconstruction strategy based on an additive regularization framework, where optimization is performed via first-order optimization. Computations of the cost function and of its gradient were performed through a dedicated implementation of the block-BiCGStab algorithm for multiple RHS systems inversion. In all our experiments, such an implementation led to a significant reduction of the computation time (up to one half) compared to BiCGStab, the gain being all the more significant that the reconstruction problem is difficult and the computation time is high. The procedure proposed in this paper is therefore particularly attractive for problems dealing with high-frequency illuminations and/or highly contrasted objects, which are of great interest: high-frequency acquisitions contain more detailed spatial information about the object under study, yielding higher resolution, and highly contrasted objects may be encountered in many application areas, such as biomedical imaging. Consequently, the proposed procedure may contribute to reducing the computational burden of microwave imaging, which still represents a crucial issue in the development of practical applications.

Improved convergence properties of block-BiCGStab over its sequential counterpart also enable block-BiCGStab to reach a low residual level for each iterative system inversion. In our experiments, the tolerance threshold on the norm of the residuals was set to ${10}^{-6}$, whereas higher values are commonly used. More accurate solutions to the forward problem yield more efficient optimization steps; this may explain why the proposed procedure achieved satisfactory reconstruction results for the most difficult objects of the Fresnel database, for which no result of equivalent quality has been reported in the literature.

Exploiting the proposed block resolution strategy may also improve the efficiency of any other microwave imaging inversion method requiring the computation of similar quantities (i.e., the total fields from the incident fields), either based on additive [13] or multiplicative [33] regularization. This strategy can also be advantageously used whatever the choice of the optimization method, such as non-linear conjugate gradient [29,64], regularized Gauss–Newton [33], Levenberg–Marquardt [67], or distorted Born iterative method [27]. It would also be worth studying the numerical performance of the block-BiCGStab algorithm in different geometrical configurations concerning, e.g., data acquisition setups or discretization schemes. Finally, it may also be efficiently used in other applications involving the computation of the total electromagnetic field accounting for scattering, e.g., radar cross-section [48,68,69], eddy-current [26,70] or magnetotelluric [71] computations. Depending on the complexity of the problem, specific tuning of the block-BiCGStab algorithm could be required, for which the strategies proposed in this paper may help.

Further improvements may be brought by parallel implementation of the proposed procedure. In our experiments, multicore architecture was taken advantage of only by performing multithreaded operations at a low level within the block-BiCGStab loop. However, massive parallelization of the resolution of independent systems (with BiCGStab) reduces the computation time by exploiting parallelization at a higher level [33]. Indeed, the two parallelization levels could be combined into a block resolution procedure, where the multiple systems would be split into several subsystems to be solved in parallel. In [53], such an idea was successfully tested on an early version of block-BiCGStab. This remains a promising possibility to accelerate microwave imaging in the context of High Performance Computing, as well as parallelization capacities brought by GPU- or multi-GPU-based implementations.