Design of Large Effective Apertures for Millimeter Wave Systems using a Sparse Array of Subarrays

TL;DR

This paper presents a method for designing large, effective millimeter wave apertures using a sparse array of subarrays optimized for beam width and sidelobe suppression, with applications in source bearing estimation.

Contribution

It introduces a multi-objective optimization framework for subarray placement in sparse arrays, balancing beam width and sidelobe levels for millimeter wave systems.

Findings

Optimized sparse array designs improve bearing estimation accuracy.

Performance approaches estimation-theoretic bounds.

Numerical examples demonstrate effectiveness with 16-element subarrays at 60 GHz.

Abstract

We investigate synthesis of a large effective aperture using a sparse array of subarrays. We employ a multi-objective optimization framework for placement of subarrays within a prescribed area dictated by form factor constraints, trading off the smaller beam width obtained by spacing out the subarrays against the grating and side lobes created by sparse placement. We assess the performance of our designs for the fundamental problem of bearing estimation for one or more sources, comparing performance against estimation-theoretic bounds. Our tiled architecture is motivated by recent progress in low-cost hardware realizations of moderately sized antenna arrays (which play the role of subarrays) in the millimeter wave band, and our numerical examples are based on 16-element (4x4) subarrays in the 60 GHz unlicensed band.

Tables

Table 1. TABLE I: Sample Array Configurations

Shape	$\alpha$	$\beta$	$\gamma$	MSLL	BW	ecc
A1	1	0.1	0.1	-8 dB	${5.7}^{\circ }$	0
A2	1	1	1	-10.2 dB	$8^{\circ }$	0
Compact	-	-	-	-12.8 dB	${12.1}^{\circ }$	0.7
Naive	-	-	-	-7.6 dB	${7.9}^{\circ }$	0.0

Equations

$$\mathtt{BW}(\mathbf{C},\mathbf{w}),\mathtt{MSLL}(\mathbf{C},\mathbf{w}),\mathtt{ecc}(\mathbf{C},\mathbf{w})$$

$$G_D(\mathbf{C},\mathbf{w})$$

$$\text{textsc}AoS(\mathbf{C})$$

$$u=\sin (\theta )\cos (\varphi ),v=\sin (\theta )\sin (\varphi )$$

$$R(u,v)=\frac{1}{N^2}{\left|\sum _{i=1}^N{e}^{jk\left(u{d}_i^x+v{d}_i^y\right)}\right|}^2$$

$$R_{(u_0,v_0)}(u,v)=R(u-u_0,v-v_0)$$

$$R_{\max }(u_0,v_0)=$$

$$=$$

$$R_{\mathrm{sidelobe}}=$$

$$\text{s.t.}(u^{\ast },v^{\ast })\ne (u_0,v_0),$$

$$R_{u_0,v_0}(u^{\ast },v^{\ast })\ge R_{u_0,v_0}(D_{ϵ}(u^{\ast },v^{\ast }))$$

$$R_{avg}=\frac{2}{4\pi }{\int }_{-1}^1{\int }_{-\sqrt{1-v^2}}^{\sqrt{1-v^2}}\frac{R_{u_0,v_0}(u,v)}{\sqrt{1-u^2-v^2}}dudv$$

$$o(\mathbf{C})=\frac{o^{\mathrm{raw}}(\mathbf{C})-\underset{\forall \mathbf{C}}{\min }\left\{o^{\mathrm{raw}}(\mathbf{C})\right\}}{\underset{\forall \mathbf{C}}{\max }\left\{o^{\mathrm{raw}}(\mathbf{C})\right\}-\underset{\forall \mathbf{C}}{\min }\left\{o^{\mathrm{raw}}(\mathbf{C})\right\}}$$

$${\mathbf{C}}^{\ast },{\mathbf{w}}^{\ast }=$$

$$\text{textsc}AoS(\mathbf{C})$$

$$f(\mathbf{C},\mathbf{w})=\alpha \mathtt{BW}(\mathbf{C},\mathbf{w})+\beta \mathtt{MSLL}(\mathbf{C},\mathbf{w})+\gamma \mathtt{ecc}(\mathbf{C},\mathbf{w})$$

$$R_{\rho }(\tilde{u},\tilde{v})$$

$$\rho$$

$${\mathbf{C}}^{\ast }=\text{arg}\underset{\mathbf{C}\in \mathcal{C}}{\min }f(\mathbf{C},\text{mathbbm}{1}_N)$$

$$\mathbf{x}=\sum _{j=1}^K{\alpha }_j\mathbf{s}({\mathbf{u}}_j)+\mathbf{z}$$

$$p(\mathbf{x}|\Theta ,\mathbf{\alpha })=\prod _{{\mathbf{u}}_j\in \Theta }\frac{1}{{\pi }^N{\sigma }^2}\exp \left(-\frac{{\Vert \mathbf{x}-{\alpha }_j\mathbf{s}({\mathbf{u}}_j)\Vert }^2}{{\sigma }^2}\right)$$

$${\mathbf{R}}_{ϵ}(\hat{\Theta })=\mathbb{E}\left[\sum _{i=1}^K(\mathbf{u}-{\hat{\mathbf{u}}}_{\mathbf{i}})(\mathbf{u}-{\hat{\mathbf{u}}}_{\mathbf{i}}{)}^T\right]$$

$${\hat{\mathbf{u}}}_{ML}=arg\underset{\mathbf{u}}{\max }{\left|\mathbf{s}(\mathbf{u}{)}^{\mathbf{H}}\mathbf{x}\right|}^2$$

$$CRB({\mathbf{R}}_{\mathbf{ϵ}})=({\mathbf{J}}_{\mathbf{F}}+{\mathbf{J}}_{\mathbf{P}}{)}^{-1}$$

$$({\mathbf{J}}_{\mathbf{F}}{)}_{ij}=-{\mathbb{E}}_{\mathbf{x},\mathbf{u}}\left[\frac{{\partial }^{2}l(\mathbf{x}|\mathbf{u})}{\partial u_i\partial u_j}\right],({\mathbf{J}}_{\mathbf{P}}{)}_{ij}=-{\mathbb{E}}_{\mathbf{u}}\left[\frac{{\partial }^{2}l(\mathbf{u})}{\partial u_i\partial u_j}\right]$$

$${\mathbb{E}}_{\mathbf{x},\mathbf{u}}\left[\frac{\partial l(\mathbf{x}|\mathbf{u})}{\partial \mathbf{u}}\right]=0$$

