Relaxed Bi-quadratic Optimization for Joint Filter-Signal Design in Signal-Dependent STAP

TL;DR

This paper introduces a convex relaxation approach to a non-convex joint filter-signal optimization problem in signal-dependent STAP, providing insights and potential solutions for improved radar signal processing.

Contribution

It applies biquadratic optimization and semidefinite relaxation techniques to transform a non-convex problem into a convex one, offering new analytical and numerical solution methods.

Findings

The original problem is generally non-convex.

The relaxed problem is convex and solvable.

Potential solutions are identified both analytically and numerically.

Abstract

We investigate an alternative solution method to the joint signal-beamformer optimization problem considered by Setlur and Rangaswamy[1]. First, we directly demonstrate that the problem, which minimizes the received noise, interference, and clutter power under a minimum variance distortionless response (MVDR) constraint, is generally non-convex and provide concrete insight into the nature of the nonconvexity. Second, we employ the theory of biquadratic optimization and semidefinite relaxations to produce a relaxed version of the problem, which we show to be convex. The optimality conditions of this relaxed problem are examined and a variety of potential solutions are found, both analytically and numerically.

Figures

Equations

$$f_d=2f_o\frac{\mathbf{\delta }{v}^T[\sin ({\varphi }_t)\sin ({\theta }_t)\sin ({\varphi }_t)\cos ({\theta }_t)\cos ({\varphi }_t)]}{c}$$

$${\mathbf{y}}_{\mathtt{t}}={\rho }_t{\mathbf{v}}_t(f_d)\otimes \mathbf{s}\otimes {\mathbf{a}}_t({\theta }_t,{\varphi }_t)$$

$$\tilde{\mathbf{y}}={\mathbf{y}}_{\mathtt{t}}+{\mathbf{y}}_{\mathtt{n}}+{\mathbf{y}}_{\mathtt{i}}+{\mathbf{y}}_{\mathtt{c}}={\mathbf{y}}_{\mathtt{t}}+{\mathbf{y}}_{\mathtt{u}},$$

$${\mathbf{y}}_{\mathtt{i},k}={\mathbf{\alpha }}_k\otimes {\mathbf{a}}_{\mathtt{i}}({\theta }_k,{\varphi }_k)$$

$${\mathbf{R}}_{\mathtt{i}}=\sum _{k=1}^K{\mathbf{R}}_{\mathbf{\alpha },k}\otimes {\mathbf{a}}_{\mathtt{i}}({\theta }_k,{\varphi }_k){\mathbf{a}}_i({\theta }_k,{\varphi }_k{)}^H.$$

$${\gamma }_{pq}\mathbf{v}(f_{\mathtt{c},pq})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_{pq},{\varphi }_{pq})$$

$$f_{\mathtt{c},pq}=2f_o\frac{{\dot{\mathbf{x}}}_{\text{r}}^T[\sin ({\varphi }_{pq})\sin ({\theta }_{pq})\sin ({\varphi }_{pq})\cos ({\theta }_{pq})\cos ({\varphi }_{pq})]}{c}.$$

$${\mathbf{y}}_{\mathtt{c},q}=\sum _{p=1}^P{\gamma }_{pq}\mathbf{v}(f_{\mathtt{c},pq})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_{pq},{\varphi }_{pq}).$$

$${\mathbf{B}}_{\mathbf{q}}=[\mathbf{v}(f_{C,1q})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_{1q},{\varphi }_{1q})\cdots \mathbf{v}(f_{C,Pq})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_{Pq},{\varphi }_{Pq})].$$

$${\mathbf{R}}_{\mathbf{\gamma }}^q={\mathbf{B}}_q{\mathbf{R}}_{\gamma }^{pq}{\mathbf{B}}_q^H$$

$${\mathbf{R}}_{\mathtt{c}}=\sum _{q=1}^Q{\mathbf{R}}_{\mathbf{\gamma }}^q$$

$${\mathbf{y}}_{\mathtt{c},q}={\gamma }_q\mathbf{v}(f_{\mathtt{c},q})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_q,{\varphi }_q)$$

$${\mathbf{B}}_{\mathbf{q}}=1_P^T\otimes \mathbf{v}(f_{\mathtt{c},q})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_q,{\varphi }_q).$$

$${\mathbf{R}}_{\gamma }^q$$

$$={\overline{R}}_{\gamma }^{pq}({\mathbf{v}}_q\otimes \mathbf{s}\otimes {\mathbf{a}}_q)({\mathbf{v}}_q\otimes \mathbf{s}\otimes {\mathbf{a}}_q{)}^H$$

$$(\mathbf{A}\otimes \mathbf{B}\otimes \mathbf{C})(\mathbf{D}\otimes \mathbf{E}\otimes \mathbf{F})=\mathbf{AD}\otimes \mathbf{BE}\otimes \mathbf{CF}.$$

$$\mathbf{v}(f_{\mathtt{c},q})\otimes \mathbf{s}\otimes \mathbf{a}({\theta }_q,{\varphi }_q)$$

$$={\mathbf{\Gamma }}_q\mathbf{s}$$

$$\tilde{\mathbf{y}}={\rho }_t\mathbf{Ts}+\sum _{q=1}^Q{\gamma }_q{\mathbf{\Gamma }}_q\mathbf{s}+{\mathbf{y}}_{\mathtt{n}}+{\mathbf{y}}_{\mathtt{i}}.$$

$$\begin{array}{rlrl}& \underset{\mathbf{w},\mathbf{s}}{\min }& & {\mathbf{w}}^H{\mathbf{R}}_{\mathtt{u}}(\mathbf{s})\mathbf{w}\\ & \text{s.t.}& & {\mathbf{w}}^H\mathbf{Ts}=\kappa \\ & & & {\mathbf{s}}^H\mathbf{s}\le P_o\end{array}.$$

$$\mathbf{w}$$

$$\mathbf{s}$$

$${\mathrm{f}}_o$$

$${\mathrm{f}}_{NI}$$

$$={\mathbf{b}}^H{\mathbf{\Psi }}_W^T{\mathbf{R}}_{\mathbf{ni}}{\mathbf{\Psi }}_W\mathbf{b}={\mathbf{b}}^H{\tilde{\mathbf{R}}}_{\mathbf{ni}}\mathbf{b}$$

$${\mathrm{f}}_C$$

$${\text{f}}_q(\mathbf{w},\mathbf{s})$$

$$=$$

$${\text{f}}_q(\mathbf{w},\mathbf{s})={\overline{R}}_{\gamma }^q|{\mathbf{b}}^H{\mathbf{\Psi }}_W^T{\mathbf{\Gamma }}_q{\mathbf{\Psi }}_S\mathbf{b}{|}^2={\overline{R}}_{\gamma }^q|{\mathbf{b}}^H{\overline{\mathbf{\Gamma }}}_q\mathbf{b}{|}^2$$

$${\text{f}}_q(\mathbf{w},\mathbf{s})$$

$$=$$

$$\mathrm{Re}({\mathbf{w}}^H{\mathbf{\Gamma }}_q\mathbf{s})$$

$$=\frac{1}{2}{\mathbf{b}}^H\left[\begin{array}{cc}{0}_{NML\times NML}& {\mathbf{\Gamma }}_q\\ {\mathbf{\Gamma }}_q^H& 0_{N\times N}\end{array}\right]\mathbf{b}$$

$$\mathrm{Im}({\mathbf{w}}^H{\mathbf{\Gamma }}_q\mathbf{s})$$

$$=\frac{1}{2}{\mathbf{b}}^H\left[\begin{array}{cc}{0}_{NML\times NML}& \frac{1}{j}{\mathbf{\Gamma }}_q\\ -\frac{1}{j}{\mathbf{\Gamma }}_q^H& 0_{N\times N}\end{array}\right]\mathbf{b}$$

$$\mathrm{Re}({\mathbf{w}}^H{\mathbf{\Gamma }}_q\mathbf{s})$$

Full text

I Introduction

In 1972, at a NATO conference in the East Midlands of England, two luminaries of signal processing agreed that an adaptive system was incapable of being simultaneously spatially and temporally optimal. Agreeing with M. Mermoz’s proposition, the Naval Underwater System Center’s Norman Owsley asserted that since “the spectrum of the signal must be known a priori” for temporal optimality, “the post beamformer temporal processor cannot be fully adaptive and fully optimum (sic) simultaneously”[2]. In many ways, this fundamental philosophy has carried forward in multichannel signal processing and control theory unabated in the last 40-plus years since that conference in Loughborough. Under the fully-adaptive radar paradigm, however, it is necessary to find transmit and receive resources that are simultaneously optimal, or at least optimal under prior knowledge assumptions. With these additional degrees of freedom, we can in fact develop a result that defies this wisdom. The purpose of this report is to document such an attempt, which is a significant extension of the work of [1] in joint waveform-filter design for radar space-time adaptive processing (STAP), that challenges the pre-existing status quo.

In the general scenario, we assume that an airborne radar equipped with an array of sensing elements observes a moving target on the ground. Furthermore, this radar can change its transmitted waveform every coherent processing interval (CPI), instead of on a per-pulse basis. In order to develop a strategy for waveform design, we consider a STAP model that includes fast-time (range) samples as well as the usual slow-time and spatial samples. This is a departure from traditional STAP, which operates on spatiodoppler responses from the radar after matched filtering [3, 4], but recent advances in the literature have considered incorporating fast-time data for more accurate clutter modeling, both transmitter- and jammer-induced [5, 6].

An obvious result of including the fast-time data in the model is that the clutter representation is now signal-dependent, since in airborne STAP, the dominant clutter source is non-target ground reflections that persist over all range bins. We assume that, if modeled as a random process, the clutter is uncorrelated with any interference or noise (which, unlike the clutter, are assumed to have no signal dependence), and the related clutter correlation matrix is also signal dependent. If we formulate our waveform-filter design in the typical minimum variance distortionless-response (MVDR) framework [7], this dependence in the correlation structure leads to what multiple authors ([8, 9, 10], among others) have empirically or intuitively identified as a non-convex optimization problem. However, the authors in [1] correctly identified that for a fixed transmit waveform, the problem is convex in the receive filter, and vice versa. This led to a collection of design algorithms based on alternating minimization, a reasonable heuristic.

While the alternating minimization method is useful, it has no claims of optimality, nor do other methods in the literature that dealt with signal-dependent interference in other contexts, like the single sensor radar and reverberant channel of [11]. Furthermore, none of these authors directly proved or cited anything that demonstrated the non-convexity of the problem or its engineering consequences. We will ameliorate some of these concerns in this paper, first by directly proving the non-convexity of the objective, then relating it to the existing literature on biquadratic programming (BQP). While the biquadratic program is demonstrably non-convex, it is possible to relax it to a convex quadratic program using semidefinite programming (SDP) (see [12] for details on SDP and [13, 14] on relaxing the BQP). This relaxation will permit us to efficiently solve the problem computationally and analytically (up to a matrix completion), as well as reveal important structural information about the solution that matches with engineering intuition.

Before continuing, we outline some mathematical notation that will be used throughout the paper. The symbol $\mathbf{x}\in\mathbb{R}^{N\times 1}(\mathbf{x}\in\mathbb{C}^{N\times 1})$ indicates a column vector of real (complex) values, while $\mathbf{X}\in\mathbb{R}^{M\times N}(\mathbf{X}\in\mathbb{C}^{M\times N})$ indicates an $N$ row, $M$ column matrix of real (complex) values. The superscripts $T$, $*$, and $H$ indicate the transpose, conjugate, and Hermitian (conjugate) transpose of a matrix or vector. The symbol $\>\otimes\>$ indicates a Kronecker product. The operator $\operatorname{vec}$ turns a matrix into a vector with the matrix stacked columnwise – that is, for $\mathbf{A}\in\mathbb{R}^{p\times q},\operatorname{vec}(\mathbf{A})\in\mathbb{R}^{pq\times 1}$. Special matrices that will recur often include:

• •

the $n\times n$ identity matrix $\mathbf{I}_{n}$;

• •

the $n\times 1$ vector of all ones $\mathbf{1}_{n}$;

• •

the $n\times m$ all-zero matrix $\mathbf{0}_{n\times m}$;

• •

and, via Magnus & Neudecker[15], the $(p,q)$ commutation matrix $\mathbf{K}_{p,q}\in\mathbb{R}^{pq\times pq}$ that is non-uniquely determined by the relation $\operatorname{vec}(\mathbf{A}^{T})=\mathbf{K}_{p,q}\operatorname{vec}(\mathbf{A})$ for $\mathbf{A}$ as above. $\mathbf{K}_{p,q}$ is an orthogonal permutation matrix (i.e., $\mathbf{K}_{p,q}^{T}=\mathbf{K}_{p,q}^{-1}$ ). Transposing the matrix swaps the indices as well ($\mathbf{K}_{p,q}^{T}=\mathbf{K}_{q,p}$), and clearly then $\mathbf{K}_{p,q}\mathbf{K}_{q,p}=\mathbf{I}_{pq}$. Additionally, $\mathbf{K}_{1,q}=\mathbf{I}_{q}$ and $\mathbf{K}_{p,1}=\mathbf{I}_{p}$. We will also frequently use the one-index form $\mathbf{K}_{p}=\mathbf{K}_{p,p}$ for the square commutation matrix.

We will often use calligraphic letters like $\mathcal{C}$ to indicate tensors or multilinear operators (like the Hessian matrix). Additionally, the vectorizations of various matrices will be indicated by the equivalent lowercase Greek symbol in bold – for example, a matrix $\mathbf{B}$ will have the vectorization $\boldsymbol{\beta}=\operatorname{vec}(\mathbf{B})$.

The rest of the paper continues as follows: First, we will outline the general signal model and relate it to the work on transfer functions by [16]. Next, we will define the joint design problem and demonstrate, after some transformation, that the problem is non-convex unless all clutter patches are known to be nulled a priori. Section IV will describe the relationship with the biquadratic problem and methods to relax the joint design into a convex quadratic semidefinite program. We then attempt to find analytic solutions and insights into this relaxed problem in Section V. These solutions are generally verified in Section VI through simulation, as well as demonstrate that the resulting rank-one solution performs similarly to the alternating minimization. Finally, we summarize our conclusions and highlight future research paths in Section VII.

II STAP Model

Consider a radar that consists of a calibrated airborne linear array of $M$ identical sensor elements, where the first element is the phase center and also acts as the transmitter. During the transmission period, the radar probes the environment with a number of pulses $s(t)$ of width $T$ seconds and bandwidth $B$ hertz at a carrier frequency $f_{o}$. Within each burst, we transmit $L$ pulses at a rate of $f_{p}$ (i.e. the pulses are transmitted every $T_{p}=1/f_{p}$ seconds) and collect them in a coherent processing interval (CPI). We also assume that the phase center is located at $\mathbf{x}_{r}$ and the platform moves at a rate that is relatively constant over the CPI.

Assume the probed environment contains a target that lies at an azimuth $\theta_{t}$ and an elevation $\phi_{t}$ relative to the array phase center, moving with a relative velocity vector $\boldsymbol{\delta}v=\begin{bmatrix}{\delta v}_{x}&{\delta v}_{y}&{\delta v}_{z}\end{bmatrix}^{T}.$ If we assume that the array interelement spacing $d$ is small relative to the distance between the platform and the target, then the target’s doppler shift is independent of the element index and given by

[TABLE]

where $c$ is the usual speed of light.

Let us now assume that we discretely sample the pulse $s(t)$ into $N$ samples, resulting in the sample vector $\mathbf{s}=[s_{1},s_{2},\cdots,s_{N}]\in\mathbb{C}^{N}$. Assuming the data is aligned to a common reference and given other assumptions from [27], we can say that at the target range gate $\tau_{t}$, the combined target response over the entire CPI can be represented by a vector $\mathbf{y}_{t}\in\mathbb{C}^{NML}$ given by

[TABLE]

where $\rho_{t}$ is the complex backscattering coefficent from the target, the vector $\mathbf{v}_{t}(f_{d})\in\mathbb{C}^{L}$ is the doppler steering vector whose $i$th element is given by $e^{-\jmath 2\pi f_{d}(i-1)T_{p}}$, and the vector $\mathbf{a}_{t}(\theta_{t},\phi_{t})\in\mathbb{C}^{M}$ is the spatial steering vector whose $i$th element is given by $e^{-\jmath 2\pi(i-1)\vartheta}$ where $\vartheta=d\sin(\theta_{t})\sin(\phi_{t})f_{o}/c$. In this form, our departure from traditional STAP is clear, given the dependence on the waveform $\mathbf{s}$.

Since no radar operates in an ideal environment, the target return $\mathbf{y}_{t}$ is corrupted by a variety of undesired returns from the environments – noise, signal-independent interference, and clutter. We can consider the overall return in the target range gate as an additive model:

[TABLE]

where the subscripts n, i, and c stand for noise, interference, and clutter, respectively, and these are all collected in the Undesired signal term. We assume that these undesired energy sources are statistically uncorrelated from each other, with unique distributions. We shall subsequently describe each of these sources, starting with the noise.

We assume the noise is zero mean and identically distributed across all sensors, pulses, and fast time samples. The covariance matrix of $\mathbf{y}_{n}$ is given by $\mathbf{R}_{\mathtt{n}}\in\mathbb{C}^{NML\times NML}$.

