Transmit MIMO Radar Beampattern Design Via Optimization on the Complex Circle Manifol

TL;DR

This paper introduces a novel optimization method on the complex circle manifold for MIMO radar beampattern design, effectively handling non-convex constant modulus constraints and improving performance over existing methods.

Contribution

It develops a new projection, descent, and retraction (PDR) algorithm for direct optimization on the complex circle manifold, with analytical guarantees and enhanced robustness.

Findings

PDR outperforms existing wideband beampattern design methods.

The algorithm guarantees monotonic cost reduction and convergence to local minima.

Orthogonality and robustness to target mismatch are successfully incorporated.

Abstract

The ability of Multiple-Input Multiple-Output (MIMO) radar systems to adapt waveforms across antennas allows flexibility in the transmit beampattern design. In cognitive radar, a popular cost function is to minimize the deviation against an idealized beampattern (which is arrived at with knowledge of the environment). The optimization of the transmit beampattern becomes particularly challenging in the presence of practical constraints on the transmit waveform. One of the hardest of such constraints is the non-convex constant modulus constraint, which has been the subject of much recent work. In a departure from most existing approaches, we develop a solution that involves direct optimization over the non-convex complex circle manifold. That is, we derive a new projection, descent, and retraction (PDR) update strategy that allows for monotonic cost function improvement while maintaining feasibility over the complex circle manifold (constant modulus set). For quadratic cost functions (as is the case with beampattern deviation), we provide analytical guarantees of monotonic cost function improvement along with proof of convergence to a local minima. We evaluate the proposed PDR algorithm against other candidate MIMO beampattern design methods and show that PDR can outperform competing wideband beampattern design methods while being computationally less expensive. Finally, orthogonality across antennas is incorporated in the PDR framework by adding a penalty term to the beampattern cost function. Enabled by orthogonal waveforms, robustness to target direction mismatch is also demonstrated.

Figures

Tables

Table 1. Table I: Computational complexity for the state of the art algorithms under the constant modulus constraint.

Method	Complexity order
WBFIT[8]	$𝒪\left(FN{M}^2\right)$
ADMM[34]	$𝒪\left(F{N}^2{M}^2\right)$
SDR[20, 21]	$𝒪\left(N^{3.5}{M}^{3.5}\right)+𝒪\left(T{N}^2{M}^2\right)$
IA-CPC [33]	$𝒪\left(F{N}^2{M}^2\right)$
PDR	$𝒪\left(F{N}^2{M}^2\right)$

Table 2. Table II: The deviation from the desired beampattern (the cost function in Eq ( 7 )) for Case 1.

Method	$10{\text{log}}_{10}\left(\rho \left(𝐱\right)\right)$	Run time (sec)	Iterations
Unconstrained	$19.93$	-	-
WBFIT[8]	$31.93$	$0.37$	$135$
SDR[20, 21]	$25.50$	$1107$	$31$
IA-CPC [33]	$28.30$	$14.39$	$172$
ADMM[34]	$24.93$	$19.20$	$200$
PDR	$22.80$	$9.27$	$121$

Table 3. Table III: The deviation from the desired beampattern (the cost function in Eq ( 7 )) for Case 2.

Method	$10{\text{log}}_{10}\left(\rho \left(𝐱\right)\right)$	Run time (sec)	Iterations
Unconstrained	$19.52$	-	-
WBFIT[8]	$33.11$	$0.39$	$199$
SDR[20, 21]	$28.41$	$1190$	$33$
IA-CPC [33]	$30.86$	$14.47$	$180$
ADMM[34]	$28.91$	$20.36$	$208$
PDR	$26.69$	$7.38$	$115$

Table 4. Table IV: The deviation from the desired beampattern (the cost function in Eq ( 7 )) for Case 3.

Method	$10{\text{log}}_{10}\left(\rho \left(𝐱\right)\right)$	Run time (sec)	Iterations
Unconstrained	$18.18$	-	-
WBFIT[8]	$34.29$	$1.55$	$102$
SDR[20, 21]	$28.21$	$1262$	$37$
IA-CPC [33]	$31.78$	$18.58$	$192$
ADMM[34]	$28.26$	$24.10$	$231$
PDR	$27.86$	$12.14$	$142$

Table 5. Table V: Cost function and Auto-correlation level for PDR vs WeCAN and SA-Method

Method	$10{\text{log}}_{10}\left(\rho \left(𝐱\right)\right)$	$\alpha$	${\text{ISL}}_0$ (dB)	Time (sec)
LFM	$32.93$	-	$-302.53$	-
WeCAN [27]	$32.71$	-	$-102.92$	$11.52$
SimulAnneal [31]	$30.72$	-	$-3.67$	$9.29$
PDR	$27.77$	$80$	$-5.27$	$10.58$
	$28.64$	$100$	$-6.92$	$12.67$
	$29.66$	$200$	$-12.12$	$15.05$

Equations

$$s_m(t)=x_m(t)e^{j2\pi f_{c}t}$$

$${\mathbf{x}}_m=[x_m(0)x_m(1)...x_m(N-1){]}^T$$

$$y_m(p)=\sum _{n=0}^{N-1}{x}_m(n)e^{-j2\pi \frac{np}{N}},p=-\frac{N}{2},...,0,...,\frac{N}{2}-1$$

$$P(\theta ,p)=|{\mathbf{a}}^H(\theta ,p){\mathbf{y}}_p{|}^2$$

$$\text{IEEEeqnarraymulticol}3l\mathbf{a}(\theta ,p)=$$

$$\displaystyle\big{[}1\quad e^{j2\pi(\frac{p}{NT_{s}}+f_{c})\frac{d\cos\theta}{c}}\,...\,e^{j2\pi(\frac{p}{NT_{s}}+f_{c})\frac{(M-1)d\cos\theta}{c}}\big{]}^{T}$$

$${\mathbf{a}}_{sp}=\mathbf{a}({\theta }_s,p),s=1,...,S$$

$$P_{sp}=|{\mathbf{a}}_{sp}^H{\mathbf{y}}_p{|}^2=|{\mathbf{a}}_{sp}^H{\mathbf{F}}_p\mathbf{x}{|}^2$$

$$\mathbf{x}=[{\mathbf{x}}_0^T{\mathbf{x}}_1^T...{\mathbf{x}}_{M-1}^T{]}^T$$

$$\begin{array}[]{cc}\displaystyle\min_{\mathbf{x}}&\sum_{s=1}^{S}\sum_{p=-\frac{N}{2}}^{\frac{N}{2}-1}[d_{sp}-|\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|]^{2}\\ \text{s.t.: }&|\mathbf{x}|=\mathbf{1}\par\end{array}$$

$$\begin{split}&\min_{\mathbf{\phi}_{sp}}|d_{sp}e^{j\phi_{sp}}-\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|^{2}\\ &=\min_{\mathbf{\phi}_{sp}}\Big{\{}d_{sp}^{2}+|\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|^{2}\\ &\ \ \ -2\text{Re}[d_{sp}|\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|\text{cos}(\phi_{sp}-\text{arg}(\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}))]\Big{\}}\\ &=[d_{sp}-|\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|]^{2}\ \ (\text{for}\ \phi_{sp}=\arg\{\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}\})\end{split}$$

$$\begin{array}[]{cc}\displaystyle\min_{\mathbf{x},\{\phi_{sp}\}_{\forall s,p}}&f(\mathbf{x})=\sum_{s=1}^{S}\sum_{p=-\frac{N}{2}}^{\frac{N}{2}-1}|d_{sp}e^{j\phi_{sp}}-\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|^{2}\\ \text{s.t.: }&|\mathbf{x}|=\mathbf{1}\end{array}$$

$$f(\mathbf{x})$$

$${\mathbf{A}}_p=\left[\begin{array}{c}{\mathbf{a}}_{1p}^H\\ ⋮\\ {\mathbf{a}}_{Sp}^H\end{array}\right],{\mathbf{d}}_p=\left[\begin{array}{c}{d}_{1p}{e}^{j{\varphi }_{1p}}\\ ⋮\\ d_{Sp}{e}^{j{\varphi }_{Sp}}\end{array}\right],$$