$${\mathbb{E}}_{\mathbf{x},\mathbf{u}}\left[jk\sum _{i=1}^N{\mathbf{d}}_i\left(x_i{e}^{jk{\mathbf{u}}^T{\mathbf{d}}_i}-x_i^{\ast }{e}^{-jk{\mathbf{u}}^T{\mathbf{d}}_i}\right)\right]=0$$

$$jk\left(\alpha \sum _{i=1}^N{\mathbf{d}}_i-{\alpha }^{\ast }\sum _{i=1}^N{\mathbf{d}}_i\right)=jk(\alpha -{\alpha }^{\ast })\sum _{i=1}^N{\mathbf{d}}_i$$

$$J_F$$

$$=2k^2\gamma {\mathbf{D}}^T\mathbf{D}$$

$${\mathbf{a}}^{\mathbf{T}}{\mathbf{R}}_{\mathbf{ϵ}}\mathbf{a}\ge {\int }_0^{\infty }\mathcal{V}\left\{\underset{\mathbf{\delta }:{\mathbf{a}}^{\mathbf{T}}\mathbf{\delta }=h}{\max }\int A(\mathbf{u},\mathbf{\delta })P_e(\mathbf{u},\mathbf{\delta })d\mathbf{u}\right\}hdh$$

$$H_0$$

$$H_1$$

