On the Design of Constant Modulus Probing Waveforms with Good Correlation Properties for MIMO Radar via Consensus-ADMM Approach

TL;DR

This paper develops a novel consensus-ADMM based method to design constant modulus MIMO radar waveforms with optimized correlation properties, ensuring convergence and computational efficiency.

Contribution

It introduces a new non-convex consensus optimization framework and customized ADMM algorithms with theoretical convergence guarantees for waveform design.

Findings

Proposed algorithms effectively optimize waveform correlation properties.

Algorithms converge under proper parameter settings.

Variants improve computational speed and efficiency.

Abstract

In this paper, we design constant modulus probing waveforms with good correlation properties for collocated multi-input multi-output (MIMO) radar systems. The main content is as follows: first, we formulate the design problem as a fourth order polynomial minimization problem with constant modulus constraints. Then, by exploiting introduced auxiliary variables and their inherent structures, the polynomial optimization model is equivalent to a non-convex consensus minimization problem. Second, a customized alternating direction method of multipliers (ADMM) algorithm is proposed to solve the non-convex problem approximately. In the algorithm, all the subproblems can be solved analytically. Moreover, all subproblems except one subproblem can be performed in parallel. Third, we prove that the customized ADMM algorithm is theoretically-guaranteed convergent if proper parameters are chosen. Fourth, two variant ADMM algorithms, based on stochastic block coordinate descent and accelerated gradient descent, are proposed to reduce computational complexity and speed up the convergence rate. Numerical examples show the effectiveness of the proposed consensus-ADMM algorithm and its variants.

Tables

Table 1. TABLE I: The proposed consensus-ADMM algorithm

Initialize: compute Lipschitz constants $\left\{L_n,n\in 𝒯\right\}$ and

$L$. Set iteration index $k=1$, initialize ${𝚽}^1$ and $\left\{{𝚲}_n^1\right\}$,

and let $\left\{{𝚽}^1={𝚽}_n^1,n\in 𝒯\right\}$.

Repeat

Step 1: compute ${\alpha }^{k+1}$ and ${𝚽}^{k+1}$ via (26a) and (26b).

Step 2: update $\left\{{𝚽}_n^{k+1},n\in 𝒯\right\}$ and $\left\{{𝚲}_n^{k+1},n\in 𝒯\right\}$

via (26c) and (23c) respectively in parallel.

Until some preset termination criterion is satisfied.

Table 2. TABLE II: Comparison of execution time (second per iteration).

$(N,M,𝒯)$	CA	D-ADMM	L-BFGS	consensus -ADMM	consensus-ADMM -SBCD-$25\%$	consensus -ADMM-AGD
(64,8,[0,16])	0.34s	316s	0.12s	0.08s	0.02s	0.09s
(128,8,[0,16])	0.67s	950s	0.14s	0.09s	0.02s	0.10s
(128,16,[0,32])	3.4s	N/A	0.24s	0.26s	0.08s	0.28s
(​1024,128,[0,256])	N/A	N/A	5.1s	9.1s	2.3s	11.7s

Equations

$${\mathbf{a}}_{\theta }={\left[1,e^{j\pi \sin \theta },\cdots ,e^{j\pi (M-1)\sin \theta }\right]}^T.$$