The interference term consists of jammers and other spurious emitters that are intentional or unintentional and are ground-based, airborne, or both. We assume that there are $K$ known interference sources, but otherwise we have no knowledge of what they transmit into the surveillance region, and (we hope) they are not dependent on our transmitted signal. Thus, we model their contributions as a zero-mean random process spread over all pulses and fast time samples. Assume that the $k$th interferer is located at the azimuth-elevation pair $(\theta_{k},\phi_{k})$. In the $l$th PRI, we assume that the waveform is a complex continuous-time signal $\alpha_{kl}(t)$. Under the same sampling scheme, this transforms into a vector $\boldsymbol{\alpha}_{kl}\in\mathbb{C}^{N}$, similar in form to $\mathbf{s}$. Stacked across all PRIs, we obtain another random vector $\boldsymbol{\alpha}_{k}=[\boldsymbol{\alpha}_{k0}^{T}\,\boldsymbol{\alpha}_{k1}^{T}\,\cdots\,\boldsymbol{\alpha}_{k(L-1)}^{T}]^{T}\in\mathbb{C}^{NL}$, whose covariance matrix we define as $\mathop{\bf E\/}\{\boldsymbol{\alpha}_{k}\boldsymbol{\alpha}_{k}^{H}\}=\mathbf{R}_{\boldsymbol{\alpha},k}$. Then, the response from the $k$th interferer can be modeled as

[TABLE]

where, as with the target, $\mathbf{a}_{i}(\theta_{k},\phi_{k})$ is the array response to the interferer. The covariance of $\mathbf{y}$ is then \mathbf{R}_{\boldsymbol{\alpha},k}\raisebox{1.0pt}{>\otimes>}\mathbf{a}_{\mathtt{i}}(\theta_{k},\phi_{k})\mathbf{a}_{i}(\theta_{k},\phi_{k})^{H}, since expectations follow through Kronecker products that don’t have a random dependence. If we assume that each of the interferers is uncorrelated with each other, then the overall covariance of the combined signal-independent interference $\mathbf{y}_{\mathtt{i}}$ is

[TABLE]

In future sections, we will use the notation $\mathbf{R_{ni}}=\mathbf{R}_{\mathtt{n}}+\mathbf{R}_{\mathtt{i}}$ to denote the combined noise and interference covariance matrix, describing the second order effects of the signal-independent corruption.

Finally, we come to the clutter. In airborne radar applications, the most significant clutter source is the ground, which produces returns persistent throughout all range gates up to the horizon. Though other clutter sources exist, like large discrete objects, vegetation, and targets not currently being surveilled, the specific stochastic model we apply here only concerns ground clutter. However, we note that a later formulation in this paper may be amenable to considering those other sources in a manner that recalls efforts in the literature on channel estimation (see, for example, [16]). For now, let us assume we have a number of clutter patches (say, $Q$) each comprising $P$ distinct scatterers. As with the target, the return from the $p$th scatterer in the $q$th patch, located spatially at the azimuth-elevation pair $(\theta_{pq},\phi_{pq})$, maintains a Kronecker structure given by

[TABLE]

where the returned complex reflectivity is $\gamma_{pq}$ and the Doppler shift observed by the platform is $f_{\mathtt{c},pq}$. This Doppler shift is solely induced by the platform motion, characterized by the aforementioned velocity vector $\dot{\mathbf{x}}_{\textmd{r}}$, and is given by

[TABLE]

Thus, the overall response from the $q$th clutter patch is

[TABLE]

In order to define the covariance matrix of this response, we require the $q$th combining matrix $\mathbf{B}_{q}\in\mathbb{C}^{NML\times P}$ as

[TABLE]

and the covariance of the reflectivity vector $[\gamma_{1q}\,\gamma_{2q}\,\cdots\,\gamma_{Pq}]^{T}$ given by $\mathbf{R}^{pq}_{\gamma}\in\mathbb{C}^{P\times P}$. Then, the overall covariance matrix for the patch is

[TABLE]

If we assume that the scatterers in one patch are uncorrelated with the scatterers in any other patch, then the total clutter response is $\mathbf{y}_{\mathtt{c}}=\sum_{q=1}^{Q}\mathbf{y}_{\mathtt{c},q}$ and its covariance is given by

[TABLE]

As above, we will denote the overall undesired response covariance matrix as $\mathbf{R}_{\mathtt{u}}=\mathbf{R_{ni}}+\mathbf{R}_{\mathtt{c}}(\mathbf{s})$.

Since this is a rather cumbersome model, we can simplify our description of the clutter as follows. Let us assume that our range resolution is large enough that we cannot resolve individual scatterers in each patch – as mentioned in [1], this is typical in STAP applications. Thus, we can regard each scatterer in the patch as lying within the same range gate and having approximately equal Doppler shifts, hence $f_{\mathtt{c},pq}\approx f_{\mathtt{c},q}$. Similarly, if we assume far-field operation, the scatterers will lie in approximately the same angular resolution cell centered at $(\theta_{q},\phi_{q})$, which means $\theta_{pq}\approx\theta_{q},\phi_{pq}\approx\phi_{q}$. Given this simplification, we can modify our representation of the patch response and its covariance.

Under this assumption, the per-patch clutter response is

[TABLE]

where $\gamma_{q}=\sum_{p=1}^{P}\gamma_{pq}$ is the combined reflectivity of all scatterers within the patch. For the covariance, the combining matrix for the $q$-th clutter patch is given by $\mathbf{B_{q}}=[\mathbf{v}(f_{\mathtt{c},q})\otimes\mathbf{s}\otimes\mathbf{a}(\theta_{q},\phi_{q}),\ldots,\mathbf{v}(f_{\mathtt{c},q})\otimes\mathbf{s}\otimes\mathbf{a}(\theta_{q},\phi_{q})]$. Since the deterministic patch response is repeated $P$ times, via standard Kronecker product properties, this is equivalent to

[TABLE]

More importantly, since $\mathbf{R}^{q}_{\boldsymbol{\gamma}}=\mathbf{B}_{q}\mathbf{R}^{pq}_{\gamma}\mathbf{B}^{H}_{q}$, we have

[TABLE]

where $\overline{R}^{q}_{\gamma}=\mathbf{1}_{P}^{T}\mathbf{R}^{pq}_{\gamma}\mathbf{1}_{P}$.

Relationship to the Channel Model construct

In [16], the authors considered a transfer function/matrix approach for simultaneous transmit and receive resource design in MIMO radar, similar to the typical literature on control theory and digital communications. Using a rather obvious mathematical fact, we can immediately reframe our model in this form. Observe the following Kronecker mixed product property: for conformable matrices/vectors/scalars $\mathbf{A}$ through $\mathbf{F}$,

[TABLE]

Clearly, we can apply this to a deterministic space-time-doppler response vector. For example, the per-patch clutter response under the simplification is \mathbf{v}(f_{\mathtt{c},q})\raisebox{1.0pt}{>\otimes>}\mathbf{s}\raisebox{1.0pt}{>\otimes>}\mathbf{a}(\theta_{q},\phi_{q}). Since we can regard the scalar 1 as a conformable matrix, we have

[TABLE]

where $\boldsymbol{\Gamma}_{q}$ is obviously defined. By the same token, the target response is equivalent to $\mathbf{T}\mathbf{s}$ where \mathbf{T}=\mathbf{v}_{t}(f_{d})\raisebox{1.0pt}{>\otimes>}\mathbf{I}_{N}\raisebox{1.0pt}{>\otimes>}\mathbf{a}_{t}(\theta_{t},\phi_{t}). With this form, the overall received vector in the range gate of interest is

[TABLE]

This is something considerably easier on the eyes and more comprehensibly relates to the transfer function approach. We will continue to use this notation in subsequent analysis, as it reveals the structure of the design problem in a much more direct fashion.

III Joint Waveform-Filter Design

With the signal model in place, we now turn to the initial purpose of [1]. Our goal is to find a STAP beamformer vector $\mathbf{w}\in\mathbb{C}^{NML}$ and a transmit signal $\mathbf{s}\in\mathbb{C}^{N}$ that minimizes the combined effect of the noise and interference represented by the covariance matrix $\mathbf{R_{ni}}\in\mathbb{C}^{NML\times NML}$) and the signal-dependent clutter. However, we also want to ensure reasonable radar operation. We describe this process below.

At the range gate we interrogate (which we assume contains the target), the return $\tilde{\mathbf{y}}$ is processed by a filter characterized by the weight vector $\mathbf{w}$, forming the output return $\mathbf{w}^{H}\tilde{\mathbf{y}}$. As mentioned above, we want to design this vector and the signal to minimize the expected undesired power $\mathop{\bf E\/}\{|\mathbf{w}^{H}\mathbf{y}_{\mathtt{u}}|^{2}\}=\mathbf{w}^{H}\mathbf{R}_{\mathtt{u}}(\mathbf{s})\mathbf{w}$. Additionally, we would like to constrain this minimization to ensure reasonable radar operation. First, for a given target space-time-doppler bin, we want a particular filter output, say, $\kappa\in\mathbb{C}$. This filter output can be represented by the Capon beamformer equation $\mathbf{w}^{H}\mathbf{T}\mathbf{s}$, where $\mathbf{T}\in\mathbb{C}^{NML\times N}$ is the target channel response above. Second, we place an upper bound on the total signal power, say, $P_{o}$. Mathematically, we can represent this optimization problem as

[TABLE]

In [1], this problem was computationally shown to be non-convex, and therefore used a variety of alternating minimization algorithms to solve it. Our purpose will first be to analytically prove that the problem as formulated is non-convex from an engineering standpoint. However,determining convexity, or lack thereof, will be difficult in the current configuration, so we will make a notational change. Let $\mathbf{b}$ be the combined beamformer-signal vector defined by $\mathbf{b}=[\mathbf{w}^{T}\;\mathbf{s}^{T}]^{T}$.

In order to recover the individual elements $\mathbf{w},\mathbf{s}$ from $\mathbf{b}$, we define two matrices:

[TABLE]

For further notational simplicity, let us also define the complete objective function $f_{o}$ as a sum of the noise-interference cost $f_{NI}$ & the total clutter cost $f_{C}$ (which is itself a sum of the $Q$ per-patch clutter costs):

[TABLE]

III-A Massaging the clutter objective

We start with the per-patch clutter objective function. Using the channel representation above, $\mathbf{R}^{q}_{\boldsymbol{\gamma}}$ can be distilled to $\mathbf{R}^{q}_{\boldsymbol{\gamma}}=\overline{R}^{q}_{\gamma}\boldsymbol{\Gamma}_{q}\mathbf{s}\mathbf{s}^{H}\boldsymbol{\Gamma}_{q}^{H}$. Therefore, the cost functional for the $q$-th clutter patch is

[TABLE]

Substituting $\mathbf{w}=\boldsymbol{\Psi}_{W}\mathbf{b}$ and $\mathbf{s}=\boldsymbol{\Psi}_{S}\mathbf{b}$ into (14b) above, the final form of this cost function in the joint vector is

[TABLE]

where $\overline{\boldsymbol{\Gamma}}_{q}=\boldsymbol{\Psi}_{W}^{T}\boldsymbol{\Gamma}_{q}\boldsymbol{\Psi}_{S}$.

III-B An alternative form of the per-patch clutter cost

Although the form derived above is quite compact & informative, we develop another equivalent formation of the per-patch clutter cost that will allow us to quickly identify the Hessian of each of the cost functions considered and, ultimately, determine the convexity of the problem.

An alternative form of the per-patch clutter cost function comes from Equation (14b). First, for any complex number $z$, its squared magnitude is $|z|^{2}=\operatorname{Re}(z)^{2}+\operatorname{Im}(z)^{2}$. Applying this to Equation (14b), we obtain:

[TABLE]

Using elementary algebra, the components of Equation (16b) are given by:

[TABLE]

Before we continue, we require a particular property of complex matrices. Any complex matrix $\mathbf{M}$ can be decomposed into the matrix sum $\mathbf{M}=\mathbf{M}_{H}+\mathbf{M}_{AH}$ where $\mathbf{M}_{H}=\frac{1}{2}(\mathbf{M}+\mathbf{M}^{H}),\mathbf{M}_{AH}=\frac{1}{2}(\mathbf{M}-\mathbf{M}^{H})$ are the Hermitian and Anti-Hermitian parts of $\mathbf{M}$. We can see by inspection that the block matrices in the real & imaginary parts listed above contain the Hermitian and anti-Hermitian parts of $\overline{\boldsymbol{\Gamma}}_{q}$. Therefore, these forms are

[TABLE]

Substituting these into the forms above and simplifying results in

[TABLE]

This is an even more compact form of $f_{q}$ that, as we will see in the next section, permits us to quickly and easily derive the Hessian of the per-patch clutter cost, the total clutter cost, and the overall objective function.

III-C Summary of Forms

Finally, we arrive at multiple equivalent final forms of the clutter function $f_{C}$, which is a sum of the per-patch clutter forms:

[TABLE]

In the joint variable $\mathbf{b}$, the overall optimization problem becomes

[TABLE]

Before we prove joint non-convexity, the form in Equation (14a) provides an immediate proof of convexity of the per-patch cost in $\mathbf{w}$ for fixed $\mathbf{s}$ and vice versa. Since the cost function $f$ is real and the matrix $\mathbf{R}^{q}_{\boldsymbol{\gamma}}$ is positive-semidefinite (by virtue of being a correlation/covariance matrix) for fixed $\mathbf{s}$, then for any non-zero $\mathbf{w}$, $f\geq 0$ and is thus convex in $\mathbf{w}$.

III-D Proving Joint (Non-)Convexity

III-D1 Methods of verifying convexity

The traditional definition of convexity is technically only defined for real-valued functions of real arguments (be they scalar, vector, or matrix). However, newly developed theory in [17] and the associated journal literature permits us to make similar claims for real-valued functions of complex arguments.

Let us define a scalar function $f$ of the complex variable $\mathbf{b}$ as

[TABLE]

Under the traditional definition (see, for example, in [12]), this function is convex if for any $\mathbf{b}_{1},\mathbf{b}_{2}\in\Omega$ (where $\Omega$ is some convex subset of $\mathbb{C}^{J\times 1}$) and $\alpha\in\left[0,1\right]\subset\mathbb{R}$,

[TABLE]

This is true because the function of interest is real-valued, despite the fact that its arguments are complex-valued. Equivalently, the Hessian matrix $\mathcal{\widetilde{H}}_{\mathbf{v}}[f]=\nabla^{2}_{\mathbf{v}}f(\mathbf{v})$ must be positive semidefinite when evaluated at any stationary point $\mathbf{v}=\mathbf{b}_{0}\in\Omega$. We will use the Hessian method, because it reveals the interesting structure of the problem in addition to showing non-convexity.

III-D2 Determining complex Hessians

Since our function is a real-valued function of complex-valued variables, we require the following from Hjørungnes’ work on complex matrix derivatives [17]. From [17, Theorem 3.2], the stationary points of a real-valued function of complex variables are the points where

[TABLE]

where, for example, the form $\mathcal{D}_{\mathbf{b}}f(\mathbf{b},\mathbf{b}^{*})$ is the complex derivative with respect to $\mathbf{b}$. As is the standard, we treat the vector $\mathbf{b}$ and its conjugate $\mathbf{b}^{*}$ as separate variables for the purposes of differentiation. This fact derives from the Wirtinger calculus and is directly proven in the reference above.

A Taylor series argument for Hessians of real-valued functions

It is well known that for a twice-differentiable real-valued scalar function, positive semi-definiteness of its Hessian matrix over a convex set is an equivalent statement of convexity. This is true even if the arguments are complex-valued. See, for example, [17, Lemma 5.2] and the consequences thereof. Namely, if the Hessian matrix in [17, Equation 5.44] is positive semidefinite at all stationary points $\mathbf{b}_{o}$ within the convex set of interest, then the zeroth-order convexity requirement in Equation 19 is automatically satisfied if we set $\mathbf{b}_{1}=\mathbf{b}_{o}+\mathop{}\!\mathrm{d}\mathbf{b},\mathbf{b}_{2}=\mathbf{b}_{o}$. Let us rearrange Equation 19 to directly prove this assertion from the reference. First, observe that

[TABLE]

Then, the convexity requirement becomes

[TABLE]

Let us set $\mathbf{b}_{1}=\mathbf{b}_{o}+\mathop{}\!\mathrm{d}\mathbf{b},\mathbf{b}_{2}=\mathbf{b}_{o}$, where $\mathbf{b}_{o}$ is a stationary point and $\mathop{}\!\mathrm{d}\mathbf{b}$ is a vector of infintesimally small perturbations. Then, the convexity requirement is

[TABLE]

At a stationary point, the first term on the left hand side is zero. Why? We can clearly see that, as $\alpha\rightarrow 0$, it is essentially a derivative, which must vanish at a stationary point by definition. Thus, the condition again changes to:

[TABLE]

Using the Taylor series expansion from [17, Lemma 5.2], we expand $f(\mathbf{b}+\mathop{}\!\mathrm{d}\mathbf{b})$,

[TABLE]

where the higher order terms function $HOT$ converges to zero in the sense provided in the reference. Certainly then, regardless of the value of $\mathbf{b}$,

[TABLE]

If we evaluate this expression at $\mathbf{b}=\mathbf{b}_{o}$, which is our aforementioned stationary point, the first order derivatives in the second term vanish by definition and therefore

[TABLE]

In order for (22) to be satisfied, we then have

[TABLE]