$$\begin{array}[]{cc}\displaystyle\min_{\mathbf{x}}&f(\mathbf{x})=\mathbf{x}^{H}\mathbf{P}\mathbf{x}-\mathbf{q}^{H}\mathbf{x}-\mathbf{x}^{H}\mathbf{q}+r\\ \text{s.t.: }&|\mathbf{x}|=\mathbf{1}\end{array}$$

$$\begin{array}[]{cc}\displaystyle\min_{\mathbf{x}\in\mathbb{C}^{L}}&\bar{f}(\mathbf{x})=\mathbf{x}^{H}(\mathbf{P}+\gamma\mathbf{I})\mathbf{x}-\mathbf{q}^{H}\mathbf{x}-\mathbf{x}^{H}\mathbf{q}\\ \text{s.t.: }&|\mathbf{x}|=\mathbf{1}\\ \end{array}$$

$$\underset{L\text{times}}{\underbrace{\mathcal{S}\times \mathcal{S}\dots \times \mathcal{S}}}$$

$${\mathcal{S}}^L\stackrel{△}{=}\{\mathbf{x}\in {\mathbb{C}}^L:|x_l|=1,l=1,2,\dots ,L\}$$

$$\underset{\mathbf{x}\in \mathcal{M}}{\min }$$

$${\mathbf{x}}_{(k+1)}={\mathbf{x}}_{(k)}+{\beta }_k{\mathbf{\eta }}_{(k)}$$

$${\mathcal{C}}^1\stackrel{△}{=}\{\mathbf{q}\in {\mathbb{R}}^2:{\mathbf{q}}^T\mathbf{q}=1\}$$

$${\text{Proj}}_{{\mathcal{T}}_{{\bar{\mathbf{z}}}_l}{\mathcal{C}}^1}({\bar{\mathbf{w}}}_l)$$

$$\text{Ret}({\bar{\mathbf{w}}}_l)$$

$$\text{IEEEeqnarraymulticol}3l{\mathbf{P}}_{{\mathcal{T}}_{\mathbf{z}}{\mathcal{S}}^L}(\mathbf{w})$$

$$=$$

$$\mathbf{R}(\mathbf{w})$$

$$\bar{\mathbf{x}}_{(k)}=\mathbf{x}_{(k)}+\beta\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}}\big{(}\boldsymbol{\eta}_{(k)}\big{)}$$

$$\mathbf{x}_{(k+1)}=\mathbf{R}\big{(}\bar{\mathbf{x}}_{(k)}\big{)}$$

$$\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}}\big{(}\boldsymbol{\eta}_{(k)}\big{)}=\boldsymbol{\eta}_{(k)}-\text{Re}\{\boldsymbol{\eta}_{(k)}^{*}\odot\mathbf{x}_{(k)}\}\odot\mathbf{x}_{(k)}$$

$$\begin{split}\bar{\mathbf{x}}_{(k)}&=\mathbf{x}_{(k)}+\beta\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}}\big{(}\boldsymbol{\eta}_{(k)}\big{)}\\ &=\mathbf{x}_{(k)}+\beta\big{(}\boldsymbol{\eta}_{(k)}-\text{Re}\{\boldsymbol{\eta}^{*}_{(k)}\odot\mathbf{x}_{(k)}\}\odot\mathbf{x}_{(k)}\Big{)}\end{split}$$

$$\mathbf{x}_{(k+1)}=\mathbf{R}\big{(}\bar{\mathbf{x}}_{(k)}\big{)}$$

Full text

I Introduction

Multiple input multiple output (MIMO) radar, in general, transmits independent waveforms from its transmitting elements and observes the backscattered signals (from the target and from interference sources). On the other hand, in standard phased array radars, many small elements are employed so that each element emits an identical signal (up to a phase shift). These phases are shifted to focus the transmit beam in a certain direction [1]. The advantage of transmitting independent waveforms in MIMO radar provides extra degrees of freedom and the waveforms may be optimized across antennas to enhance the performance of radar systems.

A central problem in MIMO radar is to design a set of waveforms such that the transmitted beampattern matches certain specifications, e.g. a desired beampattern, either for narrowband [2, 3, 4, 5, 6] or wideband [7, 8, 9, 10, 11] setups. Although the transmit beampattern is used to focus the transmitted power in a certain directions of interest, a well-designed beampattern also helps enhance the Signal-to-Noise-Ratio (SNR)[2, 12, 13, 14]. While unconstrained design is straightforward, but this is a highly challenging problem in the presence of practical constraints on the waveform [15].

A principally important constraint on the transmit waveforms is the constant modulus constraint (CMC). The CMC is crucial in the design process due to the presence of non-linear amplifiers in radar systems [16] which must operate in saturation mode. Existing approaches that deal with beampattern design under CMC can be classified into two categories: indirect and direct approaches. The first category consists of methods that approximate or relax the CMC, i.e., the design process is conducted under approximated constraints and then the produced solution is converted to the constant modulus set. Examples in this category include: the peak-to-average ratio (PAR) waveform constraint [17], [8] and the energy constraint [18]. In general, since these approaches deal with an approximation of CMC, performance may degrade considerably in an actual real-world scenario when the constraint is strict [15]. In the second category are methods that directly enforce CMC and hence lead to better performance compared to the indirect approaches, but with relatively expensive computational procedures, such as a quasi Newton iterative method in [19], Semidefinite Relaxation (SDR) with randomization [20, 21] and the sequence of convex programs approach in [11].

In this work, our goal is to break the trade-off between performance measures such as faithfulness to the desired beampattern and computational complexity of designing/optimizing the waveform code. We show that this is possible by invoking principles of optimization over manifolds. That is, we derive a new projection, descent, and retraction (PDR) based numerical algorithm that in each iteration allows for monotonic cost function improvement while maintaining feasibility over the complex circle manifold (constant modulus set). Optimization over such a manifold has been investigated for problems in communications for example [22]. However, the work in [22] deals with an entirely different physical problem (and hence cost function) from ours and does not investigate analytical properties of the solution or convergence of the algorithm, which is a key focus of our work for quadratic loss functions.

Besides CMC, imposing orthogonality across antennas has been shown to be particularly meritorious. Orthogonal MIMO waveforms enable the radar system to achieve an increased virtual array [23, 14] and, hence leads to many practical benefits [24, 25, 26]. From a beampattern design standpoint, a compelling practical challenge is that the “directional knowledge” of target and interference sources utilized in specifying the desired beampattern may not be perfect. In such scenarios, it has been shown in [23, 14] that the gain loss in the transmit-receive patterns for orthogonal waveform transmission is very small under target direction mismatch.

Some work has been done towards the joint incorporation of CMC and orthogonality constraints [27, 28, 29, 30] under the umbrella of waveforms with desired auto-and-cross correlation properties. In these works, however, the beampattern is not designed but an outcome. Recently, in [31], a MIMO beampattern design that incorporates both CMC and orthogonality was investigated via a numerical approach based on the simulated annealing algorithm. To incorporate CMC, phase of the waveform vector is optimized numerically but analytical properties/guarantees of the solution are not investigated.

Our work addresses the aforementioned challenges in transmit MIMO beampattern design focusing particularly on tractable and scalable approaches in the presence of CMC. Specifically, our contributions include:

• •

Projection, Descent and Retraction algorithm (PDR): A new approach is developed that works directly on the complex circle, i.e. descent is achieved while maintaining feasibility in the sense of CMC. The proposed numerical update consists of three steps: (1) Projection of the gradient of the cost function onto the tangent space of the complex circle manifold, i.e. the CMC set, (2) Descent on the tangent space (affine set), and (3) Retraction back to the complex circle.

• •