$$\mathbf{X}=[{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_M].$$

$${\mathbf{s}}_{\theta }=\mathbf{X}{\mathbf{a}}_{\theta }.$$

$$P_{\theta }={\mathbf{a}}_{\theta }^H{\mathbf{X}}^H\mathbf{X}{\mathbf{a}}_{\theta }.$$

$$\begin{split}&\qquad\qquad\quad n~{}\text{\rm zeros}\\ &{\bf{S}}_{n}=\left[\begin{array}[]{cccccc}{}\hfil&\overbrace{0~{}\cdots~{}0}&1&{}\hfil&\mathbf{\scalebox{1.5}{0}}&{}\hfil\\ {}\hfil&{}\hfil&{}\hfil&\ddots&{}\hfil&{}\hfil\\ {}\hfil&{}\hfil&{}\hfil&{}\hfil&1&{}\hfil\\ {}\hfil&\mathbf{\scalebox{2.5}{0}}&{}\hfil&{}\hfil&{}\hfil&{}\hfil\\ \end{array}\right].\end{split}$$

$$P_{{\theta }_i,{\theta }_j,n}={\mathbf{a}}_{{\theta }_i}^H{\mathbf{X}}^H{\mathbf{S}}_n\mathbf{X}{\mathbf{a}}_{{\theta }_j},$$

$$\begin{array}{rl}& \underset{\alpha ,\mathbf{X}}{\min }e(\alpha ,\mathbf{X})+P_c(\mathbf{X}),\\ & \mathrm{subject}\mathrm{to}|x_{i,m}|=1,i=1,\cdots ,N,m=1,\cdots ,M,\\ & \alpha \in (0,{\alpha }_{\max }],\end{array}$$

$$\begin{array}{rl}& e(\alpha ,\mathbf{X})=\sum _{\theta \in \Theta }|\alpha {\bar{P}}_{\theta }-{\mathbf{a}}_{\theta }^H{\mathbf{X}}^H\mathbf{X}{\mathbf{a}}_{\theta }{|}^2,\end{array}$$

$$\begin{array}{r}{P}_c(\mathbf{X})=\sum _{n\in \mathcal{T}\backslash 0}\sum _{{\theta }_i\in \hat{\Theta }}{w}_{\mathrm{ac}}^2|P_{{\theta }_i,{\theta }_i,n}{|}^2+\sum _{n\in \mathcal{T}}\sum _{\genfrac{}{}{0ex}{}{{\theta }_i\ne {\theta }_j}{{\theta }_i,{\theta }_j\in \hat{\Theta }}}{w}_{\mathrm{cc}}^2|P_{{\theta }_i,{\theta }_j,n}{|}^2,\end{array}$$

$$\begin{array}{rl}& \underset{\alpha ,\mathbf{\Phi }}{\min }e\left(\alpha ,\mathbf{X}(\mathbf{\Phi })\right)+P_c\left(\mathbf{X}(\mathbf{\Phi })\right),\\ & \mathrm{subject}\mathrm{to}\alpha \in (0,{\alpha }_{\max }],0⪯\mathbf{\Phi }\prec 2\pi ,\end{array}$$

$$\begin{array}{rl}& {\mathbf{a}}_{\theta ,\theta }=\mathrm{vec}({\mathbf{a}}_{\theta }{\mathbf{a}}_{\theta }^H),p=\sum _{\theta \in \Theta }{\bar{P}}_{\theta },\\ & \mathbf{q}=-\sum _{\theta \in \Theta }{\bar{P}}_{\theta }{\mathbf{a}}_{\theta ,\theta },\mathbf{A}=\sum _{\theta \in \Theta }{\mathbf{a}}_{\theta ,\theta }{\mathbf{a}}_{\theta ,\theta }^H.\end{array}$$

$$e(\alpha ,\mathbf{X}(\mathbf{\Phi }))={\mathbf{v}}^H(\alpha ,\mathbf{\Phi })\mathbf{Qv}(\alpha ,\mathbf{\Phi }),$$

$$\mathbf{v}(\alpha ,\mathbf{\Phi })$$

$$\mathbf{Q}$$

$${\mathbf{B}}_n(\mathbf{\Phi })=\left[\begin{array}{ccccc}& 0& w_{\mathrm{cc}}{P}_{{\theta }_1,{\theta }_2,n}& \cdots & w_{\mathrm{cc}}{P}_{{\theta }_1,{\theta }_K,n}\\ & w_{\mathrm{cc}}{P}_{{\theta }_2,{\theta }_1,n}& 0& \cdots & w_{\mathrm{cc}}{P}_{{\theta }_2,{\theta }_K,n}\\ & ⋮& ⋮& \ddots & ⋮\\ & w_{\mathrm{cc}}{P}_{{\theta }_K,{\theta }_1,n}& w_{\mathrm{cc}}{P}_{{\theta }_K,{\theta }_2,n}& \cdots & 0\end{array}\right],$$

$${\mathbf{B}}_n(\mathbf{\Phi })=\left[\begin{array}{ccccc}& w_{\mathrm{ac}}{P}_{{\theta }_1,{\theta }_1,n}& w_{\mathrm{cc}}{P}_{{\theta }_1,{\theta }_2,n}& \cdots & w_{\mathrm{cc}}{P}_{{\theta }_1,{\theta }_K,n}\\ & ⋮& ⋮& ⋮& ⋮\\ & w_{\mathrm{cc}}{P}_{{\theta }_K,{\theta }_1,n}& w_{\mathrm{cc}}{P}_{{\theta }_K,{\theta }_2,n}& \cdots & w_{\mathrm{ac}}{P}_{{\theta }_K,{\theta }_K,n}\end{array}\right].$$

$$P_c(\mathbf{X}(\mathbf{\Phi }))=\sum _{n\in \mathcal{T}}\Vert {\mathbf{B}}_n(\mathbf{\Phi }){\Vert }_F^2.$$

$$h(\alpha ,\mathbf{\Phi })={\mathbf{v}}^H(\alpha ,\mathbf{\Phi })\mathbf{Qv}(\alpha ,\mathbf{\Phi }),$$

$$f_n(\mathbf{\Phi })=\Vert {\mathbf{B}}_n(\mathbf{\Phi }){\Vert }_F^2,$$

$$\begin{array}{rl}& \underset{\alpha \in \text{mathds}R,\left\{\mathbf{\Phi },{\mathbf{\Phi }}_n\right\}\in \text{mathds}{R}^{N\times M}}{\min }h(\alpha ,\mathbf{\Phi })+\sum _{n\in \mathcal{T}}{f}_n({\mathbf{\Phi }}_n),\\ & \mathrm{subject}\mathrm{to}\mathbf{\Phi }={\mathbf{\Phi }}_n,n\in \mathcal{T},\\ & \alpha \in (0,{\alpha }_{\max }],0⪯\mathbf{\Phi }\prec 2\pi .\end{array}$$

$$\begin{array}{rl}& \mathcal{L}(\alpha ,\mathbf{\Phi },\{{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n,n\in \mathcal{T}\})\\ =& h(\alpha ,\mathbf{\Phi })+\sum _{n\in \mathcal{T}}\left(f_n({\mathbf{\Phi }}_n)+⟨{\mathbf{\Lambda }}_n,{\mathbf{\Phi }}_n-\mathbf{\Phi }⟩+\frac{{\rho }_n}{2}\Vert {\mathbf{\Phi }}_n-\mathbf{\Phi }{\Vert }_F^2\right),\end{array}$$

$$\begin{array}{rl}& {\mathcal{L}}_n\left(\mathbf{\Phi },{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n\right)\\ =& f_n({\mathbf{\Phi }}_n)+⟨{\mathbf{\Lambda }}_n,{\mathbf{\Phi }}_n-\mathbf{\Phi }⟩+\frac{{\rho }_n}{2}\Vert {\mathbf{\Phi }}_n-\mathbf{\Phi }{\Vert }_F^2,n\in \mathcal{T}.\end{array}$$

$$\{{\alpha }^{k+1},{\mathbf{\Phi }}^{k+1}\}=\underset{\genfrac{}{}{0ex}{}{\alpha \in (0,{\alpha }_{\max }],}{0⪯\mathbf{\Phi }\prec 2\pi }}{\mathrm{argmin}}\mathcal{L}\left(\alpha ,\mathbf{\Phi },\{{\mathbf{\Phi }}_n^k,{\mathbf{\Lambda }}_n^k,n\in \mathcal{T}\}\right),$$

$${\mathbf{\Phi }}_n^{k+1}=\underset{{\mathbf{\Phi }}_n}{\mathrm{argmin}}{\mathcal{L}}_n\left({\mathbf{\Phi }}^{k+1},{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n^k\right),n\in \mathcal{T},$$

$${\mathbf{\Lambda }}_n^{k+1}={\mathbf{\Lambda }}_n^k+{\rho }_n({\mathbf{\Phi }}_n^{k+1}-{\mathbf{\Phi }}^{k+1}),n\in \mathcal{T},$$

$$\Vert {\nabla }_{\alpha }h(\alpha ,\mathbf{\Phi })-{\nabla }_{\alpha }h(\hat{\alpha },\mathbf{\Phi }){\Vert }_F\le L_{\alpha }|\alpha -\hat{\alpha }|,$$

$$\Vert {\nabla }_{\mathbf{\Phi }}h(\alpha ,\mathbf{\Phi })-{\nabla }_{\mathbf{\Phi }}h(\alpha ,\hat{\mathbf{\Phi }}){\Vert }_F\le L\Vert \mathbf{\Phi }-\hat{\mathbf{\Phi }}{\Vert }_F,$$

$$\Vert \nabla f_n(\mathbf{\Phi })-\nabla f_n(\hat{\mathbf{\Phi }}){\Vert }_F\le L_n\Vert \mathbf{\Phi }-\hat{\mathbf{\Phi }}{\Vert }_F,n\in \mathcal{T},$$

$$L_{\alpha }\ge 2p,$$

$$L\ge 4(M-1)({\alpha }_{\max }{\bar{P}}_{\max }+M^{2}N+2M-2)|\Theta |,$$

$$L_n\ge 2w_{\mathrm{c}}^2(2M-1)(M^{2}N+2M-1)K^2.$$

$$\mathcal{L}(\alpha ,\mathbf{\Phi },\{{\mathbf{\Phi }}_n^k,{\mathbf{\Lambda }}_n^k,n\in \mathcal{T}\})\le \mathcal{U}(\alpha ,\mathbf{\Phi },\{{\mathbf{\Phi }}_n^k,{\mathbf{\Lambda }}_n^k,n\in \mathcal{T}\}),$$

$${\mathcal{L}}_n\left({\mathbf{\Phi }}^{k+1},{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n^k\right)\le {\mathcal{U}}_n({\mathbf{\Phi }}^{k+1},{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n^k),n\in \mathcal{T},$$

$$\begin{array}{rl}& \mathcal{U}(\alpha ,\mathbf{\Phi },\{{\mathbf{\Phi }}_n^k,{\mathbf{\Lambda }}_n^k,n\in \mathcal{T}\})\\ & \triangleq h({\alpha }^k,{\mathbf{\Phi }}^k)+⟨{\nabla }_{\mathbf{\Phi }}h({\alpha }^k,{\mathbf{\Phi }}^k),\mathbf{\Phi }-{\mathbf{\Phi }}^k⟩\\ & +⟨{\nabla }_{\alpha }h({\alpha }^k,{\mathbf{\Phi }}^k),\alpha -{\alpha }^k⟩+\frac{L}{2}\Vert \mathbf{\Phi }-{\mathbf{\Phi }}^k{\Vert }_F^2+\frac{L_{\alpha }}{2}|\alpha -{\alpha }^k{|}^2\\ & +\sum _{n\in \mathcal{T}}{\mathcal{L}}_n(\mathbf{\Phi },{\mathbf{\Phi }}_n^k,{\mathbf{\Lambda }}_n^k),\end{array}$$

$$\begin{array}{rl}& {\mathcal{U}}_n({\mathbf{\Phi }}^{k+1},{\mathbf{\Phi }}_n,{\mathbf{\Lambda }}_n^k)\\ & \triangleq f_n({\mathbf{\Phi }}^{k+1})+⟨\nabla f_n({\mathbf{\Phi }}^{k+1})+{\mathbf{\Lambda }}_n^k,{\mathbf{\Phi }}_n-{\mathbf{\Phi }}^{k+1}⟩\\ & +\frac{{\rho }_n+L_n}{2}\Vert {\mathbf{\Phi }}_n-{\mathbf{\Phi }}^{k+1}{\Vert }_F^2.\end{array}$$

Full text

I Introduction

Multiple-input multiple-output (MIMO) radar system is regarded as a promising paradigm for the next generation radar systems. Unlike the standard phased-array radar to transmit scaled versions of a single waveform, probing signals, transmitted via different antennas in the MIMO radar system, are independent. Through this additional waveform diversity, MIMO radar owns superior capabilities compared with the traditional phased-array radar, such as higher spatial resolution, more flexible beampattern, and better detection performance [1][2]. MIMO radar system can be classified into two categories: distributed and collocated. In the former, transmitters are widely separated in space and each of them can provide an independent view of the target, which can improve detection performance [3][4]. In the latter, antennas in the transmitter are placed in close proximity and different probing signals from various collocated antennas can generate various desired beampatterns, leading to an improved directional resolution and interference rejection capability [5]–[7].

Probing signal waveforms play a central role in the signal processing performance of a MIMO radar system. Specifically, since matching the desired spacial beampatterns and lowering spacial correlations levels can increase spacial directional gain and eliminate clutter interference from other directions, a lot of researchers have been attracted to designing probing signal waveforms to meet these goals in recent years. Authors in [8] and [9] matched the waveform covariance matrix to the desired beampattern through a semidefinite programming method, then exploited the cyclic algorithm to synthesize the constant modulus waveform and pursued good auto-/cross-correlation properties. In [10], authors formulated the waveform design problem as a fourth order polynomial minimization problem with constant modulus constraints, then proposed a quasi-Newton solving algorithm to approximate the model’s optimal solution. Moreover, the approach can be applied to the scenario of desired low correlation sidelobe levels within certain lag intervals. The authors in [11] focused on the direct or indirect control of mainlobe ripples in the beampattern design problem. They reformulated the design as a feasibility problem with the lowest system cost. To achieve a high signal to interference plus noise ratio and low sidelobe levels performance, a fixed waveform covariance matrix was proposed in [12]. However, the matrix does not exploit the full waveform diversity. In [13], the authors proposed a novel transmit beampattern matching design one-step method, which obtains the transmit signal matrix by unconstrained optimization. The drawback of the waveforms generated by this method is that their envelope is not constant modulus. To reduce the computational complexity, a closed-form covariance matrix design method was proposed in [14] based on discrete Fourier transform (DFT). The authors in [15] and [16] also applied the DFT-based technique to a planar-antenna-array, and developed a finite-alphabet constant-envelope waveforms design algorithm for the desired beampattern. However, the performance of the DFT-based method is slightly worse for a small number of antennas. The authors in [17] studied the robust transmit beampattern design problem and exploited the semidefinite relaxation technique to treat non-convex optimization problems. In [18] [19], the authors exploited successive convex relaxation techniques to handle non-convex quadratic equality constraints in the constant modulus waveform design problem. In [20], the authors proposed a double cyclic alternating direction method of multipliers (D-ADMM) algorithm to solve the non-convex beampattern design problem and in [21], they considered the joint optimization problem of the covariance matrix and antenna position. In [22], the authors applied the majorization-minimization technique to match the desired transmit beampattern, which enjoys faster convergence than D-ADMM. The authors in [23] focused on MIMO radar waveform design under the constant modulus and similarity constraints. They proposed a sequential iterative algorithm based on the block coordinate descent (BCD) framework, which has shown its superiority compared with the CA approach in [9]. In [24], the authors considered the constant modulus waveform design to achieve a desired wideband MIMO radar beampattern with space-frequency nulling. In each algorithm iteration, the authors optimized the original non-convex problem’s approximation version meaning that the proposed algorithm can be executed in parallel. Besides the above works, some researchers synthesized transmit waveforms under some practical constraints, such as mainlobe ripple constraints [25], spectral shape constraints [26], constant modulus constraints [27], similarity constraints [28] [29], and transmitted power constraints [30]. However, these works only focus on the synthesized beampattern design problem and pay little attention to the correlation properties of the waveforms.

In this paper, we extend our previous work in [31] and propose a consensus-ADMM approach to design constant modulus probing waveforms, which can match the desired spacial beampatterns while suppressing the spacial auto-correlation and cross-correlation levels in the collocated MIMO radar system. Its main contributions are as follows.

• •

Consensus problem formulation: the design problem is formulated as a fourth order polynomial minimization problem with constant modulus constraints. Then, by introducing auxiliary variables, it is further equivalent to a non-convex consensus minimization problem.

• •

Parallel solving algorithm: consensus-ADMM is customized to solve the non-convex consensus problem approximately. In the implementation, all the subproblems can be solved analytically. Moreover, except one subproblem, all subproblems can be performed in parallel. This favourable execution architecture is the main advantage of the proposed consensus-ADMM over state-of-the-art techniques, which is very suitable for practical implementation.

• •

Theoretically-guaranteed performance: we prove that the solving algorithm is guaranteed convergent to a stationary point of the non-convex optimization problem if proper parameters are chosen.

• •

Improvement strategies: two variant ADMM algorithms, based on stochastic block coordinate descent (SBCD) and accelerated gradient descent (AGD), are proposed to reduce computational complexity and speed up the convergence rate.

The rest of the paper is organized as follows. In Section II, we formulate the beampattern design problem to a non-convex consensus minimization problem. In Section III, consensus-ADMM is customized to solve the non-convex minimization problem. The performance analysis, including convergence and computational complexity of the proposed consensus-ADMM algorithm, are presented in Section IV. Two variant algorithms, named by consensus-ADMM-SBCD and consensus-ADMM-AGD, are given to improve computational complexity and convergence performance of the solving algorithm respectively in Section V. Finally, Section VI demonstrates the effectiveness of the proposed consensus-ADMM algorithms and the conclusions are given in Section VII.

Notation: bold lowercase and uppercase letters denote column vectors and matrices and italics denote scalars. $\mathds{R}$ and $\mathds{C}$ denote the real field and the complex field respectively. The superscripts $(\cdot)^{*}$, $(\cdot)^{T}$ and $(\cdot)^{H}$ denote conjugate operator, transpose operator and conjugate transpose operator respectively. $x_{i}$ denotes the $i$-th element of vector $\mathbf{x}$. $|\cdot|$ denotes the absolute value. The subscript $\|\cdot\|_{2}$ denotes the Euclidean vector norm and $\|\cdot\|_{F}$ denotes the Frobenius matrix norm. $\nabla(\cdot)$ represents the gradient of a function. $\rm{Re}(\cdot)$ takes the real part of the complex variable and ${\rm Tr}(\cdot)$ denotes the trace of a matrix. $\langle\cdot,\cdot\rangle$ and $\otimes$ are the dot product operator and convolution operator respectively. ${\rm vec}(\cdot)$ vectorizes a matrix by stacking its columns on top of one another and ${\rm mat}(\cdot,N,M)$ reshapes a vector to an $N\times M$ matrix. $\Pi(\cdot)$ denotes the projection operator. $\mathbb{E}[\cdot]$ performs the expectation of the variables and ${\mathbf{I}}$ denotes an identity matrix.

II System Model And Problem Formulation

II-A System Model

Consider a MIMO radar system equipped with $M$ antennas in a uniform linear array as shown in Figure 1. In the system, we set the inter-element spacing $d=\frac{\lambda}{2}$, where $\lambda$ is the signal wavelength. The spacial direction $\theta$ belongs to angle set $\Theta=(-90^{\circ},90^{\circ})$, which represents the antenna scanning scope. The steering vector $\mathbf{a}\in\mathds{C}^{M}$ at direction $\theta$ is given by

[TABLE]

The probing waveform transmitted by the $m$-th antenna is denoted by $\mathbf{x}_{m}=[x_{1,m},\cdots,x_{N,m}]^{T},m=1,\cdots,M$. Then, the waveforms transmitted by the MIMO radar system can be expressed by the following $N$-by-$M$ matrix

[TABLE]

The synthesized signal at direction $\theta$ (far field) is

[TABLE]

The beampattern, which describes the power distribution at direction $\theta$, is defined as

[TABLE]

To describe the correlation properties of the probing waveforms at time slot $n$, we define a $N$-by-$N$ off-line diagonal matrix $\mathbf{S}_{n}$ as follows

[TABLE]

Through $\mathbf{S}_{n}$, the time-delayed signal can be expressed by $\mathbf{S}_{n}\mathbf{X}\mathbf{a}_{\theta}$. Then, the spacial correlation of the probing waveforms and its delayed version can be obtained by

[TABLE]

where $\theta_{i},\theta_{j}\in\hat{\Theta}\subset\Theta$ and $\hat{\Theta}=\{\theta_{1},\cdots,\theta_{K}\}$ is the considered angle set of spacial directions. Specifically, when $\theta_{i}=\theta_{j}$, $P_{\theta_{i},\theta_{i},n}$ denotes the spacial auto-correlation, otherwise $P_{\theta_{i},\theta_{j},n}$ means the spacial cross-correlation.

II-B Problem Formulation

We optimize MIMO radar probing waveforms based on the following considerations: first, as mentioned in (4), since the beampattern describes the spacial power distribution, we desire that it can match the directions of interest, which can decrease clutter components and extend the probing distance; second, since low auto-correlation sidelobes can increase spacial resolution and low cross-correlation levels can reduce interferences from other directions, we desire that the optimized probing waveforms have low auto-correlation sidelobes and low cross-correlation levels; third, in order to maximize the efficiency of the power amplifier in the MIMO radar transmitter, the probing waveforms should be constant modulus, i.e., $|x_{i,m}|=1,\ i=1,\cdots,N,\ m=1,\cdots M$.

Based on the above considerations, we formulate the following optimization model to design MIMO radar probing waveforms

[TABLE]

where

[TABLE]

and $w_{\rm ac}$ and $w_{\rm cc}$ are preset positive weights and $\mathcal{T}$ is the time delay parameter set of interest. In the objective function of model (6), the first term $e(\alpha,\mathbf{X})$ represents the mismatching square error between the designed beampattern and the desired beampattern $\bar{P}_{\theta}$ and $\alpha$ is a scaling factor that needs to be optimized. The second term $P_{c}(\mathbf{X})$ relates to the auto-correlation sidelobes and cross-correlation levels at the considered spacial directions. Because $P_{\theta_{i},\theta_{j},-n}^{*}=P_{\theta_{i},\theta_{j},n}$, correlation levels for $n<0$ are not included. It is difficult to solve (6) directly since its quartic objective function and constant modulus constraints are non-convex. In the following, by exploiting its inherent structure, we show how to design an efficient solving algorithm to pursue theoretically-guaranteed solutions.

First, let $\mathbf{X}$’s phase be new variable. Since $x_{i,m}=e^{j\phi_{i,m}}$, we can drop constant modulus constraints and rewrite (6) as the following minimization problem

[TABLE]

where the constraint $0\preceq\mathbf{\Phi}\prec 2\pi$ means all the elements in $\mathbf{\Phi}$ belong to $[0,2\pi)$.

Second, we define the following quantities

[TABLE]

Then, the first term $e(\alpha,\mathbf{X}(\mathbf{\Phi}))$ in (8) can be rewritten as

[TABLE]

where

[TABLE]

Third, to let $P_{c}(\mathbf{X})$ be in a compact expression, we define $K$-by-$K$ matrices set $\{\mathbf{B}_{n}(\mathbf{\Phi})|n\!\in\!\mathcal{T}\}$, where $K$ is spacial directions of interest, i.e., set $\hat{\Theta}$’s size. Specifically, when $n=0$,

[TABLE]

and when $n\neq 0$,

[TABLE]

Then, $P_{c}(\mathbf{X}(\mathbf{\Phi}))$ in (8) can be rewritten as

[TABLE]

To facilitate the subsequent derivations, we further define

[TABLE]

and introduce a set of auxiliary variables $\{\mathbf{\Phi}_{n}|n\in\mathcal{T}\}$. Then, problem (8) can be formulated as the following consensus-like problem [32]

[TABLE]

In comparison with (8), model (14) allows subfunction $h(\alpha,\mathbf{\Phi})$ or $f_{n}(\mathbf{\Phi}_{n})$ to handle its local variable independently when $\mathbf{\Phi}$ or $\mathbf{\Phi}_{n}$ is fixed. In the next section, an efficient algorithm, named by consensus-ADMM, is proposed to solve (14) approximately. Moreover, we prove that consensus-ADMM converges to a stationary point of model (8)111Several state-of-the-art algorithms [33]–[35], can be customized to handle problem (6) (not direct). The algorithms keep objective function’s value monotonic decreasing or increasing in the iteration procedure. This monotonic property along with some mild conditions is exploited to prove the convergence of the algorithms. Specifically, in [35], the presented algorithm, similar to our proposed consensus-ADMM approach, can also be implemented in parallel. Applying these algorithms to solving (6) can be interesting research directions..

To the best of our knowledge, it is the first time that a parallel algorithm structure is introduced to match the desired beampattern for the MIMO radar system, which means that the proposed consensus-ADMM algorithm is more suitable for the large scale MIMO radar waveforms design problem. Moreover, convergence analysis and improved variants of the proposed consensus-ADMM algorithm are also considered.

III Consensus-ADMM SOLVING ALGORITHM

The augmented Lagrangian function of problem (14) can be written as

[TABLE]

where $\mathbf{\Lambda}_{n}$ and $\rho_{n}$ are the Lagrangian multiplier and penalty parameters respectively. To facilitate discussions later, we define the following functions

[TABLE]

Based on (15) and (16), the proposed consensus-ADMM algorithm [32] can be described as

[TABLE]

where $k$ is the iteration number.

Remarks on (17): since $\mathcal{L}\left(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\}\right)$ and $\mathcal{L}_{n}\left(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k}\right)$ are non-convex, it is difficult to implement (17a) and (17b). However, we have the following lemma to characterize Lipschitz properties of $\nabla h(\alpha,\mathbf{\Phi})$ and $\nabla f_{n}(\mathbf{\Phi})$ (see proof in Appendix A).

Lemma 1

gradients $\nabla h(\alpha,\mathbf{\Phi})$ and $\nabla f_{n}(\mathbf{\Phi})$ are Lipschitz continuous, i.e.,

[TABLE]

where constants

[TABLE]

Here, $\bar{P}_{\rm max}=\max\limits_{\theta\in\Theta}\{\bar{P}_{\theta}\}$, $w_{\rm c}=\max\{w_{\rm ac},w_{\rm cc}\}$, and $|\Theta|$ denotes the size of set $\Theta$.

Based on the above Lemma, $\mathcal{L}\left(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\}\right)$ and $\mathcal{L}_{n}\left(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k}\right)$ can be upper-bounded by the following strongly convex functions [36]

[TABLE]

where

[TABLE]

and

[TABLE]

Then, instead of solving (17a) and (17b) directly, we propose the following customized consensus-ADMM algorithm (23).

[TABLE]

Notice that $\mathcal{U}(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\!\in\!\mathcal{T}\})$ is a strongly convex quadratic function with respect to $\alpha$ and $\mathbf{\Phi}$ respectively and $\mathcal{U}_{n}(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k})$ is also some strongly convex quadratic function with respect to $\mathbf{\Phi}_{n}$. Therefore, the minimizers of $\mathcal{U}(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\!\in\!\mathcal{T}\})$ and $\mathcal{U}_{n}(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k})$ can be determined through the following procedures:

Set the gradients of the functions $\mathcal{U}(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\!\in\!\mathcal{T}\})$ and $\mathcal{U}_{n}(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k})$ to be zeros

[TABLE]

which lead to the following linear equations

[TABLE]

Then, by solving the above linear equations and projecting the solutions onto the corresponding feasible regions, we can obtain

[TABLE]

Combining (23c) and (26), we summarize the proposed consensus-ADMM algorithm in Table I.

IV Analysis

IV-A Convergence Issue

We have the following theorem to show convergence properties of the proposed consensus-ADMM algorithm in Table I.

Theorem 1

Let $(\alpha^{*},\mathbf{\Phi}^{*},\{\mathbf{\Phi}_{n}^{*},\mathbf{\Lambda}_{n}^{*},n\in\mathcal{T}\})$ denote some limit point of the sequence $\left(\alpha^{k},\mathbf{\Phi}^{k},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\}\right)$. $\forall n\in\mathcal{T}$, if the penalty parameters $\rho_{n}$ and Lipschitz constants $L_{n}$ satisfy $\rho_{n}\geq 9L_{n}$, the proposed consensus-ADMM algorithm is convergent, i.e.,

[TABLE]

Moreover, $(\alpha^{*},\mathbf{\Phi}^{*})$ is a stationary point of problem (8), i.e., it satisfies the following inequalities