which is clearly the definition of positive semidefiniteness for a Hermitian complex matrix. Therefore, in order for the function $f$ to be convex, the matrix $\begin{bmatrix}\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}&\mathcal{H}_{\mathbf{b}^{*}\mathbf{b}^{*}}\\ \mathcal{H}_{\mathbf{b}\mathbf{b}}&\mathcal{H}_{\mathbf{b}^{*}\mathbf{b}}\end{bmatrix}$ must be positive semidefinite at all possible $\mathbf{b}_{o}$ inside the convex domain of $f$.

Remark 1**.**

Traditionally, convexity is only defined for real functions of real variables, which may make the analysis above seem incomplete. However, as we will show below, the definition of convexity for a complex vector variable is equivalent to the definition of joint convexity for its real and imaginary parts. To show this, we will prove that the definiteness condition remains if one replaces the differentials with respect to the complex variables with those for the real variables.

First, let us define the real and imaginary parts of the complex vector as $\mathbf{x}=\operatorname{Re}(\mathbf{b}),\mathbf{y}=\operatorname{Im}(\mathbf{b})$ (which we assume are linearly independent real vectors and will treat as such). Clearly, then, we have $\mathbf{b}=\mathbf{x}+j\mathbf{y}$ and $\mathbf{b}^{*}=\mathbf{x}-j\mathbf{y}$. Applying the differential rules of [17] to these expressions, we obtain:

[TABLE]

which clearly means that we can replace the vector of complex differentials we see throughout the above derivation with

[TABLE]

where $\mathbf{T}_{RC}$ is a transformation matrix that maps the real components to the complex vectors. Making this substitution in the Taylor series expansion from above results in

[TABLE]

where ${HOT}(\mathbf{x},\mathbf{y})$ must still be positive & converge to zero in the same sense as before. As an aside, we note that the first derivatives are clearly those w.r.t. the components $\mathbf{x},\mathbf{y}$.

Continuing on, it is still true that the following inequality holds, even if we substitute for the real and imaginary components:

[TABLE]

Evaluating the inequality at the new stationary points $\mathbf{x}_{o},\mathbf{y}_{o}$ eliminates the first derivatives, yielding

[TABLE]

Again, (22) is satisfied only if the right-hand side of the above expression is greater than the value of the function at the stationary point, ergo

[TABLE]

Now, the convexity requirement is that $\mathbf{T}_{RC}^{H}\begin{bmatrix}\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}&\mathcal{H}_{\mathbf{b}^{*}\mathbf{b}^{*}}\\ \mathcal{H}_{\mathbf{b}\mathbf{b}}&\mathcal{H}_{\mathbf{b}^{*}\mathbf{b}}\end{bmatrix}\mathbf{T}_{RC}$ must be positive semidefinite at all possible $\mathbf{x}_{o},\mathbf{y}_{o}$ inside the convex domain of $f$. However, this is exactly identical to the complex condition above, since according to [18, Lemma 7.63(b)], these are biconditional statements! For real-valued functions, therefore, convexity in the complex variable is identical to joint convexity in its real and imaginary parts.

For notational simplicity, let us define this matrix – which we will call the Hessian of the function $f$ with respect to the Taylor series argument above – as:

[TABLE]

The augmented Hessian

While the Hessian $\mathcal{\widetilde{H}}[f]$ is useful for determining convexity, it can be quite difficult to find expressions in a simple and systematic way for a given function $f$. If we use a change of variables, however, we can quickly find a related Hessian using a chain rule. As in [17], we first define the augmented variable $\mathcal{Z}=\begin{bmatrix}\mathbf{b}&\mathbf{b}^{*}\end{bmatrix}$. Derivatives taken with respect to $\mathcal{Z}$ are often just stackings of the derivatives with respect to $\mathbf{b}$ and $\mathbf{b}^{*}$. The matrix of second derivatives with respect to the augmented variable $\mathcal{H}_{\mathcal{Z}\mathcal{Z}}(f)$, henceforth called the augmented Hessian of the function $f$, is given by:

[TABLE]

The Chain Rule for augmented Hessians

A benefit of the augmented Hessian is that a true chain rule in complex arguments exists for them, as opposed to the desired Hessian form. While there is a more general form of the chain rule listed in the reference [17], we present below a stripped down version more appropriate for our purposes.

If $f=h\circ\mathbf{g}$ is a scalar function of a $J\times 1$ vector $\mathbf{b}$ (a composition of the scalar function $h$ and the vector function $\mathbf{g}$), then the augmented Hessian of $f$ is given by

[TABLE]

As usual, the operator $\mathcal{D}_{\mathbf{Z}}(\mathbf{F})=\dfrac{\partial\operatorname{vec}(\mathbf{F})}{\partial\operatorname{vec}^{T}(\mathbf{Z})}$ can be considered the derivative of the matrix function $\mathbf{F}$ with respect to the complex matrix variable $\mathbf{Z}$, and similarly for vector & scalar functions.

Mapping between augmented Hessian and the actual Hessian

While the augmented Hessian is excellent for calculation purposes, it is the actual Hessian that is needed for determining convexity. Thankfully, there is a straightforward relationship between the two:

[TABLE]

and vice versa. Functionally, this means that we exchange the (1,1) block with the (2,1) block and the (1,2) block with the (2,2) block.

III-E STAP Step 0: Preliminaries

III-E1 Procedural outline

Given the various functional forms and preliminary proofs, we will now outline the rest of the procedure to determine convexity of the objective. First, we will determine the derivatives and Hessians of each component of the per-patch clutter cost $\textsf{f}_{q}$ and use the chain rule to establish its Hessian. Using elementary techniques, we will demonstrate that this matrix is always indefinite unless a certain condition is satisfied, in which case it is always zero. Next, we will use the additivity of the combined clutter cost function to compute its Hessian and prove a similar definiteness result. Finally, we will compute the Hessian of the overall objective and show its indefiniteness by construction.

III-E2 Functions of interest

Recall that in (17), we defined the per-patch clutter function is $f_{q}=\mathbf{g}^{H}_{q}(\mathbf{b})\mathbf{M}_{q}\mathbf{g}_{q}(\mathbf{b})$, where

[TABLE]

This form of the per-patch clutter function is quite compact and simplifies finding its augmented and actual Hessians, and therefore the augmented and actual Hessians of the total clutter and overall objective functions. Additionally, from this point forward, let $J=N(ML+1)$ denote the joint number of parameters (and thus the length of the combined beamformer $\mathbf{b}$).

III-F STAP Step 1: The differentials of $\mathbf{g}_{q}$

III-F1 First derivative

Using the rules of complex matrix derivatives/differentials outlined above and in Hjørungnes [17], we can find the first derivatives of the inner vector function to be

[TABLE]

III-F2 Hessian

Let us first observe that $\mathbf{g}_{q}$ is a stacking of scalar functions of $\mathbf{b}$:

[TABLE]

Using [17, Eq. 5.69, 5.70] , we see that, as a result of the above form, the augmented Hessian of $\mathbf{g}_{q}$ is just a stacking of the per-element augmented Hessians:

[TABLE]

In order to find the per-element augmented Hessians, we use the following fact. From [17, Example 5.2], for any compatible matrix $\boldsymbol{\Phi}$,

[TABLE]

Continuing with this in mind, we find the per-element augmented Hessians to be

[TABLE]

Therefore, the augmented Hessian of the vector function $\mathbf{g}_{q}$ is

[TABLE]

III-G STAP Step 2: The differentials of $h$ w.r.t. $\mathbf{g}_{q}$

III-G1 First derivative

Recall $h=\mathbf{g}^{H}_{q}\mathbf{M}_{q}\mathbf{g}_{q}$. Thus,

[TABLE]

This means that

[TABLE]

III-G2 Hessian

From [17, Example 5.2] , we know that

[TABLE]

for any compatible matrix $\boldsymbol{\Phi}$. That means that, by inspection, $\mathcal{H}_{\mathbf{g}\mathbf{g}}(h)=\mathbf{0}_{J\times J}$!

III-H STAP Step 3: The per-patch clutter Hessian

Applying the chain-rule for the augmented Hessian and the derivatives found above, the augmented Hessian of the per-patch clutter cost function is

[TABLE]

In order to reduce this to a simpler form, we need to know what the scalars in each term are:

[TABLE]

Thus, the per-patch augmented Hessian reduces to

[TABLE]

Based on the form, we observe that the only non-zero unique block of this matrix is the (2,1) block (since the (1,2) block is its transpose). What is the (2,1) block? First, recall that

[TABLE]

Temporarily omitting the $\overline{R}^{q}_{\gamma}$ term, the (2,1) block is

[TABLE]

Therefore, the final form of the augmented per-patch Hessian is

[TABLE]

Using the second fact in Section 1, the final form of the true Hessian is

[TABLE]

where we have defined $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{q})=\overline{R}^{q}_{\gamma}\left((\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\overline{\boldsymbol{\Gamma}}_{q}\right)_{H}$.

III-I STAP Step 4: Definiteness of $\mathcal{\widetilde{H}}[\textsf{f}_{q}]$

In order for $\textsf{f}_{q}$ to be convex in $\mathbf{b}$ , then $\mathcal{\widetilde{H}}[\textsf{f}_{q}]$ must be positive-semidefinite for all valid $\mathbf{b}_{o}$. Since this Hessian is block diagonal, it is PSD if and only if each of the diagonal blocks are PSD. Furthermore, since the (2,2) block is the transpose of the (1,1) block, the (2,2) block’s definiteness is the same as the (1,1) block’s definiteness. Therefore, $\mathcal{\widetilde{H}}[\textsf{f}_{q}]$ is PSD iff $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{q})=\overline{R}^{q}_{\gamma}\left((\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\overline{\boldsymbol{\Gamma}}_{q}\right)_{H}$ is PSD. We will show definitively that $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{q})$ cannot be PSD solely by its form.

First, $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{q})$ is given by

[TABLE]

It is well known (see [18, Lemma 6.32, p. 124]) that a complex matrix with this form is indefinite. Why? Due to the structure, its eigenvalues are:

• •

the $N$ singular values of $\frac{\overline{R}^{q}_{\gamma}}{2}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\boldsymbol{\Gamma}_{q}$ (the upper right corner matrix), which are positive;

• •

the negatives of the singular values of $\frac{\overline{R}^{q}_{\gamma}}{2}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\boldsymbol{\Gamma}_{q}$, which are negative;

• •

and zeroes for the rest.

What are these singular values of $\frac{\overline{R}^{q}_{\gamma}}{2}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\boldsymbol{\Gamma}_{q}$? First,

[TABLE]

It is obvious that the eigenvalues of this matrix are $\left(\frac{\overline{R}^{q}_{\gamma}}{2}\right)^{2}|\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}|^{2}{\|\mathbf{c}_{q}\|^{2}}$ with multiplicity $N$. The singular values of $\frac{\overline{R}^{q}_{\gamma}}{2}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\boldsymbol{\Gamma}_{q}$ are then the positive square roots of this eigenvalue, or $\frac{\overline{R}^{q}_{\gamma}}{2}|\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}|\|\mathbf{c}_{q}\|$.

If a matrix has negative eigenvalues, it cannot be PSD; similarly, if it has positive eigenvalues, it cannot be negative semidefinite. The only way this matrix has no negative eigenvalues is if they are all zero. Since two of the terms in the singular value expression are always positive, this eigenvalue is zero only if the optimal signal/beamformer pair nulls the $q$th clutter patch. Therefore, the per-patch clutter Hessian is indefinite and the per-patch clutter cost is nonconvex. This is true for any choice of signal or beamformer on any choice of set.

III-J STAP Step 5: Definiteness of $\mathcal{\widetilde{H}}[\textsf{f}_{C}]$

Since the overall clutter cost is a sum of the per-patch costs, its Hessian is the sum of the per-patch clutter Hessians. In other words,

[TABLE]

Similar to the per-patch clutter Hessian, the definiteness of the overall clutter Hessian is determined by the definiteness of $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{C})$. Therefore, let us take a closer look at $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{C})$:

[TABLE]

Define the $NML\times N$ complex matrix $\boldsymbol{\Pi}_{Q}=\sum_{q=1}^{Q}\frac{\overline{R}^{q}_{\gamma}}{2}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s})^{*}\boldsymbol{\Gamma}_{q}$. Hence,

[TABLE]

Once again, the previous structure occurs. By the lemma, the eigenvalues of $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{C})$ are the $N$ singular values of $\boldsymbol{\Pi}_{Q}$ (i.e., the square roots of the eigenvalues of $\boldsymbol{\Pi}_{Q}^{H}\boldsymbol{\Pi}_{Q}$), the negatives of those singular values, and $N(ML-1)$ zeros.

What are these singular values? First, let us examine $\boldsymbol{\Pi}_{Q}^{H}\boldsymbol{\Pi}_{Q}$:

[TABLE]

To simplify things slightly, we can find a different form of $\boldsymbol{\Gamma}_{q}^{H}\boldsymbol{\Gamma}_{l}$:

[TABLE]

Returning to the original matrix $\boldsymbol{\Pi}_{Q}^{H}\boldsymbol{\Pi}_{Q}$, we find

[TABLE]

If we define the matrix $\mathbf{D}=\begin{bmatrix}(\mathbf{w}^{H}\boldsymbol{\Gamma}_{1}\mathbf{s})^{*}\mathbf{c}_{1}\cdots(\mathbf{w}^{H}\boldsymbol{\Gamma}_{Q}\mathbf{s})^{*}\mathbf{c}_{Q}\end{bmatrix}$, then the scalar inside the double sum above is the $(q,l)$th element of the matrix $\mathbf{D}^{H}\mathbf{D}$, which is a Gramian matrix and is thus positive semidefinite. The double sum, then, is the scalar $\mathbf{1}_{Q}^{T}\mathbf{D}^{H}\mathbf{D}\mathbf{1}_{Q}$, which is clearly positive. Therefore, a final form of $\boldsymbol{\Pi}_{Q}^{H}\boldsymbol{\Pi}_{Q}$ is

[TABLE]

It is immediately obvious that the matrix has only one eigenvalue – $\tfrac{\mathbf{1}_{Q}^{T}\mathbf{D}^{H}\mathbf{D}\mathbf{1}_{Q}}{4}$ – with multiplicity $N$. This eigenvalue is either always positive or zero. From its form, we can surmise that it is zero only in the case when $\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}=0$ for all clutter patches $q$. Why is this so?

First, if $\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}\neq 0$ for all $q$, then $\mathbf{1}_{Q}^{T}\mathbf{D}^{H}\mathbf{D}\mathbf{1}_{Q}=0$ only if the inner products $\mathbf{c}_{q}^{H}\mathbf{c}_{l}=0$ for all $q,l$. However, if this is true, then $\mathbf{c}_{q}^{H}\mathbf{c}_{q}={\|\mathbf{c}_{q}\|^{2}}=0$ for all $q$; in other words, there is no clutter, which is a ridiculous requirement. What if we relax the inner product condition to $\mathbf{c}_{q}^{H}\mathbf{c}_{l}=0$ for $q\neq l$? This doesn’t change things much, because then $\mathbf{1}_{Q}^{T}\mathbf{D}^{H}\mathbf{D}\mathbf{1}_{Q}=\sum_{q=1}^{Q}\frac{(\overline{R}^{q}_{\gamma})^{2}}{2}|\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}|^{2}{\|\mathbf{c}_{q}\|^{2}}$ . This is absolutely positive unless $\mathbf{w}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}=0$ for all $q$, which is the condition we were trying to avoid.

The singular values of $\boldsymbol{\Pi}_{Q}$ are therefore clearly only positive, implying $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{C})$ has positive, negative, and zero eigenvalues. Hence, the clutter Hessian is indefinite, and the clutter cost function is non-convex in $\mathbf{b}$.

III-K STAP Step 6: Definiteness of $\mathcal{\widetilde{H}}[\textsf{f}_{o}]$

At this point, one might think adding a full rank PSD matrix (in this case, the combined noise-and-interference correlation matrix) somewhere will break up the useful structure from the previous sections. However, this is not the case, and a similarly useful structure appears in the Hessian of the complete cost function. We begin by noting that, since $\textsf{f}_{o}$ is a sum of $\textsf{f}_{NI}$ and $f_{C}$, the Hessian is also a sum: $\mathcal{\widetilde{H}}[\textsf{f}_{o}]=\mathcal{\widetilde{H}}[f_{NI}]+\mathcal{\widetilde{H}}[\textsf{f}_{C}]$. Since $f_{NI}=\mathbf{b}^{H}\begin{bmatrix}\mathbf{R_{ni}}&\mathbf{0}_{NML\times N}\\ \mathbf{0}_{N\times NML}&\mathbf{0}_{N\times N}\end{bmatrix}\mathbf{b}=\mathbf{b}^{H}\mathbf{\widetilde{R}_{ni}}\mathbf{b}$ is a real quadratic form, we once again turn to Hjørungnes’ Example 5.2 as cited above, which means the Hessian $\mathcal{\widetilde{H}}[f_{NI}]$ is

[TABLE]

or, more compactly, $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(f_{NI})=\mathbf{\widetilde{R}_{ni}}$.

Hence, the total Hessian of the cost function $\mathcal{\widetilde{H}}[\textsf{f}_{o}]$ is

[TABLE]

where $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{o})=\mathbf{\widetilde{R}_{ni}}+\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{C})$ . Again, due to the structure, the definiteness of $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{o})$ determines the definiteness of the complete Hessian. In order to ascertain the definiteness of $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{o})$, we need to observe its structure:

[TABLE]

Before continuing, we introduce two essential, yet equivalent, theorems from Kreindler & Jameson [19] for the definiteness of a partitioned matrix.

Theorem 1** (via [19],Theorems $I_{a^{\prime\prime}}$ & $I_{b^{\prime\prime}}$ ).**

*A matrix $\mathbf{M}$ partitioned as $\mathbf{M}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\\ \mathbf{B}^{H}&\mathbf{C}\end{bmatrix}$ is non-negative definite if and only if the following conditions are satisfied:

1. A1

$\mathbf{C}\succeq 0$** 2. A2

$\mathbf{B}=\mathbf{B}\mathbf{C}^{\dagger}\mathbf{C}$** 3. A3

$\mathbf{A}-\mathbf{B}\mathbf{C}^{\dagger}\mathbf{B}^{H}\succeq 0$**

or

1. B1

$\mathbf{A}\succeq 0$** 2. B2

$\mathbf{B}=\mathbf{A}\mathbf{A}^{\dagger}\mathbf{B}$** 3. B3

$\mathbf{C}-\mathbf{B}^{H}\mathbf{A}^{\dagger}\mathbf{B}\succeq 0$**

where $\mathbf{A}^{\dagger}$ indicates the Moore-Penrose pseudoinverse of the matrix $\mathbf{A}$ and $\mathbf{A}\succeq 0$ means that $\mathbf{A}$ is non-negative definite.

In the case of the matrix $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(f_{o})$, $\mathbf{A}=\mathbf{R_{ni}},\mathbf{B}=\boldsymbol{\Pi}_{Q}$, and $\mathbf{C}=\mathbf{0}_{N\times N}$. The proof of definiteness follows by checking necessary conditions & observing how and if they are violated. We propose that this matrix is not PSD unless we know, a priori, that all stationary points null every clutter patch. We present two equivalent proofs of our proposition from the sets of A and B conditions in the above theorem.

Proof A.

Any quadratic form involving the all-zero matrix $\mathbf{0}_{N\times N}$ is zero, implying the zero matrix is simultaneously positive semidefinite, negative semidefinite, and indefinite. Thus, Condition A1 is immediately satisfied.

Conditions A2 & A3 require the Moore-Penrose pseudoinverse of the all-zero matrix, which is itself: $\mathbf{0}_{N\times N}^{\dagger}=\mathbf{0}_{N\times N}^{\phantom{*}}$. Then, we can immediately restate condition A2 as $\boldsymbol{\Pi}_{Q}=\mathbf{0}_{NML\times N}$. As has been previously established, this only occurs if $\mathbf{w}_{o}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}_{o}=0$ for all clutter patches $q$ and all stationary points $\mathbf{w}_{o},\mathbf{s}_{o}$ in the feasible set.

After simplification, Condition A3 is equivalent to $\mathbf{R_{ni}}\succeq 0$. Since $\mathbf{R_{ni}}$ is a covariance matrix and all covariance matrices are positive semidefinite by construction, Condition A3 is immediately satisfied. ∎

Proof B.

Since $\mathbf{R_{ni}}$ is a covariance matrix and all covariance matrices are positive semidefinite by construction, Condition B1 is immediately satisfied.

As a consequence of this result, we can use [18, Lemma 7.63(b)] to immediately state that because $\mathbf{R_{ni}}$ is positive semidefinite and Hermitian, its Moore-Penrose pseudoinverse $\mathbf{R_{ni}}^{\dagger}$ is also positive semidefinite.

Furthermore, using [18, Lemma 10.46(b)(i)], we can say that for any matrix $\mathbf{X}\in\mathbb{C}^{N\times NML}$, $\mathbf{X}^{H}\mathbf{R_{ni}}^{\dagger}\mathbf{X}\succeq 0$. Let us set $\mathbf{X}=\boldsymbol{\Pi}_{Q}$. Then, this lemma allows us to state $\boldsymbol{\Pi}_{Q}^{H}\mathbf{R_{ni}}^{\dagger}\boldsymbol{\Pi}_{Q}\succeq 0$, regardless of the choice of $\mathbf{b}$ that composes $\boldsymbol{\Pi}_{Q}$.

However, this statement is, on its face, an immediate contradiction to Condition B3, unless $\boldsymbol{\Pi}_{Q}^{H}\mathbf{R_{ni}}^{\dagger}\boldsymbol{\Pi}_{Q}=\mathbf{0}_{N\times N}$. By [18, Lemma 10.12(b)], this is only true if $\boldsymbol{\Pi}_{Q}^{H}\mathbf{R_{ni}}^{\dagger}=\mathbf{0}_{N\times NML}$ or, similarly, $\mathbf{R_{ni}}^{\dagger}\boldsymbol{\Pi}_{Q}=\mathbf{0}_{NML\times N}$

If, indeed, $\mathbf{R_{ni}}^{\dagger}\boldsymbol{\Pi}_{Q}=\mathbf{0}_{N\times NML}$, then that would mean Condition B2 would read

[TABLE]

But, as we established before, this is only possible if, for each stationary point $\mathbf{w}_{o},\mathbf{s}_{o}$ in the feasible set, $\mathbf{w}_{o}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}_{o}=0$ for all clutter patches $q$. ∎

A third proof relies on a similar construction of Theorem 1 from [20], which we paraphrase in Theorem 2 below:

Theorem 2** (via [20], Theorem 7.7.9(a, b)).**

Let the matrix $\mathbf{M}$ be partitioned as in Theorem 1. The following statements are equivalent:

$\mathbf{M}$* is positive semidefinite* 2. 2.

$\mathbf{A}$* and $\mathbf{C}$ are positive semidefinite and there is a contraction $\mathbf{X}$ shaped identically to $\mathbf{B}$ such that $\mathbf{B}=\mathbf{A}^{1/2}\mathbf{X}\mathbf{C}^{1/2}$.*

In the above theorem, a contraction is defined any matrix whose largest singular value is less than or equal to one, and the superscript $1/2$ denotes the unique matrix square root defined for all PSD matrices. Observe also that this theorem implies that $\mathcal{R}(\mathbf{B})\subseteq\mathcal{R}(\mathbf{A}),\mathcal{R}(\mathbf{B}^{H})\subseteq\mathcal{R}(\mathbf{C})$. Our proof continues as follows:

Proof C.

As established previously, $\mathbf{R_{ni}}\succeq 0$ and $\mathbf{0}_{N\times N}\succeq 0$. Clearly, the matrix square root of the all zeros matrix is itself. Thus, for $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(f_{o})$ to be PSD, $\boldsymbol{\Pi}_{Q}=\mathbf{R_{ni}}^{1/2}\mathbf{X}\mathbf{0}_{N\times N}^{1/2}=\mathbf{0}_{NML\times N}$, which will be true for any contraction $\mathbf{X}$. As established previously, $\boldsymbol{\Pi}_{Q}=\mathbf{0}_{NML\times N}$ only if, for each stationary point $\mathbf{w}_{o},\mathbf{s}_{o}$ in the feasible set, $\mathbf{w}_{o}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}_{o}=0$ for all clutter patches $q$. Therefore, this is the necessary and sufficient condition for $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{o})$ to be PSD. ∎

In any case, $\mathbf{w}_{o}^{H}\boldsymbol{\Gamma}_{q}\mathbf{s}_{o}=0$ for all clutter patches $q$ and stationary points $\mathbf{w}_{o},\mathbf{s}_{o}$ in order for the Hessian to be positive-semidefinite and the objective to be convex. This is a clearly illogical condition – if such a beamformer-signal pair existed a priori, there is effectively no clutter in the region of the desired target response, the signal design is arbitrary and we would only design a beamformer! Additionally, this condition also requires $\mathbf{R_{ni}}\mathbf{w}_{o}=\mathbf{0}_{NML\times 1}$, via the stationary point definitions. If $\mathbf{R_{ni}}$ is full rank, then $\mathbf{w}_{o}=\mathbf{0}_{NML\times 1}$, and we do no processing whatsoever.

Since the only possible ways to not obtain a contradiction are themselves contradictions, we must conclude that, regardless of the beamformer-signal pair, the overall cost function Hessian is indefinite (since the same contradiction would be reached if we desired a negative definite $\mathcal{H}_{\mathbf{b}\mathbf{b}^{*}}(\textsf{f}_{o})$ instead) and the problem is not jointly convex.

IV The Biquadratic Program & Relaxations

In this section, we discuss the relatively unexplored area of biquadratic programming – that is, joint optimization of two multidimensional variables over a cost function that is quadratic in each variable – and its application to fully adaptive radar, which arises from joint signal-beamformer design schemes.

First, we will review the existing literature on biquadratic programming (henceforth, BQP) and related works on tensor approximations, in order to provide a basic mathematical framework. We will then recall the work of Setlur, et al. in order to construct a relevant BQP and its potentially solvable relaxations, with some additional insights into the other unique challenges this problem presents that are unaddressed by the literature.

Interestingly, coverage of the BQP in the literature, from the optimization community or otherwise, can be charitably described as sparse; in fact, there are only seven papers directly addressing the BQP produced by a total of ten authors, all of which have been released within the last six years. This is somewhat puzzling given the depth and breadth of related problems both cited in the works mentioned below and conceivable by an ordinary engineer – applications range from economics to quantum physics to the fully-adaptive radar concept we wish to pursue here. Nevertheless, this limitation allows a reasonable overview of the state-of-the-art without too much difficulty. Other research exists on similar problems, mostly in the areas of nearest tensor approximation, but this is extensive and highly unspecialized. Therefore, we refer the reader to the multiple overlapping sources within the papers mentioned below on these matters.

IV-A The General BQP

We begin by describing the general biquadratic program analyzed by Ling, et al. in [13]. Consider the vectors $\mathbf{x}\in\mathbb{R}^{n}$ and $\mathbf{y}\in\mathbb{R}^{m}$. The homogenous biquadratic optimization problem is

[TABLE]

where $b_{ijkl}$ is the $(i,j,k,l)$th element of the 4th-order tensor $\mathcal{B}\in\mathbb{R}^{n\times m\times n\times m}$, the subscripted $x$ and $y$ are the appropriate element of the corresponding vector, and $\|\cdot\|$ is the standard Euclidean 2-norm. As the authors note, for fixed $\mathbf{x}$ (resp. $\mathbf{y}$), this problem is quadratic in $\mathbf{y}$ (resp. $\mathbf{x}$), can be solved quite easily and, in fact, is convex if the appropriate inner matrix is positive semidefinite. This property also leads to the name of this class of problems (cf. bilinear optimization). Before we continue, it is worth noting that the objective function $b(\mathbf{x},\mathbf{y})$ can be rewritten using a tensor operator and the matrix inner product, viz.

[TABLE]

where $\mathbf{X}\bullet\mathbf{Y}=\operatorname{Tr}(\mathbf{X}^{T}\mathbf{Y})$, and the tensor-operator matrices are $\mathcal{B}\mathbf{x}\mathbf{x}^{T}=\sum_{i,k=1}^{n}b_{ijkl}x_{i}x_{k}$ and $\mathbf{y}\mathbf{y}^{T}\mathcal{B}=\sum_{j,l=1}^{m}b_{ijkl}y_{j}y_{l}$.

The above problem has been modified elsewhere in the literature, mostly in terms of the constraints. Bomze, et al. [21] constrain the objective in Equation 24 to lie on the simplex for each variable. Ling, Zhang, and Qi have described strategies when the problem is quadratically constrained in [14, 22].

Regardless of the constraints, [13] proves that the BQP is nonconvex, that both the BQP and its naive semidefinite relaxation (which we will more explicitly discuss in the next section) are NP-hard problems, and that the initial general problem does not even admit a polynomial time approximation algorithm. This would seem to bode ill for any hope of reasonable solutions, but the authors demonstrate that polynomial-time approximation algorithms exist (with bounds approximately inversely proportional to the dimension of the variables) when the objective is square-free or has squared terms in only one of the variables. Finally, they demonstrate the creation of such algorithms, showing a tradeoff between speed and accuracy for sum-of-squares based solvers (better accuracy) versus convex semidefinite relaxation solvers (faster). Additionally, if the structure is sufficiently sparse, then existing solvers can attack the semidefinite relaxation even more efficiently.

IV-B The Problem at Hand

Before continuing, we recall the basics of our problem. We operate on a radar datacube collected over $L$ pulses with $N$ fast time samples per pulse in $M$ spatial bins. Our goal is to find a STAP beamformer vector $\mathbf{w}\in\mathbb{C}^{NML}$ and a transmit signal $\mathbf{s}\in\mathbb{C}^{N}$ that minimizes the combined effect of the noise and interference represented by the covariance matrix $\mathbf{R_{ni}}\in\mathbb{C}^{NML\times NML}$) and the signal-dependent clutter. The signal-dependent clutter is modeled as $Q$ independent clutter patches. Each patch has an individual spatiodoppler response matrix $\boldsymbol{\Gamma}_{q}\in\mathbb{C}^{NML\times N}$, which ties into the per-patch clutter covariance $\mathbf{R}_{\boldsymbol{\gamma}}^{q}(\mathbf{s})=\boldsymbol{\Gamma}_{q}\mathbf{s}\mathbf{s}^{H}\boldsymbol{\Gamma}_{q}^{H}$. The overall clutter covariance is then $\sum\limits_{q=1}^{Q}\boldsymbol{\Gamma}_{q}\mathbf{s}\mathbf{s}^{H}\boldsymbol{\Gamma}_{q}^{H}$.

We constrain this minimization in two ways. First, for a given target space-time-doppler bin, we want a particular filter output, say, $\kappa\in\mathbb{C}$. This filter output can be represented by the Capon beamformer equation $\mathbf{w}^{H}\mathbf{T}\mathbf{s}$, where $\mathbf{T}\in\mathbb{C}^{NML\times N}$ is the spatiodoppler response of the target bin. Second, we place an upper bound on the total signal power, say, $P_{o}$. With these constraints in place, the overall STAP problem is given by

[TABLE]

If we treat the signal and beamformer as a single stacked variable, say $\mathbf{b}=\begin{bmatrix}\mathbf{w}^{T}&\mathbf{s}^{T}\end{bmatrix}^{T}\in\mathbb{C}^{N(ML+1)}$, we can find an equivalent form of the above optimization problem in the new variable. Define $\boldsymbol{\Psi}_{W},\boldsymbol{\Psi}_{S}$ as the matrices that recover $\mathbf{w},\mathbf{s}$ from $\mathbf{b}$, i.e. $\mathbf{w}=\boldsymbol{\Psi}_{W}\mathbf{b}$ and $\mathbf{s}=\boldsymbol{\Psi}_{S}\mathbf{b}$. Then, we can define the optimization problem as

[TABLE]

where $\mathbf{\widetilde{R}_{ni}}=\boldsymbol{\Psi}_{W}^{T}\mathbf{R_{ni}}\boldsymbol{\Psi}_{W}$, $\overline{\boldsymbol{\Gamma}}_{q}=\boldsymbol{\Psi}_{W}^{T}\boldsymbol{\Gamma}_{q}\boldsymbol{\Psi}_{S}$, and $\mathbf{\widetilde{T}}=\boldsymbol{\Psi}_{W}^{T}\mathbf{T}\boldsymbol{\Psi}_{S}$ are just “expanded” versions of the matrices seen in the previous problem.

In either case, these are clearly biquadratic programs, which we’ve shown are non-convex. The question now is: can we find a solvable representation, relaxed or otherwise?

IV-C Getting to the SDR

Our first goal is using the notation/mechanisms of the literature to find the semidefinite relaxation of the equivalent programs BQP 1 and BQP 2. First, assume we have defined a tensor $\mathcal{C}$ such that the following operator relation holds: $\mathcal{C}\mathbf{s}\mathbf{s}^{H}=\sum_{q=1}^{Q}\mathbf{R}_{\boldsymbol{\gamma}}^{q}(\mathbf{s})$. Similiarly, assume an expanded form of this tensor exists, say $\widetilde{\mathcal{C}}$, such that $\widetilde{\mathcal{C}}\mathbf{b}\mathbf{b}^{H}=\sum\limits_{q=1}^{Q}\overline{\boldsymbol{\Gamma}}_{q}\mathbf{b}\mathbf{b}^{H}\overline{\boldsymbol{\Gamma}}_{q}^{H}$. What this tensor looks like, exactly, will be seen in a later section.

On complex matrices, the inner product becomes $\mathbf{X}\bullet\mathbf{Y}=\operatorname{Re}\{\operatorname{Tr}(\mathbf{X}^{H}\mathbf{Y})\}$. We’d like to get everything into a real form so we can use this inner product somehow, as in the literature.

Note the first constraint (since $\kappa$ is complex) implies its conjugate must also be constrained, i.e.

[TABLE]

Hence, the target-based constraints are

[TABLE]

We can recast this constraint as two real constraints. First, recall the obvious things about complex numbers:

[TABLE]

If we substitute in the Capon constraints, we then obtain

[TABLE]

where the subscripts $H$ and $A$ indicate the hermitian and anti-hermitian part of the expanded matrix $\mathbf{\widetilde{T}}$, respectively. (As an aside, recall that $\mathbf{\widetilde{T}}_{H}^{H}=\mathbf{\widetilde{T}}_{H}$ and $\mathbf{\widetilde{T}}_{A}^{H}=-\mathbf{\widetilde{T}}_{A}$ by definition) This leads us to a new version of the second form of the optimization problem, as a function of the combined beamformer-signal vector:

[TABLE]

or, using the tensor form,

[TABLE]

Using the properties of trace and the fact that all quantities in this optimization problem are now real, we can recast it using the inner product from above:

[TABLE]

Note that we’ve used the Hermitian (resp. anti-Hermitian) nature of $\mathbf{\widetilde{R}_{ni}}$ and $\mathbf{\widetilde{T}}_{H}$ (resp. $\mathbf{\widetilde{T}}_{A}$) to avoid the Hermitian superscript everywhere.

A common path to semidefinite relaxations defines a matrix variable, say, $\mathbf{B}=\mathbf{b}\mathbf{b}^{H}$. It is clear that $\mathbf{Z}$ is Hermitian and positive semidefinite, and $\operatorname{rank}(\mathbf{B})=1$. This means our optimization problem above is equivalent to the matrix quadratic program

[TABLE]

Depending on the properties of $\widetilde{\mathcal{C}}$, this could be a convex problem in $\mathbf{Z}$ if not for the rank condition, which is the primary hurdle. This is where the ”relaxation” in semidefinite relaxation comes in. If we permit $\mathbf{Z}$ to take any rank, we have relaxed the problem to the (possibly) convex one:

[TABLE]

Notice we did not, at any point, progress through a matrix bilinear form in finding this semidefinite relaxation, as seen in [13], etc. This is because the constraint set in this problem prevents such a representation – namely, the Capon constraint, which is bilinear in the vector arguments already.

If we can obtain a solution to this problem, a natural question is how to apply this solution to the original problem. If we are lucky, the rank of the relaxed solution will be one, and there will be nothing left to do. If the rank of the solution is greater than one, then we can use the nearest rank-one approximation that still satisfies the constraints. Some authors have We will see in the simulations below that for this particular scenario, we can obtain solutions that are nearly rank-one in a numerical sense.

V Solutions of the Relaxed BQP

Recall that, when standardized, the semidefinite relaxation of the BQP is:

[TABLE]

Before we begin, we reformulate the problem to use the variable $\operatorname{vec}(\mathbf{B})$. First, note that $\operatorname{Tr}(\mathbf{A}\mathbf{B})=\operatorname{vec}^{T}(\mathbf{A}^{T})\operatorname{vec}(\mathbf{B})$. For Hermitian $\mathbf{A}$, this becomes $\operatorname{Tr}(\mathbf{A}\mathbf{B})=\operatorname{vec}^{H}(\mathbf{A})\operatorname{vec}(\mathbf{B})$ (since $\mathbf{A}^{T}=\mathbf{A}^{*}$ and the conjugation operator is linear). This allows us to say that, for Hermitian $\mathbf{A}$, the relevant inner product on $\mathbb{C}$ is $\mathbf{A}\bullet\mathbf{B}=\operatorname{vec}^{H}(\mathbf{A})\operatorname{vec}(\mathbf{B})$.

Next, recall that the tensor $\widetilde{\mathcal{C}}$ can be unfolded into the vector matrix $\widetilde{\mathbf{C}}_{V}=\sum_{q=1}^{Q}\operatorname{vec}(\overline{\boldsymbol{\Gamma}}_{q})\operatorname{vec}(\overline{\boldsymbol{\Gamma}}_{q})^{H}$. Using these, we recast SSDR into a vectorized form:

[TABLE]

In a further “development”, let us explicitly put the equality constraint back as a complex number and break out the positive semi-definiteness as a separate constraint. This makes VSSDR turn into

[TABLE]

According to Hjørungnes, there are certain methods we need to account for the Hermitian constraint in the following analysis. For notational ease, let us define the following vectors:

[TABLE]

Furthermore, if we define the following matrices,

[TABLE]

then these can be rewritten as

[TABLE]

V-A Proof of Convexity

Before continuing, we demonstrate directly that this relaxed optimization problem is indeed convex. Recall that the relaxed problem is

[TABLE]

A convex optimization problem minimizes a convex objective over a convex set. A well-known restriction of this concept is given in [12, pp. 136-137] as follows: for a given optimization problem, if

the objective function is convex, 2. 2.

the inequality constraint functions are convex, and, 3. 3.

the equality constraint functions are affine (i.e. it looks like $\mathbf{a}^{H}\mathbf{x}-b$.)

then the problem is convex. This mostly stems from the fact that affine equality constraints define a polytope, which is convex, and so its intersection with the existing convex set is also convex.

We start with the objective. A quadratic function is convex if the associated matrix is positive semidefinite. Thus, since $\tilde{\mathbf{C}}_{V}\succeq 0$ by definition (as a sum of rank one Hermitian matrices), then $\operatorname{vec}^{H}(\mathbf{B})\widetilde{\mathbf{C}}_{V}\operatorname{vec}(\mathbf{B})$ is a convex function. Adding an affine function to it does not change the convexity, so the overall objective is convex. Continuing to the constraints, both the power constraint and the semidefiniteness constraints are convex – the former because it is affine, the latter because it describes a convex cone. Finally, the equality constraint is convex, because it can be broken up into two real affine constraints on the real and imaginary parts. Therefore, the relaxed problem is convex. However, since $\tilde{\mathbf{C}}_{V}$ is typically not full rank (indeed, the rank is bounded above by $ML$ in this case), it is not usually strictly convex, and thus there will exist a multitude of solutions.

V-B Slater’s Condition

Recall that strong duality is said to hold for a given optimization problem is the primal problem is convex and it satisfies Slater’s condition (aka the problem is strictly feasible). A consequence of this property is that there is no duality gap between solutions of the primal and dual problems. We will show that our problem satisfies Slater’s condition given a relationship between the constraint values.

Given an optimization problem

[TABLE]

(where $\Omega$ is a convex set), let us assume $f_{o}(\mathbf{x})$ is convex and define the set $\mathcal{D}=\Omega\cap\,\operatorname{dom}(f_{o})\,\cap\,(\cap_{i=1}^{m}\operatorname{dom}(g_{i}))$ to be the total feasible domain of the problem. Slater’s condition is satisfied if there exists at least one point in the relative interior (i.e. not on the boundary) of the problem’s feasibility set that satisfy all of the equality constraints and strictly satsifies the inequality constraints. Mathematically, we can represent the condition as follows:

Lemma 1** (Slater’s Condition).**

For the standard optimization problem above, if each $g_{i}$ is convex and there exists a point $\mathbf{x}_{o}\in\operatorname{relint}(\mathcal{D})$ such that $g_{i}(\mathbf{x}_{o})<0\;\forall i$ and $h_{j}(\mathbf{x}_{o})=0\;\forall j$, then strong duality holds and there is zero duality gap.

For our problem, the feasibility set is positive semidefinite matrices, which means its relative interior is positive-definite matrices. Our constraint functions are

[TABLE]

where $\boldsymbol{\beta}=\operatorname{vec}(\mathbf{B})$ for a given Hermitian matrix $\mathbf{B}$. Thus, we need to find a matrix $\mathbf{B}\succ 0$ that satisfies the above equations.

We start with solving for $\boldsymbol{\beta}$ in our equality constraints, then ensure that the resulting matrix is both positive definite and strictly satisfies the power inequality. As a combined matrix-vector equation, we have:

[TABLE]

If solutions of this equation exist (and they should exist), then they are given by

[TABLE]

where $\dagger$ indicates the pseudoinverse and $\mathbf{z}$ is the vectorization of a matrix $\mathbf{Z}$ similarly shaped to $\mathbf{B}$.

A logical question is, then, what is that pseudoinverse (and consequently, what is the existence condition)? Since the matrix $\begin{bmatrix}\tilde{\boldsymbol{\tau}}^{T}\mathbf{K}_{J,J}\\ \tilde{\boldsymbol{\tau}}^{H}\end{bmatrix}$ is fat, we might have some hope in the form $\mathbf{C}^{\dagger}=\mathbf{C}^{H}(\mathbf{C}\mathbf{C}^{H})^{-1}$. First, we can find the matrix we hope to invert:

[TABLE]

So long as we have an actual target, this is always invertible (and hence, a solution always exists). This means that the pseudoinverse of interest is

[TABLE]

Additionally, the projection matrix in the solution is

[TABLE]

where $\boldsymbol{\tau}_{H}=\operatorname{vec}(\mathbf{T}^{H})$.

We can break our solution into the “min-norm” and “nullspace” parts. Clearly, the min-norm part is

[TABLE]

It is clear from the structure that this is not a positive definite matrix (or even a positive semidefinite one), so we require the nullspace component to place us into the relative interior. If we break the vector $\mathbf{z}$ into components as $\mathbf{z}=\mathbf{P}_{WW}^{T}\mathbf{z}_{1}+\mathbf{P}_{SW}^{T}\mathbf{z}_{2}+\mathbf{P}_{WS}^{T}\mathbf{z}_{2H}+\mathbf{P}_{SS}^{T}\mathbf{z}_{3}$, then the nullspace solution is:

[TABLE]

Hence, the overall matrix that satisfies the equality constraints is

[TABLE]

Our problem now becomes finding a matrix $\mathbf{Z}$ such that $\mathbf{B}\succ 0$ and $\operatorname{tr}(\mathbf{Z}_{3})<P_{o}$. One method of attack is setting $\mathbf{Z}_{2}=\mathbf{0}_{NML\times N}$ and attempting to reach a valid result this way. Under this assumption, our solution matrix is

[TABLE]

The positive definiteness requirement can then be expressed in one of two ways:

[TABLE]

Here, we rely on another judicious guess, setting $\mathbf{Z}_{1}=\mathbf{I}_{NML}$, which is clearly positive definite.111We can make this guess because $NML>N$. If, for whatever reason, the number of transmit resources were greater than the number of receive resources, then we could start with $\mathbf{Z}_{3}$ and continue from there.Then, we only need to find a $\mathbf{Z}_{3}$ such that $\mathbf{Z}_{3}\succ\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{4}_{F}}}\mathbf{T}^{H}\mathbf{T}$ and $\operatorname{tr}(\mathbf{Z}_{3})<P_{o}$.

Let us assume that such a matrix exists. Following [20, Corollary 7.7.4(d)], since $\mathbf{Z}_{3}\succ\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{2}_{F}}}\mathbf{T}^{H}\mathbf{T}$, then $\operatorname{tr}(\mathbf{Z}_{3})>\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{4}_{F}}}\operatorname{tr}(\mathbf{T}^{H}\mathbf{T})=\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{2}_{F}}}$. But we already know that $\operatorname{tr}(\mathbf{Z}_{3})<P_{o}$. Hence, we have the chained inequality

[TABLE]

which obviously requires that $P_{o}>\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{2}_{F}}}$ or perhaps ${\|\mathbf{T}\|^{2}_{F}}>\frac{|\kappa|^{2}}{P_{o}}$.

With this in hand, we can now construct our matrix. Since $NML>N$ and (generally) $\operatorname{rank}(\mathbf{T})=N$, then $\mathbf{T}^{H}\mathbf{T}\succ 0$.222Again, if our available resources are swapped, we can use the outer product instead and get a similar answer. For any matrix $\mathbf{A}\succ 0$ and real numbers $\alpha,\beta$, $\alpha\mathbf{A}\succ\beta\mathbf{A}$ if and only if $\alpha>\beta$. Let $\mathbf{Z}_{3}=\frac{P_{o}-\varepsilon}{{\|\mathbf{T}\|^{2}_{F}}}\,\mathbf{T}^{H}\mathbf{T}$, where $\varepsilon>0$ is some small real number. This guarantees that the trace feasibility is strictly satisfied, since $\operatorname{tr}(\mathbf{Z}_{3})=P_{o}-\varepsilon<P_{o}$. It also satisfies the chained inequality above, so long as $\epsilon<P_{o}-\frac{|\kappa|^{2}}{{\|\mathbf{T}\|^{2}_{F}}}$. Thus, the chain inequality is the only controlling element, and thus there is no duality gap so long as ${\|\mathbf{T}\|^{2}_{F}}>\frac{|\kappa|^{2}}{P_{o}}$.

In the side-looking STAP case, the target matrix has a Kronecker structure – that is, \mathbf{T}=\mathbf{v}_{t}\raisebox{1.0pt}{>\otimes>}\mathbf{I}_{N}\raisebox{1.0pt}{>\otimes>}\mathbf{a}_{t}. This means that $\mathbf{T}^{H}\mathbf{T}={\|\mathbf{v}_{t}\|^{2}}{\|\mathbf{a}_{t}\|^{2}}\mathbf{I}_{N}$. If we assume that the spatial response comes from a uniform linear array and that the Doppler response is similarly uniform, then ${\|\mathbf{a}_{t}\|^{2}}=M$ and ${\|\mathbf{v}_{t}\|^{2}}=L$. This further implies that $\mathbf{T}^{H}\mathbf{T}=ML\mathbf{I}_{N}$ and ${\|\mathbf{T}\|^{2}_{F}}=NML$. If we apply this knowledge to the positive definiteness condition, we need to find a $\mathbf{Z}_{3}$ such that $\mathbf{Z}_{3}\succ\frac{|\kappa|^{2}}{N^{2}ML}\mathbf{I}_{N}$. This clearly exists if we set $\mathbf{Z}_{3}=\alpha\mathbf{I}_{N}$ where $\alpha>\frac{|\kappa|^{2}}{NML}$. Following the path from the general case, this means that the Slater condition is satisfied so long as $NML>\frac{|\kappa|^{2}}{P_{o}}$.

V-C Obtaining the KKTs

First, let us adopt the convention following conventions: each inequality constraint (multiplier) is given by $g_{j}(\cdot)$ ($\lambda_{j}$), each equality constraint (multiplier) is given by $h_{i}(\cdot)$ ($\mu_{i}$), and an optimal value of a variable is denoted by a superscript o (i.e., $\mathbf{B}^{o}$). The KKT conditions can then be written as

$\nabla_{\mathbf{B}}\mathcal{L}(\mathbf{B}^{o},\boldsymbol{\mu}^{o},\boldsymbol{\lambda}^{o})=\mathbf{0}_{J\times J}$ 2. 2.

$\lambda_{j}^{o}\geq 0\;\forall j$ 3. 3.

$\lambda_{j}^{o}g_{j}(\mathbf{B}^{o})=0\;\forall j$ 4. 4.

$g_{j}(\mathbf{B}^{o})\leq 0\;\forall j$ 5. 5.

$h_{i}(\mathbf{B}^{o})=0\;\forall i$ 6. 6.

$\mathcal{H}[\mathcal{L}(\mathbf{B}^{o},\boldsymbol{\mu}^{o},\boldsymbol{\lambda}^{o})]\succeq 0$

That is, for a matrix $\mathbf{B}^{o}$ to be a regular minimizer of the related optimization problem, it is necessary that it satisfies these conditions. If the problem is convex, these are necessary and sufficient conditions for the optimal minimizer and associated Lagrange multipliers. For future notational simplicity, we will drop the superscript $o$ until absolutely necessary (i.e., a final statement of the optimal solution). According to [17], the gradient in the first KKT condition is $\nabla_{\mathbf{B}}=\frac{\partial}{\partial B^{*}}$ or, in other words, $\mathcal{D}_{\mathbf{B}^{*}}\mathcal{L}=\operatorname{vec}^{T}(\frac{\partial\mathcal{L}}{\partial B^{*}})$. Thus, this vectorized form can take the place of the gradient.

V-C1 KKT Condition 1

Let $\boldsymbol{\Sigma}$ be the slackness variable associated with the PSD condition on $\mathbf{B}$, and $\boldsymbol{\sigma}$ be its vectorization. The Lagrangian is

[TABLE]

Condition 1 is, after taking the derivatives,

[TABLE]

In general, the full set of solutions to the second form of Condition 1 is given by

[TABLE]

where $\mathbf{z}$ is an arbitrary vector isomorphic to a Hermitian matrix $\mathbf{Z}$, partitioned identically to $\mathbf{B}$ and $\boldsymbol{\Sigma}$. This solution exists provided the following existence condition holds:

[TABLE]

This may seem daunting, but we can unpack parts of this quite easily. Recall that the vectors $\boldsymbol{\beta},\boldsymbol{\sigma}$ can be partitioned into the sums:

[TABLE]

and that these partitions are disjoint. Due to the disjointness, a quick examination of the second form of Condition 1 gives us a few results:

[TABLE]

The first two equations result from the clutter matrix on the left hand side annihilating the $\mathbf{P}_{WW}^{T}$ and $\mathbf{P}_{SS}^{T}$ components of $\boldsymbol{\beta}$, and imply $\boldsymbol{\Sigma}_{1}=\mathbf{R_{ni}}$ & $\boldsymbol{\Sigma}_{3}=\lambda\mathbf{I}_{N}$. The third equation covers both the $\mathbf{P}_{SW}^{T}$ and $\mathbf{P}_{WS}^{T}$ components, because the $\mathbf{P}_{WS}^{T}$ component is a conjugate permutation of the $\mathbf{P}_{SW}^{T}$ component. The solution (if it exists) of that third component, generally, is

[TABLE]

where $\mathbf{z}_{2}$ is the $\mathbf{P}_{SW}^{T}$ component of the arbitrary vector $\mathbf{z}$ above. The existence condition for this solution is, unsurprisingly,

[TABLE]

or, after rearrangement,

[TABLE]