$$\mathrm{Decide}(\mathbf{u})=\{\begin{array}{ll}\mathbf{u}& \text{if }{\rho }_1>{\rho }_2\\ \mathbf{u}+\mathbf{\delta }& \text{if }{\rho }_1<{\rho }_2\end{array}$$

$$P_{nc}(\mathbf{u},\mathbf{\delta })$$

$$=Pr\left({\rho }_1<{\rho }_2|\mathbf{u}\right)$$

$$=Q_1(a,b)-\frac{1}{2}{e}^{-\frac{a^2+b^2}{2}}{I}_0(ab)$$

Full text

I Introduction

Many sensing and situational awareness applications (e.g., radar imaging for vehicles and drones) require highly directional, electronically steerable beams. Reducing beam width requires expansion of antenna aperture. This is typically accomplished by filling the aperture with antenna elements spaced at half the carrier wavelength or less, in order to avoid grating lobes. However, this approach does not scale well with aperture size since the cost, power consumption and design complexity increases with number of antenna elements.

In this paper, we investigate the problem of synthesizing narrow beams using a tiled architecture, with a sparse set of subarrays spread over a large physical aperture. Each subarray is a relatively compact antenna array with a moderate number of elements. The total number of antenna elements, summing over subarrays, is much smaller than that for a classical design spanning the entire physical aperture. The sidelobes and grating lobes resulting from such spatial undersampling must therefore be controlled in order for our proposed “array of subarrays” to be useful. Our goal here is to determine the placement of a given number of subarrays over a physical aperture in order to optimize multiple beam attributes, including beam width, maximum sidelobe level, and directivity.

While our framework is general, the design of millimeter wave arrays is of particular interest to us, because the small carrier wavelength enables synthesis of narrow beams using relatively compact apertures. As a running example throughout this paper, we consider the design of a 60 GHz array of subarrays created by placing 8 subarrays over an aperture size of $10$ cm by $10$ cm ($20\lambda\times 20\lambda$ for $\lambda=0.5$ cm), where each subarray has $4\times 4$ elements arranged in uniform rectangular grid with $0.5\lambda$ horizontal spacing and $0.6\lambda$ vertical spacing. The subarrays need to be aligned along their axes, assuming that all elements have unidirectional linear polarization. Each subarray tile occupies extra physical area on the plane, which must also be accounted for in the placement procedure. These design restrictions are consistent with existing prototype subarrays.

I-A Contributions

Our contributions are summarized as follows:

$\bullet$ We formulate the problem of subarray placement as multi-objective optimization of key performance measures such as beam width, maximum sidelobe level, eccentricity and directivity,

[TABLE]

where $\bm{C}$ is $N_{s}\times 2$ Subarray center position matrix, $\bm{w}$ is $N\times 1$ beamsteering weight vector and $\textsc{AoS}(\bm{C})$ are physical constraints to avoid overlapping subarrays. Note that orientation of subarrays is not an optimization variable in this architecture, since the polarization of the elements has to be aligned for beamforming. The problem (1) is non-convex and combinatorially explosive. We use geometric heuristics to eliminate similar configurations in the first stage of our algorithm, and then employ a second stage of refinement using small perturbations around the first stage solution.

$\bullet$ We evaluate our designs using estimation-theoretic benchmarks for two-dimensional (2D) direction of arrival (DoA) estimation. At low signal-to-noise ratio (SNR), large sidelobes can lead to large errors in the DoA estimate. At high SNR, on the other hand, the DoA estimation error is governed by beam width. We derive a Ziv-Zakai bound (ZZB), which captures the effect of both large and small estimation errors, for DoA estimation for specular paths. The ZZB exhibits a distinct transition in its behavior from low to high SNR, tending at high SNR to the Cramer-Rao bound (CRB), which captures the effect of small errors around the true parameter value. Thus, we use the ZZB transition SNR as a measure of efficacy of sidelobe reduction, and the CRB as a measure of efficacy of beam width reduction.

$\bullet$ The beam attributes of arrays designed using multi-objective optimization depend on its parameterization. We design two sparse arrays, $A1$ with primary emphasis on reducing beamwidth and $A2$ based on joint optimization of beamwidth and maximum sidelobe level. These designs are compared against two benchmark arrays. The first is termed a “compact array,” with subarrays packed closely together: this is expected to have worse beam width but smaller sidelobes than our sparse designs. The second is termed a “naive array,” obtained by placing subarrays in diamond pattern to obtain beamwidth equivalent to that of sparse array A2. Some illustrative numerical results for our running example are as follows. The sparse design A1 is $11$ dB better than the compact array in terms of CRB, while degrading less than $1$ dB in terms of ZZB threshold. The sparse design A2 is $4$ dB better in terms of CRB than the “naive array,” while also having a better ZZB threshold.

$\bullet$ We investigate DoA estimation performance numerically using a state of the art super-resolution algorithm. The impact of the higher sidelobes due to sparse placement, and hence that of our optimization procedure, is more evident when estimating DoA in the presence of multiple interfering targets. We show that, depending on the strength of interferers, our optimized arrays achieve better estimation accuracy than the “compact” and “naive” benchmark arrays at moderate to high SNR, due to a combination of sharper beamwidth and lower sidelobes. We also show that the efficacy of DoA estimation using our sparse designs, and the associated benchmarks, is maintained when we employ compressive measurements.

I-B Related Work

There is a rich body of work on sparsifying linear arrays, including minimum redundancy arrays [1], genetic optimization[2], joint Cramér Rao Bound and sidelobe level optimization [3], and simulated annealing [4]. Most popular design strategies try to find an element pattern which minimizes beamwidth, along with some notion of DoA ambiguities such as sidelobe level or probability of DoA outlier. Recent approaches like Nested 2D arrays [5] and H-arrays [6] utilize the idea of “difference co-array” to reduce the number of redundant spacings and maximize the randomness of element positions, so that the number of spatial frequencies being sampled by the array is maximized. Most of these techniques, however, assume that antenna elements can be placed freely. Hence, they do not apply in our setting, where element placement within subarrays is constrained.

The prior work most similar to our is [7], which investigates design of linear arrays with two and three subarrays. However, the focus there is on performance criteria for comparing a number of sensible designs in a far smaller design space, rather than searching over a large space of possibilities as we do here.

Our performance evaluation requires implementation of DoA estimation algorithms. Classical subspace-based algorithms such as MUSIC [8] and ESPRIT [9], as well as their extensions to arrays of subarrays such as [10, 11, 12, 13], rely on regular array geometries for efficient computation. Recently developed super-resolution algorithms such as Basis Pursuit Denoising (BPDN) [14] and Newtonized Orthogonal Matching Pursuit (NOMP) [15] are both more general and have better performance. It is worth noting that [16] shows that, for large arrays, BPDN and other sparse estimation techniques with compressive measurements outperform subspace-based methods. In this paper, we employ NOMP in our numerical experiments, since we have found it to provide better performance than BPDN at lower complexity. We also show that the performance trends are unchanged under compressive measurements, consistent with recent general theory [17].

I-C Outline

We first describe the beam attributes to be optimized while designing the sparse arrays and discuss constraints for the optimization in Section II. The geometric heuristics and design approach is described in detail in Section III. We then provide a brief review of estimation bounds for 2D bearing estimation and discuss their utility for analysing the sparse arrays in Section IV. Finally, we compare the example array configurations against benchmarks and evaluate the performance of super-resolution algorithms in Section V.

II Sparse Subarray Design

We formulate the array design problem in terms of jointly optimizing multiple beam parameters that are expected to affect DoA estimation performance.

II-A Beam Pattern Basics

We use the directional cosines

[TABLE]

to represent the DoA of target. The elevation $\theta$ and azimuth $\phi$ angles are measured from the broadside direction (perpendicular to the baseline array plane). The 2D beampattern $R(u,v)$ in direction $(u,v)$ when the beam is steered towards the broadside is given by

[TABLE]

where $N$ is the number of array elements; $[d_{i}^{x},d_{i}^{y}]^{T}\triangleq\bm{d}_{i}$ are the 2D co-ordinates of arrays elements, and $k=\frac{2\pi}{\lambda}$ is the wavenumber. We assume isotropic antenna elements with ideal steering weights and far-field sources with normalized response. The term “subarray” refers to the subset of elements with uniform half-wavelength spacing, while “super-array” refers to the placement of these subarrays, which is described by the subarray centers. Since the elements in a subarray are fixed, the array element locations, $\bm{D}$ can be expressed in terms of the subarray centers, $\bm{C}$, as $\bm{D}=\bm{C}\otimes\mathbbm{1}_{N_{e}}+\bm{D}_{e}\otimes\mathbbm{1}_{N_{s}}$, where $\bm{D}_{e}$ is the fixed $2\times N_{e}$ matrix containing the subarray element coordinates with respect to its center, $N_{s}$ is the number of subarrays, $N_{e}$ is number of elements in individual subarrays, $\mathbbm{1}_{n}$ is an $n\times 1$ column vector of ones, and $\otimes$ denotes the Kronecker product.

When beamforming in a general direction $(u_{0},v_{0})$ (broadside corresponds to $(u_{0},v_{0})=(0,0)$), the beam pattern is simply given by

[TABLE]

For Direction of Arrival (DoA) estimation, the ideal beam should have small beamwidth with minimal sidelobes and high directivity. We consider the following beam attributes, some of which depend on the steering direction $(u_{0},v_{0})$, as key performance metrics to be optimized:

• •

2D beamwidth (BW): Although the mainbeam of non-uniform array has non-trivial shape in 2D, we approximate it as an ellipse to define beamwidth. We evaluate the 2D beamwidth in terms of the 3-dB beamwidths along the major and minor axes of this ellipse, denoted by $\mathrm{BW_{max}}$ and $\mathrm{BW_{min}}$, respectively. The mean squared error of DoA estimation depends on the sum of these beamwidths (see Appendix A), hence we define beamwidth as $\mathrm{BW^{DoA}=\sqrt{BW_{max}^{2}+BW_{min}^{2}}}$.

• •

Maximum sidelobe level (MSLL): is the relative level of the strongest sidelobe in the beampattern with respect to mainlobe. $\mathtt{MSLL}=10\log\left(\mathrm{R_{max}/R_{sidelobe}}\right)$. $R_{\mathrm{max}},R_{\mathrm{sidelobe}}$ are the largest and second largest magnitude local maxima’s of beampattern $R(u,v)$ given by,

[TABLE]

where $D_{\epsilon}(u^{*},v^{*})=\left\{(u,v):|u-u^{*}|<\epsilon,|v-v^{*}|<\epsilon\right\}$ denotes $\epsilon$-neighborhood. Note that $R_{\mathrm{max}}$ does not depend on steering direction $(u_{0},v_{0})$, but $R_{\mathrm{sidelobe}}$ might.

• •

Directivity ($\mathtt{G_{D}}$): The directivity is the ratio of mainlobe power to average power, $\mathtt{G_{D}}=10\log\left(\mathrm{\frac{R_{max}}{R_{avg}}}\right)$ The average power does not have closed form expression for general planar arrays and is evaluated in $(u,v)$ domain by the integral [18],

[TABLE]

• •

Eccentricity (ecc): is a measure of the asymmetry of the mainbeam. We add this additional parameter to suppress the trivial linear placement solution, $\mathtt{ecc}=\sqrt{1-\left(\mathrm{BW_{min}/BW_{max}}\right)^{2}}$

For non-uniform planar arrays, none of these beam parameters have a closed form expression, hence they must be computed numerically. In our simulations, we compute these beam attributes using a beampattern over a $512\times 512$ grid in UV space as shown in Figure 2.

II-B Problem Formulation

In order to develop geometric heuristics for optimization, we first analyze the effect of increasing aperture width, (keeping the number of antenna elements fixed, for linear and planar arrays, with uniform and subarray-based architectures as shown in the rightmost section of Figure 3. In the plots of beam attributes in Figure 3, the dashed line represents the aperture width when inter-element spacing equals a half-wavelength, when the uniform and subarray-based configurations match.

• •

The MSLL for subarrayed array increases much faster than uniform configuration due to grating lobe appearing close to mainbeam. This attribute is sensitive to element distribution and behaves unpredictably for non-uniform arrays.

• •

The 3dB Beamwidth ($\mathtt{BW}^{max}$ for planar array) for both array types reduces congruently reaffirming that it is inversely proportional to aperture width independent of the distribution of elements.

• •

Directivity increases as we increase the inter-element spacing, but only upto certain limit and then becomes constant [19]. This generalizes well to planar arrays as shown in the Figure 3. As one can see the directivity for subarrayed configurations remains approximately constant with increasing aperture width beyond standard spacing. Hence we do not include this metric in our cost function.

The objectives that we wish to trade off against each other do not have the same units: for example, MSLL is measured in dB from max sidelobe level, whereas BW is measured in $\mathrm{deg}$. We therefore normalize each raw objective value, $o^{\mathrm{raw}}$ by its range as follows:

[TABLE]

The range of each objective is computed numerically while constructing the dictionary of all configurations, $\bm{C}\in\mathscr{C}$, described later in III-A.

The constrained multi-objective optimization can now be formulated as follows:

[TABLE]

where

[TABLE]

is the weighted cost function in terms of the normalized objective functions and [$\alpha,\beta,\gamma$] are weights that can be used to sweep through the optimal surface for this optimization. We show some example arrays obtained for different choices of weights in Table I.

This is a combinatorial optimization problem which is non-convex, with non-differentiable geometric constraints. The objective function can be evaluated exactly given the subarray centers, but exploring the entire solution space is computationally infeasible. Furthermore, beam characteristics in general depend on the steering weights, which in turn depend on the direction $(u_{0},v_{0})$ in which we are steering. We therefore employ two key simplifications:

• •

We remove the dependence of the cost function on beamforming direction, and hence on steering weights, by computing the objectives based on an expanded beam pattern, as discussed in Section II-C.

• •

We employ geometric heuristics to cut down the solution space to a reasonable size, as described in Section III.

II-C Invariance to Beamforming Direction

The cost function in (6) is evaluated using the beampattern $R(u-u_{0},v-v_{0})$, which depends on the beamsteering direction ($u_{0},v_{0}$). It would be prohibitively expensive to evaluate the beam attributes over all such beampatterns for finding the optimal ($\bm{C,w}$). However, for arrays with isotropic elements, we can define an Expanded Beam pattern (EBP) which subsumes beampatterns of all steering direction in a Region of Interest (ROI) [20]. Suppose that our maximum steering angle in the ROI is $\theta_{max}$. From (2), we see that $(u_{0},v_{0})$ lies within a circle of radius $\sin\theta_{max}$. On the other hand, sidelobes can appear at any $(u,v)$ within a circle of radius $1$. It is easy to see, therefore, that $(u-u_{0},v-v_{0})$ is guaranteed to lie within a circle of radius $1+\sin\theta_{max}$. We can therefore compute beam attributes using the following EBP:

[TABLE]

Figure 4 shows the EBP, $R_{1.5}(\tilde{u},\tilde{v})$ for ROI with $\theta_{max}=30^{\circ}$, and the beampattern for the steering angle ($(u_{0},v_{0})=(0.3,0.4)$). The shape of the main beam is preserved under the transformation (7), hence beam width and eccentricity can be directly evaluated from (7). The MSLL evaluated from EBP is a worst-case value, corresponding to an argument $\rho(\tilde{u},\tilde{v})$ which, in principle, might not correspond to a feasible value of $(u-u_{0},v-v_{0})$ in (4). However, the maximum sidelobe always lies within the main lobe of the subarray beam pattern, so that physically implausible values of $\rho(\tilde{u},\tilde{v})$ do not correspond to large local maxima of the EBP.

With the introduction of the EBP, we can, without loss of generality, assume that the main beam is being steered towards broadside, setting $\bm{w}=\mathbbm{1}_{N}$. Our problem now reduces to finding the optimal configuration $C^{*}$.

III Placement Optimization

In order to optimize the placement, we need to evaluate the cost function over all array configurations. The number of configurations depends on the allowed form factor, and the size of the subarray module, and an exhaustive search over all configurations is computationally infeasible: for example, the number of configurations for a discrete grid of size $20\times 20$ is of the order $10^{20}$. We therefore propose a two-stage approach, first performing a combinatorial search on a reduced search space, and then obtaining the final solution by searching over perturbations around the solution from the first stage.

III-A Combinatorial search

We reduce the solution space by removing geometrically “similar” arrays. We employ the covariance of the element positions, and pairwise element separations, as measures of similarity. The choice of covariance of element positions $\bm{\Sigma}(\bm{C})$ as similarity metric is motivated by its inverse proportionality to Cramér Rao Bound on accuracy of DoA estimation (see Section IV-A1). However, array configurations with similar array covariance but diverse beam attributes also exist: Figure 5 shows an example of two array configurations with different shapes but the same covariance. The arrays have similar beamwidth but their MSLL levels are different. We observe that the variance of pairwise element distances $\psi(C)$ is different for these arrays, and use it as an indicator for these large scale deviations.

This allows us to reduce the dimensionality of the solution space from $2\times N_{sub}$ down to $3$ array shape parameters: the eigenvalues ($\lambda_{1},\lambda_{2}$) of $\bm{\Sigma}(\bm{C})$ and $\psi(C)$.

III-A1 Subarray Placement Algorithm

We construct a prefix tree dictionary to find feasible solutions using a breadth first search based enumeration technique. The element position covariance for an array of subarray can be uniquely represented by the covariance of its subarray centers, $\bm{\Sigma}_{D}=\bm{\Sigma}_{C}+\bm{\Sigma}_{D_{e}}$. Hence the super-array center covariance can be used instead of that of the full array in the Dictionary search algorithm. Each node in the tree stores a subarray center position, and the path from root to a node at the $n^{th}$ layer of prefix tree represents a unique configuration of $n$ subarrays. The subarray centers are constrained to lie on a fixed set of discrete grid points $G$.

The algorithm is described in Algorithm 1. We briefly discuss the key steps below.

• •

Initialize: In order to allow for sufficient exploration, we employ multiple random initializations $C_{i}^{1}$ of the root node being placed on $N_{\mathrm{init}}$ different locations on the grid. (For example, circular configurations cannot be obtained if the root subarray is fixed at the center.)

• •

List: Define the operator $\mathscr{T}_{G}:G^{|C_{i}^{n}|}\rightarrow G^{|V_{i}|}$ which maps the set of subarrays centers $C_{i}^{n}\in\mathscr{C}^{n}$ to the set of $|V_{i}|$ vacant gridpoints in $G$ available for placement of next subarray which are not blocked by the subarrays already placed at $C_{i}^{n}$. This operator also accounts for additional surface area occupied by the subarray module apart from the physical antenna elements (see Appendix C for details).

• •

Append: The $(n+1)^{th}$ subarray configuration is constructed from the vacancies, $C_{i}^{n}\times V_{i}$ where $\times$ denotes cartesian product of sets. A temporary dictionary $\hat{\mathscr{C}}$ is formed by inserting $|V_{i}|=\kappa$ child nodes for each node in the $n^{th}$ layer.

• •

Prune: Nodes corresponding to “similar” configurations are deleted based on the array shape parameters

Find eigenvalues ($\lambda_{1},\lambda_{2}$) of the subarray center covariance matrix, $\Sigma(C)$ and variance of array separations, $\psi(C)=\mathbb{E}\left[(l_{ij}-\mathbb{E}[l_{ij}])^{2}\right]$, where, $l_{ij}$ denotes the distances betweeen $i^{th}$ and $j^{th}$ elements. 2. 2.

Enumerate unique configurations by binning the ($\lambda_{1},\lambda_{2},\psi$) triplets over a 3-D grid with resolution $\tau$ and randomly picking one configuration from each bin (see Appendix D for criteria to choose $\tau$).

The procedure is repeated until dictionary atoms reach the desired number of subarrays i.e $n=N_{S}$. All arrays in the dictionary $\mathscr{C}$ obtained from this algorithm satisfy the AoS$(\bm{C})$ constraint by construction. This simplifies the constrained Multi-objective optimization in Eq. (5) to

[TABLE]

III-B Iterative Placement Refinement

In the second stage, we try to improve the cost function (8) by applying small local perturbations (within a bin of the grid $G$) to the subarray positions obtained from the combinatorial search in the first stage, as described in Algorithm 2.

Figure 6 shows a sample of how costs are minimized using sequential refinement. After running few iterations, the final array has $0.8^{\circ}$ lower beamwidth, while keeping other beam attributes relatively unchanged.

IV Estimation-Theoretic Benchmarks

We now seek to evaluate the efficacy of our sparse designs for the canonical application of 2D DoA estimation. We compare different array designs in terms of estimation-theoretic bounds as well as simulated performance using a super-resolution algorithm. For clarity in exposition, from here onwards we overload $\bm{u}\triangleq[u,v]$ to denote the DoA.

IV-A Signal Model

We model the received signal from $K$ sources in the scene with distinct DoAs $\Theta=[\bm{u_{1},u_{2},\cdots,u_{k}}]$ as

[TABLE]

where $\bm{s(u_{j})}=\left[e^{jk\bm{u_{j}^{T}d_{1}}},\ldots,e^{jk\bm{u_{j}^{T}d_{N}}}\right]^{T}$ is the array response, $\bm{z}=[z_{1},\ldots,z_{N}]^{T}$ is complex white noise such that $\mathbb{E}(zz^{H})=\sigma^{2}I_{N}$, and $\{\alpha\}_{j=1}^{K}$ are complex gains which are unknown deterministic constants. The joint probability density of received signal conditioned on $(\Theta,\{\alpha\}_{j=1}^{K})$ is given by,

[TABLE]

For any DoA estimator $\hat{\Theta}$, the covariance of estimation error is defined as,

[TABLE]

$\bm{R}_{\epsilon}$ can be geometrically interpreted by its trace $\sqrt{tr(R_{\epsilon})}$ which represents the expected overall Root mean square error (RMSE) in DoA estimation (see Appendix A). We use this measure to compare the performance of array designs in Section V. For single source case $(K=1)$, the joint maximum likelihood estimator of $\bm{u}$ and $\alpha$ yields a noncoherent estimator for $\bm{u}$ as follows:

[TABLE]

For this case, we derive the Cramer Rao (CRB) and Ziv-Zakai (ZZB) bounds on $\bm{R}_{\epsilon}$ to assess the best possible estimation accuracy of different designs. Although derived for single source case, we use these bounds for multiple source case as well to compare DoA estimation performance.

IV-A1 Cramér Rao Bound

The Bayesian Cramér Rao Bound for this signal model is given by [20]:

[TABLE]

where $\bm{J_{F},J_{P}}$ denote the Fisher Information Matrix (FIM) contributions from the observation and the prior distribution of DoA respectively.

[TABLE]

where $l(\bm{x|u}),l(\bm{u})$ are the conditional log likelihood and prior log likelihoods respectively. In addition the following regularity condition needs to be satisfied,

[TABLE]

In order to always satisfy this condition, we enforce the array element positions to be centered i.e. $\sum_{i=1}^{N}{\bm{d}_{i}}=\bm{0}$

For single source, the FIM is given by,

[TABLE]

which depends only on the element positions, $D$ and Signal to Noise ratio (SNR) ($\gamma=|\alpha|^{2}/\sigma^{2}$). Assuming the DoA prior to be uniformly distributed in the ROI ($\theta\leq 30^{\circ}$), the prior FIM simplifies to $\bm{J_{P}}=1.343\bm{I}_{2}$.

IV-A2 Ziv-Zakai Bound

The CRB is a local bound, which accounts for estimation performance dependent on mainbeam, hence it is only useful at high SNR. In order to better characterize the estimation performance of Sparse arrays at low SNR, we calculate the Ziv-Zakai Bound (ZZB) which incorporates the effect of sidelobes and predicts the threshold behavior. For any directional vector $\bm{a}=[\cos\xi,\sin\xi]^{T}$, the ZZB is given by [21]

[TABLE]

where, $A(\bm{u,\delta})=\min{\left\{p(\bm{u}),p(\bm{u+\delta})\right\}}$, $\mathcal{V}(.)$ is the valley filling function and $P_{e}(\bm{u,\delta})$ is error probability of the following vector parameter binary detection problem,

[TABLE]

This error probability can be lower bounded by the minimum probability of error of the following optimal non-coherent detector:

[TABLE]

where $\rho_{1}=|\bm{x^{H}s(u)}|$, $\rho_{2}=|\bm{x^{H}s(u+\delta)}|$. Given $\bm{u=u_{0}}$, $\rho_{1},\rho_{2}$ are rician distributed with scale parameter $s=\sigma^{2}/M$ and non-centrality parameter $\nu=|\alpha|N,|\alpha R(\bm{\delta})|N$ respectively where $R(\bm{\delta})=R(\delta_{x},\delta_{y})$ is the beampattern from (3). The error probability is given by [22]

[TABLE]

which is not a function of $\bm{u}$. For ROI in our case, the maximum error $h^{\mathrm{(max)}}=\left(\bm{a}^{T}\bm{\delta}\right)^{\mathrm{(max)}}=1$. Note that for a uniformly distributed DoA in spherical coordinates ($\theta,\phi$), the distribution of $\bm{u}$ is not uniform. However for simplicity of analysis, we make the assumption that $\bm{u}$ is uniformly distributed on a circular disc. Hence, the ZZB expression simplifies to

[TABLE]

The maximum error probability over all directions $\bm{\delta}$ cannot be expressed as closed form expression. However, due to the monotonicity of Marcum’s Q function, $Q_{1}(.)$ and Bessel function of $0^{th}$ order, $I_{0}(.)$, the error probability in (15) is maximized only when $R(\bm{\delta})$ is maximized. Therefore, for each values of $h$ we compute the $\max_{\bm{\delta}:\bm{a^{T}\delta}=h}|R(\bm{\delta})|$ numerically by search over a discrete set of points on the line segment $\bm{a^{T}\delta}=h$ and substitue in (16).

IV-B DoA estimation algorithm

Grid-based sparse estimation for a set of DoAs models the the received signal (9) as follows:

[TABLE]

where $\bm{S(\Psi)}=\left[\bm{s(u_{1})}\cdots\bm{s(u_{|\Psi|})}\right]$ contains the array response at discretized set of DoAs $u_{i}\in\bm{\Psi}$ as columns. The nonzero entries in $\bm{b}$ point to presence of target in the corresponding DoA in $\Psi$. The DoA and gain pair $(\bm{\hat{u}}_{i},\hat{\alpha}_{i})_{i=1}^{K}$ can be estimated by jointly minimizing the residual power,

[TABLE]

The NOMP algorithm summarized below provides a two stage estimator:

Detection: Using precomputed $\bm{S(\Psi)}$, coarse estimates of DoA and complex gain are obtained

[TABLE] 2. 2.

Refinement: The estimates are refined using the Newton method:

[TABLE]

where $H_{\nabla}T$ and $\nabla T$ denote the Hessian and gradient of $T(\bm{u},\alpha)$ with respect to $\bm{u}$ at current estimate $(\bm{\hat{u}},\hat{\alpha})$ (see [23] for details).

The algorithm is repeated with the residual signal, $\bm{r}=\bm{x}-\hat{\alpha}^{\prime}\bm{s(\hat{\bm{u}}^{\prime})}$ to estimate other DoAs. The refinement steps are repeated after each new detection for all DoAs in a cyclic manner for few rounds to improve accuracy.

The algorithm yields $K_{est}$ DoA estimates, with estimation performance degrading when $K_{est}$ does not match the true number of DoAs, $K_{DoA}$. Hence, in order to evaluate the arrays independent of such errors, we implement both algorithm where $K=K_{DoA}$ is known.

We also run extensive simulations with another state of the art algorithm, BPDN [14], with default parameters and 5 refinement stages. The computational complexity of BPDN is significantly higher than that of NOMP. The grid $\bm{\Psi}$ needs to be adapted at each iteration for BPDN, depending on the DoAs, whereas it remains fixed for NOMP $\bm{S(\Psi)}$, and can be precomputed, which makes it suitable for faster implementation in large arrays. Since the NOMP algorithm also yields somewhat better estimation accuracy than BPDN, we only present results obtained with NOMP here.

V Numerical Results

V-A Design of arrays

Using the combinatorial search algorithm, we create a search space of array configurations, $\bm{C}\in\mathscr{C}$ of size $|\mathscr{C}|=657,000$ for $N_{sub}=8$ subarrays. We explore this search space to obtain arrays using the placement optimization with various choices of weights given to the cost function in Eq. (8). The following key observations were made about the effect of these weights on subarray placement:

• •

When primary weight is given towards minimizing beamwidth (larger $\alpha$), the subarrays are widely distributed over the available aperture area.

• •

When we increase $\beta$ to suppress MSLL, the array becomes restricted to smaller area.

• •

The weight $\gamma$ controls both the shape of main lobe as well as the position of sidelobes.

Based on these observations, we set the weights to obtain following sample array configurations:

A1: Primary emphasis is given towards minimizing beamwidth by setting $\alpha=1$. If multiple arrays are found near the minimum beamwidth, we give secondary priority to other attributes by giving small weight to them ($\beta,\gamma=0.1$). 2. 2.

A2: In this case, we emphasize all beam attributes by setting all weights equal to 1.

We compare these array designs against two simple array configurations, 1)“compact” array where subarrays are placed together such that overall element pattern becomes a uniform rectangular array, 2) “naive” array where subarrays are spread along a diamond shape such that its resultant beamwidth is equal to that of A2. Table I lists the weights and resulting beam attributes of these arrays.