[TABLE]

Remarks: Theorem 1 indicates that the proposed consensus-ADMM algorithm is theoretically-guaranteed to be convergent to a stationary point of model (8) under the conditions $\rho_{n}\geq 9L_{n},n\in\mathcal{T}$. Here, we should note that these conditions are easily satisfied since we can choose $L_{n}$’s value according to (19c) in Lemma 1 and penalty parameters $\rho_{n}$ can be set accordingly to satisfy $\rho_{n}\geq 9L_{n}$. The key idea of proving Theorem 1 is to find out that potential function $\mathcal{L}(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n},n\in\mathcal{T}\})$ decreases sufficiently in every ADMM iteration and is lower-bounded. To reach this goal, we first prove several related lemmas in Appendix B. Then, we give the detailed proof of Theorem 1 in Appendix C.

IV-B Implementation Analysis

In the following, we show how to compute $\nabla_{\alpha}h(\alpha,{\mathbf{\Phi}})$, $\nabla_{\mathbf{\Phi}}h(\alpha,{\mathbf{\Phi}})$ and $\nabla f_{n}({\mathbf{\Phi}})$ efficiently by exploiting their inside structures.

1. $\nabla_{\alpha}h(\alpha,{\mathbf{\Phi}})$ and $\nabla_{\mathbf{\Phi}}h(\alpha,{\mathbf{\Phi}})$. They can be expressed as follows respectively

[TABLE]

where $\frac{\partial\mathbf{v}(\alpha,\mathbf{\Phi})}{\partial\alpha}=\left[1;\mathbf{0}\right]$ and

[TABLE]

In (30), $\frac{\partial\mathbf{v}(\alpha,\mathbf{\Phi})}{\partial\phi_{i,m}}$ can be calculated through

[TABLE]

where $i=1,\dotsb,N$, $m=1,\dotsb,M$, and $\frac{\partial(\mathbf{X}^{H}(\mathbf{\Phi})\mathbf{X}(\mathbf{\Phi}))}{\partial\phi_{i,m}}$ can be computed through (32).