Algorithm analysis and convergence guarantees: For quadratic cost functions, we formally prove that the cost function is monotonically decreasing while updating from one point to another on the complex circle and further convergence is guaranteed to a local minima.

• •

Incorporating orthogonality: We show that the aforementioned PDR technique can be applied to enforce orthogonality across antennas by introducing an orthogonality encouraging penalty term with the cost function.

• •

Numerical simulations insights and validation: We compare the PDR algorithm against the state-of-the-art approaches in MIMO beampattern design which address the CMC constraint. Results show that we can achieve better fidelity against a desired beampattern, at a remarkably lower computational cost. We also show that when orthogonality is incorporated, the PDR algorithm can achieve a beampattern design that exhibits robustness to target direction mismatch, which is hugely desirable in real-world scenarios where the specification of an idealized beampattern may not be exact.

The rest of the paper is organised as follows. The problem formulation is presented in Section II. In Section III, we provide a brief background on optimization over manifolds and develop the proposed PDR algorithm for beampattern design in the presence of the constant modulus constraint, equivalent to optimization over the complex circle manifold. New analytical results are provided in this setting that prove monotonic cost function improvement as well as convergence. Also presented in this section an extension towards incorporating the orthogonality constraint. Section IV compares and contrasts the proposed PDR algorithm against the state-of-the-art approaches in wideband transmit MIMO beampattern design. Concluding remarks and possible future research directions are discussed in Section V.

Notation: $\mathbb{C}^{N}$ and $\mathbb{R}^{N}$ denote the $N$-dimensional complex and real vector spaces, respectively. We use boldface upper case for matrices and boldface lower case for vectors. $\|\mathbf{x}\|_{2}$ is the 2-norm of the vector $\mathbf{x}$ and $\|\mathbf{X}\|_{F}$ is the Frobenius norm of the matrix $\mathbf{X}$. The transpose, the conjugate, and the conjugate transpose (Hermitian) operators are denoted by $(.)^{T}$, $(.)^{*}$, and $(.)^{H}$ respectively. $\odot$ and $\otimes$ denote the Hadamard and Kronecker product respectively. $\text{Re}(.)$ and $\text{Im}(.)$ denote extraction of the real part, and imaginary part of a complex number (or vector), respectively. $|x|$ denotes modulus of the complex number $x$ and $|\mathbf{x}|$ is a vector of element wise absolute values of $\mathbf{x}$, i.e., |\mathbf{x}|=\big{[}|x_{1}|\ |x_{2}|\ \ldots\ |x_{L}|\big{]}^{T}. $\nabla_{\mathbf{s}}(f(\mathbf{s}))$ denotes the gradient of the function $f$ w.r.t. the vector $\mathbf{s}$. $\text{vec}(.)$ denotes the vectorization operator. $\mathbf{I}_{L}$ is an $L\times L$ identity matrix. For matrices $\mathbf{A}$ and $\mathbf{B}$, $\mathbf{A}\geq\mathbf{B}\Rightarrow\mathbf{A}-\mathbf{B}\geq 0$, i.e., the matrix $\mathbf{A}-\mathbf{B}$ is positive-semi definite. $\mathbbm{1}_{M}$ denotes an $M\times M$ matrix with all ones. $\text{ddiag}(\mathbf{A})$ sets all off-diagonal entries of the matrix $\mathbf{A}$ to zero.

II Problem formulation

Consider a Uniform Linear Array (ULA) MIMO radar system that employs $M$ transmit antennas with inter-element spacing $d$, as shown in Fig. 1. The transmitted bandpass signal transmitted by the $m^{\text{th}}$ antenna is given by

[TABLE]

where $f_{c}$ is the carrier frequency and $x_{m}(t)$ is the baseband signal.

The baseband signal $x_{m}(t)$ is sampled to $N$ samples with sampling rate $T_{s}=1/B$ and the samples are collected in the following vector

[TABLE]

where $x_{m}(n)\triangleq x_{m}(t=nT_{s}),\ n=0,1,\ldots,N-1$, and $B$ is the bandwidth in Hz.

Let $y_{m}(p)$ be the discrete Fourier Transform (DFT) of $x_{m}(n)$ and expressed as:

[TABLE]

The discrete frequency beampattern in the far-field at spatial angle $\theta\in[0^{\circ},180^{\circ}]$ is [8]:

[TABLE]

where $\mathbf{a}(\theta,p)\in\mathbb{C}^{M\times 1}$ is the steering vector at frequency $\frac{p}{NT_{s}}+f_{c}$ defined as

[TABLE]

and $\mathbf{y}_{p}=[y_{0}(p)\quad y_{1}(p)\quad...\quad y_{M-1}(p)]^{T}$.

Furthermore, the spatial angle $\theta$ can be discretized by dividing the interval $[0^{\circ},180^{\circ}]$ into $S$ sub-intervals, i.e., $\{\theta_{s}\}_{s=1}^{S}$ and hence the steering vector $\mathbf{a}(\theta,p)$ can be written in terms of $\theta_{s}$ and $p$ as

[TABLE]

Using these notations, the beampattern in Eq (4) can be expressed in discrete angle-frequency as

[TABLE]

where $\mathbf{x}\in\mathbb{C}^{MN\times 1}$ is the concatenation of the waveforms vectors defined in Eq (2), i.e.,

[TABLE]

and $\mathbf{F}_{p}=\mathbf{e}_{p}\otimes\mathbf{I}_{M}$ where $\mathbf{e}_{p}=[1\quad e^{-j2\pi\frac{p}{N}}\quad...\quad e^{-j2\pi\frac{(N-1)p}{N}}]$. The beampattern design problem under constant modulus constraint can be formulated as

[TABLE]

where $d_{sp}\in\mathbb{R}$ is the desired beampattern and the constant modulus constraint ($|\mathbf{x}|=\mathbf{1}$) implies that $|x_{m}(n)|=1$ for $m=0,1,\ldots,M$ and $n=0,1,\ldots,N$. We note that the a cost function in Eq (7) has been the focus of much past work [8, 9, 10, 11] with different ways of specifying $d_{sp}$.

Note that the cost function in Eq (7) is not tractable due to the $|.|$ (absolute value) operator, this operator makes the function non-differentiable w.r.t. the variable $\mathbf{x}$. To overcome this issue, it has been shown in [8] that for a generic term $[d_{sp}-|\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}|]^{2}$ of Eq (7), the following holds

[TABLE]

In view of the above, the authors in [8] formulate the following problem over $\phi_{sp}$ and $\mathbf{x}$ jointly

[TABLE]

With this form, the beampattern design problem will be carried out over two minimization stages: one w.r.t. $\phi_{sp}$ for fixed $\mathbf{x}$ with a minimizer $\phi_{sp}=\arg\{\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}\}$ and the second one will be w.r.t. $\mathbf{x}$ for fixed $\phi_{sp}$. It is worthwhile observing that the phase variable $\phi_{sp}$ is not inherent to the problem but introduced to make the problem tractable, i.e., when the phases of $d_{sp}e^{j\phi_{sp}}$ and $\mathbf{a}^{H}_{sp}\mathbf{F}_{p}\mathbf{x}$ agree, the cost function is a quadratic w.r.t. $\mathbf{x}$. Therefore, in this work, we are seeking to optimize $f(\mathbf{x})$ w.r.t. $\mathbf{x}$ for fixed $\phi_{sp}$ under CMC.

Then, the cost function in problem (9) can be rewritten compactly, with fixed $\phi_{sp}$, as:

[TABLE]

where

[TABLE]

and $\mathbf{P}=\sum_{p}\mathbf{F}^{H}_{p}\mathbf{A}^{H}_{p}\mathbf{A}_{p}\mathbf{F}_{p}$, $\mathbf{q}=\sum_{p}\mathbf{F}^{H}_{p}\mathbf{A}^{H}_{p}\mathbf{d}_{p}$ and $r=\sum_{p}\mathbf{d}_{p}^{H}\mathbf{d}_{p}$. Consequently, the problem in (9) is reframed as

