Sparse Reconstruction for Near-Field MIMO Radar Imaging using Fast Multipole Method

Radar imaging using multiple input multiple output systems are becoming popular recently. These applications typically contain a sparse scene and the imaging system is challenged by the requirement of high quality real-time image reconstruction from under-sampled measurements via compressive sensing. In this paper, we deal with obtaining sparse solution to near- field radar imaging problems by developing efficient sparse reconstruction, which avoid storing and using large-scale sensing matrices. We demonstrate that the “fast multipole method” can be employed within sparse reconstruction algorithms to efficiently compute the sensing operator and its adjoint (backward) operator, hence improving the computation speed and memory usage, especially for large-scale 3-D imaging problems. For several near-field imaging scenarios including point scatterers and 2-D/3-D extended targets, the performances of sparse reconstruction algorithms are numerically tested in comparison with a classical solver. Furthermore, effectiveness of the fast multipole method and efficient reconstruction are illustrated in terms of memory requirement and processing time.


I. INTRODUCTION
Radar imaging has many applications such as subsurface or behind wall imaging, improvised explosive device detection, and collision avoidance and has been of interest in the literature, recently [1][2][3][4]. For high resolution imaging systems, planar arrays with a large number of antenna elements is used and the antenna is generally a multiple-input-multiple-output (MIMO) array composed of spatially distributed transmitting and receiving sub-arrays, operating sequentially or simultaneously. In many of these applications the target to be imaged lies in the near field of the antenna array. The spatial diversity of the sub-arrays, in association with wideband operation, provides high resolution and image quality with few number of antennas [5][6].
A near-field MIMO imaging problem can be numerically modeled as a linear inverse problem whose solution contains the position, shape, and reflectivity distribution of the target. However, the solution is usually not unique due to underdetermined nature of the problem. Besides, ill-posed structure of the linear system makes it difficult to reconstruct high quality images from noisy measurements. In order to generate useful images, additional regularization approaches must be considered [7].
In the last two decades, researchers have developed many regularization methods to overcome shortcomings of the inverse problem and yield improvement at the quality of its solution. One fundamental approach that has been suggested as a regularizer, particularly for compressive sensing (CS)based imaging settings [8][9][10][11]27,28,39], is sparsity constraint. The sparsity constraint provides an approximate solution that contains few non-zero entries when compared to its dimension. This approach is also known as sparse approximation [12][13] and mathematically expressed as where is referred to as regularization (or cost) function. Depending on the choice of , reconstruction characteristic of the imaged space varies. The most basic approach is to generate a maximally sparse representation by selecting ( ) = ‖ ‖ 0 where ‖ • ‖ 0 ∶ ℝ → ℝ is ℓ0 "norm", which gives the number of non-zero entries of a vector in ℝ .
Orthogonal matching pursuit (OMP) is a greedy method based on ℓ0 "norm". It approximates a sparse solution by iteratively selecting a column of the sensing matrix, which contribute to the sparsity most. The iterations are continued until a predetermined sparsity level is reached. Stagewise orthogonal matching pursuit (StOMP) is another type of greedy method based on OMP. It simply differs from OMP by selecting multiple columns at each iteration, namely stage, which makes it converge faster than OMP. Regularized orthogonal matching pursuit (ROMP) is a modified form of OMP that does not possess a threshold value for sparsity level; instead it selects the columns having similar dot products with the solution vector. Despite the fact that these methods are able to converge in very short runtime, they generally do not offer any guarantee for a sparse solution [14][15][16]. Moreover, the complexity of the search for all possible sparse subsets is generally exponential in the number of columns, which makes the solution of (1) with ℓ0 "norm" NP-hard [17]. Greedy methods work under specific conditions and does not require exhaustive search.
Another useful approach for sparse approximation is to replace the ℓ0 "norm" with ℓ1-norm, which converts the inverse problem into a convex optimization problem. Alternating direction method of multipliers (ADMM) is based on ℓ1-norm and has been used for various imaging problems [18][19][20][21]. It is a form of augmented Lagrangian method (ALM), which handles solution of the optimization problem more efficiently by dividing it into smaller components. ADMM splits primal variables, augments the Lagrangian of the optimization problem (as in method of multipliers) and carries out iterative variable minimization steps. There also exist split augmented Lagrangian shrinkage algorithm (SALSA) [21,22] and constrained-SALSAs (C-SALSA-1/C-SALSA-2) [19]. These techniques transform the unconstrained expression of the problem into a constrained one by performing variable splitting and using an ALM, specifically ADMM.
Tikhonov regularization, which is based on ℓ2-norm, is a simple tool specialized for regularizing ill-posed problems [23]. Although it regularizes conditioning of problems, it generally does not offer a sparse solution. In [24], two sparse reconstruction algorithms are proposed based on generalized Arnoldi-Tikhonov regularization, approximating ℓ1-norm and total variation (TV) in terms of ℓ2-norm. These algorithms can provide a sparse solution for the problem within fewer iterations than classical Tikhonov regularization algorithm at the cost of higher relative error.
All the sparse approximation methods mentioned above attempt to solve the inverse problem through iterative matrixvector multiplications, which are equivalent to computing the (forward) sensing operator and its (backward) adjoint operator. As the imaging problem gets electrically larger, solution of the forward problem becomes difficult to store and requires immense computational resources with complexities of ( 2 ) and ( 2 ) for processing time and memory, respectively, where is the number of unknowns and is the iteration count. In the literature, therefore, much effort has been made, searching for methods to reduce computational complexity of the 3-D imaging problems. For instance, in [25], authors propose an accelerated algorithm based on Bayesian learning and approximate message passing to solve large-scale electrical impedance tomography problem. Commonly, diagonalizability of the sensing operator and fast Fourier transforms are exploited for efficient large-scale sparse reconstruction, for example in optical volumetric imaging [26] and compressive spectral imaging [27]. Similarly, for wideband near-field radar imaging, in [28], efficient CS reconstruction is achieved by decomposing the sensing operator into Fourier transform and sampling operations, and in [29], sensing operator is computed using an interpolationfree holographic imaging algorithm that also involves fast Fourier transforms. These approaches reduce memory usage and computation time for a near-field radar imaging with a monostatic configuration.
Here we demonstrate that the "fast multipole method" can be employed within sparse reconstruction algorithms to efficiently compute the sensing operator and its adjoint (backward) operator, for a general multistatic imaging setting. The forward part of the radar imaging problem is an electromagnetic scattering problem and the fast multipole method (FMM) is a powerful tool to efficiently calculate the large matrix-vector products for the forward (sensing) operator. In fact, FMM can be used to calculate the matrixvector product without forming the sensing matrix explicitly and has the ability to decrease the computation time to ( 3/2 ) and reduce memory requirement to ( 3/2 ) [30][31][32][33].
In this paper, we seek accelerated sparse solution to nearfield multiple-input-multiple-output (MIMO) radar imaging problem. For this purpose, we initially construct the problem as a convex optimization problem with sparsity (ℓ1-norm and TV regularization) and solve it by using the augmented Lagrangian framework. Then, for large-scale imaging problems, we propose employing the FMM formulation to efficiently compute the matrix-vector multiplications within the reconstruction algorithm. Hence, the novelty of this study is to develop efficient sparsity-based reconstruction methods that exploit FMM formulation for matrix-vector multiplications involved, enabling significant reduction in computation time and memory requirement. Several numerical simulations are carried out for 2D and 3D objects as well as point scatterers in order to validate the study. The effectiveness of the FMM is also demonstrated in terms of memory usage and computation time by comparing it with classical direct and iterative linear system solvers.
This paper is organized as follows. Section II introduces imaging configuration and its representation as a linear system. In section III, application of the FMM to solution of inverse problem is described briefly. Section IV presents mathematical basis of inverse problem. Section V deals with the algorithmic approaches for numerical solution of the inverse problem. In section VI, numerical simulations carried out for various near-field imaging scenarios are given and discussed. Section VII concludes the paper.