Since $\mathbf{Q}\in\mathds{C}^{(M^{2}+1)\times(M^{2}+1)}$, $\mathbf{v}\in\mathds{C}^{M^{2}+1}$, and $\frac{\partial\mathbf{v}(\alpha,\mathbf{\Phi})}{\partial\alpha}=\left[1;\mathbf{0}\right]$, we can obtain $\frac{\partial\mathbf{v}^{H}(\alpha,\mathbf{\Phi})}{\partial\alpha)}$ through no more than $(M^{2}+1)^{2}$ complex multiplications. Moreover, since $\frac{\partial\mathbf{X}^{H}(\mathbf{\Phi})\mathbf{X}(\mathbf{\Phi})}{\partial\phi_{i,m}}$ involves $2(M-1)$ nonzero elements, then there are $2(M-1)MN$ nonzero elements in $\frac{\partial\mathbf{v}^{H}(\alpha,\mathbf{\Phi})}{\partial{\rm vec}(\mathbf{\Phi})}$. It means that it takes $(M^{2}+1)^{2}+(2M-1)MN$ complex multiplications to compute $\nabla_{\mathbf{\Phi}}h(\alpha,\mathbf{\Phi})$. Then, we can see that the total computation cost on $\nabla_{\alpha}h(\alpha,{\mathbf{\Phi}})$ and $\nabla_{\mathbf{\Phi}}h(\alpha,{\mathbf{\Phi}})$ is roughly $\mathcal{O}(M^{4}+2M^{2}N)$.

2. $\nabla f_{n}({\mathbf{\Phi}})$. The elements in $\nabla f_{n}({\mathbf{\Phi}})$ can be obtained through

[TABLE]

To compute $P_{\theta_{i},\theta_{j},n}$ and $\frac{\partial P_{\theta_{i},\theta_{j},n}}{\partial\phi_{i,m}}$ for every $\frac{\partial f_{n}({\mathbf{\Phi}})}{\partial\phi_{i,m}}$ efficiently, we define $\mathbf{s}_{\theta_{i}}=\mathbf{X}\mathbf{a}_{\theta_{i}}$ and $\mathbf{\bar{s}}_{\theta_{i}}$ denoting $\mathbf{s}_{\theta_{i}}$’s reversing vector. Then, $\forall n\in\mathcal{T}$, since $P_{\theta_{i},\theta_{j},n}=\mathbf{a}_{\theta_{i}}^{H}\mathbf{X}^{H}\mathbf{S}_{n}\mathbf{X}\mathbf{a}_{\theta_{j}}$, it can be obtained through convolution operation $\mathbf{s}_{\theta_{i}}^{*}\otimes\mathbf{\bar{s}}_{\theta_{j}}.$ It means that the cost to obtain all the $P_{\theta_{i},\theta_{j},n}$ is roughly $|\mathcal{T}|K^{2}MN$ complex multiplications. Moreover, corresponding gradients $\frac{\partial P_{\theta_{i},\theta_{j},n}}{\partial\phi_{i,m}}$ can be obtained through $\mathbf{s}_{\theta_{i}}^{*}\otimes\frac{\partial\bar{\mathbf{s}}_{\theta_{j}}}{\partial\phi_{i,m}}+\bar{\mathbf{s}}_{\theta_{j}}\otimes\frac{\partial\mathbf{s}_{\theta_{i}}^{*}}{\partial\phi_{i,m}}.$ Since there is only one nonzero element in either $\frac{\partial\bar{\mathbf{s}}_{\theta_{j}}}{\partial\phi_{i,m}}$ or $\frac{\partial\mathbf{s}_{\theta_{i}}^{*}}{\partial\phi_{i,m}}$, it takes only two complex multiplications to compute $\frac{\partial P_{\theta_{i},\theta_{j},n}}{\partial\phi_{i,m}}$. Then, we can obtain all of them through $2|\mathcal{T}|K^{2}MN$ complex multiplications. Taking complex multiplication operations to obtain $P_{\theta_{i},\theta_{j},n}^{*}\frac{\partial P_{\theta_{i},\theta_{j},n}}{\partial\phi_{i,m}}$ into account, we can see that the total cost of computing $\frac{\partial f_{n}({\mathbf{\Phi}})}{\partial\phi_{i,m}}$ $\forall n\in\mathcal{T}$ is roughly $\mathcal{O}(3|\mathcal{T}|K^{2}MN)$.

Observing the proposed consensus-ADMM algorithm in Table I and the corresponding (26) and (23c), we can see that the main computational cost lies in computing $\nabla_{\alpha}h(\alpha,{\mathbf{\Phi}})$, $\nabla_{\mathbf{\Phi}}h(\alpha,{\mathbf{\Phi}})$ and $\nabla f_{n}({\mathbf{\Phi}})$, which are much larger than other terms. Therefore, we conclude that the total cost in each ADMM iteration is roughly $\mathcal{O}(M^{4}+2M^{2}N+3|\mathcal{T}|K^{2}MN)$.

V Improvements

V-A Reduce Complexity

In the proposed consensus-ADMM algorithm, $\mathbf{\Phi}_{n}$, $\forall n\in\mathcal{T}$, are updated independently (or in parallel) as are the Lagrangian multipliers $\mathbf{\Lambda}_{n}$. This fact admits us to update only a part of the variables $\{\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n},n\in\mathcal{T}\}$ in each ADMM iteration to reduce computational complexity.

Specifically, consider a randomized updating strategy called stochastic block coordinate descent (SBCD) [37]. In the $k$-th iteration, let $\mathcal{N}^{k}$ denote some $\mathcal{T}$’s subset. We choose elements from $\mathcal{T}$ to construct $\mathcal{N}^{k}$ with the probability

[TABLE]

If some $n\in\mathcal{N}^{k}$, the corresponding variables $\mathbf{\Phi}_{n}^{k}$ and $\mathbf{\Lambda}_{n}^{k}$ are updated using (26c) and (23c) respectively. Otherwise, we just set $\mathbf{\Phi}_{n}^{k+1}=\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k+1}=\mathbf{\Lambda}_{n}^{k}.$

In this way, it is obvious that computational complexity in each ADMM iteration can be reduced significantly. The approach provides an option for some practical systems when their computation resources are very limited. Moreover, this kind of implementation strategy can still guarantee the algorithm converges with high probability to a stationary point of problem (8) under some wild conditions (we provide a sketch of the proof in Appendix D).

V-B Speed Up Convergence

Besides the computational complexity in each iteration, convergence speed is another concern from a practical viewpoint. In this paper, inspired by Nesterov’s accelerated gradient descent method (AGD) [38], which is originally applied to a convex problem, we develop its variant in the following

[TABLE]

where $\gamma^{k}=\frac{k-1}{k+t-1}$ and $t\geq 3$ is some preset constant. The algorithm starts from $\mathbf{\hat{\Phi}}_{n}^{1}=\mathbf{\Phi}_{n}^{1}$. Here, we should note that AGD method’s convergence can be proved under strong assumptions, such as solving convex or strongly convex optimization problem [39]–[41]. Here, the direct combination of consensus-ADMM and AGD can be cast as a heuristic method. It is difficult to prove that it can improve convergence rate theoretically for the considered non-convex probing waveform design problem [41]–[44]. However, simulation results, presented in the following section, show that its practical performance is superior to the original one in Table I.

VI Simulation results

In this section, numerical results are presented to illustrate the performance of the proposed MIMO radar beampattern design algorithm. We consider the number of antennas as $M=8,16,128$ with the length of each sequence $N=64,128,1024$ respectively. The set of the spacial angles covers $(-90^{\circ},90^{\circ})$ with spacing $0.1^{\circ}$. The residuals of the proposed consensus-ADMM algorithm in the $k$-th iteration are defined as $\displaystyle\sum_{n\in\mathcal{T}}\left\lVert\mathbf{\Phi}_{n}^{k+1}\!-\!\mathbf{\Phi}^{k+1}\right\rVert_{F}$ and $\displaystyle\sum_{n\in\mathcal{T}}\left\lVert\mathbf{\Phi}_{n}^{k+1}\!-\!\mathbf{\Phi}_{n}^{k}\right\rVert_{F}.$ The termination criteria is set as both of the residuals are less than $10^{-4}$ or the maximum iteration number 60000 is reached. The weights $(w_{\rm ac},w_{\rm cc})$ are $(10,10)$. The desired beampattern is

[TABLE]

where $\theta_{1}=-40^{\circ}$ and $\theta_{2}=30^{\circ}$. The parameter $t$ in the AGD method is $3$. The penalty parameter $\rho_{n}$ can affect ADMM algorithm’s convergence rate. For example, much larger $\rho_{n}$ will let the optimization problem become singular and slow it down. Here, we recommend its value as $\rho_{n}=\|{\rm vec}(\nabla f_{n}(\mathbf{\Phi}^{1}))\|_{\infty}$333Here, $\mathbf{\Phi}^{1}$ is initial value of the phase variable. Using this setting, simulation results are pretty good and convergence can always be observed.. For D-ADMM, the penalty parameters are set as 20. The random phase sequence is chosen to initialize all approaches. All experiments are performed in MATLAB 2016b/Windows 7 environment on a computer with 2.1GHz Intel 4110$\times$2 CPU and 64GB RAM.