[TABLE]

This optimization problem is known to be a hard non-convex NP-hard problem [32]. Some of the best known methods that develop a solution for this form are: Wideband Beampattern Formation via Iterative Techniques (WBFIT) [8], Semi-Definite relaxation with Randomization (SDR) [20, 21], the Monotonically Error-bound Improving Technique (MERIT) [32], the Iterative Algorithm for Continuous Phase Case (IA-CPC) [33], design algorithm based on Alternating Direction Method of Multipliers (ADMM) [34]. Some of these approaches suffer from low performance accuracy (in terms of deviation from the desired beampattern) while others involve relatively expensive computational procedures. Importantly, CMC is extracted in different parts of the optimization or in some other methods approached asymptotically [11], but a direct optimization over the non-convex CMC remains elusive.

We take a drastically different approach by invoking principles of optimization over non-convex manifolds. Our focus is on developing a gradient based method, which can enable descent on the complex circle manifold (formal manifold terminology for the CMC) while maintaining feasibility.

Before describing our solution in the next Section, we make an alteration to the cost function by adding $\gamma\mathbf{x}^{H}\mathbf{x}$:

[TABLE]

where $L=MN$, and $\gamma\geq 0$ (it will be used later in Lemma III.2 to control convergence). Since the problem in (12) enforces CMC, the term $\gamma\mathbf{x}^{H}\mathbf{x}$ is constant (i.e. $\gamma\mathbf{x}^{H}\mathbf{x}=\gamma L$). Hence, the optimal solution of the problem in (11) and the optimal solution of the problem in (12) are identical for any $\gamma\geq 0$.

III Constant modulus constraint and optimization over manifolds

The search space or the feasible set of the problem in (12) can be thought of as the product of $L$ (complex) circles, i.e.,

[TABLE]

where $\mathcal{S}$ is one (complex) circle which is defined as $\mathcal{S}\overset{\triangle}{=}\{x\in\mathbb{C}:x^{*}x=\text{Re}\{x\}^{2}+\text{Im}\{x\}^{2}=1\}$. This set ($\mathcal{S}$) can be seen as a sub-manifold of $\mathbb{C}$ [35] and hence, the product of such $L$ circles is a sub-manifold of $\mathbb{C}^{L}$ [35]. This manifold is known as the complex circle manifold and defined as

[TABLE]

Before we proceed with the solution of the optimization problem in (12), we first provide a background on optimization over manifolds [36] and subsequently develop a new technique to optimize directly over the complex circle manifold.

III-A Optimization over manifolds

The term optimization over manifolds refers to a class of problems of the form

[TABLE]

where $g(\mathbf{x})$ is a smooth real-valued function and $\mathcal{M}$ is some manifold. Many classical line-search algorithms from unconstrained nonlinear optimization in $\mathbb{C}^{L}$ such as gradient descent can be used in optimization over manifolds but with some modifications. In general, line-search methods in $\mathbb{C}^{L}$ are based on the following update formula

[TABLE]

where $\boldsymbol{\eta}_{(k)}\in\mathbb{C}^{L}$ is the search direction and $\beta_{k}\in\mathbb{R}$ is the step size. The most obvious choice for the search direction is the steepest descent direction which is the negative gradient of $g(\mathbf{x})$, i.e., $\boldsymbol{\eta}_{(k)}=-\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)})$. In the literature [36, 37], the following high level structure is suggested:

• •

The descent will be performed on the manifold itself rather than in the Euclidean space by means of the intrinsic gradient. The intrinsic gradient $\nabla_{\mathcal{M}}g(\mathbf{x}_{(k)})$ of $g(\mathbf{x})$ at point $\mathbf{x}_{(k)}\in\mathcal{M}$ is a vector in the tangent space $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}$ (for the definition of $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}$, see [36], Section 3.5.7). This intrinsic gradient can be obtained by projecting the standard (Euclidean) gradient $\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)})$ onto $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}$ by means of a projection operator \mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}}\big{(}\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)})\big{)}.

• •

The update is performed on the tangent space along the direction of \mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}}\big{(}\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)})\big{)} with a step $\beta$, i.e., $\bar{\mathbf{x}}_{(k)}=\mathbf{x}_{(k)}-\beta\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}}(\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)}))$ where $\bar{\mathbf{x}}_{(k)}\in\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}$.

• •

Since $\bar{\mathbf{x}}_{(k)}\notin\mathcal{M}$, it will be mapped back to the manifold by the means of a retraction operator, \mathbf{x}_{(k+1)}=\mathbf{R}\big{(}\bar{\mathbf{x}}_{(k)}\big{)}.

For some manifolds, the projection $\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{M}}(.)$ and retraction $\mathbf{R}(.)$ operators admit a closed form. Interested readers may refer to [36] for details. In the following subsections, we develop new results that employ the ideas articulated above in optimization over manifolds to the complex circle manifold (constant modulus set), i.e., $\mathcal{M}=\mathcal{S}^{L}$. In particular, 1) we derive expressions for the projection and retraction operators for the complex circle manifold $\mathcal{S}^{L}$, 2) establish new analytical results which include proof of monotonic cost function improvement while maintaining feasibility in $\mathcal{S}^{L}$, and 3) provide convergence guarantees.

III-B Complex circle manifold

The projection and retraction operators for the complex circle manifold are inspired by those for the unit circle manifold [36], which is defined as

[TABLE]

That is true because each entry of any feasible vector $\mathbf{x}\in\mathcal{S}^{L}$ can be viewed as a point in $\mathbb{R}^{2}$ with its components on the (real) unit-circle ($\mathcal{C}^{1}$). Thus, in order to derive these operators for the complex circle manifold we need to start with the operators for the unit circle manifold.

First, let us define the following two vectors $\mathbf{z}\in\mathcal{S}^{L}$ and $\mathbf{w}\in\mathbb{C}^{L}$, and let $\bar{\mathbf{z}}_{l}$ and $\bar{\mathbf{w}}_{l}$ be 2-dimensional vectors represent the components of the $l^{\text{th}}$ element of $\mathbf{z}$ and $\mathbf{w}$ respectively, i.e.,

• •

$\bar{\mathbf{z}}_{l}=[\text{Re}\{z_{l}\}\ \text{Im}\{z_{l}\}]^{T}=[z_{lr}\ z_{li}]^{T}\in\mathbb{R}^{2}$

• •

$\bar{\mathbf{w}}_{l}=[\text{Re}\{w_{l}\}\ \text{Im}\{w_{l}\}]^{T}=[w_{lr}\ w_{li}]^{T}\in\mathbb{R}^{2}$

where $l=1,2,\ldots,L$. Note that since the vector $\mathbf{z}\in\mathcal{S}^{L}$, then the vector $\bar{\mathbf{z}}_{l}\in\mathcal{C}^{1}$, i.e., $\bar{\mathbf{z}}_{l}^{T}\bar{\mathbf{z}}_{l}=1$.

The projection of $\bar{\mathbf{w}}_{l}\in\mathbb{R}^{2}$ to the tangent space of the unit circle manifold111The tangent space for unit circle manifold at $\mathbf{q}\in\mathcal{C}^{1}$ is defined as $\mathcal{T}_{\mathbf{q}}\mathcal{C}^{1}=\{\mathbf{h}\in\mathbb{R}^{2}:\mathbf{q}^{T}\mathbf{h}=0\}$[36]. $\mathcal{T}_{\bar{\mathbf{z}}_{l}}\mathcal{C}^{1}$ at $\bar{\mathbf{z}}_{l}$ is given by [36]:

[TABLE]

and the retraction of any $\bar{\mathbf{w}}_{l}\in\mathbb{R}^{2}$ to $\mathcal{C}^{1}$ is given by [36]

[TABLE]

These operators for projecting $\bar{\mathbf{w}}_{l}$ onto $\mathcal{T}_{\bar{\mathbf{z}}_{l}}\mathcal{C}^{1}$ and retracting it to $\mathcal{C}^{1}$ are illustrated in Fig. 2.

