International Journal of Antennas and Propagation

Volume 2015 (2015), Article ID 924067, 21 pages

http://dx.doi.org/10.1155/2015/924067

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

Department of Information Technology (INTEC), Ghent University, Sint-Pietersnieuwstraat 41, 9000 Gent, Belgium

Received 19 December 2014; Accepted 24 March 2015

Academic Editor: Amelie Litman

Copyright © 2015 Jürgen De Zaeytijd and Ann Franchois. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. Introduction

During the last decade significant progress has been reported in fully 3D quantitative microwave imaging algorithms, for example, [1–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–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–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 the “modified gradient” approach avoids these computationally demanding forward problem solutions, the “Newton-type” 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 piecewise-homogeneous 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 , with the Jacobian matrix, is obtained from linearizing the scattered field 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–35]. The resulting “least squares” update system with a dense (approximate) Hessian matrix 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 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 . 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 single-frequency 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.

#### 2. Formulation

The object with unknown relative complex permittivityis immersed in a homogeneous background medium with relative complex permittivity , where is the relative permittivity, the conductivity, the permittivity of vacuum, the 3D position vector, and 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 (see Figure 1(a)), scattering data are collected by successively illuminating the object with different known incident electric fields , which are radiated by elementary electrical dipoles positioned in points exterior to and with polarizations along unit vectors ; for each illumination the resulting scattered electric field components are measured in positions along unit vectors [8]. All these (complex scalar) scattered field components are stored in one -dimensional measured field vector , with . Note that (resp., ) refers to a unique combination of a position (resp., ) and a polarization (resp., an orientation ).