With some algebraic fiddling, one can verify that Equations 29a, 29b, 30, and 31 are equivalent to most of the much uglier solution above. The remaining necessary mathematical spackle is setting $\mathbf{B}_{1}=\mathbf{Z}_{1}$ and $\mathbf{B}_{3}=\mathbf{Z}_{3}$, again due to the annihilating nature of the clutter matrix on the block-diagonal terms. This more or less says the diagonal blocks of the matrix are arbitrary, subject to satisfying the other KKT conditions.

We conclude this section by repeating the necessary and sufficient conditions for the first KKT condition to be satisfied, as found above:

[TABLE]

where $\mathbf{P}_{\mathbf{C}_{V}}^{\perp}=\mathbf{I}_{N^{2}ML}-\mathbf{C}_{V}\mathbf{C}_{V}^{\dagger}=\mathbf{I}_{N^{2}ML}-\mathbf{C}_{V}^{\dagger}\mathbf{C}_{V}$ is the orthogonal projector onto the nullspace of $\mathbf{C}_{V}$.

V-C2 KKT Conditions 2-4: The Inequality Constraints

With the gradient condition exhausted, we turn to the inequality constraints (power bound, positive-semidefiniteness of the solution) & their related conditions. For convenience, we shall attack these somewhat independently in separate subsections, though they will interact.

The power constraint

The KKT conditions related to the power constraint are as follows:

[TABLE]

Rearranging with our partitioned variables, we get

[TABLE]

A nice result of the slackness condition is $\lambda\operatorname{Tr}(\mathbf{B}_{3})=\lambda P_{o}$, which we will use later. The other “result” is that $\lambda=0$ when the solution does not reach the power bound, and $\lambda>0$ when it does. This will inform interpretations of the matrix situation below.

Positive semidefiniteness of $\mathbf{B}$

The conditions for semidefiniteness of the relaxed beamformer-signal basis are slightly more complex, but reveal a significant amount of structure to the solution. First, the direct form of these conditions are

[TABLE]

Of course, in this form, they are not especially useful. However, recall that from the first KKT condition, we know $\boldsymbol{\Sigma}$ to some extent:

[TABLE]

Via Theorem 1, we have conditions for the positive-semidefiniteness of the basis matrix and its slackness variable. To wit, the basis matrix $\mathbf{B}=\begin{bmatrix}\mathbf{B}_{1}&\mathbf{B}_{2}\\ \mathbf{B}_{2}^{H}&\mathbf{B}_{3}\end{bmatrix}$ is PSD if and only if

[TABLE]

and the slackness matrix $\boldsymbol{\Sigma}$ is PSD if and only if

[TABLE]

where $\lambda^{\dagger}$ is a scalar pseudoinverse, equaling $\lambda^{-1}$ when $\lambda>0$ and zero otherwise. Using Theorem 2, we can also say that Conditions $\Sigma$P1–$\Sigma$P3 are satisfied if a contraction $\mathbf{X}\in\mathbb{C}^{NML\times N}$ exists such that $\boldsymbol{\Sigma}_{2}=\sqrt{\lambda}\mathbf{R}_{\mathbf{ni}}^{1/2}\mathbf{X}$.

The complementary slackness condition for the matrix case reduces to 4 equalities, given below

[TABLE]

We can use these conditions to find equivalent forms for $\lambda$. First, taking the trace of Condition CS4, we have

[TABLE]

Since $\lambda$ is both real and nonnegative, this means that $\operatorname{tr}(\mathbf{B}_{2}^{H}\boldsymbol{\Sigma}_{2})$ is real and nonpositive! We also have $\operatorname{tr}(\mathbf{B}_{2}^{H}\boldsymbol{\Sigma}_{2})=\operatorname{tr}(\boldsymbol{\Sigma}_{2}^{H}\mathbf{B}_{2})$, which we can apply to the trace of Condition CS1 to obtain another form of $\lambda$:

[TABLE]

We can also apply this logic to Equation (29c) to reveal an interesting consequence about the cost function. When vectorized, $\operatorname{tr}(\mathbf{B}_{2}^{H}\boldsymbol{\Sigma}_{2})=\boldsymbol{\beta}_{2}^{H}\boldsymbol{\sigma}_{2}$. Then, we have

[TABLE]

(Incidentally, this means that $\tilde{\mu}^{*}\kappa$ is real, and hence the optimal phase of $\tilde{\mu}$ is that of $\kappa$.) From the above, however, we can see that $-\boldsymbol{\beta}_{2}^{H}\boldsymbol{\sigma}_{2}=\lambda P_{o}$, and so

[TABLE]

The right hand side of the final equation is immediately recognizable as our objective function, which in some sense implies that this would be the only remaining aspect of the dual.

V-D KKT Condition 5: Equality Constraints

This is the final major KKT condition left to examine, because KKT Condition 6 is trivially satisfied by $\mathbf{C}_{V}$ being positive semidefinite. Here, our primary concern is the equality constraint $\boldsymbol{\beta}^{H}_{2}\boldsymbol{\tau}=\kappa$. According to the first KKT condition, we know that

[TABLE]

Substituting this into the equality constraint gives us

[TABLE]

Additionally, recall that, for this solution to exist, the following condition on $\boldsymbol{\sigma}_{2}$ must hold:

[TABLE]

VI Consequences of the KKTs

The optimality conditions shown above are, at first blush, a complicated set of matrix equations to solve. However, we can derive some insight into the nature of the relaxed solution by attacking small pieces of it. In this section, we will first demonstrate generic properties of every solution to the KKTs, then show a potential power-bounded solution under certain conditions. Finally, we will provide a roadmap for non-power-bounded solutions and their feasibility.

VI-A General properties of the relaxed solution

First, we show that the solution does not reach the power bound $P_{o}$ iff the slackness matrix $\boldsymbol{\Sigma}_{2}=\mathbf{0}_{NML\times N}$.

Lemma 2**.**

$\lambda=0\iff\boldsymbol{\Sigma}_{2}=\mathbf{0}_{NML\times N}$**

Proof:

First, we proceed in the forward direction. If $\lambda=0$, then the slackness matrix becomes

[TABLE]

To be part of a feasible solution, this must be positive semidefinite, which is only possible if $\boldsymbol{\Sigma}_{2}=\boldsymbol{\Sigma}_{2}(\mathbf{0}_{N\times N})^{\dagger}(\mathbf{0}_{N\times N})=\mathbf{0}_{NML\times N}$ (see Theorem 1). Hence, the forward direction is proved.

In the reverse direction, we prove through contradiction. Assume that $\boldsymbol{\Sigma}_{2}=\mathbf{0}_{NML\times N}$ and $\lambda>0$. The matrix slackness condition CS2 dictates that $\lambda\mathbf{B}_{2}=-\mathbf{B}_{1}\boldsymbol{\Sigma}_{2}$. Under our first assumption, this becomes $\lambda\mathbf{B}_{2}=\mathbf{0}_{NML\times N}$, which simplifies to $\mathbf{B}_{2}=\mathbf{0}_{NML\times N}$ under the second assumption. However, any feasible solution must also satisfy the Capon constraint $\operatorname{tr}(\mathbf{B}_{2}^{H}\mathbf{T})=\kappa\neq 0$. Clearly, $\mathbf{B}_{2}=\mathbf{0}_{NML\times N}$ violates this constraint, which leads to our contradiction and completes the proof. ∎

Next, we show that every feasible solution reaches the power bound if and only if the noise-and-interference correlation matrix $\mathbf{R_{ni}}$ is full rank. This is a considerably longer and more complex proof.

Proposition 1**.**

$\lambda>0\iff\operatorname{rank}(\mathbf{R_{ni}})=NML$**

Proof:

We will prove this by contradiction as well. Assume that $\mathbf{R_{ni}}$ is full rank and $\lambda=0$. Given the second condition, Lemma 2 requires that $\boldsymbol{\Sigma}_{2}=\mathbf{0}_{NML\times N}$. If we apply this to the slackness condition CS1, then $\mathbf{B}_{1}\mathbf{R_{ni}}=\mathbf{0}_{NML\times NML}$. However, since $\mathbf{R_{ni}}$ is full rank, this implies $\mathbf{B}_{1}=\mathbf{0}_{NML\times NML}$. This violates our non-triviality, but we will continue with the proof to show we reach a further contradiction. Since the overall solution matrix must be PSD, $\mathbf{B}_{2}=\mathbf{0}_{NML\times N}$ as well. As in the proof of Lemma 2, we have reached a contradiction because the Capon constraint is violated, which completes the proof. ∎

Proof:

First, observe that for a non-trivial PSD solution matrix, $\mathcal{R}(\mathbf{B}_{2}^{H})\subseteq\mathcal{R}(\mathbf{B}_{3})$ by Theorem 2. Next, we turn to slackness condition CS4, which states $\lambda\mathbf{B}_{3}=-\mathbf{B}_{2}^{H}\boldsymbol{\Sigma}_{2}$. If $\lambda>0$, then clearly $\mathcal{R}(\mathbf{B}_{3})\subseteq\mathcal{R}(\mathbf{B}_{2}^{H})$. By the standard rules of subset inclusion, then, $\mathcal{R}(\mathbf{B}_{2}^{H})=\mathcal{R}(\mathbf{B}_{3})$ and $\operatorname{rank}(\mathbf{B}_{2}^{H})=\operatorname{rank}(\mathbf{B}_{3})$. This fact will become useful later.

We will now use the results of [23] on another slackness condition, CS3, and use a substitution of CS4 to get to our destination. Using [23, Theorem 2.2] on CS3 to solve for $\mathbf{R_{ni}}$, we have the following requirements for $\mathbf{R_{ni}}$ to be at least PSD (which it is):

$-\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H}\mathbf{B}_{2}\succeq 0$: If we substitute CS4 into this, we obtain $\lambda\mathbf{B}_{3}^{2}\succeq 0$, which is satisfied because $\lambda>0$ and $\mathbf{B}_{3}$ is PSD. 2. 2.

$\mathcal{R}(\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H})\subseteq\mathcal{R}(\mathbf{B}_{2}^{H})$: This is satisfied because $\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H}=-\frac{1}{\lambda}\mathbf{B}_{2}^{H}\boldsymbol{\Sigma}_{2}\boldsymbol{\Sigma}_{2}^{H}$ via CS4, and the range inclusion follows directly. 3. 3.

$\operatorname{rank}(-\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H}\mathbf{B}_{2})=\operatorname{rank}(\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H})$: This is not immediately satisfied by the other conditions, but it does imply that $\operatorname{rank}(\mathbf{B}_{3})=\operatorname{rank}(\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H})$.

The more interesting requirement comes from [24], which adds the following: $\mathbf{R_{ni}}$ is positive definite (and thus full rank) if and only if $\operatorname{rank}(-\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H}\mathbf{B}_{2})=\operatorname{rank}(\mathbf{B}_{2}^{H})$. We know from the third PSD requirement above that $\operatorname{rank}(-\mathbf{B}_{3}\boldsymbol{\Sigma}_{2}^{H}\mathbf{B}_{2})=\operatorname{rank}(\mathbf{B}_{3})$, and thus $\mathbf{R_{ni}}$ is full rank if and only if $\operatorname{rank}(\mathbf{B}_{3})=\operatorname{rank}(\mathbf{B}_{2}^{H})$. However, as seen from above, this condition is already satisfied if $\lambda>0$, and so our proof is complete. ∎

Since the proof of Proposition 1 provides us with the evidence to show that power-bounded solutions occur only in non-singular noise-and-interference environments, we can also demonstrate an additional property of any power-bounded solution: namely, the rank of the overall solution matrix. Recall that for $\mathbf{B}$ to be PSD, $\mathcal{R}(\mathbf{B}_{2})\subseteq\mathcal{R}(\mathbf{B}_{1})$. However, slackness condition CS1 provides that $\mathbf{B}_{1}\mathbf{R_{ni}}=-\mathbf{B}_{2}\boldsymbol{\Sigma}_{2}^{H}$. Since $\mathbf{R_{ni}}$ is full rank, CS1 becomes $\mathbf{B}_{1}=-\mathbf{B}_{2}\boldsymbol{\Sigma}_{2}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}$, which implies $\mathcal{R}(\mathbf{B}_{1})\subseteq\mathcal{R}(\mathbf{B}_{2})$. Thus, $\mathcal{R}(\mathbf{B}_{1})=\mathcal{R}(\mathbf{B}_{2})$, and $\operatorname{rank}(\mathbf{B}_{1})=\operatorname{rank}(\mathbf{B}_{2})=\operatorname{rank}(\mathbf{B}_{3})\leq N$, where the last equality is implied by Proposition 1’s proof and the inequality is obvious.

In the sidelooking STAP case, the noise & interference covariance is always full-rank, so we conclude that most practical solutions will achieve the power bound. However, for numerical purposes, these are worthwhile observations, because computation might result in a rank-deficient estimate of $\mathbf{R_{ni}}$. We will explore this in more detail in Section VI-D below.

Lemma 3**.**

If $\lambda>0$, $\tilde{\mu}\neq 0$.

Proof:

We prove by contradiction. Assume $\tilde{\mu}=0$. When applied to (29c), we then have $\boldsymbol{\sigma}_{2}=\mathbf{C}_{V}\boldsymbol{\beta}_{2}$. Premultiplying with $\boldsymbol{\beta}_{2}^{H}$, we have $\boldsymbol{\beta}_{2}^{H}\boldsymbol{\sigma}_{2}=\boldsymbol{\beta}_{2}^{H}\mathbf{C}_{V}\boldsymbol{\beta}_{2}\geq 0$ as the quadratic form of a PSD matrix. But we’ve already established that for $\lambda>0$, $\boldsymbol{\beta}_{2}^{H}\boldsymbol{\sigma}_{2}<0$, which is a contradiction and our proof is complete. ∎

Lemma 4**.**

In a power-bounded solution, $\operatorname{rank}(\mathbf{B})=N-\operatorname{rank}(\lambda\mathbf{I}_{N}-\boldsymbol{\Sigma}_{2}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\boldsymbol{\Sigma}_{2})$.

Proof:

Our above note on Proposition 1 is our starting point. Since the matrix product $\mathbf{B}\boldsymbol{\Sigma}$ is zero, so is its rank. This further implies $\operatorname{rank}(\mathbf{B})=J-\operatorname{rank}(\boldsymbol{\Sigma})$. In general, since $\boldsymbol{\Sigma}$ is PSD, $\operatorname{rank}(\boldsymbol{\Sigma})=\operatorname{rank}(\mathbf{R_{ni}})+\operatorname{rank}(\lambda\mathbf{I}_{N}-\boldsymbol{\Sigma}_{2}^{H}\mathbf{R}_{\mathbf{ni}}^{\dagger}\boldsymbol{\Sigma}_{2})$. Therefore, generally, $\operatorname{rank}(\mathbf{B})=(NML-\operatorname{rank}(\mathbf{R_{ni}}))+(N-\operatorname{rank}(\lambda\mathbf{I}_{N}-\boldsymbol{\Sigma}_{2}^{H}\mathbf{R}_{\mathbf{ni}}^{\dagger}\boldsymbol{\Sigma}_{2}))$. In a power bounded solution, $\operatorname{rank}(\mathbf{R_{ni}})=NML$ and the inverse exists, thus $\operatorname{rank}(\mathbf{B})=N-\operatorname{rank}(\lambda\mathbf{I}_{N}-\boldsymbol{\Sigma}_{2}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\boldsymbol{\Sigma}_{2})$. ∎

VI-B A connection with waterfilling

Using some of these general results from above, we can show that the optimal solution to the KKTs at the power bound satisfies equations that resemble the well-known ”waterfilling” concept.

First, we begin with a decomposition of the clutter matrix $\mathbf{C}_{V}$. Recall that the rank of the clutter ”subspace” (i.e., the rank of $\mathbf{R}_{\mathtt{c}}(\mathbf{s})$ and thus $\mathbf{C}_{V}$)is limited by both the physical extent and nature of non-target scatterers and our overall ability to observe this state of nature. For side-looking airborne arrays, the well-known Brennan rule [25, 4] is a reasonable approximation if certain conditions hold, but as [26] showed, a more robust result is obtained by applying the Landau-Pollak theorem. In any case, let us assume that the ”true” rank of the clutter is $Q_{\mathrm{eff}}\leq Q\leq N^{2}ML$ (the subscript denoting the effective number of clutter patches). Then, we can use the economy eigendecomposition to find two equivalent representations of $\mathbf{C}_{V}$, namely:

[TABLE]