According to our discussion at the beginning of this subsection, the relation between the projection of $w_{l}$ onto the tangent space222The tangent space of $\mathcal{S}$ at the point $z_{l}\in\mathcal{S}$ is defined as [35] $\mathcal{T}_{z_{l}}\mathcal{S}=\{y\in\mathbb{C}:\text{Re}\{y^{*}z_{l}\}=0\}$ $\mathcal{T}_{z_{l}}\mathcal{S}$ and the projection of $\bar{\mathbf{w}}_{l}$ onto the tangent space $\mathcal{T}_{\bar{\mathbf{z}}_{l}}\mathcal{C}^{1}$ is similar to the relation between $w_{l}$ and $\bar{\mathbf{w}}_{l}$ which is $\bar{\mathbf{w}}_{l}=[\text{Re}\{w_{l}\}\ \text{Im}\{w_{l}\}]^{T}$. Given this relation, the projection defined in Eq (17) can be used to project $\mathbf{w}$ element-wise onto the tangent space333This tangent space is the product of $L$ tangent spaces of those for the manifold $\mathcal{S}$, i.e., $\mathcal{T}_{\mathbf{z}}\mathcal{S}^{L}=\mathcal{T}_{z_{1}}\mathcal{S}\times\mathcal{T}_{z_{2}}\mathcal{S}\ldots\times\mathcal{T}_{z_{L}}\mathcal{S}$. $\mathcal{T}_{\mathbf{z}}\mathcal{S}^{L}$, denoted as $\mathbf{P}_{\mathcal{T}_{\mathbf{z}}\mathcal{S}^{L}}(\mathbf{w})$. Rearranging the entries of the 2-dimensional vector in Eq (17) as real and imaginary components for the $l^{\text{th}}$ entry of $\mathbf{P}_{\mathcal{T}_{\mathbf{z}}\mathcal{S}^{L}}(\mathbf{w})$ yields

[TABLE]

Similarly, for the retraction operator in $\mathbb{C}^{L}$, rearranging the entries of the vector in Eq (18) for each entry of $\mathbf{R}(\mathbf{w})$ yields

[TABLE]

III-C Projection, Descent and Retraction (PDR) algorithm

Given the projection in Eq (19) and the retraction in Eq (20), the optimization steps over manifolds (described in the subsection III-A) can be computed to solve the problem in (14) over the complex circle manifold ($\mathcal{M}=\mathcal{S}^{L}$). Precisely, this problem can be solved iteratively by preforming the following steps at each iteration $k$:

A Projection of the search direction $\boldsymbol{\eta}_{(k)}=-\nabla_{\mathbf{x}}g(\mathbf{x}_{(k)})$ onto the tangent space $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}$ using Eq (19). 2. 2.

A Descent on this tangent space to update the current value of $\mathbf{x}_{(k)}$ on the tangent space $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}$ as

[TABLE] 3. 3.

A Retraction of this update to $\mathcal{S}^{L}$ by using Eq (20) as

[TABLE]

Note that $\mathbf{x}_{(k)},\mathbf{x}_{(k+1)}$ are both on the complex circle manifold, i.e. CMC points, while $\bar{\mathbf{x}}_{(k)}$ is generally a non-CMC point with magnitude $\geq 1$. The proposed algorithm from these steps is named as Projection-Descent-Retraction (PDR) and it is visually illustrated in Fig. 3.

PDR for beampattern design: The cost function of the beampattern design problem defined in (12) is a quadratic w.r.t. the complex variable $\mathbf{x}$ and its gradient is given by $\nabla_{\mathbf{x}}\bar{f}(\mathbf{x})=2({\mathbf{P}}+\gamma{\mathbf{I}}){\mathbf{x}}-2{\mathbf{q}}$. The procedure of minimizing the cost function $\bar{f}(\mathbf{x})$ using the PDR approach is formally described in Algorithm 1. The convergence of Algorithm 1 is studied in the following sub-section. Assuming that Algorithm 1 converges, the convergence of Algorithm 2 (which alternates between $\phi_{sp}$ and $\mathbf{x}$ and achieves practical beampattern design) is guaranteed and already established in past works [38, 39, 8, 11].

Computational complexity of PDR (Algorithm 1): As can be inferred from the step-wise description of Algorithm 1, there are two key steps of complexity $\mathcal{O}(L^{2})$ and $\mathcal{O}(L)$ respectively. For large $L$, PDR’s complexity per iteration when optimizing quadratic cost functions is dominantly $\mathcal{O}(L^{2})$, where $L=NM$. Table I shows computational complexity for the following state-of-the-art approaches along with PDR: Wideband Beampattern Formation via Iterative Techniques (WBFIT) [8], Semi-Definite relaxation with Randomization (SDR) [20, 21], Iterative Algorithm for Continuous Phase Case (IA-CPC) [33] and design algorithm based on Alternating Direction Method of Multipliers (ADMM) [34], where $F$ is the total number of iterations and $T$ denotes the number of randomization trails for SDR. From this table, it can be seen that PDR exhibits lower complexity compared to SDR and similar complexity to ADMM and IA-CPC (per iteration). The ADMM and IA-CPC approaches however need more iterations (larger $F$) to achieve the same performance as PDR (as demonstrated in Section IV). Note that the complexity of competing methods is reported as derived in their respective/past work.

III-D Convergence Analysis

The update \bar{\mathbf{x}}_{(k)}=\mathbf{x}_{(k)}+\beta\mathbf{P}_{\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}}\big{(}\boldsymbol{\eta}_{(k)}\big{)} in Step 4 of Algorithm 1 will produce a point on the tangent space $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}$ with a decrease in the cost if the step size $\beta$ is chosen carefully. A condition on the step size that ensures a decrease in the cost function is provided in the following lemma.

Lemma III.1

Let $\lambda_{\mathbf{P}+\gamma\mathbf{I}}$ denote the largest eigenvalue of the matrix ($\mathbf{P}+\gamma\mathbf{I}$). If the step size $\beta$ satisfies

[TABLE]

then $\bar{f}(\mathbf{x}_{(k)})\geq\bar{f}(\bar{\mathbf{x}}_{(k)})$.

Proof:

See subsection A-A of the Appendix. ∎

In the next lemma, we show that the cost function $\bar{f}(\mathbf{x})$ is non-increasing through the retraction step given that the positive scalar $\gamma$ satisfies a certain condition.

Lemma III.2

Let $\lambda_{\mathbf{P}}$ denote the largest eigenvalue of the matrix $\mathbf{P}$. If

[TABLE]

then $\bar{f}(\bar{\mathbf{x}}_{(k)})\geq\bar{f}(\mathbf{x}_{(k+1)})$.

Proof:

See subsection A-B of the Appendix. ∎

In the following lemma, we show that the original cost function $f(\mathbf{x})$ defined in (11) is non-increasing and the iterative procedure converges.

Lemma III.3

Given $\gamma\geq\frac{L}{8}\lambda_{\mathbf{P}}+\|\mathbf{q}\|_{2}$ and $0<\beta<1/{\lambda_{\mathbf{P}+\gamma\mathbf{I}}}$ the sequence $\{\bar{f}(\mathbf{x}_{(k)})\}_{k=0}^{\infty}$ generated by Algorithm 1 is non-increasing (from Lemmas III.1 and III.2) and hence the sequence $\{f(\mathbf{x}_{(k)})\}_{k=0}^{\infty}$ is also non-increasing. Moreover, since $f(\mathbf{x})\geq 0$ (sum of norms), $\forall\ \mathbf{x}$ , it is bounded below and converges to a finite value $f^{*}$.

Proof:

See subsection A-C of the Appendix. ∎

Remark: The proposed PDR has essentially enabled a gradient based update while maintaining feasibility on the (non-convex) complex circle manifold with guarantees of monotonic cost function decrease and convergence. We note that the guarantees provided here may not necessarily generalize to other non-convex manifolds. It is indeed structure specific to this problem that enables our construction as shown in Fig. 3.