Figure 7 shows the array designs obtained using our optimization approach and their beam patterns. The A2 array has a sharp beamwidth and only $2.6$ dB worse MSLL compared to the compact array. On the other hand, a naive sparse array with circular arrangement of subarrays yields $2.6$ dB higher MSLL compared to A2 for similar beamwidth. Our designs A1, A2 exhibit several small sidelobes (the highest sidelobe for A2 is -10.2 dB), whereas the naive array exhibits fewer but more pronounced sidelobes. Since large sidelobes and grating lobes can cause large errors in DoA estimation, we expect our designs to yield better estimation performance, which is borne out by the results presented in the next section.

V-B Comparison of Estimation Performance

We evaluate the arrays based on their DoA estimation accuracy at different SNRs for both single and multiple source cases. We use the RMSE in estimating DoA for comparison which is given by $\bar{\epsilon}=\sqrt{\mathbb{E}[|\bm{\hat{u}-u_{0})|^{2}}]}=\sqrt{tr(\bm{R_{\epsilon}})/2}$.

V-B1 Estimation bounds

The Cramér Rao Bound is evaluated using (12), $CRB(\bar{\epsilon})=\sqrt{tr(CRB(\bm{R_{\epsilon}}))/2}$ . The Ziv Zakai Bound is evaluated using (16),