Figure 2 plots the performance curves of objective function versus the iteration number for our proposed consensus-ADMM algorithms and other three state-of-the-art approaches: CA approach [9], L-BFGS method [10], and D-ADMM [20]. Besides the objective function in (6), its two parts: $e(\alpha,\mathbf{X}(\mathbf{\Phi}))$ and $P_{c}(\mathbf{X}(\mathbf{\Phi}))$ are also presented in the figures since both of them have obvious physical significance. It should be noted that the D-ADMM method only focuses on beampattern matching problem. Its spacial correlation characteristics are not included in the corresponding figure. From Figure 2, we can see that the resulting $e(\alpha,\mathbf{X}(\mathbf{\Phi}))$ is much larger than $P_{c}(\mathbf{X}(\mathbf{\Phi}))$. Different approaches for $e(\alpha,\mathbf{X}(\mathbf{\Phi}))$ have similar convergence performance. However, convergence results for $P_{c}(\mathbf{X}(\mathbf{\Phi}))$ are quite different. Specifically, the AGD strategy (35) can speed up convergence very well. Other algorithms tend to achieve the similar value to AGD method with relatively large number of iterations. In comparison, consensus-ADMM-SBCD-$25\%$’s convergence is a little bit slow. However, we should note that it has lower computational complexity. In practice, parameter $25\%$ can be adjusted to make a tradeoff between convergence rate and computational complexity.

Figure 3 shows synthesized spacial beampatterns by our proposed consensus-ADMM approaches and three other approaches. From the figure, it can be observed that all of the approaches can match the desired spacial beampattern very well at different antenna numbers and waveform lengths.

Figures 4 – 6 show the normalized spacial correlation level $C_{\theta_{i},\theta_{j},n}$ with different simulation parameters. Here, the normalized spacial correlation function $C_{\theta_{i},\theta_{j},n}$ for a certain interval in dB is defined as

[TABLE]

where $n\in\mathcal{T}$ and $\theta_{i},\theta_{j}\in\hat{\Theta}$. It is obvious that for $i=j$ and $i\neq j$, $C_{\theta_{i},\theta_{i},n}$ indicates spacial auto-/cross-correlation characteristics of the designed MIMO radar probing waveforms respectively. From the figures, we can see that the normalized spacial auto-correlation functions are symmetric and the normalized cross-correlation functions of $(-40^{\circ},30^{\circ})$ are symmetric to that of $(30^{\circ},-40^{\circ})$. The figures also indicate that either increasing $N$ and $M$ or decreasing $\mathcal{T}$ can lower the correlation levels. These facts are reasonable since larger $N$ and $M$ or smaller $\mathcal{T}$ indicate more degrees of freedom in designing probing waveforms. Moreover, we can also see that the consensus-ADMM-AGD approach enjoys the best auto-/cross-correlation characteristics. This fact is in accordance with the simulation result in Figure 2.

Table II shows averaged execution time (per iteration) of the proposed consensus-ADMM approach and three other state-of-the-art approaches [9][10][20]. Here, N/A means that iteration operations cannot be finished in reasonable time. From the table, we can see that when $M$ and $N$ are small, the proposed consensus-ADMM algorithm and its variants (SBCD-$25\%$, AGD) have less execution time. When $N$ and $M$ are increased, for example, from (64,8) to (1024,128), execution time of the consensus-ADMM approach becomes comparable to L-BFGS approach. Combining these execution time per iteration with the curves in Figure 2, the total computational time can be computed easily. Moreover, we should note that in the proposed consensus-ADMM algorithm or its variants, parallel execution architecture plays an essential role leading to better implementation efficiency than state-of-the-art methods including L-BFGS, D-ADMM, and MM, which means that it is more suitable for large-scale applications from a practical viewpoint of implementation444In the real radar system, the algorithm is usually implemented using Field Programmable Gate Array (FPGA). This kind of integrated chip is very suitable for implementing an algorithm with a parallel structure..

VII Conclusion

In this paper, we focus on designing constant modulus probing waveforms with good correlation properties for the collocated MIMO radar system. By introducing auxiliary variables and exploiting the designing problem’s inherent structures, we formulate a consensus-like optimization model. Then, the ADMM technique is customized to solve the corresponding non-convex problem approximately. We prove that the proposed ADMM approach is theoretically-guaranteed convergent if proper parameters are chosen. Simulation results show the effectiveness of the proposed consensus-ADMM algorithm and its variants, especially suitable for large-scale MIMO radar systems.

Acknowledgment

The authors appreciate Dr. Ian Clarkson, associate editor of this paper, and the anonymous reviewers who help the authors to elevate the quality of this paper.

Appendix A Proof of Lemma 1

In the following, we prove that both $\nabla h({\alpha,\mathbf{\Phi}})$ and $\nabla f_{n}(\mathbf{\Phi})$ are Lipschitz continuous via the definition of Lipschitz continuity. To state the proof clearly, we rewrite (9) and (11) in the following

[TABLE]

To facilitate the subsequent derivations, we denote

[TABLE]

Then, $\mathbf{Qv}(\alpha,\mathbf{\Phi})$ can be expressed by

[TABLE]

First, we can obtain $\frac{\partial\mathbf{v}(\alpha,\mathbf{\Phi})}{\partial\alpha}=\left[1;\mathbf{0}\right]$. Plugging it and (39) into (29a), we can have

[TABLE]

where $\alpha,\hat{\alpha}\in(0,\alpha_{\rm max}]$. From (40), we can see that $\nabla_{\alpha}h({\alpha,\mathbf{\Phi}})$ is Lipschitz continuous with constant $L_{\alpha}\geq 2p$.

Second, for $\nabla_{\mathbf{\Phi}}h({\alpha,\mathbf{\Phi}})$, we have the following derivations

[TABLE]

According to Lagrange’s mean value theorem, since $h(\alpha,\mathbf{\Phi})$ is continuous and differentiable, there exists some point $\bar{\phi}_{i,m}$ between $\phi_{i,m}$ and $\hat{\phi}_{i,m}$ which satisfies

[TABLE]

Combining (41) and (42), we obtain

[TABLE]

From (29b), we can obtain

[TABLE]

which can be further derived as

[TABLE]

Since $\frac{\partial\mathbf{v}(\alpha,\bar{\mathbf{\Phi}})}{\partial{\bar{\phi}_{i,m}}}=\left[0;\frac{\partial\mathbf{z}(\bar{\mathbf{\Phi}})}{\partial\bar{\phi}_{i,m}}\right]$, the first term in (44)’s right side can be derived as (see (37))

[TABLE]

Since $\mathbf{z}(\bar{\mathbf{\Phi}})={\rm vec}(\mathbf{X(\bar{\mathbf{\Phi}})}^{H}\mathbf{X}(\bar{\mathbf{\Phi}}))$ and either $\frac{\partial\mathbf{X}^{\!H}(\bar{\mathbf{\Phi}})\mathbf{X}({\bar{\mathbf{\Phi}}})}{\partial{\bar{\phi}_{i,m}}}$ or $\frac{\partial^{2}\mathbf{X}^{\!H}(\bar{\mathbf{\Phi}})\mathbf{X}({\bar{\mathbf{\Phi}}})}{\partial{\bar{\phi}_{i,m}}^{2}}$ is an $M$-by-$M$ matrix (see them in (32) and (46) respectively), we can see that $\frac{\partial\mathbf{z}(\mathbf{\Phi})}{\partial\phi_{i,m}}$ and $\frac{\partial^{2}\mathbf{z}(\mathbf{\Phi})}{\partial\phi_{i,m}^{2}}$ has $2(M-1)$ nonzero constant modulus elements respectively. Since the maximum modulus of elements (MME) in $\bf{A}$ is $|\Theta|$, (45) can be further derived as

[TABLE]

Similarly, since $\frac{\partial^{2}\mathbf{v}{(\alpha,\bar{\mathbf{\Phi}})}}{\partial\phi_{i,m}^{2}}=\left[0;\frac{\partial^{2}\mathbf{z}{(\bar{\mathbf{\Phi}})}}{\partial\bar{\phi}_{i,m}^{2}}\right]$ and (39) holds, the second term in (44)’s right side can be denoted as

[TABLE]

Since MME in $\mathbf{a}_{\theta,\theta}$ is 1, we can see that MMEs in $\mathbf{q}$ and $\mathbf{A}$ are $\bar{P}_{\rm max}|\Theta|$ and $|\Theta|$ respectively, where $\bar{P}_{\rm max}=\max\limits_{\theta\in\Theta}\{\bar{P}_{\theta}\}$, i.e., $\bar{P}_{\rm max}$ is the maximum element of the desired beampattern. Since $\mathbf{X}(\mathbf{\Phi})\in\mathds{C}^{N\times M}$ and $x_{i,m}=e^{j\phi_{i,m}}$, MME in $\mathbf{z}$ is $N$. Then, (48) can be derived as