II. FORWARD PROBLEM
The construction of the forward problem is the key to establish the mathematical relation between the reflectivity distribution of the imaging scene and the measurement data, and hence is essential for the development of the imaging algorithm. The solution of the inverse problem can then be carried out by iterations and this requires solving the forward problem at each iteration step. Therefore, a fast and robust algorithm that computes the forward (sensing) operator is required for efficient and accurate image reconstruction.
In this study, an ultra-wideband (UWB) near field MIMO array imaging structure, as introduced in [34,35], is considered. We will assume that the data is collected by transmitting sequentially from each element of the transmit array and receiving simultaneously by the receiving array elements of the MIMO array. The imaging structure is depicted in Fig. 1. A two dimensional plus-shaped MIMO array located in = 0 plane is used. The receiving antennas are positioned along the x-axis and their locations are denoted by ( , 0, ) while the transmitting antennas are placed along the z-axis and their locations are denoted by ( , 0, ). Despite the vector nature of the scattering phenomenon, we will use a simple scalar model that is commonly used in the literature [36]. Applying the Born approximation, the scattered field at the corresponding receiving antenna position due to any scatterer in the imaging volume ( , , ) is mathematically expressed as [30,31] ( , , , , where ( ) is the transmitted pulse in time domain, is speed of light, and ( , , ) is the three dimensional reflectivity distribution function of the target. and denote the distances to the point ( , , ) from the transmitting and the receiving antennas, respectively, and can be written as By taking temporal Fourier transform of (2), the received signal can be written in the frequency domain as where ( ) is the Fourier transform of the transmitted pulse with being the wavenumber. Considering a computerized image reconstruction process, continuous expressions in (2) and (5) can be discretized by expressing the three dimensional reflectivity distribution function in terms of voxels, thus, (5) can be written as: where is number of voxels and is reflectivity of the th voxel. The transmitting antennas radiate discrete frequencies, which are equally spaced by frequency step of ∆ , in the operational bandwidth.
, denotes the distance from the th transmitting antenna to the center of the th voxel and , denotes the distance from the center of the th voxel to the th receiving antenna. Furthermore, the number of transmitting and receiving antennas are and , respectively.
Note that the scattered field expressed in (6) is obtained by ignoring multiple reflections among the voxels. The transmitted pulse is assumed to be directly reflected by each voxel to the corresponding receiving antenna without any contribution from the rest of the voxels.
The discrete model defines a linear system as given in (7). The reflectivity values of the voxels are organized in the vector in a lexicographic order and the measurements are listed in the same order in the right-hand-side vector . The matrix ∈ ℂ × is the sensing (system) matrix and its total number of rows , is equal to × × , whereas the number of columns , is equal to .
The open form of this linear system is given in (8).
For the solution of such a linear system using an iterative algorithm, × × × multiplications must be performed for one matrix-vector product. This means that as the dimension of the unknown vector increases, the solution becomes quite inefficient due to the requirements for excessive amount of computational resources and memory.

III. FAST MULTIPOLE METHOD
In this study, we solve the imaging problem by iterative methods, where each iteration requires the computation of the forward operator (multiplication of a vector with the sensing matrix ) and its adjoint (multiplication of a vector with the matrix ). As expressed in (8), entries of the sensing matrix are product of two Green's functions, which is the fundamental factor that allows us to use FMM for the solution of the forward problem through a two stage process. The first Green's function represents the path that the transmitted signal travels from the transmitting antenna to the target and it is treated by the first stage of FMM. The second Green's function represents the propagation of the reflected signal back to the receiving antenna and the second stage of the FMM is applied to this part. Details of both stages are given below.

A. THE FIRST STAGE
As depicted in Fig. 2(a), we assume that only one transmitting antenna operates at a time. The aggregation step, therefore, is not applied in this stage. The transmitted signal is directly translated from th transmitting antenna to the geometric center of the imaging domain and then distributed to each voxel by disaggregation step. This phenomenon is mathematically expressed as ( ,̂, ) = ̂⋅ (11) The quantities in (10) and (11) are referred to as the translation function and the disaggregation function, respectively, where is the order of truncation for the translation function, ℎ (1) ( ) is the spherical Hankel function of the first kind and th order, and ( ) is the Legendre polynomial of order . The excess bandwidth formula (EBF) is used to determine the truncation order [37,38]. Integration is over a unit sphere and it is implemented numerically by using Gauss-quadrature rule. The sample points for the integration along elevation and azimuth axes are determined by using the truncation order. Note that as many translation functions as the number of transmitting antennas are required at this stage.
The number of operations carried out in this stage is × × since (9) must be repeated for each voxel and receiving antenna positions as well as each frequency steps. Fig. 2(b) demonstrates that the scattered signals (excluding the multiple reflections) from every voxel of the target is recollected at geometric center of the image volume in aggregation step. The total signal is then translated back to geometric center of the receiving antennas, and, the received signal is redistributed among the receiving antennas through disaggregation step. Mathematical expression of this stage is given as follows

B. THE SECOND STAGE
is the aggregation function, is the distance from the center of the image volume to the geometric center of the receiving antennas, is the distance from th voxel to the center of the image volume, and is the distance from geometric center of the receiving antennas to the th receiving antenna. This stage lasts shorter than the first stage since it contains only one translation. Note that (12) must be repeated for each voxel position, receiving antenna position and frequency steps. Hence, the total number of operations required in this stage is × × . The above explanations describe how the forward (sensing) operator is computed using FMM formulation. The same procedure can also be applied for the computation of the adjoint (backward) operator by simply exchanging order of the stages, i.e., in the 1 st stage, aggregation, translation and disaggregation operations are carried out, while translation and disaggregation operations are performed in the 2 nd stage.
Consequently, FMM reduces the overall operations needed for the evaluation of a single matrix-vector product to ( + ) × × as compared to classical direct muliplication, which requires × × × operations. This reduction leads to an efficient solution in terms of processing time and memory for the solution of imaging problems where relatively large antenna arrays are used.

IV. INVERSE PROBLEM
This section deals with estimation of reflectivity distribution of the imaged scene from measurements using the MIMO array. The reflectivity of a physical object is generally a smooth function and its values at neighboring points are highly correlated along down-and cross-range directions, which allows it to be represented in a sparse form when an appropriate regularization function ( ) is introduced [40]. The sparsity-inducing regularization problem is formulated in unconstrained form as where > 0 is regularization parameter, which controls the sparsity and determines the trade-off between the regularization and data fidelity [14]. In particular, should be increased as signal-to-noise ratio (SNR) decreases. An alternative form is the constrained problem of where is error tolerance, which is introduced when the measurement data is noisy and its value is estimated using SNR of the received signal. In (14) and (15), (•) acts as a transform operator, which can be selected with respect to the physical characteristic of the imaged scene, e.g., it can be selected as ℓ1-norm (‖ ‖ 1 ) for point scatterers with weak background, whereas common choice for two-and threedimensional extended targets is discrete gradient operator of the following form where , , and corresponds to difference operators along -, -, and -axis, respectively. This operation is called total variation ( ) of the unknown and leads to a reconstruction where sharp edges and rapidly changing structures are preserved on the reconstructed image [40]. The mathematical expression of three-dimensional is [18] ( where [ , , ] denotes the reflectivity belonging to the , , th voxel of the imaging scene.

V. RECONSTRUCTION ALGORITHMS
We have considered two convex optimization based algorithms for sparse reconstruction, namely, ADMM and C-SALSA-2 and they are explained in the following subsections.

A. ADMM
One approach to solve the problem in (14) is to employ the ADMM whose algorithmic steps are given above. In each iteration of the algorithm, and vectors are updated in order by applying alternating minimization of augmented Lagrangian.
In a large-scale imaging problem, 3 rd step is computationally most demanding part of the algorithm, since it solves a minimization problem. Here we solve this minimization problem iteratively with an inner iteration (e.g. using conjugate gradient least squares (CGLS)), where multiplication of a vector with the sensing matrix and its Hermitian is efficiently carried out through FMM. In the literature, the optimal value of is found by a search method (e.g. cross-validation), which requires solution of (14) multiple times to obtain its best value [41,42]. This increases the computation time by the number of trials for different values. This is another reason to employ the FMM formulation to accelerate the algorithm.

B. C-SALSA-2
In [19], C-SALSA is proposed to solve the constrained sparsity-regularized problem in (15) by adding an indicator function to its objective function, which basically transforms it into an unconstrained problem. Then, resulting problem is transformed into another problem by variable splitting and it is solved by ADMM. The parameter selection is simpler in the constrained problem since (15) does not possess the regularization parameter .
The algorithmic steps for C-SALSA-2 are given above. As can be seen, 3 rd , 4 th , 6 th , and 8 th steps of C-SALSA-2 involve the product of a vector with the sensing matrix and its Hermitian . Here we compute these operations efficiently using FMM as described in Section III. The computationally most demanding step of the algorithm is the 4 th step where the inverse of ( + ) is calculated. This step can be rewritten as a matrix equation in +1 and since is known from the 3 rd step, the linear system of equations in (20) can be solved iteratively again, enabling to incorporate efficient implementations of multiplications with or . For an imaging problem with an unknown vector of high dimension, solving (20) iteratively is still not feasible by forming the large-scale sensing matrix and computing the necessary matrix-vector products. Here, we accelerate this computation by applying two consecutive FMMs to obtain ( +1 ), i.e., 1 st FMM is for +1 and 2 nd FMM is for (•), which provide a great reduction in computation of 4 th step. Furthermore, in practice, the discrete gradient operators ( , , and ) are also not required to be formed explicitly since they can be computed by filtering the unknown vector with appropriate derivative kernels, e.g. [−1 0 1], along all directions.

VI. NUMERICAL RESULTS
Several near-field imaging problems are considered to investigate the performance of augmented Lagrangian approach and the related algorithms described above. As sketched in Fig. 3(a), a plus-shaped MIMO array is used in our imaging scenarios, containing equally spaced 30 transmitting antennas and 30 receiving antennas. Aperture of the array is 0.725 m × 0.725 m and spacing between the antennas is 0.025 m. The target is located 0.55 m away from the center of the array (which is assumed to be the nine point scatterers depicted in Fig. 3(b)). Operational frequency bandwidth ranges from 7 GHz to 13 GHz with 7 frequency steps of 1 GHz. It must be noted that number of frequency steps are deliberately kept low so that we can analyze how the algorithms and the sparsity constraint perform for underdetermined nature of the problems.
bandwidth. Vertical (along -axis) and horizontal (alongaxis) cross-range "Rayleigh" resolutions are given by and respectively. , and , are widths of the transmitting and receiving arrays along -axis, whereas , and , are widths along -axis [34]. The down-range "Rayleigh" resolution is given by The down-range resolution of the setup shown in Fig. 4(a) is 0.025 m. The cross-range resolutions calculated by (28) and (22) are equal to 0.025 m and they correspond to the definition of "Rayleigh resolution", which can be written as Rayleigh Resolution = × Null-to-Null Width of Point Spread Function (24) where the point spread function (PSF) is response of imaging setup to a point scatterer and is a subjective constant which is typically 0.5 or 1.0. Since the MIMO array has a rectangular aperture, it is selected as 1.0 in this study. Note that this definition is valid for high SNR cases. On the other hand, for very low SNR cases, Cramér-Rao bound must be considered, however, it is not included in this study. We considered three target scenarios and numerically simulated them in MATLAB using a workstation with Intel(R) Xeon(R) CPU E5-2650 2.60 GHz processor and 128 GB RAM. For comparing the performance of the sparsityinducing algorithms (ADMM and C-SALSA-2), conventional CGLS is also implemented.
In our simulations, the signal power is defined as where the term in the parenthesis is the average return power from the image scene. The total power is given by times this quantity, because we are using transmitters, receivers, and frequency steps. Note that this definition strongly depends on the shape and reflectivity distribution of the target, since averaging is done over non-empty voxels.
Quality of the results are estimated with respect to peak signal-to-noise ratio (pSNR) of reconstructed reflectivity images, defined by pSNR = 10log 10 where, MAX is maximum voxel value of reconstructed reflectivity image and MSE is mean squared error which can be simply defined as the mean of the square of the difference between the noise-free × reference image and ̂, which is the image reconstructed from the noisy measurement as Note that all simulations given below are carried out with 30 dB input SNR.

A. POINT SCATTERERS
In the first scenario, the scene contains an array of 9 point scatterers located on a grid of 0.025 m on each side on the = 0.55 m plane. The reflectivity level of the point scatterer at the center is set to 0.5, whereas all other scatterers have a reflectivity of 1. Fig. 4(a) shows actual reflectivities of the scatterers for this scenario. The images are reconstructed by solving (14) using the developed efficient ADMM approach with FMM. Here ‖ ‖ 1 is used as regularization function. Note that this problem can also be solved by the developed C-SALSA-2 approach with the same regularization function, however, we preferred using ADMM in order to demonstrate applicability of FMM in different sparse reconstruction algorithms. The images are formed on = 0.55 m plane with a resolution of 2 mm in both directions, resulting in a total of 2601 pixels. The reconstructions are normalized by the largest value and corresponding images obtained by CGLS and ADMM are given in Fig. 4(b) and Fig. 4(c), respectively.
The results show that the theoretical cross range "Rayleigh" resolutions can be achieved by solving the problem with ADMM. On the other hand, CGLS cannot provide a focused image under the same conditions, which emphasize the The quality of the reconstructed images is quantitatively evaluated by calculating their pSNRs, which are 12.36 dB and 42.17 dB for CGLS and ADMM, respectively. In computational image formation, pSNR lower than 20 dB is accepted to be low-quality while pSNR higher than 40 dB is considered to be almost excellent quality [39]. Hence, the calculated pSNRs and difference between them indicate that ADMM performs more effectively for the reconstruction of a sparse scene compared to the conventional CGLS.

B. EXTENDED TARGETS
For the reconstruction of extended targets, two different cases are considered: (i) a 2-D hexagonal target and (ii) 3-D concentric cylinders. The developed efficient C-SALSA-2 algorithm with FMM is used to reconstruct the images that are given in this subsection. TV of the magnitude is used as the regularization function.
(i) Hexagonal Target The actual reflectivity of the hexagonal target is illustrated in Fig. 5(a). The target is a planar structure of size 0.4 m × 0.4 m located at = 0.55 m plane, and the image is constructed using 4 mm by 4 mm pixels, which makes 10201 pixels in total. Reflectivity level distribution of the target is nonuniform, i.e., reflectivity of the outer part is 1, whereas it is 0.6 for the inner parts except for the U-shaped section, which has a reflectivity of 0.4. The measurements are simulated using the discrete model given in (7) and the reflectivity images are reconstructed by CGLS and C-SALSA-2. Fig. 5(b) and Fig. 5(c) respectively shows the normalized reconstructions obtained by CGLS and C-SALSA-2, where the color scale shows the magnitude of the reflectivity distribution of the imaging scene along x-and z-directions. The results demonstrate that the characteristics of the hexagonal target such as the edges where the reflectivity rapidly changes and the U-shaped section is blurry, when the image is reconstructed using CGLS. On the other hand, solving with C-SALSA-2 provides a well-focused image, containing all details of the target. Furthermore, the background that does not contain any scatterer is clearer in the reconstructed image obtained by C-SALSA-2. The pSNR values for CGLS and C-SALSA-2 are calculated as 22.95 dB and 36.23 dB, respectively, which mathematically support our observations. For one matrix-vector product, efficient computation time is measured as ~54.9 seconds, while memory usage is ~1.07 GBytes.

(ii) Concentric Cylinders
In the last scenario, we considered a 3-D object that is formed by three concentric cylinders, lying along the range direction SALSA-2 are depicted in the middle and right columns of Fig.  6, respectively. The slices in Fig. 6(b), Fig. 6(e), and Fig. 6(h) show that it is not possible to attain well-focused reconstruction by CGLS e.g., the concentric cylinders are blurry and reflectivity changes are not exactly how they should be. Conversely, as can be clearly seen in Fig. 6(c), Fig. 6(f), and Fig. 6(i), the reconstructions by C-SALSA-2 results visibly more accurate images in terms of reflectivity distribution and contrast transitions between the cylinders. Besides, the calculated pSNRs for CGLS and C-SALSA-2 reconstructions are 22.87 dB and 36.50 dB, respectively, which support our observations. Besides, efficient computation time and memory requirement per matrix-vector product respectively take ~162.7 seconds and 3.1 GBytes.
We further investigate the performance of FMM-based sparse reconstruction method for different input SNR values. To do this, we repeated the simulations of the concentric cylinders with C-SALSA-2 and CGLS, each time decreasing the input SNR gradually from 50 dB to 0 dB and measuring the pSNR values of the reconstructed reflectivity images. Fig. 7 plots the pSNRs as a function of input SNR value for both C-SALSA-2 and CGLS. As can be seen, C-SALSA-2 can still provide pSNR above 20 dB although the input SNR is about 5 dB. On the other hand, conventional CGLS cannot handle noise presence as well as C-SALSA-2 due to the lack of regularization term and its pSNR continuously decreases after the input SNR gets lower than 25 dB. Fig. 8 shows the reconstructions obtained by CGLS and C-SALSA-2 for the hexagonal target when different input SNR values are applied.
Finally, we analyzed the computational performance of FMM. As stated earlier, the sparsity-based reconstruction algorithms require solving multiple forward problems at each iteration. Direct computation of the forward problems with unknowns leads to immense computational burden for large imaging problems due to the per-iteration complexity of ( 2 ). On the other hand, FMM calculates the matrix-vector product in a group-by-group manner (see Section III for our grouping method), hence, reduces the computational  complexity to ( 3/2 ) for both processing time and memory at each iteration.
The efficiency of FMM is also demonstrated by simulating the same target scenario with different number of antennas. Fig. 9 plots processing time (in seconds) and memory requirement (in Mbytes) as a function of total number of antennas (i.e., = ). As proposed, FMM outperforms the direct matrix-vector product above a specific number of antennas, which we call crossover point, in terms of processing time and memory. Note that this crossover point depends on several things that must be considered, i.e., type of application, how the forward (sensing) operator is implemented, truncation order of the FMM, etc.

VII. CONCLUSION
In this paper, we solve near-field imaging problems and seek sparse solution for them. A MIMO radar configuration is considered for the imaging and corresponding forward problem is implemented using Born approximation. Then, the problem is discretized and turned into a linear system of equations, the structure of which is exploited to apply a two stage FMM, thus, large-scale problems can be solved in an efficient way. The problems are solved by two sparsityinducing algorithms (ADMM and C-SALSA-2) in conjunction with appropriate regularization functions (ℓ1norm and ) according to the properties of the imaged scene and reflectivity characteristics of the targets. We validated our approach by a series of numerical tests conducted in comparison with a classical solver. Results demonstrate that well-focused images can be constructed for different type of near-field imaging scenarios, ranging from point scatterers to 3-D continuous targets, despite low input SNR. When FMM is applied to the solution of the problem, quality of the reconstructions can still be preserved. Furthermore, as the dimensions of the problem increases, which is number of antennas in our case, FMM makes the solution very efficient in terms of processing time and memory requirement.