[TABLE]

where $\bm{a_{1},a_{2}}$ denote the directions of maximum and minimum beamwidths of the array. We also computed the ML estimation (MLE) error by Monte Carlo simulation using (11) with an overcomplete dictionary of array responses. Figure 8 shows the CRB, ZZB and MLE curves for all the arrays. CRB is proportional to beamwidth (see Appendix B for details). The ZZB bound converges to CRB at the so-called “ZZB threshold” SNR: when the SNR is below this threshold, far-ambiguities in DoA estimation caused by large sidelobes dominate the MSE. The tradeoff between beamwidth and MSLL is thus expected to translate to one between CRB (better with smaller beamwidth) and ZZB threshold (worse with larger MSLL). Thus, as expected, “A1” array achieves the lowest CRB, followed by “naive” and “A2” with equal CRB, while “Compact” array has the largest beamwidth and hence highest CRB. The trend in MSLL is weakly reflected in the ZZB thresholds (for a single target, sidelobes do limited damage): the degradation in ZZB threshold, relative to that of the compact array for the optimized arrays (A1, A2) is less than $1$ dB, while the gain in CRB due to smaller beamwidth is $4$ dB and $2$ dB, respectively. The MLE error curve also agrees with the threshold behavior predicted by ZZB. We see in the next set of results, however, that the size of the sidelobes becomes much more important when we consider multiple targets.