[TABLE]

Plugging (47) and (49) into (44), we have

[TABLE]

Moreover, plugging (50) into (43), we can obtain

[TABLE]

Therefore,$\nabla_{\mathbf{\Phi}}h(\alpha,\mathbf{\Phi})$ is Lipschitz continuous with the constant $L\geq 4(M\!-\!1)(\alpha_{\rm max}\bar{P}_{\rm max}\!+\!M^{2}N\!+\!2M\!-\!2)|\Theta|$.

Third, for $\nabla f_{n}(\mathbf{\Phi})$, there exists

[TABLE]

For $\frac{\partial^{2}f_{n}({\bar{\mathbf{\Phi}}})}{\partial\bar{\phi}_{i,m}^{2}}$, we have

[TABLE]

The elements in $\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})$ are $P_{\theta_{i},\theta_{j},n}$ or [math], where $P_{\theta_{i},\theta_{j},n}=\mathbf{a}_{\theta_{i}}^{H}\mathbf{X}^{H}\mathbf{S}_{n}\mathbf{X}\mathbf{a}_{\theta_{j}}$. Since $\mathbf{a}_{\theta}=\left[1,e^{j\pi\sin\theta},\cdots,e^{j\pi(M-1)\sin\theta}\right]^{T}$ and (53) holds, we can get

[TABLE]

From (54), it is easy to see that MME in all the $P_{\theta_{i},\theta_{j},n}$ is no larger than $M^{2}N$. Therefore, MME in $\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})$ is no larger than $w_{\rm c}M^{2}N$, where $w_{\rm c}=\max\{w_{\rm ac},w_{\rm cc}\}$.

Next, we consider MMEs in $\frac{\partial\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})}{\partial\bar{\phi}_{i,m}}$ and $\frac{\partial^{2}\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})}{\partial\bar{\phi}_{i,m}^{2}}$. From (53), we can get

[TABLE]

which indicate

[TABLE]

Therefore, we can conclude that MME in either $\frac{\partial\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})}{\partial\bar{\phi}_{i,m}}$ or $\frac{\partial^{2}\mathbf{B}_{n}({\bar{\mathbf{\Phi}}})}{\partial\bar{\phi}_{i,m}^{2}}$ is no larger than $w_{\rm c}(2M-1)$. Plugging these results into (52), we can obtain

[TABLE]

where $K$ is the number of considered spacial directions. Plugging (55) into (51), we can obtain

[TABLE]

Therefore, we can see that functions $\nabla f_{n}(\mathbf{\Phi}),n\in\mathcal{T}$ are Lipschitz continuous with the constant $L_{n}\!\!\geq\!\!2w_{\rm c}^{2}(\!2M\!-\!1\!)(\!M^{2}N\!\!+\!\!2M\!\!-\!\!1)K^{2}$. This concludes the proof. $\hfill\blacksquare$

Appendix B Proof of Lemmas 2–5

Lemma 2

The following two inequalities exist

[TABLE]

[TABLE]

Proof:

For (57), its left side can be derived as

[TABLE]

In Lemma 1, we show that gradient $\nabla h(\alpha,\mathbf{\Phi})$ is Lipschitz continuous. According to the Decent Lemma [36], we have

[TABLE]

which can be further derived to the following inequality

[TABLE]

Plugging (61) into (59), we have the following derivations

[TABLE]

Furthermore, since

[TABLE]

we can get (57) from (62).

Exploiting the Lipschitz continuous property of $\nabla f_{n}(\mathbf{\Phi})$ and similar derivations, (58) can also be obtained. $\hfill\blacksquare$ ∎

Lemma 3

in each iteration, $\forall~{}n\in\mathcal{T}$, $\|{\bf\Lambda}_{n}^{k+1}-{\bf\Lambda}_{n}^{k}\|_{F}^{2}$ can be bounded, i.e.,

[TABLE]

Proof:

The optimal solutions of the problems (23b) can be obtained by solving $\nabla_{\mathbf{\Phi}_{n}}\mathcal{U}_{n}\left(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k}\right)=0$, $\forall~{}n\in\mathcal{T}$, i.e.,

[TABLE]

Combining (64) and (23c), we can obtain

[TABLE]

Plugging (65) into $\|\mathbf{\Lambda}_{n}^{k+1}-\mathbf{\Lambda}_{n}^{k}\|_{F}^{2}$, we have the following derivations

[TABLE]

where the second inequality comes from Lemma 1. This completes the proof. $\hfill\blacksquare$ ∎

Lemma 4

let $\bar{c}_{n}\!=\!\rho_{n}^{3}\!-\!7\rho_{n}^{2}L_{n}\!-\!{8\rho_{n}L_{n}^{2}\!-\!32L_{n}^{3}}$ and $\tilde{c}_{n}\!=\!\rho_{n}^{3}\!-\!12\rho_{n}L_{n}^{2}\!-\!48L_{n}^{3}$. If $\bar{c}_{n}>0$ and $\tilde{c}_{n}>0$, in each consensus-ADMM iteration, the augmented Lagrangian function $\mathcal{L}(\cdot)$ decreases sufficiently, i.e.,

[TABLE]

Proof:

Define the following quantities

[TABLE]

Then, we can get

[TABLE]

For $\Delta_{\alpha,\mathbf{\Phi}}^{k}$, we have the following inequality

[TABLE]

Letting $\alpha=\alpha^{k}$ and $\mathbf{\Phi}=\mathbf{\Phi}^{k}$ and plugging them into (57), we can get

[TABLE]

Moreover, since $\mathcal{U}(\alpha^{k}\!,\!\mathbf{\Phi}^{k}\!,\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\})$ is up-bound of function $\mathcal{L}(\alpha^{k}\!,\!\mathbf{\Phi}^{k}\!,\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\})$, i.e., ,

[TABLE]

Combining the above two inequalities, we can conclude

[TABLE]

Then, inequality (68) can be further derived as follows

[TABLE]

Since $\mathcal{U}(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\!\in\!\mathcal{T}\})$ is strongly quadratic convex function with respect to $\alpha$ and $\mathbf{\Phi}$ (see (21)) [45], we have

[TABLE]

Moreover, since $\{\!\alpha^{k\!+\!1}\!,\!\mathbf{\Phi}^{k\!+\!1}\!\}\!\!\!=\!\!\!\!\underset{\alpha\!\in\!(0,\alpha_{\!\max}],\atop{{0\preceq\mathbf{\Phi}\prec 2\pi}}}{\arg\min}\!\mathcal{U}\!\!\left(\!\alpha,\!\mathbf{\Phi}\!,\!\{\!\mathbf{\Phi}_{n}^{k}\!,\!\mathbf{\Lambda}_{n}^{k}\!,\!n\!\in\!\mathcal{T}\}\!\right)$, we can get

[TABLE]

Plugging (70) and (71) into (69), we can obtain

[TABLE]

For $\Delta_{\mathbf{\Phi}_{n}}^{k}$, it can be rewritten as

[TABLE]

Plugging $\mathbf{\Phi}_{n}=\mathbf{\Phi}_{n}^{k}$ into (58), we can obtain

[TABLE]

Because the functions $\mathcal{U}_{n}\left(\mathbf{\Phi}^{k+1},\mathbf{\Phi}_{n},\mathbf{\Lambda}_{n}^{k}\right)$, $n\in\mathcal{T}$, are strongly convex with respect to $\mathbf{\Phi}_{n}$. We have

[TABLE]

Plugging the above two inequations into (73), we have

[TABLE]

Plugging (23c) into the above inequality, we can get

[TABLE]

Furthermore, plugging (63) into (75), we can see that $\Delta_{\mathbf{\Phi}_{n}}^{k}$ is lower bounded by

[TABLE]

For $\Delta_{\mathbf{\Lambda}_{n}}^{k}$, through similar derivations, we have

[TABLE]

Plugging (69), (76), and (77) into (67), we can obtain

[TABLE]

where $\bar{c}_{n}\!=\!\rho_{n}^{3}\!-\!7\rho_{n}^{2}L_{n}\!-\!{8\rho_{n}L_{n}^{2}\!-\!32L_{n}^{3}}$ and $\tilde{c}_{n}\!=\!\rho_{n}^{3}\!-\!12\rho_{n}L_{n}^{2}\!-\!48L_{n}^{3}$. So, $\forall n\in\mathcal{T}$, if $\bar{c}_{n}>0$ and $\tilde{c}_{n}>0$, then the augmented Lagrangian function $\mathcal{L}(\cdot)$ decreases sufficiently. This completes the proof. $\hfill\blacksquare$ ∎