III-E Orthogonal waveform design across antennas

The feasible set of the optimization problem that incorporating CMC and orthogonality to design the beampattern can be understood as the intersection of the aforementioned complex circle manifold and the complex Stiefel manifold [36, 40, 41]. Working directly on the intersection of these manifolds is difficult (or impossible) because the intersection of two manifolds is not always a manifold, and even when it is, it may not be easy to describe. Our strategy to deal with this problem is to optimize over the complex circle manifold while modifying the beampattern design cost function with the addition of a penalty term that emphasizes orthogonality. Specifically, the cost function in (12) can be altered by adding the following penalty term $\alpha\|\mathbf{X}^{H}\mathbf{X}-N\mathbf{I}_{M}\|^{2}_{F},\ \alpha>0$, where $\mathbf{X}\in\mathbb{C}^{N\times M}$ is the transmit waveform matrix $\mathbf{X}=[\mathbf{x}_{0}\quad\mathbf{x}_{1}\quad...\quad\mathbf{x}_{M-1}]$ and related to the vector $\mathbf{x}$ defined in Eq (6) through the vectorization operator, i.e., $\mathbf{x}=\text{vec}(\mathbf{X})$. This way, $\mathbf{X}$ will be “encouraged” to be orthogonal/unitary, and with this penalty term, the optimization problem that assimilates both constraints can be written as

[TABLE]

The gradient of the cost function $h(\mathbf{x})$ w.r.t. $\mathbf{x}$ can be computed by utilizing the relation between the variables $\mathbf{x}$ and $\mathbf{X}$. The penalty term is a scalar function of the matrix $\mathbf{X}$, and the gradient of a scalar function with respect to $\mathbf{X}$ and $\mathbf{x}$ can be also related through the vec operator [42], i.e.,

[TABLE]

using this relation, the gradient of the penalty term w.r.t. $\mathbf{x}$ is given by

[TABLE]

and hence the gradient of the cost function $h(\mathbf{x})$ is given by

[TABLE]

The optimization problem in (26) can now be solved by executing Algorithm 1 but with a different cost function, i.e. $h(\mathbf{x})$. As before, practical beampattern design in the presence of CMC while encouraging orthogonality can be obtained by Algorithm 2, which invokes Algorithm 1 in Step 6.

IV Numerical Results

Various numerical simulations are provided to assess the performance of PDR based beampattern design and compare it to state-of-the-art approaches. Results from the following three simulations are reported next: 1) beampattern design under the constant modulus constraint, 2) beampattern design under both constant modulus and orthogonality constraints, and 3) an investigation to examine the robustness of the produced waveforms under the two constraints to direction mismatch.

Unless otherwise specified, the following settings are used in this section. Consistent with past work [8, 11], we assume a ULA MIMO radar system with the following parameters: the number of transmit antennas $M=10$, the number of time samples $N=32$, carrier frequency $f_{c}=1$ GHz, bandwidth $B=200$ MHz, sampling rate $T_{s}=1/B$, inter-element spacing $d=c/2(f_{c}+B/2)$, and the spatial angle $\theta$ is divided into $S=180$ points.

IV-A Beampattern design under CMC

We examine the beampattern design problem under CMC and compare the performance of PDR to the following state-of-the-art approaches: 1) Wideband Beampattern Formation via Iterative Techniques (WBFIT) [8], 2) Semi-Definite relaxation with Randomization (SDR) [20, 21], 3) Iterative Algorithm for Continuous Phase Case (IA-CPC) [33], and 4) Design algorithm based on Alternating Direction Method of Multipliers (ADMM) [34]. Similar to our work, SDR and IA-CPC are devised as approaches that optimize quadratic cost functions while enforcing CMC. That is, IA-CPC and SDR will be applied to optimize the same cost function that PDR also addresses. Finally, we also report results for the case where no constraints are posed on the waveform code $\mathbf{x}$. This unconstrained design is impractical but provides a bound on the performance of all constrained methods. For a fair comparison, PDR and the competing methods are initialized with the same waveform; a pseudo-random vector of unit magnitude complex entries.

We consider three distinct specifications of the desired beampattern. Case 1 (based on [8]) only has angular dependence and has been specified to uniformly illuminate a broadside region. Case 2 (based on [11]) has both angle and frequency dependence. Case 3 has more specifications such as restricting the transmission to be in a certain frequency band for spectrally crowded scenarios [43, 44]. The step size for PDR was chosen as $\beta=0.00005$, $\beta=0.00004$, and $\beta=0.00004$ for Case 1, Case 2, and Case 3, respectively. The parameters for all competing methods were set as prescribed in their respective papers or by using code given directly by the authors.

Case 1: The desired beampattern is given by

[TABLE]

In Table II, the values of the deviation from the desired beampattern $\rho(\mathbf{x})$ (where $\rho(\mathbf{x})$ is the cost function in Eq (7)) are reported. Table II confirms that PDR provides substantial gains, about 2.13 dB, 2.7 dB, and 5.5 dB over ADMM, SDR, and IA-CPC, respectively.

In Fig. 4 a 2D visualization of the designed beampattern is shown for each of the competing methods. Clearly, PDR achieves a beampattern that is closest to the unconstrained case, which naturally serves as a bound on the performance.

Case 2: The desired beampattern is given by

[TABLE]

Similar to Case 1, the values of the deviation from the desired beampattern are listed in Table III. Clearly, the PDR algorithm gives the closest value to the unconstrained case with a gap of 2 dB over the second best method. In Fig. 5, the designed beampattern is visualized for all the competing methods. Clearly, the beampattern that results from PDR is closer to the desired one than those resulting form competing methods.

Case 3: For this scenario, the beampattern will be suppressed in two angular-frequency regions as follows

[TABLE]

The transmission is restricted in certain practical frequency bands in accordance with [43, 44]. Precisely, these frequencies are: $\ f=[1.025\ \text{GHz},1.0625\ \text{GHz}]$ (see Fig. 6).

As Table IV reveals, also in this case, PDR outperforms all competing methods in terms of deviation from the desired beampattern.

Remark: Note that Tables II, III, and IV also report run time (in seconds) as the time taken to optimize the waveform code $\mathbf{x}$. For fairness of comparison, we used the same platform for all implementations: MATLAB R16, CPU Core i5, 3.1 GHz and 8 GB of RAM. Overall, Tables II, III, and IV show that the different methods exhibit complementary merits. WBFIT, one of the earliest methods proposed to solve this problem is the fastest, but its deviation from the idealized beampattern is highest. With more sophisticated optimization techniques, IA-CPC and SDR offer performance gains but also increase complexity. It is readily apparent from these tables that PDR is highly efficient computationally bettered only by WBFIT. And PDR offers nearly 7 dB of gain in performance over WBFIT. Hence the results corroborate our assertion that PDR provides the most favorable complexity-performance trade-off.

IV-B Joint CMC and orthogonality constraints

In this numerical simulation, the beampattern design under the CMC and orthogonality constraints using PDR (run with a step size $\beta=0.00003$) will be examined. The level of orthogonality will be measured by using the following quantity:

[TABLE]

where ISL is the Integrated Sidelobe Level between the transmitted signal defined in [27]. Orthogonality across antennas is equivalent to the auto correlation case (ISL at time-lag [math]), i.e., $\text{ISL}_{0}$. In terms of the desired beampattern, we consider the same scenario as Case 1 in the previous subsection. The values of the deviation from the desired beampattern are reported in Table V. PDR is now compared against approaches that also directly or approximately capture both CMC and orthogonality; this includes: 1) the well-known and widely used linear frequency modulated (LFM ) waveform code [45, 46], 2) Weighted-cyclic algorithm-new (WeCAN) [27], and 3) the recent simulated annealing based approach (SimulAnneal) [31], which is one of the few known techniques that performs explicit beampattern design under both CMC and orthogonality constraints.