Here, $\breve{\mathbf{U}}_{C}\in\mathbb{C}^{N^{2}ML\times Q_{\mathrm{eff}}}$ is the matrix that forms the basis for the $Q_{\mathrm{eff}}$-dimensional vectorized clutter subspace, whose $i$th column is the eigenvector $\breve{\mathbf{u}}_{i}\in\mathbb{C}^{N^{2}ML}$, $i\in\{1,\,\dots,\,Q_{\mathrm{eff}}\}$. $\mathbf{D}_{C}\in\mathbb{C}^{Q_{\mathrm{eff}}\times Q_{\mathrm{eff}}}$ is a diagonal matrix whose $(i,i)$th element is the nonzero eigenvalue $\nu_{i}$. We can also consider the full eigendecomposition $\mathbf{C}_{V}=\breve{\mathbf{U}}\mathbf{D}\breve{\mathbf{U}}$. Here, the unitary matrix $\breve{\mathbf{U}}=[\breve{\mathbf{U}}_{C}\;\breve{\mathbf{U}}_{N}]$, where $\breve{\mathbf{U}}_{C}$ is as above and $\breve{\mathbf{U}}_{N}$ collects the $N^{2}ML-Q_{\textrm{eff}}$ eigenvectors in the nullspace of $\mathbf{C}_{V}$, which we can regard as the vectors $\breve{\mathbf{u}}_{i},i\in\{Q_{\mathrm{eff}}+1,\,\dots,\,N^{2}ML\}$. $\mathbf{D}$ is just the direct sum of $\mathbf{D}_{C}$ and a $N^{2}ML-Q_{\mathrm{eff}}\times N^{2}ML-Q_{\mathrm{eff}}$ all-zeros matrix. This decomposition also provides us with an alternative representation of the signal-dependent clutter covariance matrix $\mathbf{R}_{\mathtt{c}}(\mathbf{s})$. Let us assume that there is a set of $N^{2}ML$ matrices $\mathbf{U}_{i}\in\mathbb{C}^{NML\times N}$ whose vectorizations are $\operatorname{vec}(\mathbf{U}_{i})=\breve{\mathbf{u}}_{i}$ – that is, they correspond to the eigenvectors of $\mathbf{C}_{V}$. Then, the signal-dependent clutter covariance can also be given by

[TABLE]

We note that $\mathbf{U}_{i}^{H}\mathbf{U}_{i}=\frac{1}{N}\mathbf{I}_{N}$ for all $i$, not just those in the clutter indices. Thus, $\sqrt{N}\mathbf{U}_{i}$ is a rank-$N$ partial isometry.

Now, given this formulation, we can show a waterfilling-like effect by combining the matrix slackness condition (CS3) and one of the Lagrangian conditions. Recall that the Lagrangian requires $\boldsymbol{\sigma}_{2}=\mathbf{C}_{V}\boldsymbol{\beta}_{2}-\tilde{\mu}^{*}\boldsymbol{\tau}$. In matrix form, and using the set of partial isometries, this expands to

[TABLE]

Plugging this into (CS3) gives us

[TABLE]

Let us assume that we are power-bounded and everything that implies from the lemmas above. Thus, we can “directly” find $\mathbf{B}_{2}$ by applying $\mathbf{R}_{\mathbf{ni}}^{-1}$ above to obtain:

[TABLE]

From here, we make a judicious guess. Premultiplying both sides with another normalized channel matrix $\mathbf{U}_{j}^{H}$ and taking the trace gives us

[TABLE]

We can extract the $j$th term from the sum and collect it on the left hand side to obtain

[TABLE]

Now, the second term on the left-hand side of the above equation could be divided out as long as we had a guarantee it was always positive. First, if $\nu_{j}=0$, this term is 1, which is positive. Otherwise, $\nu_{j}>0$ since they are the non-zero eigenvalues of a positive semidefinite matrix. Next, we need examine the trace statement. $\mathbf{B}_{3}$ is positive semidefinite by construction, and $\mathbf{U}_{j}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{U}_{j}$ is positive definite because each $\mathbf{U}_{j}$ is a full rank scaled partial isometry and $\mathbf{R_{ni}}$ is positive definite (see [20, Theorem 7.7.2]). Thus, the trace of this matrix product is always positive and real. With this in hand, we can now say

[TABLE]

Since $\mathbf{B}_{3}$ is a relaxed version of $\mathbf{s}\mathbf{s}^{H}$, we can regard the expression $\nu_{j}\,\operatorname{tr}(\mathbf{U}_{j}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{U}_{j}\mathbf{B}_{3})$ to be a relaxed form of $\nu_{j}\,\mathbf{s}^{H}\mathbf{U}_{j}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{U}_{j}\mathbf{s}$, which is (effectively) a clutter-to-noise-and-interference ratio for the $j$th basis matrix. Similarly, $\operatorname{tr}(\mathbf{U}_{j}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{T}\mathbf{B}_{3})$ is a joint target-and-clutter to noise-and-interference ratio, and the cross term $\nu_{i}\,\operatorname{tr}(\mathbf{U}_{j}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{U}_{i}\mathbf{B}_{3})$ captures the ratio of patch-to-patch interaction and the noise-and-interference level. Thus, the first term on the right hand side is effectively a normalized measure of the clutter-matched target spectrum given a signal basis $\mathbf{B}_{3}$, and the right hand side is a normalized measure of the intraclutter spectrum given that same basis.

Traditional waterfilling dictates that power is preferentially injected to bands where the overall signal-to-noise ratio is high, proceeding to “worse”-off bands until the available power is exhausted. This is somewhat inverted in (41) because the left-hand side is an unscaled version of the subspace alignment between the solution $\mathbf{B}_{2}$ and the $j$th canonical clutter transfer matrix $\mathbf{U}_{j}$. Naively, we would like to minimize this alignment for $j\in\{1,\dots,Q_{\mathrm{eff}}\}$, since aligning with the clutter would nominally degrade our matching of the target. However, the coupled nature of (41) requires more nuance than that. First, observe that for the non-clutter directions, i.e., $j\in\{Q_{\mathrm{eff}}+1,...,N^{2}ML\}$, $\nu_{j}=0$ and (41) becomes

[TABLE]

This means that the solution’s alignment in the non-clutter directions follows that of the whitened target’s, less the combined crossover between non-clutter and clutter in the whitened spectrum. If the target is strong in these directions, we are all set, because the solution should align towards them. Otherwise, the solution must match the target in-spectrum as near as possible in directions both where the whitened clutter power is the lowest and the alignment with other directions is minimized. This corresponds to the findings in [27], who showed a similar behavior in two-step mutual information waveform design over consecutive transmit epochs. In this case, the “two steps” can be regarded as the trans-receive pair instead of sequential temporal designs.

But what of the value $\tilde{\mu}$? This clearly relates to the total available resources – in this case, $\kappa$ and $P_{o}$ – and provides us with the “water” in waterfilling. We can find a form of $\tilde{\mu}$ by reexamining (38) as follows. If we premultiply by the target $\mathbf{T}^{H}$, we obtain

[TABLE]

Taking the trace of both sides, and recognizing that a feasible solution satisfies $\operatorname{tr}(\mathbf{T}^{H}\mathbf{B}_{2})=\kappa^{*}$, we now have

[TABLE]

This is where the traditional waterfilling appears, since we are saying that the gain across the target (which we know in this case to be $\kappa^{*}$) is the upper bound on the available resources (given by the first term) minus the overall impact of the clutter weighted by its alignment (the second term). Continuing on, we can rearrange this form to be

[TABLE]

Since according to the above argument, $\operatorname{tr}(\mathbf{T}^{H}\mathbf{R}_{\mathbf{ni}}^{-1}\mathbf{T}\mathbf{B}_{3})$ is never zero, $\tilde{\mu}^{*}$ is therefore:

[TABLE]

We can insert this into (41), collect terms, and solve for a matrix equation that we will leave for later study. Thus, the optimal solution to the relaxed problem describes a generalized whiten-and-match trans-receive filter process that exhibits waterfilling behavior, shaping the transmit process to match the target’s response in clutter as much as possible.

VI-C A Potential Power-Bounded Solution

If $\mathbf{C}_{V}$ is less than full rank (which, in all practical scenarios, is always true), there is a potential path to a solution. First, assume that $\mathbf{R_{ni}}\succ 0$, so $\lambda>0$. Let us propose that the optimal vectorized corner matrix $\boldsymbol{\beta}_{2}$ is completely orthogonal to the clutter, i.e. $\mathbf{C}_{V}\boldsymbol{\beta}_{2}=\mathbf{0}_{N^{2}ML\times 1}$. A natural requirement for this to be a feasible solution is that $\boldsymbol{\tau}$ is not completely embedded in the clutter spectrum, i.e. $\boldsymbol{\tau}\notin\mathcal{R}(\mathbf{C}_{V})$. First, this means that $\boldsymbol{\sigma}_{2}=-\tilde{\mu}^{*}\boldsymbol{\tau}$, or $\boldsymbol{\Sigma}_{2}=-\tilde{\mu}^{*}\mathbf{T}$. Since we know that $\boldsymbol{\beta}_{2}^{H}\mathbf{C}_{V}\boldsymbol{\beta}_{2}=\tilde{\mu}^{*}\kappa-\lambda P_{o}$, a nullspace solution implies $\lambda=\frac{\tilde{\mu}^{*}\kappa}{P_{o}}$. Solving directly for $\tilde{\mu}$, we have $\tilde{\mu}=\frac{\operatorname{tr}(\mathbf{B}_{1}\mathbf{R_{ni}})}{\kappa^{*}}$.

Applied to Condition CS2, we then obtain $\mathbf{B}_{2}=\frac{\tilde{\mu}^{*}}{\lambda}\mathbf{B}_{1}\mathbf{T}$, which, given the substitution above implies $\mathbf{B}_{2}=\frac{P_{o}}{\kappa}\mathbf{B}_{1}\mathbf{T}$. In order for this to be completely orthogonal to the clutter, \mathbf{C}_{V}\operatorname{vec}(\mathbf{B}_{1}\mathbf{T})=\mathbf{C}_{V}(\mathbf{T}^{T}\raisebox{1.0pt}{>\otimes>}\mathbf{I}_{NML})\operatorname{vec}(\mathbf{B}_{1})=\mathbf{0}_{N^{2}ML\times 1}.

Following on to Condition CS4, we also have $\mathbf{B}_{3}=\frac{P_{o}^{2}}{\lvert\kappa\rvert^{2}}\mathbf{T}^{H}\mathbf{B}_{1}\mathbf{T}$. Since scalar feasibility requires that $\operatorname{tr}(\mathbf{B}_{3})=P_{o}$, this implies that $\operatorname{tr}(\mathbf{T}^{H}\mathbf{B}_{1}\mathbf{T})=\frac{\lvert\kappa\rvert^{2}}{P_{o}}$, which also satisfies the equality constraints. These results can also be shown to satisfy Conditions CS1 and CS3, as well as all of the positive semidefiniteness conditions, so long as $\operatorname{tr}(\mathbf{B}_{1}\mathbf{R_{ni}})\leq\frac{\lvert\kappa\rvert^{2}}{{\|\mathbf{T}\|^{2}_{F}}}\frac{\operatorname{tr}(\mathbf{R_{ni}})}{P_{o}}$. If the problem satisfies Slater’s condition, we also have, directly, $\operatorname{tr}(\mathbf{B}_{1}\mathbf{R_{ni}})<\operatorname{tr}(\mathbf{R_{ni}})$, so the previous necessary condition is immediate. Hence, we have as a possible incomplete solution

[TABLE]

where $\mathbf{B}_{1}$ minimizes $\operatorname{tr}(\mathbf{B}_{1}\mathbf{R_{ni}})$ subject to the requirements

[TABLE]

The eagle-eyed reader will note correctly that this sequence of equations resembles another optimization problem. Indeed, we can reframe the implicit matrix completion in (47) as the solution to

[TABLE]

or, perhaps, after collapsing the nullspace requirement into its necessary constituent parts,

[TABLE]

This is a standard semidefinite program so long as we use the Hermitian part of $\mathbf{T}\boldsymbol{\Gamma}_{q}^{H}$, which should be a semidefinite matrix, for all $q$. For now, we will not directly solve this subproblem, but we will note that some preliminary numerical analysis has shown that this is a reasonable path forward for the future.

VI-D Non-Power-Bounded Solutions

As mentioned above, there may be situations when $r_{\mathtt{NI}}=\operatorname{rank}(\mathbf{R_{ni}})<NML$ – for example, if we actually use an estimate $\mathbf{\hat{R}_{ni}}$ of the noise-and-interference covariance, or if we consider a noise-free case for analysis. In these cases, we know from Proposition 1 that $\lambda=0$. This implies that the solution is not power bounded, and $\operatorname{tr}(\mathbf{B}_{3})<P_{o}$. Furthermore, due to Lemma 2, $\boldsymbol{\Sigma}_{2}=\mathbf{0}_{NML\times N}$. With these in mind, we can produce, at the very least, a flowchart of solution properties in the non-power-bounded case.

As a simple beginning, we note that the matrix complementary slackness conditions reduce to two:

[TABLE]

This implies that if the solution is not power-bounded, it must null the entire noise-and-interference spectrum. Hence, we know that the matrices $\mathbf{B}_{1}$ and $\mathbf{B}_{2}$ have the general form:

[TABLE]

where $\mathbf{V}_{1}$ is an arbitrary $NML\times NML$ PSD matrix, $\mathbf{V}_{2}$ is an arbitrary $NML\times N$ matrix, and $\mathbf{P}_{\mathbf{R_{ni}}}^{\perp}=\mathbf{I}_{NML}-\mathbf{R_{ni}}\mathbf{R_{ni}}^{\dagger}$ is the orthogonal projection matrix onto the nullspace of $\mathbf{R_{ni}}$. Additionally, because of the first slackness condition, we know that $\operatorname{rank}(\mathbf{B}_{1})\leq NML-r_{\mathtt{NI}}$.

This form of $\mathbf{B}_{2}$ has a feasibility consequence. Recall that a solution is feasible only if the Capon constraint $\operatorname{tr}(\mathbf{B}_{2}^{H}\mathbf{T})=\kappa\neq 0$ is satisfied. If we substitute our new form of $\mathbf{B}_{2}$ into the constraint, it becomes

[TABLE]

Since $\mathbf{V}_{2}$ is arbitrary (but not trivial), feasibility is violated when $\mathbf{P}_{\mathbf{R_{ni}}}^{\perp}\mathbf{T}=\mathbf{0}_{NML\times N}$, which only occurs when $\mathcal{R}(\mathbf{T})\subseteq\mathcal{R}(\mathbf{R_{ni}})$ (or, when vectorized, $\boldsymbol{\tau}\notin\mathcal{R}(\breve{\mathbf{R}}_{\mathbf{ni}})$). Hence, there is no feasible solution if the target is embedded in the noise & interference spectrum. We state this directly in the following Lemma.

Lemma 5**.**

If $\lambda=0$ and $\mathcal{R}(\mathbf{T})\subseteq\mathcal{R}(\mathbf{R_{ni}})$, there is no feasible solution to (VSSDR) or its equivalent problems.

Continuing on, our Lagrangian minimization becomes

[TABLE]

where $\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}^{\perp}=\mathbf{I}_{N}\otimes\mathbf{P}_{\mathbf{R_{ni}}}^{\perp}$ and $\mathbf{v}_{2}=\operatorname{vec}(\mathbf{V}_{2})$.

Let us momentarily consider the generally unrealistic case when $\mathbf{C}_{V}$ is full rank (i.e. $Q_{eff}=N^{2}ML$). Under the $\lambda=0$ hypothesis, we can directly solve for $\boldsymbol{\beta}_{2}$:

[TABLE]

However, since $\boldsymbol{\beta}_{2}\in\mathcal{N}(\breve{\mathbf{R}}_{\mathbf{ni}})$, $\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}\boldsymbol{\beta}_{2}=\mathbf{0}_{N^{2}ML\times 1}$, so we can also say that

[TABLE]

If we want this to be a feasible result (i.e. $\tilde{\mu}\neq 0$), then the whitened target $\mathbf{C}_{V}^{-1}\boldsymbol{\tau}\in\mathcal{N}(\breve{\mathbf{R}}_{\mathbf{ni}})$ (and by extension, $\boldsymbol{\tau}\in\mathcal{N}(\breve{\mathbf{R}}_{\mathbf{ni}})$). Hence, if the clutter is full rank, the target must be clear of the noise and interference (and the dimensionality of the available resources must be such that this is possible) for a feasible solution to exist. This is effectively a further restriction of Lemma 5.

Assuming $\mathbf{C}_{V}$ is less than full rank, we can find another condition on the target response for the Lagrange multipliers by premultiplying the above equation by the clutter nullspace projection matrix $\mathbf{P}_{\mathbf{C}_{V}}^{\perp}$, which becomes:

[TABLE]

This means that either $\tilde{\mu}=0$ or $\boldsymbol{\tau}\in\mathcal{R}(\mathbf{C}_{V})$. Thus, if $\boldsymbol{\tau}\notin\mathcal{R}(\mathbf{C}_{V})$ (that is, it has a non-zero component outside of the clutter), then $\tilde{\mu}=0$.

The latter scenario, which we treat first, simplifies things considerably, because now:

[TABLE]

Hence, either $\mathbf{v}_{2}$ (and thus $\boldsymbol{\beta}_{2}$) lies in $\mathcal{N}(\mathbf{C}_{V}\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}^{\perp})$ or $\mathbf{v}_{2}=\boldsymbol{\beta}_{2}=\mathbf{0}_{N^{2}ML\times 1}$. The second possibility violates the Capon constraint, therefore $\mathbf{v}_{2},\;\boldsymbol{\beta}_{2}\in\mathcal{N}(\mathbf{C}_{V}\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}^{\perp})$. However, if (for whatever reason), the target does not lie in this space, then we have reached another infeasiblity result via the Capon constraint.