V-B2 Estimation algorithm performance

We obtain DoA estimates using the NOMP algorithm [15, 23] with a known number of sources to compare the best case performance of these arrays. (The NOMP algorithm also performs as well as the brute force MLE for a single target discussed earlier–results omitted here.) The RMSE is evaluated across $N=1024$ DoAs uniformly sampled over the ROI (spherical cap of half angle $30^{\circ}$). For evaluating the estimation performance in presence of multiple targets, ($K=5$) we compute the RMS error in the DoA estimate for a primary target fixed at broadside, while interfering targets are distributed uniformly in ROI at separation of $\Delta\bm{u}\geq 0.16$ or $\Delta\theta\geq 9.2^{\circ}$ away from primary target. This separation is imposed because the estimation problem is ill-posed for DoAs in close proximity. For uniform arrays, the minimum separation is typically defined with respect to the DFT bin size (e.g. $\mathrm{\Delta_{DFT}}=2\pi/L$ for an $L$-element linear array). Since this quantity cannot be defined for non-uniform planar arrays, we choose a minimum separation halfway between the RMS beamwidths of the “compact” and sparse arrays, to capture the effect of both local errors and far ambiguity errors due to sidelobes.

In addition to RMSE vs SNR curves, we also analyze the distribution of error magnitudes. The complementary cumulative distribution (CCDF) of the estimation errors is used to compare the “outage probability” corresponding to too large an error, which captures the impact of large sidelobes.