First, compared to Table II, the deviation as reported in Table V is higher. This is of course to be expected because we are not just enforcing CMC, but also orthogonality which means a smaller feasible set of waveform codes to optimize over. The results in Table V reveal that LFM and WeCAN lead to somewhat high deviation - this is unsurprising since neither approach explicitly designs the beampattern. When $\alpha=100$, PDR leads to waveform codes that are orthogonal for all practical purposes as evidenced by the ISL measure. As Table V reveals, PDR achieves the closest beampattern to the desired one (lowest deviation in dB) with SimulAnneal as the second best. Remarkably, even with $\alpha=200$, i.e. when the emphasis on orthogonality is particularly strong (ISL of $-12.12$ dB for PDR vs. $-3.67$ for SimulAnneal), PDR provides a gain of 1 dB. The gain of PDR over SimulAnneal is about 3 dB for comparable ISL values.

IV-C Benefits of orthogonality: robustness to direction mismatch

The transmit-receive pattern $G_{TR}(\theta,\theta_{0})$ measures the beamformer response when the beam is digitally steered in the direction $\theta$ and when the true location of the target is at angle $\theta_{0}$. It has been shown in [23] and [14] (as an advantage of orthogonal over coherent signals) that for orthogonal signalling, the effect of the direction mismatch (when the target is not located in the center of the transmit beam) is minimal.

The pattern $G_{TR}(\theta,\theta_{0})$ can be expressed as

[TABLE]

where $M_{R}$ is the number of receiver antennas, $\mathbf{a}_{T}(\theta)$ and $\mathbf{a}_{R}(\theta)$ are the steering vectors on the transmitter and receiver sides, respectively. These steering vectors are defined in Equations (4.6) and (4.7) in [14]. $\mathbf{R}_{s}$ is the transmit signal correlation matrix defined as

[TABLE]

where $\mathbf{X}\in\mathbb{C}^{N\times M}$ is the transmit waveform matrix defined in Section III-E. For the two extreme cases 1) coherent transmission444In coherent transmission, the transmit signals from all antennas are phase-shifted versions from one reference signal [14]. ($\mathbf{R}_{s}={\mathbbm{1}_{M}}$) and 2) orthogonal transmission ($\mathbf{R}_{s}=\mathbf{I}_{M}$ ) [23], the transmit-receive patterns for these cases will be

[TABLE]

[TABLE]

Fig. 7 shows the pattern $G_{TR}(\theta,\theta_{0})$ for (a) LFM ($\mathbf{R}_{s}=N\mathbf{I}_{M}$), (b) Coherent Transmission, (c) WBFIT, (d) IA-CPC, and (e) PDR ($\alpha=200$) signals where the beam in the transmit mode is directed to $\theta=125^{\circ}$ and the true location of the target is such that $\Delta\theta=(\theta-\theta_{0})=0^{\circ},\ 10^{\circ},\ 20^{\circ}$. The desired beampattern for this simulation is the same as Case 1 in Section IV-A, except that PDR was now optimized to generate orthogonal waveforms.

The advantage of using orthogonal signals (LFM) over a coherent scheme is very noticeable: the gain loss is nonexistent for LFM when the target deviates from the transmission direction (i.e., target is not located in the center of the transmit beam). Remarkably, PDR signals with $\alpha=200$ achieve results comparable to the LFM case, whereas in the absence of orthogonal processing, WBFIT in Fig. 7 (c) and IA-CPC in Fig. 7 (d) like the coherent case in Fig. 7 (b) suffer significant loss in mainlobe strength under target direction mismatch.

The focus of our work is on transmit waveform design but in future investigations, receive processing may also be optimized to obtain the most desirable transmit-receive beampattern.

V Discussions and Conclusion

We consider the problem of transmit beampattern design for MIMO radar under compelling practical constraints. The non-convex constant modulus constraint (CMC) is our main focus whose presence is known to yield a hard optimization problem. For tractability, CMC is addressed in the literature often by relaxations and approximations or by approaches that are computationally burdensome. Our proposed PDR algorithm invokes the principles of optimization over manifolds to address the non-convex CMC and we demonstrate via simulations that the said PDR offers a favorable performance-complexity trade-off. Analytical guarantees of monotonic cost function decrease and convergence are provided for quadratic cost functions that arise in beampattern design. Finally, a tractable extension was developed to incorporate orthogonality of waveforms across antennas. Experimentally, the benefit of orthogonality combined with CMC is the synthesis of practical, real-world waveforms that demonstrate robustness to target direction mismatch.

A viable future work direction is to exploit for CMC constrained beampattern design new algorithms that use the framework of the sequential approximation such as [47] and recent advances in [48], [28]. Of particular interest is investigating KKT optimality of the resulting solution.

Acknowledgement

The authors would like to thank the authors of the work [33] for providing us with the MATLAB code for the Iterative Algorithm for Continuous Phase Case (IA-CPC) algorithm.

Appendix A Appendix A

In this section we provide the proof of the three lemmas that presented in this work and the equivalence between Eq (42) and the retraction step.

A-A Proof of Lemma III.1

Before we establish the proof, the projection operator in Eq (19) will be reformulated as follows. The projection of a vector $\mathbf{w}\in\mathbb{C}^{L}$ onto the tangent space $\mathcal{T}_{\mathbf{z}}\mathcal{S}^{L}$ at a point $\mathbf{z}\in\mathcal{S}^{L}$ can be rewritten using the elementary properties the Hadamard product and the fact that \text{Re}\{\mathbf{w}^{*}\odot\mathbf{z}\}=\frac{1}{2}\big{[}\mathbf{w}^{*}\odot\mathbf{z}+\mathbf{w}\odot\mathbf{z}^{*}\big{]} as

[TABLE]

where $\mathbf{D}_{\mathbf{z}}=\text{ddiag}(\mathbf{z}\mathbf{z}^{T})$, and $\tilde{\mathbf{w}}=\mathbf{D}_{\mathbf{z}}\mathbf{w}^{*}$. Note that since $\mathbf{z}\in\mathcal{S}^{L}$, then $\mathbf{D}_{\mathbf{z}}^{H}\mathbf{D}_{\mathbf{z}}=\mathbf{D}_{\mathbf{z}}\mathbf{D}_{\mathbf{z}}^{H}=\mathbf{I}$, and hence

[TABLE]

For notation simplicity, will use $\mathbf{w}$ instead of $\nabla_{\mathbf{x}}\bar{f}(\mathbf{x}_{(k)})$ to represent the gradient of $\bar{f}(\mathbf{x})$ defined in Eq (12) at iteration $k$, i.e., $\mathbf{w}=\nabla_{\mathbf{x}}\bar{f}(\mathbf{x}_{(k)})=2\mathbf{R}\mathbf{x}_{(k)}-2{\mathbf{q}}$, and hence the search direction $\boldsymbol{\eta}_{(k)}$ will be

[TABLE]

The projection of the search direction $\boldsymbol{\eta}_{(k)}$ onto the tangent space $\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}$ at a point $\mathbf{x}_{(k)}$ will be computed using the new projection formula Eq (34), and with a step $\beta$ along this direction starting from $\mathbf{x}_{(k)}$, the update on the tangent space will be

[TABLE]

Proof:

Let $\mathbf{R}=\mathbf{P}+\gamma\mathbf{I}$, the value of the cost function $\bar{f}(\mathbf{x})=\mathbf{x}^{H}\mathbf{R}\mathbf{x}-\mathbf{q}^{H}\mathbf{x}-\mathbf{x}^{H}\mathbf{q}$ on the tangent space at the point $\bar{\mathbf{x}}_{(k)}$ (using the value of $\bar{\mathbf{x}}_{(k)}$ in Eq (36)) will be

[TABLE]

Now, the difference $\bar{f}(\mathbf{x}_{(k)})-\bar{f}(\bar{\mathbf{x}}_{(k)})$ will be

[TABLE]

where in Eq (37) the observation in Eq (35) is used.