Lemma 5

if $\rho_{n}>5L_{n}$, augmented Lagrangian function $\mathcal{L}(\alpha^{k+1},\mathbf{\Phi}^{k+1},\{\mathbf{\Phi}_{n}^{k+1},\mathbf{\Lambda}_{n}^{k+1},n\in\mathcal{T}\})\geq 0$, $\forall k$.

Proof:

First, we consider $\mathcal{L}_{n}({\mathbf{\Phi}}^{k+1},\mathbf{\Phi}_{n}^{k+1},\mathbf{\Lambda}_{n}^{k+1})$. Plugging (65) into $\mathcal{L}_{n}({\mathbf{\Phi}}^{k+1},\mathbf{\Phi}_{n}^{k+1},\mathbf{\Lambda}_{n}^{k+1})$, it can be equivalent to

[TABLE]

Since $\nabla f_{n}(\mathbf{\Phi})$ is Lipschitz continuous, the last term in (78) satisfies the following inequality

[TABLE]

Plugging it into (78), we can get

[TABLE]

Moreover, we can also exploit $\nabla f_{n}(\mathbf{\Phi}_{n})$’s Lipschitz continuous property and get the following inequality

[TABLE]

Plugging it into (79), we can obtain

[TABLE]

Second, plugging (80) into $\mathcal{L}(\alpha^{k+1},\mathbf{\Phi}^{k+1},\!\{\mathbf{\Phi}_{n}^{k+1},\!\mathbf{\Lambda}_{n}^{k+1}\!,\!n\!\in\!\mathcal{T}\}\!)$, we can get

[TABLE]

Since $h(\alpha^{k+1},\mathbf{\Phi}^{k+1})\geq 0$ and $f_{n}(\mathbf{\Phi}^{k+1})\geq 0$, we can conclude that, if $\rho_{n}>5L_{n}$, $\mathcal{L}(\alpha^{k+1},\mathbf{\Phi}^{k+1},\{\mathbf{\Phi}_{n}^{k+1},\mathbf{\Lambda}_{n}^{k+1},n\in\mathcal{T}\})\geq 0$. This completes the proof. $\hfill\blacksquare$ ∎

Appendix C Proof of theorem 1

First, we prove (27) in Theorem 1.

The pre-conditions that satisfying Lemma 4 and Lemma 5 hold are

[TABLE]

Letting $\beta=\frac{\rho_{n}}{L_{n}}$ and plugging $\beta$ into the above inequalities, we can obtain

[TABLE]

The roots of the cubic function in the left side of the first inequality can be determined through the famous Cardano formula [46], which are $4.72,-2.36+2.15i,-2.36-2.15i$. Similarly, we can obtain the roots of the cubic function in the second inequality, which are $8.41,-0.70+1.82i,-0.70-1.82i$. Combining the above results with $\beta>5$, we can find that when $\beta\geq 8.41$, all the inequalities in (82) hold simultaneously. To simplify the description, we choose $\rho_{n}\geq 9L_{n}$. Therefore, Lemma 4 and Lemma 5 are tenable when $\rho_{n}\geq 9L_{n},\forall n\in\mathcal{T}$.

Since Lemma 4 holds, we sum both sides of the inequality (66) when $k=1,2,\dotsb,+\infty$ and obtain

[TABLE]

Since $\displaystyle\lim_{k\rightarrow+\infty}\!\!\mathcal{L}(\alpha^{k\!+\!1},\!\mathbf{\Phi}^{k\!+\!1},\{\mathbf{\Phi}_{n}^{k\!+\!1},\mathbf{\Lambda}_{n}^{k\!+\!1},\!n\!\in\!\mathcal{T}\})\!\geq\!0$, we can get the following inequality

[TABLE]

Since $\bar{c}_{n}\!>\!0$, $\tilde{c}_{n}\!>\!0$, $L>0$, and $L_{\alpha}>0$, the above inequality indicates that summation of infinite positive terms is less than some constant. Therefore, we can obtain (83), (84), and (85).

[TABLE]

[TABLE]

[TABLE]

Plugging (83), (84) into (63)’s right side, we can get

[TABLE]

Combining (86) and (23c), we further have

[TABLE]

Since $\alpha\in[0,\alpha_{\max})$ and $0\preceq\mathbf{\Phi}\prec 2\pi$, we can obtain the following convergence results from (84) and (85).

[TABLE]

Plugging (88b) into (87), we can conclude

[TABLE]

Plugging (87) into (65), we can derive

[TABLE]

Since $\|\nabla f_{n}({\mathbf{\Phi}})\!-\!\nabla f_{n}(\hat{\mathbf{\Phi}})\|_{F}\!\leq\!L_{n}\|{\mathbf{\Phi}}\!-\!\hat{\mathbf{\Phi}}\|_{F},\ n\in\mathcal{T}$ and $0\preceq\mathbf{\Phi},\hat{\mathbf{\Phi}}\prec 2\pi$, we can conclude that all the elements in $\nabla f_{n}({\mathbf{\Phi}})$ are bounded. From (90), it indicates that $\mathbf{\Lambda}_{n}^{k}$ is also bounded. Combining this result with (86), we can get

[TABLE]

Second, we prove $(\alpha^{*},\mathbf{\Phi}^{*})$ is a stationary point of problem (8).

Since $\{\!\alpha^{k+1}\!,\mathbf{\Phi}^{k+1}\!\}\!\!=\!\!\!\!\underset{\alpha\in(0,\alpha_{\max}],\atop{{0\preceq\mathbf{\Phi}\prec 2\pi}}}{\arg\min}\!\mathcal{U}\!\left(\alpha,\!\mathbf{\Phi},\!\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\}\right)$ and $\mathcal{U}\left(\alpha,\mathbf{\Phi},\{\mathbf{\Phi}_{n}^{k},\mathbf{\Lambda}_{n}^{k},n\in\mathcal{T}\}\right)$ is convex quadratic function with respect to $\alpha$ and $\mathbf{\Phi}$, we have the following stationary conditions.

[TABLE]

Plugging

[TABLE]

into the above stationary conditions, we can get

[TABLE]

When $k\rightarrow+\infty$, plugging convergence results (88), (89), and (91) into (92), we can obtain

[TABLE]

Since $\nabla f_{n}({\mathbf{\Phi}}^{*})=-\mathbf{\Lambda}_{n}^{*},\forall n\in\mathcal{T}$, then (93b) can be further derived as

[TABLE]

Moreover, since $e\left(\alpha,\mathbf{X}(\mathbf{\Phi})\right)=h(\alpha,\mathbf{\Phi})$ and $P_{c}\left(\mathbf{X}(\mathbf{\Phi})\right)=\sum\limits_{n\in\mathcal{T}}f_{n}({\mathbf{\Phi}})$, (93a) and (94) can be rewritten as

[TABLE]

which completes the proof. $\hfill\blacksquare$

Appendix D Proof of the convergence of the proposed algorithm with SBCD method

Performing expectation on both sides of (66), we can get the following inequality.

[TABLE]

where $p_{\min}\geq 0$ is the probability and $\bar{c}_{n}\!=\!\rho_{n}^{3}\!-\!7\rho_{n}^{2}L_{n}\!-\!8\rho_{n}L_{n}^{2}\!-\!32L_{n}^{3}$ and $\tilde{c}_{n}\!=\!\rho_{n}^{3}\!-\!12\rho_{n}L_{n}^{2}\!-\!48L_{n}^{3}$. We set $\rho_{n}\geq 9L_{n},\forall n\in\mathcal{T}$ to guarantee $\bar{c}_{n}>0$ and $\tilde{c}_{n}>0$. Then, the augmented Lagrangian function decreases sufficiently in each consensus-ADMM-SBCD iteration.

Performing expectation on both sides of (81), we can get the following inequality.

[TABLE]

where $p_{n}>0$ is the probability. Since $h(\alpha^{k+1},\mathbf{\Phi}^{k+1})\geq 0$ and $\rho_{n}\geq 5L_{n},\forall n\in\mathcal{T}$, from (96) we can conclude

[TABLE]

Summing both sides of the inequality (95) for $k=1,2,\dotsb,+\infty$, we can obtain

[TABLE]

Plugging (97) into above inequality, it can be simplified as

[TABLE]

Since $\bar{c}_{n}\!>\!0$, $\tilde{c}_{n}\!>\!0$, $L>0$, and $L_{\alpha}>0$, the above inequality indicates that summation of infinite positive terms is less than some constant. Then, we can conclude the following results

[TABLE]

Since (23c) and (63) hold, we can further get

[TABLE]

Through similar discussions (88) to (91), we can obtain the following convergence results

[TABLE]

which concludes the proof. $\hfill\blacksquare$