• •

With multiple sources, the estimation accuracy is degraded by interference from other sources, and RMSE does not converge to the single-target CRB. Fig. 9 shows the estimation performance for strong and weak interference.

Weak interference: When the interfering sources are 6 dB weaker than primary target, sparse arrays offer more than $5$ dB SNR gain compared to compact arrays for $SNR>-5$ dB. Also, the difference between A2, naive array widens to about $1$ dB in the threshold region indicating the benefit of suppressing sidelobes. 2. 2.

Strong interference: When interfering sources have same magnitude, RMSE severely degrades for both arrays with high sidelobes (A1, Naive) as well as arrays with high beamwidth (Compact). On the other hand, A2 has lowest RMSE at $SNR>-5$ dB because of the dual benefit of small beamwidth and lower sidelobes.

• •

The increase in estimation errors at high SNR is attributed to ambiguity errors from sidelobes, hence the overall sidelobe suppression for the arrays can be compared using the distribution of these error magnitudes. Fig 10 shows the CCDF curves of all arrays at SNR=$-5$ dB. The initial curvature of these curves (RMSE upto -22 dB) is expected to depend on local errors, hence the rate of change follows same order as CRB which is $\emph{A1}>\emph{A2 = naive}>\emph{compact}$. But the curvature reverses order at higher RMSE indicating the tradeoff with far-errors. We can see that A2 achieves lower outage probability in both scenarios (e.g. for RMSE threshold set to $-15$ dB) as it strikes a balance between near and far errors. In contrast, both A1 and naive exhibit high outage probability due to frequent far ambiguity errors caused by higher sidelobes.

Therefore, depending on the expected magnitude of interferers either one of the designed arrays with suitable sidelobe suppression can be selected. For a desired beamwidth reduction our design algorithm yields an array superior to a naively designed sparse array.

VI Compressive Estimation

We now evaluate the arrays for sparse estimation using compressive measurements at each subarray given by:

[TABLE]

where $\bm{x}$ is the full measurement from (17) and $\bm{\Phi}=\text{diag}(\bm{\Phi}_{1},\cdots,\bm{\Phi}_{N_{s}})$ is the $M\times N$ measurement matrix consisting of the subarrays measurement matrices as its block diagonals. Each subarray takes $M_{i}$ compressive measurements with an independent $\bm{\Phi}_{i}\in\mathbb{C}^{M_{i}\times N_{e}}$ whoose elements are chosen uniformly and independently from QPSK samples $\frac{1}{\sqrt{M_{i}}}\left\{\pm 1,\pm j\right\}$. In addition, columns of $\bm{\Phi}_{i}$ have unit norm to preserve signal norm on average ($\mathbb{E}\left[||\Phi\bm{S(u)}||^{2}\right]=||\bm{S(u)}||^{2}$) while scaling noise variance by $N/M$. The underlying DoA, $\bm{u}$ can be extracted by minimizing the ML cost function:

[TABLE]

The efficacy of compressive parameter estimation in AWGN depends on preserving the geometric structure of the parametrized signals [17]. Specifically, if $\Phi$ satisfies the $2K$ isometry property for discretized basis $\bm{S(\Psi)}$ [17],

[TABLE]

where $C$ is a constant for any arbitrarily chosen $2K$ sparse vector $\bm{b}$, the performance of the compressive system follows that for the original system, except for an SNR penalty of $M/N$.

Figure 11 shows the minimum and maximum values of this ratio over $10^{6}$ random realization of $8$ sparse $\bm{b}$ for sparse array. The ratio is within $[-5,3]$ dB for $M>32$ signifying that $32$ compressive measurements are sufficient to estimate $K=4$ DoAs.

Figure 12 shows estimation performance with $M=32$ compressive measurements collected across eight subarrays $(M_{i}=4,i\in\{1..8\})$. Comparing with Fig 8, we observe that the estimation algorithms preserve the same characteristics as with full measurements with approximately $6$ dB SNR penalty as expected ($N/M=4$).

VII Conclusions