If $\boldsymbol{\tau}$ does lie entirely within the clutter, then there are a few complex scenarios. First, consider a situation where $\mathcal{N}(\breve{\mathbf{R}}_{\mathbf{ni}})\subseteq\mathcal{N}(\mathbf{C}_{V})$, which means the clutter is subsumed entirely into interference and noise (and thus $\mathcal{R}(\mathbf{C}_{V})\subseteq\mathcal{R}(\breve{\mathbf{R}}_{\mathbf{ni}})$ & $Q_{eff}\leq Nr_{\mathtt{NI}}$). In this case, $\mathbf{C}_{V}\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}^{\perp}=\mathbf{0}_{N^{2}ML\times N^{2}ML}$. However, this also means that the target is now also embedded in the noise and interference. By Lemma 5, there is no feasible solution! In any other case, a solution exists if and only if $\tilde{\mu}^{*}\boldsymbol{\tau}\in\mathcal{R}(\mathbf{C}_{V}\mathbf{P}_{\breve{\mathbf{R}}_{\mathbf{ni}}}^{\perp})$.

VII Simulations and Results

VII-A The Solvers

While we have preliminarily reduced solving (SSDR) to a matrix completion problem, it is possible to solve (in some sense) the problem numerically with commercial solvers. Most of the analysis presented here is enabled by the modeling package CVX[28, 29], which permitted us to construct two equivalent representations of the problem with rather different results. In all cases, we considered only the solver package SDPT3 [30, 31].

The major difference occurs in how the quadratic form in the problem is presented to CVX. The first method (hereinafter the “QuadSolver”) uses the quad_form() function, which preserves the form $\operatorname{vec}(\mathbf{B})^{H}\tilde{\mathbf{C}}_{V}\operatorname{vec}(\mathbf{B})$. The second method (hereinafter the “NormSolver”) uses a composition of the pow_pos() and norm() functions, which implements the equivalent form ${\|\tilde{\mathbf{C}}_{V}^{1/2}\operatorname{vec}(\mathbf{B})\|^{2}}$. Though these forms are symbolically equivalent, the latter is supposedly more amenable to optimization by a conic representation which is more efficiently processed by the default solvers available to CVX. We will demonstrate that while this might be true, the solvers behave very differently under certain conditions.

A note on dimensionality

Regardless of the solver particulars, it is important to note that all presented scenarios and solutions are necessarily constrained by available computing power and memory constraints. While CVX and MATLAB combine to make a powerful tool, they are limited in the size of problems that can be solved on a regular workstation. To wit, we have found that when $NML>400$, our current workstation (AMD AthlonTM II X2 B24 processor at 2.3GHz with 8GB memory) throws out-of-memory errors during CVX’s setup phase. We anticipate performing this analysis on more robust systems in the future, which will allow us to extend it to more realistic radar scenarios.

VII-B Computational Analysis of the KKTs

In the following simulations, we used certain common parameters for ease of comparison. As noted above, computational limitations constrain our overall problem size; therefore, in this case, we assume $N=5$ fast time samples, $M=5$ array elements, and $L=16$ pulses. The radar operates on a carrier frequency of 1 GHz and transmits pulses with a bandwidth of 50 MHz. The receive array has elements spaced at half a wavelength, i.e. $d=\lambda_{o}/2$. We recognize that these are not particularly realistic parameters (at least for the dimensionality), but as stated above, computational limitations currently conspire against a more representative simulation.

Unless otherwise stated, the scenarios presented to the solver were as follows: the noise covariance matrix was a scaled correlation matrix with correlation function $\exp(-0.05|n|)$ for $n\in\{1,\dots,NML\}$. Interferers were placed at the azimuth-elevation pairs (0.3941, 0.3) radians and (-0.4951,0.3), with correlation matrices given by the Toeplitz matrix associated with the correlation function $\exp(0.2|l|)$ for $l\in\{1,\dots,NL\}$. The clutter was simulated by placing $Q=25$ patches of $P=5$ scatterers each, equally spaced in azimuth over the interval $(-\pi/2,\pi/2)$ at an elevation angle of $\tfrac{\pi}{4}$ radians. Adjacent scatterers within each patch are correlated with coefficient -0.2. In all cases, we have made the Capon constraint real.

To obtain signals and beamformers from the relaxed problem, we generated the best ”rank-1” approximation by summing the first $B$ singular vectors of the solution matrix $\mathbf{B}$ weighted by their associated singular value. In the Frobenius sense, only the largest singular vector/value pair is the best approximation, but we have found that this approximation does, in fact, sketch out the nature of the unrelaxed problem. In the results below, we will often sweep the value $B$ from 1 to $J$, capturing, in some sense, the entire space. We recognize, however, that this is merely a heuristic that in this case happens to work well.

VII-B1 Variation in $\kappa$

First, we examine varying the Capon constraint with a fixed power constraint. In these scenarios, $\kappa$ ranges from $10^{-2}$ to $NML\times 10^{4}$ while $P_{o}=10^{7}$. Before continuing, we note that Slater’s condition dictates that the duality gap should be zero so long as $|\kappa|^{2}<NMLP_{o}=4\times 10^{9}$, which means that our parameter sweep should lie within this region. For reasonable comparison, we assumed no interference in this scenario. Plotted scenarios exclude solver results where either the solver failed or the solution violated the KKTs significantly despite the solver claiming success. In practice, this means we excluded “solutions” where $\lambda<0$.

We begin with the eigenspread of the output solution matrix $\mathbf{B}$ for each solver and each value of $\kappa$. Figure 1 shows the eigenspreads produced by Quad Solver, with Figure 2 showing finer detail of the first $2N$ eigenvalues. Clearly, all scenarios produce functionally low-rank solutions – the ”true” rank never exceeds $N$. Indeed, the effective numeric rank only exceeds one in the $\kappa=1$ scenario. Since the power is constant, the peak eigenvalue is nearly same in all scenarios, but increasing $\kappa$ drives the non-peak eigenvalues lower.

This dependence on $\kappa$ essentially repeats for the NormSolver-produced solutions, as can be seen in Figure 3 (the overall eigenspread) and Figure 4 (the first $2N$ eigenvalues). The primary difference here is that all valid solutions produced have an effective rank of one, and the overall dropoff from the $N$th eigenvalue to the $N-1$th eigenvalue is not nearly as steep. This can potentially be attributed to less accumulated numerical error in NormSolver’s process, leading to a cleaner solution.

As predicted by the KKTs, increasing $\kappa$ increases the Lagrange multiplier $\lambda$ in equal measures for a fixed power level. This effect is demonstrated directly in Figures 5 and 6. We see that the returned $\lambda$ values appear to be similar for both solvers across the entire $\kappa$ parameter space which indicates they converge to similar solutions, at least in terms of noise-and-interference suppression.

Next, we examine the effectiveness of our low-rank approximation by analyzing the cost function of (12)/(BQP 1). In the following figures, we use the term ”basis length” to indicate how many of the $J$ singular vectors of $\mathbf{B}_{opt}$ are combined to form our estimate of $\mathbf{b}$. That is, if $\sigma_{\mathbf{B},i}$ is the $i$th largest singular value of the solution matrix $\mathbf{B}$ and $\mathbf{u}_{\mathbf{B},i}$ is the associated singular vector, a basis length of $B$ produces an approximate vector of $\mathbf{b}_{\mathrm{appx},B}=\sum_{i=1}^{B}\sigma_{\mathbf{B},i}\mathbf{u}_{\mathbf{B},i}$, with the associated beamformer $\mathbf{w}_{\mathrm{appx},B}$ and signal $\mathbf{s}_{\mathrm{appx},B}$ recovered from the stacked vector as in (13).

Figure 7 shows the overall sweep of basis length for the approximate solutions recovered from QuadSolver as a function of $\kappa$, with Figure 8 focusing only on the first $2N$ values. With the notable exception of $\kappa=1$, as $\kappa$ increases, the minimum original cost increases and the difference over the first $N$ basis vectors becomes less pronounced. We can regard the result at $\kappa=1$ to be a threshhold scenario where, since the effective numerical rank of the relaxed solution is closer to $N$, more of the basis is needed for better approximation.

In contrast, NormSolver produces solutions that generally perform poorly in the original cost function. As seen in Figures 9 and 10, only the $\kappa=1$ scenario replicates the general behavior seen from QuadSolver, with the NormSolver-produced minimum value exceeding its QuadSolver counterpart by nearly 20 dB. For higher values of $\kappa$, the previously observed flattening effect is in full force, indicating that this solver produces identically structured solutions that merely scale with $\kappa$.

Finally, our investigation of the effect of $\kappa$ concludes by examining the alignment between part of the relaxed solution and the clutter patches. Since we established in Proposition 1 that the submatrix $\mathbf{B}_{2}$’s left eigenspace spans the columns of $\mathbf{B}_{1}$ and its right eigenspace spans those of $\mathbf{B}_{3}$, this effectively means that the subspace alignment between $\mathbf{B}_{2}$ and each clutter transfer matrix $\boldsymbol{\Gamma}_{q}$ is a proxy for the combined ability of the trans-recieve pair to suppress clutter. We can measure subspace alignment through the cosine of their principal angles as follows. For given matrices $\mathbf{X},\mathbf{Y}\in\mathbb{C}^{p\times q}$, the cosine of the principal angle $\psi_{\mathbf{X},\mathbf{Y}}$ between the subspaces $\operatorname{span}\{\mathbf{X}\}$ and $\operatorname{span}\{\mathbf{Y}\}$ is

[TABLE]

The closer $\cos(\psi_{\mathbf{X},\mathbf{Y}})$ is to 1, the more aligned the subspaces. In Figures 11 and 12, the resulting subspace alignments for both target and clutter channel matrices are plotted for the solution from QuadSolver and NormSolver, respectively. For reference, we denote the existing alignment between the target channel matrix $T$ and each clutter channel matrix. For both solvers, variation in $\kappa$ seems to have little effect on the overall alignment. However, QuadSolver’s solution for $\kappa=1$, the solution’s subspace is slightly better aligned with the target at the cost of marginally higher alignment to the clutter. In any case, we see a significant difference between the two solvers: QuadSolver attempts to avoid any alignment with the clutter, while NormSolver attempts to match the target’s alignment spectrum (in shape, if not magnitude). This behavior in NormSolver appears to come at a cost to the overall solution-target alignment.

VII-B2 Variation in $P_{o}$

Next, we examine the scenario opposite the above: varying the power constraint while keeping the Capon constraint fixed. In these scenarios, $P_{o}$ ranges from $10^{3}$ to $10^{9}$ while $\kappa=100$. The Slater condition in all cases is easily satisfied; however, we do note that the solver failed to satisfy the KKTs for the final two values. This is generally attributable to scaling issues. For reasonable comparison, we assumed no interference in this scenario.

First, we return to the eigenspread of the solution matrix $\mathbf{B}$. Figures 13 & 14 show the total eigenspread and its first $2N$ eigenvalues for QuadSolver. As $P_{o}$ increases, we see that the peak eigenvalue increases to maintain equality at the power bound; this also corresponds with an overall ”lifting” of the eigenspectrum, which is most likely due to numerical precision limitations. In any case, the effective numerical rank of the solution is still one, while the overall effective rank of the solution is, at most, $N$.

Similarly, the total eigenspread and first $2N$ eigenvalues for the NormSolver solution are provided in Figures 15 and 16, respectively. We see similar behavior to the QuadSolver, but note the tighter spread in effectively non-zero eigenvalues with increasing $P_{o}$. This indicates that NormSolver is converging to the same solution at each $P_{o}$, scaled to match the power constraint.

As above, we also explore the effect of the constraint on the Lagrange multiplier $\lambda$, shown in Figures LABEL:fig:PCS_Quad_Lamba (for QuadSolver) and 18 (for NormSolver). In both cases, the multiplier sharply declines as the constraint becomes looser, indicating that the overall filtered noise power remains nearly the same and is mostly a structural issue.

We also reexamine the approximate rank-1 cost function described in the above section. For QuadSolver, we can see from Figures 19 (the entire sweep of $J$ basis vectors) and 20 (the first $2N$ basis vectors) that increasing the power lowers the overall cost when approximating using the first $N$ basis vectors. Outside these vectors, the cost dramatically increases for higher values of $P_{o}$. This indicates that low-rank approximation under weak power constraints is effectively unnecessary, while basis selection must be more careful in high-power scenarios.

In contrast, NormSolver’s solution remains functionally fixed in cost over the sweep (as seen in Figures 21 and 22), which agrees with the behavior seen in 16 – namely, the overall solution is minimally sensitive to the power constraint. This is most obvious when $P_{o}={10}^{6}$ , as even its “dramatic” increase is within a tenth of a decibel of its minimum.

Finally, we investigate the subspace alignment as a function of the power constraint. As shown in Figure 23, QuadSolver’s solution attempts to null the clutter spectrum entirely, with higher power corresponding to better nulling. That said, all solutions provide alignment figures less than ${10}^{-6}$, which is well below the target’s alignment and any reasonable interpretation of “aligned”. Power does not seem to affect the solution-target alignment to a significant degree.

Again, in contrast, NormSolver’s solution, seen in Figure 24, attempts to align with the target’s representation in the clutter and remains insensitive to the power constraint. This provides a final confirmation of the waterfilling interpretation above.

VII-C Interference Effects

In this scenario, we consider the impact of interference on the relaxed solution and the effect each solver type has on this impact. Using the general scenario from above with both no interference and an interferer located at $(\theta_{I},\phi_{I})=(0.3941,0.3)$ radians, we examine the impact on both subspace alignment and the overall adapted pattern.

First, we consider the subspace alignment as in Section VI. As above, the upper right corner of the solution matrix, $\mathbf{B}_{2}$, is compared against each clutter patch and the target. The subspace alignment cosine, $\cos(\psi)$, ranges from 0 to 1, with higher values indicating that the spans of the two matrices are more aligned. Figure 25 shows the alignment with no interference, while Figure 26 shows the alignment under the interference. In both cases, as in the previous sections, QuadSolver produces a solution that is effectively unaligned with the clutter, whereas NormSolver’s solution attempts to follow the preexisting target-clutter alignment. The major difference under this metric is a clear broadening and shifting of the ”mainlobe” of NormSolver’s solution-clutter alignment spectrum near the direction of the interferer and beyond. QuadSolver’s solution-clutter alignment appears generally unaffected, but we note that this does not mean the interferer has no impact, as we will soon see.

Next, we consider the traditional adapted pattern for STAP, which plots the following function

[TABLE]

That is, the adapted pattern for a given beamformer-signal pair $\mathbf{w}_{o},\mathbf{s}_{o}$ is a function of the Doppler frequency $f_{d}$ and the azimuth $\theta$ at a given elevation $\phi$. We consider the same scenarios as above, with and without interference, and provide the overall adapted pattern as well as Doppler cuts at the target azimuth and azimuth cuts at the target Doppler frequency.

With no interference, the naive implementation performs relatively well, as illustrated in Figures 27 (for QuadSolver) and 28 (for NormSolver). The aforementioned clutter nulling peculiarities specific to each solver reappear, but the target is well localized in each case, with the peak of the pattern at its location.

This behavior is confirmed by the cuts along the target Doppler (Figures 29 and 30) and target azimuth (Figures 31 and 32), with the target location in each denoted by the dashed line. While the azimuth plots have large sidelobes, this is mostly attributable to the limited number of antennas in this simulation ($M=5$) and not the optimization process.

Now suppose we inject a broadband interferer ”close” to the target – as a reminder, the target is at $(\theta_{t},\phi_{t})=(0.3,\tfrac{\pi}{3})$ radians and the interferer is at $(\theta_{I},\phi_{I})=(0.3941,0.3)$ radians. Figures 33 and 34 (for QuadSolver and NormSolver, respectively) show the adapted pattern under these conditions. As predicted by the subspace alignment in Figure 26, the peak in NormSolver’s pattern is shifted away from the target quite significantly. What may be surprising, however, is that QuadSolver’s adapted pattern shows the exact same error, despite the subspace alignment appearing to suggest otherwise. There are a variety of factors that contribute to this, but the most salient is that the elevation angles of clutter patches are not identical to that of the interferer. Hence, the subspaces being nulled by QuadSolver do project into the adapted pattern as nulls, but are not identical to the subspace spanned by the interference. This lends credence to our proposition that the relaxed optimization process is a trans-recieve generalized whiten-and-match filter. In both cases, the resources available permit us to only whiten the interference, especially because it is not signal-dependent.

For completeness, we also show the adapted pattern cuts under the interference scenario. Again, in all figures, the dashed line represents the target location. Figures 35 and 36 show the cut along the target Doppler. The significant bias seen above is more clear here, which would result in a rather large angle estimation error.

The cuts along the target azimuth, shown in Figures 37 and 38, still appropriately localize the target in some sense. However, note that in both cases, the peak gain at the target Doppler has dropped about 5 dB from the gain in Figures 31 and 32. Thus, there is a clear loss in both dimensions due to the interferer.

We note that this interference impact (in all cases) is lessened by increasing $M$ in our simulations. This is not surprising, since more antenna elements over the same aperture increases the degrees of freedom available to null interference and localize the target.

VIII Conclusions

In this report, we reconsidered the problem of [1] under the relaxed biquadratic program framework. We demonstrated definitively that the problem is non-convex, then showed the relaxation process. We then showed that the KKTs require power-bounded solutions when the noise-and-interference matrix is full rank, and that such a solution admits a waterfilling interpretation. Simulations demonstrated that numerical solvers can provide divergent solution paths, depending on scaling – one closer to the traditional clutter-nulling process, the other relying on a matched target-in-clutter response. Future work will attempt to generalize the findings made in this report and demonstrate its utility for other radar system models that admit a channel response representation.