Recall the definition of the matrix $\mathbf{P}=\sum_{p}\mathbf{F}^{H}_{p}\mathbf{A}^{H}_{p}\mathbf{A}_{p}\mathbf{F}_{p}$, for each $p$ the matrix $(\mathbf{A}_{p}\mathbf{F}_{p})^{H}(\mathbf{A}_{p}\mathbf{F}_{p})$ is positive semi definite, i.e., $\mathbf{y}^{H}(\mathbf{A}_{p}\mathbf{F}_{p})^{H}(\mathbf{A}_{p}\mathbf{F}_{p})\mathbf{y}=\|\mathbf{A}_{p}\mathbf{F}_{p}\mathbf{y}\|^{2}_{2}\geq 0\ \forall\ \mathbf{y}\in\mathbb{C}^{L}$ and any non-negative linear combination of positive semidefinite matrices is positive semidefinite (Observation 7.1.3 [49]), then $\mathbf{P}$ is positive semidefinite. Since $\mathbf{P}$ is positive semidefinite, then $\mathbf{R}$ is positive definite. Now, using Theorem 7.1 in [50], the matrix $\mathbf{R}$ can be diagonalized as

[TABLE]

where $\Lambda$ is diagonal matrix with the eigenvalues of $\mathbf{R}$ in the main diagonal, and $\mathbf{U}$ is unitary, i.e., $\mathbf{U}^{H}\mathbf{U}=\mathbf{U}\mathbf{U}^{H}=\mathbf{I}$.

Using Eq (39) in Eq (38) yields

[TABLE]

where $\mathbf{h}=\mathbf{U}^{H}(\mathbf{w}-\tilde{\mathbf{w}})\in\mathbb{C}^{L}$. If the step $\beta$ is chosen in a way such that $\frac{\beta}{4}\mathbf{h}^{H}(\mathbf{I}-\beta\Lambda)\mathbf{h}\geq 0\ \forall\ \mathbf{h}$ or the matrix $(\mathbf{I}-\beta\Lambda)$ is positive semi-definite and $\beta\geq 0$, then $\bar{f}(\mathbf{x})-\bar{f}(\bar{\mathbf{x}})$ will be non-negative. Then, the following condition on $\beta$ must be satisfied

[TABLE]

where $\lambda_{\mathbf{R}}$ is the largest eigenvalue of the matrix $\mathbf{R}$ ($\lambda_{\mathbf{R}}>0$ since $\mathbf{R}$ is positive definite). ∎

A-B Proof of Lemma III.2

Proof:

: Recall that the input to retraction operator , $\bar{\mathbf{x}}_{(k)}$, is the update of $\mathbf{x}_{(k)}$ on the tangent space, i.e., $\bar{\mathbf{x}}_{(k)}\not\in\mathcal{S}^{L}$ and hence $|\bar{{x}}_{l(k)}|\geq 1\ \forall\ l$. On the other hand, the output from the retraction step is a constant modulus point, i.e, $\mathbf{x}_{(k+1)}=\mathbf{R}(\bar{\mathbf{x}}_{(k)})\in\mathcal{S}^{L}$. Then, $\bar{\mathbf{x}}_{(k)}$ can be written in terms of $\mathbf{x}_{(k+1)}$ as

[TABLE]

where $\Psi=\text{diag}(\psi_{1},\psi_{2},\ldots,\psi_{L})$ is non-negative diagonal matrix. Note that Eq (42) is derived from the retraction step, see subsection A-D of this Appendix. Using Eq (42) the difference between the values of the cost function at $\bar{\mathbf{x}}_{(k)}$ and $\mathbf{x}_{(k+1)}$ will be

[TABLE]

The inequality Eq (44) holds since $(\mathbf{P}+\gamma\mathbf{I})$ is positive semidefinite and the equality Eq (44) holds since $\mathbf{x}_{(k+1)}^{H}\Psi\mathbf{x}_{(k+1)}=\|\Psi\mathbf{x}_{(k+1)}\|_{1}$ and this is true since $\mathbf{x}_{(k+1)}$ is constant modulus vector. To go further in this simplification, and since $\mathbf{P}\geq 0$ and $\Psi\geq 0$, the following theorem about \big{(}\Psi\mathbf{P}+\mathbf{P}\Psi\big{)} (Theorem 7.5 in [50]) can be utilized

[TABLE]

where $\lambda_{\Psi}$ and $\lambda_{\mathbf{P}}$ are the largest eigenvalue of $\Psi$ and $\mathbf{P}$, respectively. Using Eq (46) in Eq (44) yields

[TABLE]

where Eq (A-B) holds from Eq (46), Eq (A-B) holds since $\|\Psi\mathbf{x}_{(k+1)}\|_{\infty}=\text{max}_{l}|\psi_{l}x_{l(k+1)}|=\lambda_{\Psi}$, Eq (A-B) holds since $\|\Psi\mathbf{x}_{(k+1)}\|_{\infty}\leq\|\Psi\mathbf{x}_{(k+1)}\|_{1}$, Eq (A-B) holds since $\mathbf{q}^{H}\Psi\mathbf{x}_{(k+1)}\leq\|\Psi\mathbf{x}_{(k+1)}\|_{2}\|\mathbf{q}\|_{2}$, Eq (A-B) holds since $\|\Psi\mathbf{x}_{(k+1)}\|_{2}\leq\|\Psi\mathbf{x}_{(k+1)}\|_{1}$, and finally Eq (53) holds if

[TABLE]

∎

A-C Proof of Lemma III.3

Proof:

: From Lemmas III.1 and III.2, we have

[TABLE]

Since $\mathbf{x}_{(k)}^{H}\mathbf{x}_{(k)}=\mathbf{x}_{(k+1)}^{H}\mathbf{x}_{(k+1)}=L$, then

[TABLE]

Then the sequence $\{f(\mathbf{x}_{(k)})\}_{k=0}^{\infty}$ is non-increasing and since $f(\mathbf{x})\geq 0,\ \forall\ \mathbf{x}$ and hence it is bounded below, then it converges to a finite value $f^{*}$. ∎

A-D The equivalence between Eq (42) and the retraction step

In this subsection, we show the equivalence between Eq (42) and the retraction step (Step 5 in Algorithm 1). This Eq is used to write $\bar{\mathbf{x}}_{(k)}$ in terms of $\mathbf{x}_{(k+1)}$. Starting from the retraction formula, this equivalence can be shown as follows:

[TABLE]

where \bar{\Psi}=\text{diag}\Big{(}\frac{1}{|\bar{x}_{1(k)}|},\frac{1}{|\bar{x}_{2(k)}|},\ldots,\frac{1}{|\bar{x}_{L(k)}|}\Big{)}. Solving for $\bar{\mathbf{x}}_{(k)}$ yields

[TABLE]

Since $\bar{\mathbf{x}}_{(k)}\in\mathcal{T}_{\mathbf{x}_{(k)}}\mathcal{S}^{L}$ and hence $|\bar{{x}}_{l(k)}|\geq 1,\ l=1,2,\ldots,L$, the quantity $|\bar{{x}}_{l(k)}|$ can be written as $|\bar{{x}}_{l(k)}|=1+\psi_{l}$ with $\psi_{l}\geq 0\ \forall\ l$. Using this, the matrix $\bar{\Psi}$ will be \bar{\Psi}=\text{diag}\Big{(}\frac{1}{1+\psi_{1}},\frac{1}{1+\psi_{2}},\ldots,\frac{1}{1+\psi_{L}}\Big{)}, the matrix $\bar{\Psi}^{-1}$ will be

[TABLE]

where \Psi=\text{diag}\big{(}\psi_{1},\psi_{2},\ldots,\psi_{L}\big{)}. Substituting this value of $\bar{\Psi}^{-1}$ in Eq (58), the vector $\bar{\mathbf{x}}_{(k)}$ will be

[TABLE]

Eq (60) hence reduces to Eq (42).