Our results demonstrate that trading off beam width versus side lobes when synthesizing a large effective aperture does indeed produce performance gains in bearing estimation. Compared to a compact placement of subarrays, an optimized sparse placement produces smaller mean squared error because of its smaller beam width. Compared to a naive sparse placement, the control of side lobes via our optimized placement produces the most significant gains when estimating the bearing for multiple sources. There are a number of interesting directions for further exploration. It may be possible to improve our optimization framework by alternative approaches for pruning solutions, as well as applying generic optimization techniques (e.g., genetic optimization) initialized by our solutions. Exploring the application of our framework for communications, where transmit and receive beamforming gains are fixed by the number of elements, but control of beam width and sidelobes affects interference, is an interesting direction. In the context of sensing, our work may be viewed as design of an individual sensor which can be placed within a more comprehensive architecture, such as a network of sensors for localization and tracking.

Appendix A Mean square error in 2D DoA estimation

For 2D DoA estimation, the error along any given angle $\xi$ is given by $\bm{a^{T}R_{\epsilon}a}$, where $\bm{a}=[\cos\xi,\sin\xi]^{T}$ is the directional cosine and $\bm{R_{\epsilon}}$ is the error covariance matrix. Assuming that $\{\nu_{i},\bm{q_{i}},i=1,2\}$, denote the eigenvalues and eigenvectors of $\bm{R_{\epsilon}}$, the MSE averaged over $\bm{a}$ (assume $\xi$ uniform over $[0,2\pi]$) is given by

[TABLE]

where we have used $\mathbb{E}[\cos^{2}\xi]=\mathbb{E}[\sin^{2}\xi]=\frac{1}{2}$.

Appendix B 2D Beamwidth & CRB

We define 2D beamwidth using the Taylor series expansion of beampattern $R_{\bm{u_{o}}}(\bm{u})$ around mainlobe $R_{\bm{u_{o}}}(\bm{0})$. Since the beampattern around the main lobe. and hence the beamwidth, is invariant to beamforming direction (see II-C), we assume $\bm{u_{o}}=\bm{0}$ without loss of generality, and drop the subscript: $R_{\bm{0}}(\bm{u})\triangleq R(\bm{u})$. By taking the derivatives of (3), theTaylor series expanision up to second order is obtained as

[TABLE]

We define Half Power Beam Contour (HPBC) as the closed contour around mainbeam with $\{\bm{u}:R(\bm{u})=0.5R(\bm{0})\}$, which is approximated as an ellipse using (21) as follows:

[TABLE]

Consider the eigendecomposition of $\bm{D}^{T}\bm{D}$ given by,

[TABLE]

The eigenvectors $\bm{p}_{2},\bm{p}_{1}$ correspond to major and minor axis of HPBC ellipse respectively, and depend only on the element positions.

Figure 13 shows the mainlobe of a beam and the dotted shaded region represents its HPBC ellipse. $\check{\bm{u}}_{o},\check{\bm{u}}^{1},\check{\bm{u}}^{2}$ correspond to unit vectors in the direction of main beam, vertex and co-vertex of the HPBC where, $\check{\bm{u}}=[u,v,\sqrt{1-u^{2}-v^{2}}]$ is the unit vector towards directional cosine $\bm{u}=[u,v]$. These can be expressed as,

[TABLE]

where $\phi_{\mathrm{max}},\phi_{\mathrm{min}}$ are perpendicular azimuthal angles and $\vartheta,\varphi$ are the maximum and minimum beamwidth angles subtended from the major and minor axis of this ellipse to the mainbeam.

[TABLE]

Substituting the major and minor axis from (23) in (22), we obtain

[TABLE]

That is, the beamwidths along extremal directions are inversely proportional to the square roots of the eigenvalues of $\bm{D}^{T}\bm{D}$.

Relation to CRB

Using (13), the error covariance matrix is lower bounded by

[TABLE]

Using Appendix A, the MSE can be lowerbounded by

[TABLE]

where $\mathrm{BW^{DoA}=\sqrt{BW_{max}^{2}+BW_{min}^{2}}=\sqrt{\vartheta^{2}+\varphi^{2}}}$ is defined as MSE beamwidth ($\sin\theta\approx\theta$ for small angles $\theta$).

Appendix C Vacancy search operator $\mathscr{T}$

Our reference subarray module shown in Fig. 14 occupies space in addition to antenna elements. In order to keep element polarizations aligned, these modules can be placed in either up ($0^{\circ}$) or down ($180^{\circ}$) pose. We outline a procedure to list the vacant gridpoints $V_{i}=\mathscr{T}(C_{i}^{n})$ where the new subarray can be placed without overlapping with already placed dormant subarrays at $C_{i}^{n}$. We define the subarray state as the center $c$ of the element pattern and its pose $\nu$, since vacancies depend on both parameters.

[TABLE]

The pose variable $\nu\in\{\nu^{u},\nu^{d},\nu^{f}\}$ denotes whether subarray can be placed in up only($\nu^{u}$), down only($\nu^{d}$) or free pose ($\nu^{f}$, either up or down) at the location $c$. For a given set of dormant subarray states, $\tilde{C}_{i}^{n}$ we identify all vacant states $\tilde{V}_{i}$ for placing the new subarray. Once a new subarray is placed, the states of all dormant subarrays are updated (e.g., a free pose may switch to an up pose if the down pose becomes infeasible).

Appendix D Perturbation of array

In order to design a bin size for pruning array configurations, we analyze the effect of perturbing the location of a single array element on the eigenvalues of the array covariance matrix. Consider a small perturbation $\bm{\upsilon}=[\upsilon_{x},\upsilon_{y}]$ added to $i^{th}$ array element position: $\bar{\bm{d}}_{i}=\bm{d}_{i}+\bm{\upsilon}$. The covariance for the perturbed array is

[TABLE]

where

[TABLE]

Using Weyl’s inequality [24] for real symmetric matrices, the eigenvalue perturbation is bounded as

[TABLE]

The frobenius norms of $\bm{\mathcal{G,H}}$ are

[TABLE]

where $R_{i}=\sqrt{d_{xi}^{2}+d_{yi}^{2}}$ is the distance of the $i^{th}$ element from the array center. Hence, the overall variation of eigenvalues with variation $\Delta_{e}=\sqrt{\left(\upsilon_{xi}^{2}+\upsilon_{yi}^{2}\right)}$ of the $i^{th}$ element is

[TABLE]

Thus, the eigenvalues are more sensitive to perturbations in the locations of elements further from the center. For perturbations within one grid size used in our placement search algorithm, the eigenvalue of the subarray center covariance can vary at most by $\frac{(2R_{i}+1)}{N_{sub}}\sqrt{2}\Delta_{g}$, and we use this as a guideline for discretizing the eigenvalues for.removing geometrically similar configurations.

Acknowledgment

This work was supported in part by Semiconductor Research Corporation (SRC), DARPA, Center for Scientific Computing at UCSB and National Science Foundation under grants CNS-1518812, CNS-1518632 and CNS-0960316.