1 Introduction↩︎

Integrated sensing and communication (ISAC) unify sensing and communication (S\(\&\)C) tasks into a single system, improving efficiency and performance by sharing various resources like spectrum and hardware [1], [2]. With its potential to support emerging applications requiring high-quality wireless connections and accurate sensing, such as autonomous driving and smart homes, ISAC is widely regarded as a key enabler for next-generation wireless networks [3], [4]. However, S\(\&\)C subsystems have distinct performance requirements [5]: sensing favors unimodular and deterministic waveforms, while communication relies on waveforms with high degrees of freedom (DoFs) and randomness for efficient information transmission. Unimodular waveform here means the modulus of the complex baseband signal is unity at all times. As a result, the interference between them is inevitable.

Since wireless systems often need to simultaneously serve multiple users and meet their various sensing and communication needs, multiple-input multiple-output (MIMO) technology [6], [7] becomes particularly important in ISAC systems. Early ISAC research focused on the MIMO radar-communication coexistence by mitigating interference through spectrum-sharing techniques like dynamic spectrum access [8] and null-space projection [9]. Although these approaches enabled S\(\&\)C coexistence, the systems were designed separately and required side-information exchange, leading to additional cooperation costs [10]. As an advancement from the spectrum-sharing scheme, some works [11], [12] focused on designing or implementing new antenna structures for ISAC purpose, such as the fluid antenna-aided ISAC systems and its combination with reconfigurable intelligent surface. Besides the newly appeared antenna structures, dual-functional waveform design has recently drawn a lot of interest [13] for its ability to sense targets and transmit information using a single device simultaneously, eliminating the need for side-information exchange in S\(\&\)C cooperation [10] and the change of antenna structures.

Many studies have focused on the design of dual-functional waveforms for ISAC; see [1] and the references therein. However, most of these works overlook the crucial aspect of controlling the sidelobes of the correlation function. In modern communication systems such as 5G NR, the transmitted signals typically occupy a certain bandwidth [14]. When these wideband signals are reused for sensing in ISAC systems, the resulting range sidelobes can significantly degrade sensing performance. High sidelobe levels are undesirable because they may interfere with or even obscure the reflections from distant targets or those with small radar cross-sections (RCS), potentially leading to missed detections [15]. Moreover, even when the radar receiver beampattern is well-optimized, excessive range sidelobes can still cause false alarms. Therefore, effective sidelobe control is essential in ISAC waveform design to ensure reliable sensing performance.

Additionally, it is well known that power amplifiers operate most efficiently when the input signals are unimodular [12], which further motivates the need for unimodular waveform design in ISAC systems. In this paper, we address the problem of dual-functional waveform design for monostatic downlink transmission in a MIMO-ISAC system. Our objective is to jointly maximize the achievable communication rate and minimize the mean square error (MSE) of the radar beampattern while simultaneously controlling the peak sidelobe level.

1.1 Literature Review↩︎

Existing dual-functional waveform design research can be roughly divided into three categories: sensing-centric, communication-centric, and joint waveform design, which will be detailed below.

Sensing-centric waveform design focuses on embedding the communication information into sensing waveforms [1]. For example, the communication information can be embedded into chirp waveforms [16], spatial beampattern [17], and ambiguity function [18]. While these approaches generally exhibit strong radar sensing capabilities, they often face low communication rates due to the limited number of embedded information bits. In contrast, communication-centric waveform design implements radar sensing using existing communication waveforms, such as orthogonal frequency division multiplexing (OFDM) [19] and orthogonal time frequency space (OTFS) [20]. However, the sensing performance of communication-centric designs is unpredictable due to the inherent randomness of communication signals and potential distortion from a high peak-to-average power ratio.

To address the limitations of separate designs and achieve trade-offs between S\(\&\)C, many works have focused on joint waveform design [1]. This approach constructs waveforms by solving optimization problems under various S\(\&\)C constraints. More specifically, in [21], the authors jointly maximized a weighted sum rate while minimizing the radar beampattern approximation MSE, constrained by per-antenna power limits, enabling rate-splitting multiple access and interference management in MIMO-ISAC systems. In [22], the authors minimized the Cramér-Rao bound (CRB) for direction-of-arrival (DoA) estimation by designing the beamforming matrix under individual signal-to-interference-plus-noise ratio (SINR) constraints at each communication receiver and the transmit power budget. Additionally, the work [23] proposed to jointly precode communication and radar waveforms to achieve maximum DoFs in waveform design. Recent progress on joint waveform design in MIMO-ISAC systems has been made in [24] and [25], in which [24] optimized a weighted combination of the sum rate and the CRB for target estimation and [25] maximized the system energy efficiency by constraining the transmit power budget, communication SINR, and target estimation CRB.

The works mentioned previously accomplished a balance between S\(\&\)C to some extent, but they did not address the issue of range sidelobe control. As for now, there have been few studies on sidelobe control for ISAC systems. The work [15] proposed a MIMO-ISAC waveform design framework that realized an integrated sidelobe level (ISL) reduction by minimizing a weighted sum of beampattern MSE, ISL, and MUI. Alternatively, the work [26] focused on maximizing the SINR at the radar output while ensuring communication performance, which also achieves the ISL reduction. However, the communication channels in the two works are assumed to be frequency-flat fading, whereas the ISAC signals capable of distinguishing symbol-level delays correspond to frequency-selective fading communication channels. In addition to the ISL metric, the peak sidelobe level (PSL) should be more important for sidelobe control. That’s because the PSL at the radar receiver dictates the false alarm probability [27]. High PSL can lead to false alarms [28], making PSL control significant for achieving a low false alarm rate (FAR). Despite its importance, the PSL minimization problem has received relatively little attention in both the radar signal [27] and ISAC waveform design. One of the main difficulties in directly minimizing the PSL is that the design metric is not differentiable, and the corresponding optimization problem is a minimax problem. To the best of our knowledge, the existing works on PSL control focus on pure radar sensing, and PSL control in ISAC scenarios remains unexplored. For instance, the work [29] proposed a two-step scheme to control the PSL in MIMO radar. The work [30] approximated PSL suppression by minimizing an \(\ell_{p}\) norm metric with large \(p\) for single-antenna radar systems.

Taking the above factors into consideration, the goal of this paper is to design a downlink MIMO-ISAC waveform that maximizes the achievable rate at the communication receivers, makes the radar receiver’s beampattern as close as possible to the desired one, and controls the PSL under the constraints of unimodular transmit sequences. However, as discussed in [29], designing waveforms with good correlation properties under the constraint of unimodular sequences is already a complicated task, and it will be more challenging further to require good communication performance and the desired beampattern. To simplify the optimization problem and obtain more DoFs on waveform design, we consider the joint design of the receive filter and transmit sequence in this paper.

1.2 Our Contributions↩︎

The main contributions of this paper are as follows.

• Unimodular ISAC waveform design with range sidelobe level control: Many factors are taken into account in the problem formulation, including beampattern MSE, communication MUI, PSL, and unimodular transmit sequence. An optimization problem is formulated to minimize a weighted sum of the radar receive beampattern MSE and the MUI at the communication receiver. To avoid solving a minimax problem by minimizing the PSL directly, we control the PSL by constraining the level of sidelobe at each range bin. The formulated problem is a large-scale optimization problem with convex and nonconvex constraints.

• Efficient solution: An inexact augmented Lagrangian method (ALM) algorithm is proposed to solve the formulated problem, where each ALM nonconvex subproblem is solved approximately to avoid the huge number of iterations for obtaining an exact stationary point. Specifically, a block successive upper-bound minimization (BSUM) scheme is custom-designed to solve the subproblems in the ALM algorithm, and the updates for all blocks of variable in the BSUM scheme admit closed-form solutions, which makes the proposed algorithm efficient.

• Convergence guarantee: We analyze the convergence of the proposed inexact ALM algorithm with an adaptive penalty parameter. We show that the proposed algorithm is guaranteed to find a feasible stationary point of the formulated problem. This is the best that one can expect for this nonconvex optimization problem (with many nonconvex constraints).

Extensive simulation results are provided to demonstrate the effectiveness of the proposed algorithm. Specifically, Monte Carlo simulations are performed to evaluate the convergence performance of the proposed algorithm. The impacts of different system parameters on the system performance are examined. Finally, the proposed algorithm is compared with the ALM algorithm with a fixed penalty parameter and the modified work of [29] to show the superiority of the proposed algorithm.

It is worth mentioning that a similar inexact-ALM framework has been proposed in [31]. In [31], radar sensing beamforming is achieved through quantized constant-envelope waveform design while satisfying the communication performance constraints. However, it does not address how to suppress the sidelobes of the correlation function. Low-range bin sidelobes and beampattern are equally important for good radar imaging. Furthermore, after considering sidelobe suppression, the convergence analysis of the optimization problem formulated in this paper becomes more complex compared to [31] since the constraint terms introduced by the auxiliary variables in [31] are linear, whereas those in this paper are nonlinear.

The rest of the paper is organized as follows. In Section 2, we introduce the system model and formulate the optimization problem. In Section 3, we propose an efficient algorithm for solving the formulated problem and analyze the convergence of the proposed algorithm. We present extensive simulation results in Section 4. Finally, we conclude the paper in Section 5.

Notations: We use \(x\), \(\mathbf{x}\), \(\mathbf{X}\), and \(\mathcal{X}\) to represent scalar, column vector, matrix, and set, respectively. The notation vec(\(\cdot\)) represents the vectorization of a matrix by stacking its columns. The notations \(\|\cdot\|_{1}\), \(\|\cdot\|_{2}\), and \(\|\cdot\|_{\text{F}}\) denote correspondingly the \(\ell_1\), \(\ell_2\), and Frobenius norms of a matrix, respectively. \(\mathbb{C}\) and \(\mathbb{R}\) denote the sets of complex and real numbers, respectively; \(\mathcal{R}\{\cdot\}\) and \(\mathcal{I}\{\cdot\}\) are the real and imaginary parts of a complex number, respectively. The superscript \((\cdot)^{*}\), \((\cdot)^{\text{T}}\) and \((\cdot)^{\text{H}}\) represent the conjugate, the transpose and the conjugate transpose operations, respectively. \(\mathbb{I}_{\mathcal{A}}(A)\) denotes an indicator function of \(\mathcal{A}\), and it takes \(0\) if \(A \in \mathcal{A}\) and \(+\infty\), otherwise. \(\otimes\) represents the Kronecker product. \(\mathbf{A}_{[i:j, :]}\) represents the sub-matrix from the \(i\)-th row to the \(j\)-th row of a matrix \(\mathbf{A}\), and \(A[i,j]\) represents the element in the \(i\)-th row and \(j\)-th column of matrix \(\mathbf{A}\). \(\mathrm{tr}(\mathbf{A})\) means the trace of a matrix \(\mathbf{A}\). \(\mathbf{A} \succeq \mathbf{B}\) means \(\mathbf{A} - \mathbf{B}\) is positive semidefinite. \([n]\) denotes \(\{1, 2, \dots, n\}\) for a positive integer \(n\).

2 System Model and Problem Formulation↩︎

As shown in Fig. 1, we consider an ISAC base station (BS) equipped with a uniform linear array (ULA) of \(N_{\text{T}}\) antennas. \(N_{\text{C}}\) single-antenna communication receivers are in the downlink transmission. We consider a colocated radar receive station with a \(N_{\text{R}}\)-antenna ULA. Below, we introduce the mathematical frameworks used in both S\(\&\)C functions and formulate the optimization problem.

2.1 Communication Signal Model↩︎

Considering the frequency-selective feature of the channel in the current communication system, [14], we implement the channel model proposed in [32], in which the channel between \(m\)-th receiver and \(n\)-th BS transmitting antenna, denoted as \(\mathbf{H}_{m,n}\), is frequency-selective and assumed to be quasi-stationary during one block transmission. A \(q\)-ray model defines the discrete-time channel gains. The channel \(\mathbf{H}_{m,n}\) can be characterized by a lower triangular Toeplitz matrix with the first column being \([h_0^{(m, n)}, h_{1}^{(m, n)}, \dots, h_{L}^{(m, n)}, \mathbf{0}_{T}]^{\text{T}}\), where \(L = \lceil \frac{\tau_{\max}}{T_s} \rceil\), \(T_s\) is the duration of one time-slot, \(\tau_{\max}\) is the maximum delay spread, and \(T\) is the length of the block.

Due to the multi-taps feature of the channel, inter-block interference (IBI) occurs between every two consecutive transmissions (blocks). To eliminate IBI, a cyclic prefix (CP) of length \(L\) is added at each transmission and will be removed at the receiver. The noiseless received signal at the \(m\)-th receiver is given by \(\mathbf{y}_{m} = \sum_{n=1}^{N_{\text{T}}} \tilde{\mathbf{H}}_{m,n} \mathbf{x}_{n}\), where \(\tilde{\mathbf{H}}_{m,n} = \boldsymbol{\Pi}_{\text{cp}} \mathbf{H}_{m,n} \boldsymbol{\Gamma}_{\text{cp}}\) is a \(T \times T\) matrix, \(\boldsymbol{\Gamma}_{\text{cp}} = [\mathbf{I}_{\text{cp}}, \mathbf{I}_{T}]^{\text{T}}\) and \(\boldsymbol{\Pi}_{\text{cp}} = [\mathbf{0}_{T \times L}, \mathbf{I}_{T}]\) denote the CP-inducing and CP-removing matrices, respectively, and \(\mathbf{I}_{\text{cp}}\) contains the last \(L\) columns of a \(T\)-dimensional identity matrix \(\mathbf{I}_{T}\). Then, the received signals at all \(N_{\text{C}}\) users can be represented as \[\label{eq:Comm95model} \mathbf{y} = \mathbf{H} \mathbf{x} + \mathbf{n},\tag{1}\] where \(\mathbf{y} = [\mathbf{y}_{1}^{\text{T}}, \mathbf{y}_{2}^{\text{T}}, \dots, \mathbf{y}_{N_{\text{C}}}^{\text{T}}]^{\text{T}}\) with \(\mathbf{y}_{i} \in \mathbb{C}^{T \times 1}\) being the received signal at the \(i\)-th user, \(\mathbf{x} = [\mathbf{x}_{1}^{\text{T}}, \mathbf{x}_{2}^{\text{T}}, \dots, \mathbf{x}_{N_{\text{T}}}^{\text{T}}]\) with \(\mathbf{x}_{i} \in \mathbb{C}^{T \times 1}\) being the \(i\)-th row vector of the transmission signal matrix \(\mathbf{X} \in \mathbb{C}^{N_T \times T}\), \(\mathbf{n} \in \mathbb{C}^{N_{\text{T}}(T+L) \times 1}\) is an additive white Gaussian noise (AWGN) vector with power \(\sigma_{n}^2\), and \(\mathbf{H} \in \mathbb{C}^{N_{\text{C}}T \times N_{\text{T}} T}\) is defined as \[\mathbf{H} = \left[ \begin{smallmatrix} \tilde{\mathbf{H}}_{1,1} & \tilde{\mathbf{H}}_{1,2} & \dots & \tilde{\mathbf{H}}_{1,N_{\text{T}}}\\ \tilde{\mathbf{H}}_{2,1} & \tilde{\mathbf{H}}_{2,2} & \dots & \tilde{\mathbf{H}}_{2,N_{\text{T}}}\\ \vdots & & \ddots & \vdots \\ \tilde{\mathbf{H}}_{N_{\text{C}},1} & \tilde{\mathbf{H}}_{N_{\text{C}},2} & \dots & \tilde{\mathbf{H}}_{N_{\text{C}},N_{\text{T}}} \end{smallmatrix} \right].\] We assume that the BS has perfect channel state information (CSI), which allows us to clearly highlight the performance gain and design principles of the proposed system under ideal conditions. The waveform design under partial or statistical CSI, as in [33], will be our future work.

Let \(\mathbf{S} \in \mathbb{C}^{N_{\text{C}} \times T}\) denote the transmitted information symbol matrix, where each entry of \(\mathbf{S}\) is randomly drawn from a given constellation. The received signal at the communication users can then be represented in vector form as \[\mathbf{y} = \mathbf{s} + (\mathbf{H} \mathbf{x} - \mathbf{s}) + \mathbf{n},\] where \(\mathbf{s} = [\mathbf{s}_{1}^{\text{T}}, \mathbf{s}_{2}^{\text{T}}, \dots, \mathbf{s}_{N_{\text{C}}}^{\text{T}}]^{\text{T}}\) with \(\mathbf{s}_{i} \in \mathbb{C}^{T \times 1}\) being the \(i\)-th row vector of \(\mathbf{S}\), the term \(\mathbf{H} \mathbf{x} - \mathbf{s}\) can be viewed as the MUI[1] that interferes with the symbol detection at the communication receiver side. The received SINR at the \(i\)-th user is defined as \[\gamma_i = \frac{\mathbb{E}_{\mathbf{s}_{i}}\{\|\mathbf{s}_{i}\|_{2}^{2}\}}{\mathbb{E}_{\mathbf{s}}}\{\|\mathbf{H}_{i} \mathbf{x} - \mathbf{s}\|_{2}^{2}\} + \sigma_n^2,\] where \(\mathbf{H}_{i} = [\tilde{\mathbf{H}}_{i,1}, \tilde{\mathbf{H}}_{i,2}, \dots, \tilde{\mathbf{H}}_{i,N_{\text{T}}}]\). The achievable downlink sum rate of the users can be given as \[R = \sum_{i=1}^{N_{\text{C}}} \log_{2}(1 + \gamma_i).\]

Given that \(\mathbb{E}_{\mathbf{s}_{i}}\{\|\mathbf{s}_{i}\|_{2}^{2}\}\) is a fixed value for a specific constellation strategy, minimizing \(P_{\text{MUI}} = \|\mathbf{H}\mathbf{x} - \mathbf{s}\|_{2}^{2}\) leads to an increase in SINR, thereby indicating a higher achievable sum rate⁶.

2.2 Signal Model for Radar Imaging↩︎

As shown in Fig. 2, after removing the CP, the received echos for radar imaging at the \(i\)-th range bin, \(\mathbf{D}_{(i)} \in \mathbb{C}^{N_{R} \times T}\), are obtained from the \(t_{i+L}\)-th time slot to the \(t_{i+L+T-1}\)-th time slot, where \(t_{i+L}\) is the starting time slot of the \(i\)-th range bin, and \(\mathbf{D}_{(i)} \in \mathbb{C}^{N_{R} \times T}\) can be expressed as [29] \[\begin{gather} \label{eq:radar95echo} \mathbf{D}_{(i)} = \sum_{p_{(i)} = 1}^{P_{(i)}} h_{p_{(i)}} \mathbf{a}(\theta_{p_{(i)}}) \mathbf{v}(\theta_{p_{(i)}})^{\text{T}} \mathbf{X} \\ + \sum_{k_{(i)} \in \mathbf{\Omega}_{K}} \sum_{p_{k_{(i)}}^{\prime} = 1}^{P_{k_{(i)}}^{\prime}} h_{p_{k_{(i)}}^{\prime}} \mathbf{a}(\theta_{p_{k_{(i)}}^{\prime}}) \mathbf{v}(\theta_{p_{k_{(i)}}^{\prime}})^{\text{T}} \mathbf{X} \mathbf{J}_{k} + \mathbf{Z}_{\text{R}}, \end{gather}\tag{2}\] where \(\mathbf{X}\) is the transmitted ISAC waveform defined above, \(\mathbf{Z}_{\text{R}}\) is the additive Gaussian noise matrix with zero mean and covariance matrix \(\mathbf{R}_{\text{s}}\), \(\{\mathbf{a}(\theta_{p_{(i)}}) \in \mathbb{C}^{N_{R} \times 1}\}_{p_{(i)}=1}^{P_{(i)}}\) with its \(n\)-th element being \(e^{-j 2 \pi n \frac{d}{\lambda} \sin(\theta_{p, i})}\), \(n=0, 1, \dots, N_R-1\), and \(\{\mathbf{v}(\theta_{p_{(i)}}) = [1, e^{-j 2 \pi \frac{d}{\lambda} \sin(\theta_{p, i})}, \dots, e^{-j 2 \pi (N_{T}-1) \frac{d}{\lambda} \sin(\theta_{p, i})}] \in \mathbb{C}^{N_{T} \times 1}\}_{p_{(i)} = 1}^{P_{(i)}}\) are the radar receive and transmit steering vectors for the \(P_{(i)}\) targets in the \(i\)-th range bin, respectively, \(d = \lambda/2\), and \(\{h_{p_{(i)}} \in \mathbb{C}\}_{p_{(i)}=1}^{P_{(i)}}\) are complex amplitudes proportional to the RCS of these \(P_{(i)}\) targets. The parameters with subscript \(p_{k_{(i)}}^{\prime}\) are defined for the \(P_{k_{(i)}}^{\prime}\) scatterers at the \(k_{(i)}\)-th range bin with the same meanings, the temporal shifting matrix \(\mathbf{J}_{k}\) is defined as \[\mathbf{J}_{k} = \mathbf{J}_{-k}^{\text{T}} = \left[\begin{smallmatrix} \mathbf{0}_{(T - k) \times k} & \mathbf{I}_{T - k} \\ \mathbf{I}_{k} & \mathbf{0}_{k\times (T - k)} \\ \end{smallmatrix}\right],\] and \(\mathbf{\Omega}_{K} = \{-K, \dots, -1, 1, 2, \dots, K\}\), where \(K\) is the maximum difference of arrival times between backscattered signals from the range bin of interest and signals from neighboring range bins.

To obtain more DoFs on the waveform design, we jointly design the radar receive filter and transmit sequence. Denote the receive filter at the radar receiver as \(\mathbf{F} \in \mathbb{C}^{N_{\text{R}} \times T}\), the radar image at angle \(\theta\) and range bin \(i\) is given by \[\label{eq:DF95radar} r_{\theta, i} = |\mathbf{a}(\theta)^{\text{H}} \mathbf{D}_{(i)} \mathbf{F}^{\text{H}} \mathbf{a}(\theta)|.\tag{3}\] Under the signal models 2 and 3 , we consider optimizing two metrics to improve radar imaging performance. One is the radar receive beampattern MSE, defined as \(\| \mathbf{X} \mathbf{F}^{\text{H}} - \mathbf{R}_{\text{d}} \|_{\text{F}}\), measuring how close the receive beampattern approximates the desired one, where \(\mathbf{R}_{\text{d}}\) is the desired spatial response determined by the prior knowledge about the targets [31] so that we can have stronger angle responses of the targets.For example, if we roughly know there are three targets at \(-10^{\circ}\), \(0^{\circ}\), and \(10^{\circ}\), we can solve the correlation matrix \(\mathbf{R}_{\text{d}}\) through the algorithm in [34]. Otherwise, we set \(\mathbf{R}_{\text{d}}\) as an identity matrix, which corresponds to an omnidirectional beampattern, if we have no knowledge about the targets. Another metric is the maximum sidelobe level \(\max_{k \in \mathbf{\Omega}_{K}}\{\|\mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}}\|_{\text{F}}\}\), which is required to be less than a preset level \(\xi\) for all \(k \in \mathbf{\Omega}_{K}\) to reduce the clutters from neighboring range bins and to ensure a clear response on the range bin of interest.

2.3 Problem Formulation↩︎

Based on the above discussions, the joint design problem of the radar receive filter and the ISAC waveform is formulated to minimize a weighted sum of the receive radar beampattern MSE and MUI at communication receivers under the constraints of unimodular transmit sequence and the full power limit of the radar receive filter⁷. The optimization problem is presented below: \[\begin{align} & \min_{\mathbf{F}, \mathbf{X}} f(\mathbf{F}, \mathbf{X}) \tag{4} \\ & \;\text{s.t. } \| \mathbf{F} \|_{\text{F}}^2 = P_{\text{F}}, \tag{5} \\ & \;\;\;\;\;\;\, \| \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} \|_{\text{F}} \leq \xi,~\forall~k \in \mathbf{\Omega}_{K}, \tag{6} \\ & \;\;\;\;\;\;\; |x_{ij}| = P_{x}, \forall~x_{ij} \in \mathbf{X}\tag{7}, \end{align}\tag{8}\] where \(f(\mathbf{F}, \mathbf{X}) = \alpha \| \mathbf{X} \mathbf{F}^{\text{H}} - \mathbf{R}_{\text{d}} \|_{\text{F}}^2 + (1 - \alpha) \|\mathbf{H} \mathbf{x} - \mathbf{s}\|_{2}^2\), and 7 refers to the unimodular constraint on the transmit sequence.

Solving problem 8 is not straightforward due to the presence of both convex PSL constraints and nonconvex constraints on the receive filter and transmit sequence. In the next section, we propose an inexact ALM algorithm for solving problem 8 efficiently. “Inexact” here refers to solving the subproblems in the ALM algorithm inexactly. By doing so, we can significantly reduce the computational cost of solving the ALM subproblems without sacrificing the solution quality.

3 Proposed Approach↩︎

In this section, we propose an inexact ALM algorithm for solving problem 8 . More specifically, we first introduce the main idea of the proposed inexact ALM algorithm in Section 3.1. Then, we propose a scheme for finding a feasible point of problem 8 in Section 3.2, which plays an important role in guaranteeing the convergence of the proposed algorithm to a feasible stationary point. We custom-design a BSUM algorithm for efficiently solving the ALM subproblem in Section 3.3. Finally, we analyze the convergence of the proposed algorithm in Section 3.4.

3.1 Framework of Proposed ALM Algorithm↩︎

The inexact ALM algorithm is proposed to solve problem 8 . Its basic idea is to decompose the original optimization problem into smaller and more manageable subproblems. At each iteration, the algorithm updates primal variables associated with each ALM subproblem with a penalty term for penalizing the violation of the constraints in the original problem, followed by an update of the dual variables. Unlike the classic ALM algorithm, which solves each subproblem exactly, we solve each subproblem inexactly at each iteration for computational efficiency.

To present the inexact ALM algorithm, we first reformulate problem 8 by introducing auxiliary variables \(\boldsymbol{C}= (\mathbf{C}_{k})_{k \in \mathbf{\Omega}_{K}}\) with \(\mathbf{C}_{k} \in \mathbb{C}^{N_{T} \times N_{R}}\) for each \(k \in \mathbf{\Omega}_{K}\). The problem 8 becomes \[\begin{align}\label{prob:FX95with95C} & \min_{\mathbf{F}, \mathbf{X}, \boldsymbol{C}} f(\mathbf{F}, \mathbf{X}) \\ & \;\; \text{s.t. } \| \mathbf{F} \|_{\text{F}}^2 = P_{\text{F}}, \\ & \;\;\;\;\;\;\;\, \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} = \mathbf{C}_{k}, \\ & \;\;\;\;\;\;\;\, \|\mathbf{C}_{k}\|_{\text{F}} \leq \xi,~\forall~k \in \mathbf{\Omega}_{k}, \\ & \;\;\;\;\;\;\;\; |x_{ij}| = P_{x}, \forall~x_{ij} \in \mathbf{X}. \end{align}\tag{9}\] By denoting the Lagrange multipliers associated with the constraints \(\mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} = \mathbf{C}_{k}\) for all \(k \in \mathbf{\Omega}_K\) as \(\boldsymbol{U} = (\mathbf{U}_{k})_{k \in \mathbf{\Omega}_{K}}\), the corresponding augmented Lagrangian function (ALF) [35] of problem 9 is \[\begin{align} \label{eq:ALF95origin} & \mathcal{L}_{\rho}(\boldsymbol{C}, \mathbf{F}, \mathbf{X}; \boldsymbol{U}) = f(\mathbf{F}, \mathbf{X}) + \sum_{k \in \mathbf{\Omega}_{K}} \Big[\frac{\rho}{2} \|\mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} - \mathbf{C}_{k}\|_{\text{F}}^2 \nonumber \\ & \qquad \quad \qquad \qquad + \mathcal{R}\{tr[{\mathbf{U}_{k}}^{\text{H}} (\mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} - \mathbf{C}_{k})]\}\Big]. \end{align}\tag{10}\]

The proposed ALM algorithm is based on the above ALF and mainly consists of four main components, which are the initialization, the solution of the ALM subproblem, the update of the Lagrange variable (also called dual variable), and the update of the penalty parameter. In the next, we introduce them one by one and highlight their important roles in the whole ALM algorithm. To simplify the notations, the point \((\boldsymbol{C}, \mathbf{F}, \mathbf{X})\) is represented by a multiplet \(\boldsymbol{z} = (\boldsymbol{C}, \mathbf{F}, \mathbf{X})\).

3.1.1 Initialization↩︎

To initialize the proposed algorithm, we choose a feasible point \(\boldsymbol{z}^{(\text{feas})} \in \mathcal{S} \cap \mathcal{S}_{0}\), where \[\begin{align} \mathcal{S}_{0} & = \left\{(\boldsymbol{C}, \mathbf{F}, \mathbf{X}) \mid \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} = \mathbf{C}_{k},~\forall~k \in \mathbf{\Omega}_{k}\right\}, \\ \mathcal{S} & = \{(\boldsymbol{C}, \mathbf{F}, \mathbf{X}) \mid \mathbf{C}_{k} \in \mathcal{S}_{\text{C}},~\forall~k \in \mathbf{\Omega}_{k}; \mathbf{F} \in \mathcal{S}_{\text{F}}; \mathbf{X} \in \mathcal{S}_{\text{X}}\}, \end{align}\] and \[\begin{align} \mathcal{S}_{\text{C}} & = \left\{\mathbf{C} \in \mathbb{C}^{N_{T} \times N_{T}} \mid \|\mathbf{C}\|_{\text{F}} \leq \xi\right\}, \\ \mathcal{S}_{\text{F}} & = \{\mathbf{F} \in \mathbb{C}^{N_{T} \times T} \mid \|\mathbf{F}\|_{\text{F}}^2 = P_{\text{F}}\}, \\ \mathcal{S}_{\text{X}} & = \{\mathbf{X} \in \mathbb{C}^{N_{T} \times T} \mid |x_{ij}| = P_{x}, \forall~x_{ij} \in \mathbf{X}\}. \end{align}\] The details of finding a feasible point of problem 9 is relegated to Section 3.2 to maintain smoothness in presenting the proposed inexact ALM algorithm.

By choosing the penalty parameter \(\rho^{(0)} > 0\), choosing an arbitrary initial point \(\boldsymbol{z}^{(0)} \in \mathcal{S}\), and setting a finite Lagrange multiplier \(\boldsymbol{U}^{(0)}\), we can choose a finite constant number \(\zeta\) satisfying \[\label{eq:zeta95the95upper95bound} \zeta \geq \max\left\{f(\mathbf{F}^{(\text{feas})}, \mathbf{X}^{(\text{feas})}), \mathcal{L}_{\rho^{(0)}}(\boldsymbol{z}^{(0)}; \boldsymbol{U}^{(0)})\right\}.\tag{11}\] The upperboundness of the objective function and the ALF are important in proving the convergence of the proposed ALM algorithm to a feasible stationary point. In our proposed ALM algorithm, the penalty parameter is adaptively updated. The adaptive update of the penalty parameter enables a fast convergence but also makes the convergence analysis difficult. The upperboundness here can guarantee that any limit point of the sequence generated by the proposed algorithm is always a feasible stationary point, even when the penalty parameter becomes positive infinity. More details on the convergence analysis will be presented in Section 3.4.

3.1.2 Solving the ALM Subproblem Inexactly↩︎

At the \(\ell\)-th iteration, the classic ALM algorithm updates the primal variables by solving the following ALM subproblem: \[\label{prob:sub95problem95ALM} \boldsymbol{z}^{(\ell+1)} \in \mathop{\mathrm{arg\,min}}_{\substack{\boldsymbol{z} \in \mathcal{S}}} \mathcal{L}_{\rho^{(\ell)}}( \boldsymbol{z}; \boldsymbol{U}^{(\ell)}).\tag{12}\] However, problem 12 is nonconvex due to the nonconvex constraints of full power radar receiving filter and unimodular transmit sequence. To improve the computational efficiency, we propose to inexactly solve problem 12 to an \(\varepsilon^{(\ell)}\)-stationary point, where the \(\varepsilon\)-stationary point is defined below:

The notations \(\partial \mathbb{I}_{\mathcal{S}_{\text{F}}}(\mathbf{F})\), \(\partial \mathbb{I}_{\mathcal{S}_{\text{X}}}(\mathbf{X})\), and \(\partial \mathbb{I}_{\mathcal{S}_{\text{C}}}(\mathbf{C}_{k})\) are subdifferentials of indicator functions; see [36] and [37] for more details. The \(\varepsilon\)-stationary point defined in Definition 1 reduces to the standard stationary point when setting \(\varepsilon=0\).

In inexactly solving the ALM subproblem in 12 , we also require that the inexact solution \(z^{(\ell+1)}\) also satisfies that the value of the ALF is bounded by \[\label{ineq:upper95bound95requirement} \mathcal{L}_{\rho^{(\ell)}}( \boldsymbol{z}^{(\ell+1)}; \boldsymbol{U}^{(\ell)}) \leq \zeta,\tag{13}\] where \(\zeta\) is defined in 11 . To find an \(\varepsilon\)-stationary point that satisfies (9), we shall develop a BSUM algorithm, which will be elaborated in detail in Section 3.3. This algorithm offers closed-form updates for all blocks of variable and is guaranteed to find an \(\varepsilon\)-stationary point within a finite number of iterations.

3.1.3 Updating Lagrange Multipliers↩︎

Once the \(\varepsilon^{(\ell)}\)-stationary point is obtained, we update the Lagrange multipliers in \(\boldsymbol{U}\) for all \(k \in \mathbf{\Omega}_{K}\) according to the following updating rule: \[\begin{align} & \tilde{\mathbf{U}}_{k} = \mathbf{U}_{k}^{(\ell)} + \rho^{(\ell)} (\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}), \tag{14} \\ & U_{k}^{(\ell+1)}[i,j] = \begin{cases} \frac{u_{\max}}{|\tilde{U}[i, j]|} \tilde{U}[i, j], & \text{if } |\tilde{U}[i, j]| > u_{\max}; \\ \tilde{U}[i, j], & \text{otherwise}. \end{cases}\tag{15} \end{align}\] In the above, \(\mathbf{U}_{k}\) is first updated through the standard updating rule [35] in 14 ; then, each element in \(\mathbf{U}_{k}\) is projected onto an interval \([-u_{\max}, u_{\max}]\) with \(u_{\max}>0\) being a preset constant. The projection onto the bounded set in 14 guarantees that the Lagrange multipliers are uniformly bounded and play a central role in guaranteeing the convergence of the proposed algorithm. The boundness of the multiplier sequence in the nonconvex setting remains an open research question [36]. Readers can refer to [36] for detailed discussions.

3.1.4 Updating the Penalty Parameter↩︎

We propose an adaptive update rule for the penalty parameter to accelerate the convergence speed. Denote the violation of the constraints after the \(\ell\)-th iteration as \(v^{(\ell+1)}\), i.e., \[\label{eq:constraints95violation} v^{(\ell+1)} = \sqrt{\sum_{k \in \mathbf{\Omega}_{K}} \| \mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {{\mathbf{F}}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\|_{\text{F}}^2}.\tag{16}\] The penalty parameter \(\rho^{(\ell+1)}\) will be updated through the following rule: \[\label{eq:rho95update} \rho^{(\ell+1)} = \begin{cases} \gamma \rho^{(\ell)}, & \text{if } v^{(\ell+1)} > \delta v^{(\ell)}; \\ \rho^{(\ell)}, & \text{otherwise}. \end{cases}\tag{17}\] where \(\delta > 0\) and \(\gamma > 1\). The update rule in 17 will increase the penalty parameter if the violation at the current iteration is not reduced sufficiently compared with the previous one. By dynamically adjusting the penalty parameter based on the current state of constraint violations and dual variable updates, this approach helps maintain an effective balance between enforcing constraint feasibility and minimizing the objective function.

The proposed inexact ALM algorithm for solving problem 9 is summarized in Algorithm 3. By setting a positive sequence \(\{\varepsilon^{(\ell)}\}_{\ell \geq 0}\) with \(\varepsilon^{(\ell)} \rightarrow 0\) as \(\ell \rightarrow +\infty\), the limit point of the sequence generated by Algorithm 3, \((\boldsymbol{z}^{(\infty)}, \boldsymbol{U}^{(\infty)})\), is a stationary point of the ALF in 10 . Besides, since \(\mathbf{X}^{(\infty)} \mathbf{J}_{k} {{\mathbf{F}}^{(\infty)}}^{\text{H}} = \mathbf{C}_{k}^{(\infty)}\) holds for all \(k \in \mathbf{\Omega}_{K}\), which will be proved further ahead, then \(\boldsymbol{z}^{(\infty)}\) is a feasible stationary point of problem 9 . According to [38], any feasible stationary point of 9 is also a feasible stationary point of 8 , i.e., \((\mathbf{F}^{(\infty)}, \mathbf{X}^{(\infty)})\) is a feasible stationary point of 8 . The detailed convergence analysis will be given in Section 3.4.

3.2 A Feasible Point of Problem 9↩︎

In this subsection, we focus on finding a feasible point of problem 9 , i.e., a pair of matrices \((\mathbf{F}, \mathbf{X}) \in \mathcal{S}_{\text{F}} \times \mathcal{S}_{\text{X}}\) such that \[\label{ineq:feasible95condition} \| \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} \|_{\text{F}} \leq \xi,~\forall~k \in \mathbf{\Omega}_{K}.\tag{18}\] The auxiliary variable \(\boldsymbol{C}\) can be easily obtained by setting \(\mathbf{C}_{k} = \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}}\) for all \(k \in \mathbf{\Omega}_{K}\).

We formulate two subproblems for finding \((\mathbf{F}, \mathbf{X}) \in \mathcal{S}_{\text{F}} \times \mathcal{S}_{\text{X}}\) that satisfies 18 . To be specific, after randomly initializing \((\mathbf{F}, \mathbf{X})\), the first subproblem is formulated to reduce the ISL, i.e., \(\sum_{k \in \mathbf{\Omega}_{K}} \| \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} \|_{\text{F}}^2\), concerning \(\mathbf{X}\). We update \(\mathbf{X}\) several times to minimize the ISL: if \(\max_{k \in \mathbf{\Omega}_{K}} \{\| \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}} \|_{\text{F}}\}\) is already less than \(\xi\), we then find a feasible point; otherwise, we formulate the second subproblem to minimize the ISL concerning \(\tilde{\mathbf{F}} = \mathbf{F}^{\text{H}} \mathbf{F}\). When solving for the variable \(\mathbf{X}\), the following lemma is useful.

3.2.1 Reducing the ISL Regarding \(\mathbf{X}\)↩︎

By choosing any \(\mathbf{F} \in \mathcal{S}_{\text{F}}\), and denoting \(\mathbf{\Theta} = \sum_{k \in \mathbf{\Omega}_{K}} \mathbf{J}_{k} \mathbf{F}^{\text{H}} \mathbf{F} \mathbf{J}_{k}^{\text{H}}\), the problem of reducing the ISL regarding \(\mathbf{X}\) is formulated as \[\label{prob:X95feasible} \begin{align} & \min_{\mathbf{X}} \mathrm{tr}(\mathbf{X} \mathbf{\Theta} \mathbf{X}^{\text{H}}) \\ & \;\text{s.t. } \mathbf{X} \in \mathcal{S}_{\text{X}}. \end{align}\tag{19}\]

According to Lemma 1, the tightest upper bound \(\lambda\) to majorize \(\mathrm{tr}(\mathbf{X} \mathbf{\Theta} \mathbf{X}^{\text{H}})\) is the maximum eigenvalue of \(\mathbf{\Theta}\), denoted by \(\lambda_{\max}(\mathbf{\Theta})\). However, \(\mathbf{\Theta}\) may change after updating \(\mathbf{F}\). To avoid the high computational cost of computing \(\lambda_{\max}(\mathbf{\Theta})\), an alternative upper bound \(\lambda_{\theta}\) is implemented, where \[\lambda_{\theta} = \|\mathbf{\Theta}\|_{1} \geq \lambda_{\max}(\mathbf{\Theta}). \nonumber\] Then, the quadratic term \(\mathrm{tr}(\mathbf{X} \mathbf{\Theta} \mathbf{X}^{\text{H}})\) can be majorized as follows: \[\begin{gather} \mathrm{tr}(\mathbf{X} \mathbf{\Theta} \mathbf{X}^{\text{H}}) \leq \lambda_{\theta} \|\mathbf{X}\|_{\text{F}}^2 \\ + 2 \mathcal{R}\left\{ tr\left[\left(\mathbf{X}^{(t)} \mathbf{\Theta} - \lambda_{\theta} \mathbf{X}^{(t)}\right)^{\text{H}} \mathbf{X} \right] \right\} + C_{\text{X}}, \end{gather}\] where \(C_{\text{X}}\) is the term irrelevant to \(\mathbf{X}\) and \(\mathbf{X}^{(t)}\) is the result of previous iteration. Thus, the majorized problem of 19 is \[\label{prob:majorized95feasible95xi} \begin{align} & \min_{\mathbf{X}} \mathcal{R}\left\{ tr\left[\left(\mathbf{X}^{(t)} \mathbf{\Theta} - \lambda_{\theta} \mathbf{X}^{(t)}\right)^{\text{H}} \mathbf{X} \right] \right\} \\ & \;\text{s.t. } \mathbf{X} \in \mathcal{S}_{\text{X}}. \end{align}\tag{20}\] It is clear that solving problem 20 is equivalent to solving the following problem: \[\label{prob:majorized95feasible95xi95simplified} \begin{align} & \min_{\mathbf{X}} \|\mathbf{X} - \mathbf{Z}^{(t)}\|_{\text{F}}^{2} \\ & \;\text{s.t. } \mathbf{X} \in \mathcal{S}_{\text{X}}, \end{align}\tag{21}\] where \(\mathbf{Z}^{(t)} = \lambda_{\theta} \mathbf{X}^{(t)} - \mathbf{X}^{(t)} \mathbf{\Theta}\), whose closed form solution is \(\mathbf{X}^{(t+1)} = e^{j \arg(\mathbf{Z}^{(t)})}\).

3.2.2 Minimizing the Sidelobe Level Regarding \(\mathbf{F}\)↩︎

Suppose that the obtained solution of problem 19 is not feasible. We further minimize the sidelobe level with respect to \(\mathbf{F}\) in this part. By denoting \(\mathbf{\Phi} = \sum_{k \in \mathbf{\Omega}_{K}} \mathbf{J}_{k}^{\text{H}} \mathbf{X}^{\text{H}} \mathbf{X} \mathbf{J}_{k}\) and its singular value decomposition (SVD) as \(\mathbf{\Phi} = \mathbf{U}_{\Phi} \mathbf{\Sigma}_{\Phi} \mathbf{U}_{\Phi}^{\text{H}}\) (assuming that the singular values are arranged in descending order), the ISL becomes \[\label{eq:trace95wt95F} \mathrm{tr}(\mathbf{Z} \mathbf{\Sigma}_{\Phi} \mathbf{Z}^{\text{H}}) = \sum_{i=1}^{N} \lambda_i(\mathbf{\Phi}) \|\mathbf{z}_i\|_2^2,\tag{22}\] where \(\mathbf{Z} = \mathbf{F} \mathbf{U}_{\Phi}\), \(\mathbf{z}_i\) is the \(i\)-th row vector of \(\mathbf{Z}\), and \(\lambda_i(\mathbf{\Phi})\) is the \(i\)-th singular value of \(\mathbf{\Phi}\). Since \(\mathbf{\Phi}\) is positive semidefinite, we obtain the minimum of 22 by setting all rows of \(\mathbf{Z}\) as zero except for the row corresponding to the minimum singular value of \(\mathbf{\Phi}\). That is to say, we obtain an \(\mathbf{F}\) that minimizes the ISL by setting \(\mathbf{F} = \mathbf{\Sigma}_{\text{F}} \mathbf{U}_{\Phi}^{\text{H}}\), where \(\mathbf{\Sigma}_{\text{F}} = \mathrm{Diag}(\mathbf{d}_{\text{F}})\) with \(\mathbf{d}_{\text{F}} = [\mathbf{0}_{N_{\text{T}} - 1}^{\text{T}}, \sqrt{P_{\text{F}}}]^{\text{T}}\), and \(\mathrm{Diag}(\mathbf{d}_{\text{F}})\) is a diagonal matrix with its main diagonal entries being \(\mathbf{d}_{\text{F}}\).

3.3 BSUM Algorithm for Inexactly Solving Subproblem 12↩︎

In this subsection, we propose a BSUM scheme to solve the ALM subproblem 12 inexactly. In particular, the initial point at the \(\ell\)-th iteration is chosen as follows: \[\label{eq:init95point95at95lth95ALg1} \boldsymbol{z}_{\text{init}}^{(\ell)} = \begin{cases} \boldsymbol{z}^{(\text{feas})}, & \text{if } f(\mathbf{F}^{(\text{feas})}, \mathbf{X}^{(\text{feas})}) < \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(\ell)}; \boldsymbol{U}^{(\ell)}); \\ \boldsymbol{z}^{(\ell)}, & \text{otherwise}, \end{cases}\tag{23}\] where \(\boldsymbol{z}^{(\ell)}\) is the approximate solution obtained at the \((\ell-1)\)-th iteration of Algorithm 3, and \(\boldsymbol{z}^{(\text{feas})}\) is a feasible point obtained from Section 3.2. It follows from 10 and the definition of \(\zeta\) in 11 that \[\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(\text{feas})}; \boldsymbol{U}^{(\ell)}) \leq f(\mathbf{F}^{(\text{feas})}, \mathbf{X}^{(\text{feas})}) \leq \zeta.\] This shows that the choice of the initial point in 23 guarantees \(\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}_{\text{init}}^{(\ell)}; \boldsymbol{U}^{(\ell)}) \leq \zeta\). Moreover, the proposed BSUM method has a sufficient descent property, which will be proved in Section 3.4. Therefore, for any \(\ell \geq 0\), we always have \[\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(\ell+1)}; \boldsymbol{U}^{(\ell)}) \\ \leq \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}_{\text{init}}^{(\ell)}; \boldsymbol{U}^{(\ell)}) \leq \zeta.\] Consequently, the uniform upper bound requirement in 13 can always be satisfied under the initial point choice strategy of 23 .

Denote \(\boldsymbol{z}^{(t)} = (\boldsymbol{C}^{(t)}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)})\). Next, we present the update of each block of variables using the BSUM method.

3.3.1 Solving for Auxiliary Variables↩︎

The update of auxiliary variables in \(\boldsymbol{C}\) can be formulated as \(|\mathbf{\Omega}_{K}|\) independent optimization problems for each \(\mathbf{C}_{k}\). Instead of solving 12 with respect to (w.r.t.) each \(\mathbf{C}_{k}\) directly, we consider the following subproblem: \[\begin{align} \label{prob:block95aux95c} & \min_{\mathbf{C}_{k}} \frac{\rho^{(\ell)}}{2} \left\| \mathbf{X}^{(t)} \mathbf{J}_{k} {\mathbf{F}^{(t)}}^{\text{H}} - \mathbf{C}_{k} + \frac{\mathbf{U}_{k}^{(\ell)}}{\rho^{(\ell)}}\right\|_{\text{F}}^2 + \frac{\beta}{2} \left\|\mathbf{C}_{k} - \mathbf{C}_{k}^{(t)}\right\|_{\text{F}}^2, \nonumber \\ & \text{ s.t. } \|\mathbf{C}_{k}\|_{\text{F}} \leq \xi, \end{align}\tag{24}\] where \(\beta > 0\). As mentioned in [36], by setting \(\beta > 0\), the distance between two consecutive iterations of each \(\mathbf{C}_{k}\) becomes controllable, and every update of \(\mathbf{C}_{k}\) achieves a sufficient decrease. It is clear that problem 24 has a closed-form solution as follows: \[\label{eq:update95C} \mathbf{C}_{k}^{(t+1)} = \begin{cases} \tilde{\mathbf{C}}_{k}, & \text{if } \|\tilde{\mathbf{C}}_{k}\|_{\text{F}} \leq \xi; \\ \frac{\xi}{\|\tilde{\mathbf{C}}_{k}\|_{\text{F}}} \tilde{\mathbf{C}}_{k}, & \text{otherwise}, \end{cases}\tag{25}\] where \(\tilde{\mathbf{C}}_{k} = \frac{1}{\rho^{(\ell)} + \beta} \left(\rho^{(\ell)}\mathbf{X}^{(t)} \mathbf{J}_{k} {\mathbf{F}^{(t)}}^{\text{H}} + \mathbf{U}_{k}^{(\ell)} + \beta \mathbf{C}_{k}^{(t)}\right)\).

3.3.2 Solving for Radar Receive Filter↩︎

The problem of solving 12 w.r.t. \(\mathbf{F}\) can be rewritten as \[\label{prob:PSL95F951} \begin{align} & \min_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}) \\ & \;\text{s.t. } \| \mathbf{F} \|_{\text{F}}^2 = P_{\text{F}}, \end{align}\tag{26}\] where \(\mathcal{L}_{\text{F}}(\mathbf{F})\) is defined as follows: \[\begin{gather} \label{eq:LAF95F95PSL} \mathcal{L}_{\text{F}}(\mathbf{F}) = \mathrm{tr}(\mathbf{F} \mathbf{Q} \mathbf{F}^{\text{H}}) - 2 \alpha \mathcal{R}\{\mathrm{tr}(\mathbf{R}_{\text{d}}^{\text{H}} \mathbf{X} \mathbf{F}^{\text{H}})\} \\ - \sum_{k \in \mathbf{\Omega}_{K}} \mathcal{R}\{tr[(\rho \mathbf{C}_{k} - \mathbf{U}_{k})^{\text{H}} \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}}]\}, \end{gather}\tag{27}\] which is reformulated from 10 by ignoring the terms irrelevant to \(\mathbf{F}\) and denoting \(\mathbf{Q} = \alpha \mathbf{X}^{\text{H}} \mathbf{X} + \sum_{k \in \mathbf{\Omega}_{K}} \frac{\rho}{2}\mathbf{J}_{k}^{\text{H}} \mathbf{X}^{\text{H}} \mathbf{X} \mathbf{J}_{k}\), and the superscripts of \(\mathbf{X}^{(t)}\), \(\boldsymbol{C}^{(t+1)}\), \(\boldsymbol{U}^{(\ell)}\), and \(\rho^{(\ell)}\) are omitted for notational simplicity in this part. The quadratic term \(\mathrm{tr}(\mathbf{F} \mathbf{Q} \mathbf{F}^{\text{H}})\) can be majorized at \(\mathbf{F}^{(t)}\) through Lemma 1 as follows: \[\label{eq:F95295MM} \mathrm{tr}(\mathbf{F} \mathbf{Q} \mathbf{F}^{\text{H}}) \leq \lambda_{q} \|\mathbf{F}\|_{\text{F}}^2 + 2 \mathcal{R} \{tr[(\mathbf{F}^{(t)}(\mathbf{Q} - \lambda_{q}\mathbf{I}))^{\text{H}}\mathbf{F}]\} + C_{\text{F}},\tag{28}\] where \(\lambda_{q} > \|\mathbf{Q}\|_{1}\), and \(C_{\text{F}}\) is the term irrelevant to \(\mathbf{F}\).

By substituting 28 into 27 , we have \(\mathcal{G}_{\text{F}}(\mathbf{F}\mid\mathbf{F}^{(t)})\) as the tight upper bound of \(\mathcal{L}_{\text{F}}(\mathbf{F})\), i.e., \(\mathcal{G}_{\text{F}}(\mathbf{F}\mid\mathbf{F}^{(t)}) \geq \mathcal{L}_{\text{F}}(\mathbf{F})\) for all \(\mathbf{F} \in \mathcal{S}_{\text{F}}\) with equality holding when \(\mathbf{F} = \mathbf{F}^{(t)}\), and \[\label{eq:neg95F95MM2} \mathcal{G}_{\text{F}}(\mathbf{F}\mid\mathbf{F}^{(t)}) = \lambda_{q} \|\mathbf{F}\|_{\text{F}}^2 - 2 \mathcal{R} \left\{tr \left( {\mathbf{\Xi}^{(t)}}^{\text{H}} \mathbf{F} \right)\right\} + C_{\text{F}},\tag{29}\] where \[\mathbf{\Xi}^{(t)} = \sum_{k \in \mathbf{\Omega}_{K}} \frac{1}{2}(\rho \mathbf{C}_{k} - \mathbf{U}_{k})^{\text{H}} \mathbf{X} \mathbf{J}_{k} + \alpha \mathbf{R}_{\text{d}}^{\text{H}} \mathbf{X} - \mathbf{F}^{(t)} \mathbf{Q} + \lambda_{q} \mathbf{F}^{(t)}.\] Then, the majorized problem of 26 is \[\label{prob:PSL95F95MMed95lg0} \begin{align} & \min_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}\mid\mathbf{F}^{(t)}) \\ & \;\text{s.t. } \| \mathbf{F} \|_{\text{F}}^2 = P_{\text{F}}. \end{align}\tag{30}\] Problem 30 has a closed-form solution as follows: \[\label{sol:F} \mathbf{F}^{(t+1)} = \frac{\sqrt{P_{\text{F}}}}{\| \mathbf{\Xi}^{(t)} \|_{\text{F}}}\mathbf{\Xi}^{(t)}.\tag{31}\]

3.3.3 Solving for Transmit Waveform↩︎

The problem of solving 12 w.r.t. \(\mathbf{X}\) can be rewritten as \[\label{prob:X} \begin{align} & \min_{\mathbf{X}} \mathcal{L}_{\text{X}} (\mathbf{X}), \\ & \; \text{s.t. } \mathbf{X} \in \mathcal{S}_{\text{X}}, \end{align}\tag{32}\] where \[\begin{gather} \label{eq:ALF95wrt95X} \mathcal{L}_{\text{X}}(\mathbf{X}) = \mathrm{tr}(\mathbf{X} \mathbf{P} \mathbf{X}^{\text{H}}) - 2 \alpha \mathcal{R}\{\mathrm{tr}(\mathbf{X} \mathbf{F}^{\text{H}} \mathbf{R}_{\text{d}}^{\text{H}})\} \\ - \sum_{k \in \mathbf{\Omega}_{K}} \mathcal{R}\{tr[(\rho \mathbf{C}_{k} - \mathbf{U}_{k})^{\text{H}} \mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}}]\}\\ + (1 - \alpha) (\mathbf{x}^{\text{H}} \mathbf{H}^{\text{H}} \mathbf{H} \mathbf{x} - 2 \mathcal{R}\{\mathbf{x}^{\text{H}} \mathbf{H}^{\text{H}} \mathbf{s}\}), \end{gather}\tag{33}\] and \(\mathbf{P} = \alpha {\mathbf{F}^{(t+1)}}^{\text{H}} \mathbf{F}^{(t+1)} + \frac{\rho}{2} \sum_{k \in \mathbf{\Omega}_{K}} \mathbf{J}_{k} {\mathbf{F}^{(t+1)}}^{\text{H}} \mathbf{F}^{(t+1)} \mathbf{J}_{k}^{\text{H}}\).

We consider majorizing the two quadratic terms, \(\mathrm{tr}(\mathbf{X} \mathbf{P} \mathbf{X}^{\text{H}})\) and \(\mathbf{x}^{\text{H}} \mathbf{H}^{\text{H}} \mathbf{H} \mathbf{x}\) separately. According to Lemma 1, we have \[\label{eq:mm95195X} \mathrm{tr}(\mathbf{X} \mathbf{P} \mathbf{X}^{\text{H}}) \leq \lambda_{p} \|\mathbf{X}\|_{\text{F}}^2 + 2 \mathcal{R}\{\mathrm{tr}(\mathbf{X}^{(t)} \mathbf{P} - \lambda_{p} \mathbf{X}^{(t)})^{\text{H}} \mathbf{X}\} + C_{\text{X}},\tag{34}\] where \(\lambda_{p} > \|\mathbf{P}\|_{1}\), and \(C_{\text{X}}\) is the term irrelevant to \(\mathbf{X}\). Similarly, \[\label{eq:mm95295x} \mathbf{x}^{\text{H}} \mathbf{H}^{\text{H}} \mathbf{H} \mathbf{x} \leq \lambda_{h} \|\mathbf{x}\|_{2}^{2} + 2 \mathcal{R} \{\mathbf{x}^{\text{H}} (\mathbf{H}^{\text{H}} \mathbf{H} - \lambda_h \mathbf{I}) \mathbf{x}^{(t)}\} + C_{\text{x}},\tag{35}\] in which \(\lambda_{h} = \lambda_{\max}(\mathbf{H}^{\text{H}} \mathbf{H})\) and \(C_{\text{x}}\) is irrelevant to \(\mathbf{x}\). By substituting 34 and 35 into 33 , the tight upper bound of \(\mathcal{L}_{\text{X}}(\mathbf{X})\) is \[\mathcal{G}_{\text{X}}(\mathbf{X} | \mathbf{X}^{(t)}) = [\lambda_{p} + (1 - \alpha) \lambda_{h}] \|\mathbf{X}\|_{\text{F}}^2 - 2 \mathcal{R}\{\mathrm{tr}[{\mathbf{\Psi}^{(t)}}^{\text{H}} \mathbf{X}]\} + \tilde{C}_{\text{x}},\] where \(\mathbf{\Psi}^{(t)} = \alpha \mathbf{R}_{\text{d}} \mathbf{F} + \sum_{k \in \mathbf{\Omega}_{K}} \frac{1}{2} (\rho \mathbf{C}_{k} - \mathbf{U}_{k}) \mathbf{F} \mathbf{J}_{k}^{\text{H}} - (1 - \alpha) {\mathbf{\Phi}^{(t)}}^{\text{T}} - \mathbf{X}^{(t)} (\mathbf{P} - \lambda_{p} \mathbf{I})\) with \(\mathbf{\Phi}^{(t)} = \mathrm{mat}\{(\mathbf{H}^{\text{H}} \mathbf{H} - \lambda_h \mathbf{I}) \mathbf{x}^{(t)} - \mathbf{H}^{\text{H}} \mathbf{s}\}\), \(\mathrm{mat}\{\cdot\}\) means reshaping a column vector into a matrix, and \(\tilde{C}_{\text{x}} = (1 - \alpha) C_{\text{x}} + C_{\text{X}}\). Then, the majorized problem of 32 is \[\begin{align} & \min_{\mathbf{X}} \mathcal{G}_{\text{X}}(\mathbf{X} | \mathbf{X}^{(t)}) \\ & \; \text{s.t. } \mathbf{X} \in \mathcal{S}_{\text{X}}, \end{align}\] whose closed-form solution is \[\label{sol:X} \mathbf{X}^{(t+1)} = e^{j \mathrm{arg}(\mathbf{\Psi}^{(t)})},\tag{36}\]

The proposed BSUM algorithm for solving the ALM problem 12 inexactly is summarized in Algorithm 4. It is worth mentioning that the SQUAREM scheme [34] is implemented when updating \(\mathbf{F}\) and each row of \(\mathbf{X}\) in Algorithm 4 to accelerate the convergence. The BSUM algorithm performs efficiently since every variable admits closed-form updates. The computational cost of updating \(\mathbf{F}\), \(\mathbf{X}\), \(\boldsymbol{C}\), and \(\boldsymbol{U}\) is dominated by calculating \(\mathbf{\Xi}\), \(\mathbf{\Psi}\), and \(\mathbf{X} \mathbf{J}_{k} \mathbf{F}^{\text{H}}\) for each \(k \in \mathbf{\Omega}_{K}\), respectively, whose computational complexities are \(\mathcal{O}(|\mathbf{\Omega}_{K}| \cdot [(N_{\text{T}} + T) N_{\text{R}} T])\), \(\mathcal{O}(|\mathbf{\Omega}_{K}| \cdot [(N_{\text{R}} + T) N_{\text{T}} T])\), \(\mathcal{O}(|\mathbf{\Omega}_{K}| \cdot [(N_{\text{R}} + T) N_{\text{T}} T])\), and \(\mathcal{O}(|\mathbf{\Omega}_{K}| \cdot [(N_{\text{R}} + T) N_{\text{T}} T])\), respectively. Therefore, the total per-iteration complexity of Algorithm 4 is \(\mathcal{O}(|\mathbf{\Omega}_{K}| \cdot [(N_{\text{T}}T + N_{\text{R}} T + N_{\text{T}} N_{\text{R}} ) T])\), which scales linearly with the number of range bins for sidelobe suppression \(|\mathbf{\Omega}_{K}|\), the number of transmit antennas \(N_{\text{T}}\), and the number of receive antennas \(N_{\text{R}}\).

3.4 Convergence Analysis↩︎

The convergence analysis of the proposed inexact ALM algorithm (i.e., Algorithm 3) consists of four parts. We first show the sequence of the variables generated by Algorithm 4 enjoys a sufficient descent property in Lemma 2; then, we show that the subgradient of \(\mathcal{L}_{\rho}(\boldsymbol{C}, \mathbf{F}, \mathbf{X}; \boldsymbol{U})\) is bounded in Lemma 3; based on the results of Lemma 2 and Lemma 3, we show that Algorithm 4 achieves an \(\varepsilon\)-stationary point within a finite number of iterations in Theorem 1; finally, we prove that any limit point generated by Algorithm 3 is a feasible stationary point in Theorem 2.

The convergence analysis relies significantly on the uniform boundness of the primal variables \((\boldsymbol{C}, \mathbf{F}, \mathbf{X})\) and the multiplier \(\boldsymbol{U}\) generated by Algorithm 4. The boundness of \(\boldsymbol{U}\) can be guaranteed by the updating rule proposed in 15 , and the boundness of \((\boldsymbol{C}, \mathbf{F}, \mathbf{X})\) can be ensured by the updating rules proposed in 25 , 31 , and 36 .

4 Numerical Results↩︎

In this section, we present simulation results to evaluate the performance of the proposed algorithms under various parameter configurations. The communication channel \(\mathbf{H}\) is generated by following [32] with Extended Pedestrian A (EPA) fading profile [14]. The information symbol matrix \(\mathbf{S}\) is modulated using a unit-power QPSK alphabet, with each entry in \(\mathbf{S}\) having unit power. The desired radar receive spatial response, \(\mathbf{R}_{\text{d}}\), is considered to be a \(3\) dB beamwidth of \(20^{\circ}\) focusing at \(0^{\circ}\), generated by the algorithm in [40]. It is worth noting that the algorithm proposed in this paper can be directly extended to the Extremely large-scale MIMO (XL-MIMO) system [7] for generating unimodular low sidelobe level ISAC sequences.

The stopping criterion for Algorithm 4 at the \(\ell\)-th outer loop is set as \(\varepsilon^{(\ell)} = \sup_{\boldsymbol{z}^{(\ell)}} \|\boldsymbol{J}^{(0)}\|_{\text{F}}/\ell\), where \(\sup_{\boldsymbol{z}^{(\ell)}} \|\boldsymbol{J}^{(0)}\|_{\text{F}}\) is the subgradient upper bound in ?? . Unless otherwise specified, the parameter configurations are listed in Table 1. Algorithm 3 terminates after 500 iterations or when the following optimality violation criterion is met \[\label{eq:criteria} \max\{e^{(\ell)}, v^{(\ell)}\} \leq \sqrt{T} \times 10^{-3},\tag{37}\] where \(v^{(\ell)}\) is defined in 16 , and \[e^{(\ell)} = \sup_{\boldsymbol{z}^{(\ell)}} \|\mathbf{J}_{\text{F}}^{(\ell)}\|_{\text{F}} + \sup_{\boldsymbol{z}^{(\ell)}} \|\mathbf{J}_{\text{X}}^{(\ell)}\|_{\text{F}} + \sum_{k \in \mathbf{\Omega}_{K}} \sup_{\boldsymbol{z}^{(\ell)}} \|\mathbf{J}_{C_{k}}^{(\ell)}\|_{\text{F}}, \nonumber\] which are the summation of the upper bounds in 43 , 48 , and 50 .

4.1 Convergence Performance↩︎

In this subsection, we evaluate the convergence behavior of the proposed inexact ALM algorithm using 1000 Monte Carlo (MC) simulations. In each simulation, the channel matrix \(\mathbf{H}\) and the information symbol matrix \(\mathbf{S}\) are generated independently and randomly, with the variables \((\boldsymbol{C}, \mathbf{F}, \mathbf{X}), \boldsymbol{U}\) also randomly initialized. The predetermined feasible points are consistent across all MC simulations. The shaded areas in the figures represent the standard deviation of the results from the 1000 simulations.

Fig. 5 shows the average behavior of Algorithm 3. Fig. 5 plots the optimality violation in 37 alongside the stopping threshold (dashed line) for Algorithm 3. On average, the algorithm reaches the stopping criterion after approximately 480 iterations. Due to the random initialization of variables, the optimized sum rate and beampattern MSE vary. Fig. 5 and Fig. 5 illustrate the average curves of sum rate, beampattern MSE, and sidelobe level, showing that the final results perform well on average, with the maximum sidelobe level approaching the constraint. Fig. 5 shows the number of inner iterations required to find the \(\varepsilon^{(\ell)}\)-stationary point. Although the maximum allowed inner iterations is 50, over 99.5\(\%\) of the inner loops converge within 10 iterations. Combined with closed-form updates and the “nearest-vector” algorithm, the proposed BSUM algorithm efficiently solves the ALM subproblems.

Fig. 6 shows the average cross-correlation level and radar receive beampattern of the sequences generated by Algorithm 3 after convergence. The results indicate that all sidelobes in the interested area are well suppressed, and the optimized radar beampattern closely matches the desired one.

4.2 Communication and Radar Performance under Different System Configurations↩︎

In this subsection, we evaluate the communication and radar performance of the proposed algorithm under different system configurations. Specifically, we examine the impact of the range sidelobe suppression area size determined by \(K\) and the maximum sidelobe level constraints \(\xi^{\prime}\) (dB).

Fig. 7 illustrates the impact of the sidelobe suppression area sizes determined by \(K\) on the optimized results. As shown in Fig. 7, even with the increment of \(K\), the optimized sidelobe levels can still approach the predetermined threshold. In Fig. 7, the achievable sum rate decreases, the beampattern MSE increases, and the mainlobe level decreases slightly (shown at the bottom-left of Fig. 7) as \(K\) grows, since a larger suppression area imposes more constraints, reducing the feasible region and degrading the performance of optimized results.

Fig. 8 shows the impact of maximum sidelobe level constraints on the optimized results. As seen in Fig. 8, the sidelobe constraints can always be met. Lower sidelobe constraints require \(\mathbf{X}\) and \(\mathbf{F}\) to be less correlated, making the problem more complex and resulting in worse performance in both S\(\&\)C. This can be observed in Fig. 8, in which a lower maximum sidelobe level constraint leads to a worse achievable sum rate and a larger beampattern MSE.

4.3 Comparison with Existing SOTA Approaches↩︎

In this subsection, we compare the radar and communication performance of the proposed approach with the modified work [29] and the ALM algorithm with a fixed penalty parameter. In [29], the joint design of the radar receive filter and transmit waveform focuses on PSL control for pure radar sensing. We adapt the algorithm in [29] by first solving \(\min_{\mathbf{X} \in \mathcal{S}_{\text{X}}} \|\mathbf{Hx} - \mathbf{s}\|_{2}^2\) using the BSUM method and then following the approach in [29] to obtain the radar receive filter \(\mathbf{F}\). We refer to the modified algorithm of [29] as the “Exact Design” scheme. Furthermore, we also implemented the ADMM algorithm and BSUM-based ALM algorithm with a fixed penalty parameter to show the superiority of the proposed adaptive penalty parameter scheme. We refer to them as “ADMM” and “Fixed Penalty” schemes, respectively. Without loss of generality, the penalty parameters for both the two algorithms are fixed to be \(1\), and the total iterations of the “ADMM” algorithm are set the same as the proposed ALM algorithm. In the legends of the result figures, “Direc” refers to the receive beampattern focused at \(-30^{\circ}\), \(0^{\circ}\), and \(30^{\circ}\) with a \(3\) dB beamwidth of \(10^{\circ}\), generated by the algorithm in [40], and “Omni” denotes the omnidirectional desired beampattern, i.e., \(\mathbf{R}_{\text{d}} = (P_{\text{X}}/N_{T})\mathbf{I}_{N_{T}}\).

Fig. 9 presents the optimized results of different algorithms. In Fig. 9, the proposed ALM algorithm achieves the highest sum rate under the directional beampattern case. For an omnidirectional beam, we need to set a larger \(\alpha\) to make the final beampattern closer to omnidirectional because an omnidirectional beam requires that \(\mathbf{X}\) and \(\mathbf{F}\) be orthogonal in each row and column, which significantly shrinks the feasible region. However, this degrades the minimization on \(P_{\text{MUI}}\) and leads to a lower achievable rate. In the Exact Design scheme, although the transmit sequence is optimized solely to minimize the MUI at the receiver, it still has a lower spectral efficiency compared to the result of the proposed ALM algorithm. This suggests that directly minimizing \(\|\mathbf{H}\mathbf{x} - \mathbf{s}\|_2^2\) is prone to suboptimal local minima. In contrast, our proposed algorithm demonstrates a strong ability to escape such local minima and attain better solutions, even under more stringent constraints. Other methods, due to the inappropriate penalty parameter, provide lower achievable rates. Fig. 9 shows that the Exact Design scheme nearly perfectly matches the desired beampatterns in both the directional and omnidirectional cases while the proposed inexact ALM algorithms have slight differences from the omnidirectional beampattern. The received beampatterns under the Fixed Penalty scheme and the ADMM scheme were far from the desired one because the inappropriate penalty parameter makes the algorithm focus too much on reducing the violation. Fig. 9 shows the proposed ALM algorithm has the best sidelobe control performance.

Fig. 10 shows the MIMO radar images formed by the results of different algorithms⁸. The proposed inexact ALM algorithm achieved the best imaging performance. Although the Exact Design scheme perfectly meets the desired beampattern, it has a high sidelobe level, which results in severe interference caused by the clutter at the range bin of interest and, thus, a poorer MIMO radar image. Due to the inappropriate penalty parameter, the results under both the Fixed Penalty scheme and the ADMM scheme failed to produce radar images.

5 Conclusion↩︎

This paper addressed the joint design of the receive filter and transmit waveform for MIMO-ISAC systems. We formulated an optimization problem to minimize a weighted sum of radar beampattern MSE and MUI at communication receivers, subject to practical constraints in convex and nonconvex forms. An inexact ALM algorithm was developed to solve the problem by iteratively minimizing tight upper bounds for each variable block. We proved the algorithm’s convergence to a feasible stationary point, which is the best that one can expect for this optimization problem with many nonconvex constraints. The trade-offs between sum rate, beampattern MSE, sidelobe suppression size, and maximum sidelobe level were analyzed. Simulation results demonstrated that the proposed algorithm outperforms others across various metrics.

6 Proof of Lemma 2↩︎

From 24 , it is clear that \[\begin{gather} \label{eq:C95ineq} \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{C}^{(t)}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{C}^{(t+1)}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)}; \boldsymbol{U}^{(\ell)}) \\ \geq \frac{\beta}{2} \sum_{k \in \mathbf{\Omega}_{K}} \|\mathbf{C}_{k}^{(t)} - \mathbf{C}_{k}^{(t+1)}\|_{\text{F}}^2. \end{gather}\tag{38}\] It is simple to verify that updating \(\mathbf{F}^{(t+1)}\) by solving problem 30 is equivalent to \[\begin{gather} \label{prob:prox95F} \mathbf{F}^{(t+1)} \in \mathop{\mathrm{arg\,min}}_{\mathbf{F} \in \mathbb{C}^{N_{\text{T}} \times L}} \{ \mathcal{R}\{\mathrm{tr}[\nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t)})^{\text{H}} (\mathbf{F} - \mathbf{F}^{(t)})]\} \\ + \frac{\lambda_{q}^{(t)}}{2} \|\mathbf{F} - \mathbf{F}^{(t)}\|_{\text{F}}^2 + \mathbb{I}_{\mathcal{S}_{\text{F}}}(\mathbf{F})\}, \end{gather}\tag{39}\] where \(\lambda_{q}^{(t)}\) is defined in 28 that relies on \(\mathbf{X}^{(t)}\). From Lemma 2 in [42], we have \[\mathcal{L}_{\text{F}}(\mathbf{F}^{(t)}) - \mathcal{L}_{\text{F}}(\mathbf{F}^{(t+1)}) \geq \frac{\lambda_{q}^{(t)} - L_{f}^{(t)}}{2}\|\mathbf{F}^{(t+1)} - \mathbf{F}^{(t)}\|_{\text{F}}^2,\nonumber\] where \(L_{f}^{(t)}\) denotes the Lipschitz constant of \(\nabla_{\mathbf{F}}^{*} \mathcal{L}_{\text{F}}(\mathbf{F})\) at the \(t\)-th iteration of Algorithm 4, and \(\nabla_{\mathbf{F}}^{*} \mathcal{L}_{\text{F}}(\mathbf{F})\) is the conjugate gradient of \(\mathcal{L}_{\text{F}}(\mathbf{F})\). We always have \(\lambda_{q}^{(t)} - L_{f}^{(t)} > 0\), since \(L_{f}^{(t)}\) is the maximum eigenvalue of \(\mathbf{Q}^{(t)}\) according to the results in [43], where \(\mathbf{Q}^{(t)}\) is defined in 28 and its superscript \(^{(t)}\) means it relies on \(\mathbf{X}^{(t)}\). Then \(\lambda_{q}^{(t)} > L_{f}^{(t)}\) for \(\forall~t > 0\) according to 28 . Therefore, there exists a \(\tau_{f} > 0\), such that the following holds for all \(t > 0\): \[\begin{gather} \label{eq:Ft95Ft195ineq} \mathcal{L}_{\rho^{(\ell)}}(\mathbf{C}^{(t+1)}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\mathbf{C}^{(t+1)}, \mathbf{F}^{(t+1)}, \\ \mathbf{X}^{(t)}; \boldsymbol{U}^{(\ell)})\geq \tau_{f} \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}^2. \end{gather}\tag{40}\] By applying the above analysis procedure to 32 , we obtain \[\begin{gather} \label{eq:Xt95Xt195ineq} \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{C}^{(t+1)}, \mathbf{F}^{(t+1)}, \mathbf{X}^{(t)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{C}^{(t+1)}, \mathbf{F}^{(t+1)}, \\ \mathbf{X}^{(t+1)}; \boldsymbol{U}^{(\ell)})\geq \tau_{x} \|\mathbf{X}^{(\ell)} - \mathbf{X}^{(t+1)}\|_{\text{F}}^{2}, \end{gather}\tag{41}\] where \(\tau_{x} > 0\). Combining 38 , 40 , and 41 , we have \[\begin{gather} \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(t)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(t+1)}; \boldsymbol{U}^{(\ell)})\geq \tau_{f} \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}^2 \\ + \tau_{x} \|\mathbf{X}^{(t)} - \mathbf{X}^{(t+1)}\|_{\text{F}}^{2} + \frac{\beta}{2} \sum_{k \in \mathbf{\Omega}_{K}} \|\mathbf{C}_{k}^{(t)} - \mathbf{C}_{k}^{(t+1)}\|_{\text{F}}^2. \end{gather}\] Therefore, Lemma 2 holds.

7 Proof of Lemma 3↩︎

7.0.1 Upper Bound of \(\mathbf{J}_{\rm F}^{(t+1)}\)↩︎

By solving problem 30 , we have \(\mathbf{0} \in \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t+1)}\mid\mathbf{F}^{(t)}) + \partial \mathbb{I}_{\mathcal{S}_{\text{F}}}(\mathbf{F}^{(t+1)})\). Hence, there exists a subgradient \(\boldsymbol{A}^{(t+1)} \in \partial \mathbb{I}_{\mathcal{S}_{\text{F}}}(\mathbf{F}^{(t+1)})\) such that \[\label{eq:subgrad95GF} \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t+1)}\mid\mathbf{F}^{(t)}) + \boldsymbol{A}^{(t+1)} = \mathbf{0},\tag{42}\] which further implies that there exists a subgradient \(\mathbf{J}_{\text{F}}^{(t+1)} \in \partial_{\mathbf{F}} [\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}; \boldsymbol{U}^{(\ell)}) + \mathbb{I}_{\mathcal{S}_{\text{F}}}(\mathbf{F})] \mid _{\mathbf{F} = \mathbf{F}^{(t)}}\) such that \[\label{eq:subgrad95LF} \mathbf{J}_{\text{F}}^{(t+1)} = \nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t+1)}) + \boldsymbol{A}^{(t+1)}.\tag{43}\] By combining 42 and 43 , we have \[\begin{gather} \label{eq:subgrad95LF951} \|\mathbf{J}_{\text{F}}^{(t+1)}\|_{\text{F}} = \|\nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t+1)}) - \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t+1)}\mid\mathbf{F}^{(t)})\|_{\text{F}} \\ \overset{(a)}{\leq} \|\nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t+1)}) - \nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t)})\|_{\text{F}} \\ + \|\nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t+1)}\mid\mathbf{F}^{(t)}) - \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t)}\mid\mathbf{F}^{(t)})\|_{\text{F}}, \end{gather}\tag{44}\] where \((a)\) holds due to \(\nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t)}) = \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t)}\mid\mathbf{F}^{(t)})\). Since \[\begin{gather} \label{ieq:L95F95f} \quad\|\nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t+1)}) - \nabla_{\mathbf{F}} \mathcal{L}_{\text{F}}(\mathbf{F}^{(t)})\|_{\text{F}} \\ = \|\mathbf{Q} (\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)})\|_{\text{F}} \leq \|\mathbf{Q}\|_{\text{F}} \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}, \end{gather}\tag{45}\] and \[\begin{gather} \label{ieq:G95F95f} \|\nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t+1)}\mid\mathbf{F}^{(t)}) - \nabla_{\mathbf{F}} \mathcal{G}_{\text{F}}(\mathbf{F}^{(t)}\mid\mathbf{F}^{(t)})\|_{\text{F}} \\ = \lambda_{q}\|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}} \leq \|\mathbf{Q}\|_{\text{F}} \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}, \end{gather}\tag{46}\] it follows that \[\label{eq:sub95grad95bound95F} \|\mathbf{J}_{\text{F}}^{(t+1)}\|_{\text{F}} \leq \bar{L}_{\text{F}} \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}},\tag{47}\] where the inequality holds due to \[\begin{gather} \|\mathbf{Q}\|_{\text{F}} \leq \alpha \|\mathbf{X}\|_{\text{F}} + \sum_{k \in \mathbf{\Omega}_{K}} \frac{\rho^{(\ell)}}{2} \|\mathbf{X} \mathbf{J}_{k}\|_{\text{F}} \\ \leq \sqrt{P_{\text{X}}} \left(\alpha + \frac{1}{2}\rho^{(\ell)} \sqrt{T} |\mathbf{\Omega}_{K}|\right) = \frac{\bar{L}_{\text{F}}}{2}. \nonumber \end{gather}\]

7.0.2 Upper Bound of \(\mathbf{J}_{\rm X}^{(t+1)}\)↩︎

By following the similar analysis procedure in the previous part, we have \[\label{eq:sub95grad95bound95X} \|\mathbf{J}_{\text{X}}^{(t+1)}\|_{\text{F}} \leq \bar{L}_{\text{X}} \|\mathbf{X}^{(t)} - \mathbf{X}^{(t+1)}\|_{\text{F}},\tag{48}\] where \(\bar{L}_{\text{X}} = 2 [\sqrt{P_{\text{F}}} (\alpha + \frac{1}{2} \rho^{(\ell)} \sqrt{T} |\mathbf{\Omega}_{K}|) + (1-\alpha)\|\mathbf{H}^{\text{H}}\mathbf{H}\|_{\text{F}}]\).

7.0.3 Upper Bound of \(\mathbf{J}_{C}^{(t+1)}\)↩︎

We now consider calculating the upper bound of the subgradient vector \[\mathbf{J}_{C}^{(t+1)} = \left[(\mathbf{J}_{C_{1}}^{(t+1)})^{\text{T}}, (\mathbf{J}_{C_{2}}^{(t+1)})^{\text{T}}, \dots, (\mathbf{J}_{C_{|\mathbf{\Omega}_{K}|}}^{(t+1)})^{\text{T}}\right]^{\text{T}}, \nonumber\] where \(\mathbf{J}_{C_{k}}^{(t+1)}\) denotes the subgradient of \(\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}; \boldsymbol{U}^{(\ell)}) + \mathbb{I}_{\mathcal{S}_{\text{C}}}(\mathbf{C}_{k})\) w.r.t. \(\mathbf{C}_{k}\) at \(\mathbf{C}_{k}^{(t+1)}\). Denoting the objective function in 24 as \(\mathcal{H}_{\rho^{(\ell)}, k}(\mathbf{C}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)})\) for each \(k \in \mathbf{\Omega}_{K}\), we have \[\label{eq:subgrad95Ck} \mathbf{0} = \nabla_{\mathbf{C}_{k}} \mathcal{H}_{\rho^{(\ell)}, k}(\mathbf{C}_{k}^{(t+1)}, \mathbf{F}^{(t)}, \mathbf{X}^{(t)}) + \boldsymbol{D}_{k}^{(t+1)},\tag{49}\] where \(\boldsymbol{D}_{k}^{(t+1)} \in \partial \mathbb{I}_{\mathcal{S}_{\text{C}}}(\mathbf{C}_{k}^{(t+1)})\). Then there exists a subgradient such that \[\mathbf{J}_{C_{k}}^{(t+1)} = \nabla_{\mathbf{C}_{k}} \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z};\boldsymbol{U}^{(\ell)}) \mid _{\boldsymbol{z} = \boldsymbol{z}^{(t+1)}} + \boldsymbol{D}_{k}^{(t+1)}.\]

Following the similar derivation procedure regarding \(\mathbf{J}_{\text{F}}^{(t+1)}\), we have the upper bound of \(\|\mathbf{J}_{C_{k}}\|_{\text{F}}\). By summing \(\mathbf{J}_{C_{k}}^{(t+1)}\) over \(k \in \mathbf{\Omega}_{K}\), we finally get \[\begin{gather} \label{eq:sub95grad95bound95C} \|\mathbf{J}_{C}^{(t+1)}\|_{\text{F}} \leq \frac{1}{2} \rho^{(\ell)} \sqrt{T} |\mathbf{\Omega}_{K}| \left(\sqrt{P_{\text{F}}} \| (\mathbf{X}^{(t)} - \mathbf{X}^{(t+1)}) \|_{\text{F}} \right. \\ \left. + \sqrt{P_{\text{X}}}\|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}\right)+ \beta \sum_{k \in \mathbf{\Omega}_{K}} \| (\mathbf{C}_{k}^{(t+1)} - \mathbf{C}_{k}^{(t)}) \|_{\text{F}}. \end{gather}\tag{50}\] Combining 47 , 48 , and 50 together, we arrive at Lemma 3.

8 Proof of Theorem 1↩︎

From the results of Lemma 2 and Lemma 3, with penalty parameter \(\rho^{(\ell)}\) at the \(\ell\)-th iteration of Algorithm 3, there exists a positive real number \(M>0\) such that \[\begin{gather} \label{eq:ALF95subgrad} M [\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(t)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(t+1)}; \boldsymbol{U}^{(\ell)})] \geq M [\|\mathbf{X}^{(t)} \\ - \mathbf{X}^{(t+1)}\|_{\text{F}}^2 + \|\mathbf{F}^{(t)} - \mathbf{F}^{(t+1)}\|_{\text{F}}^2 + \frac{\beta}{2} \sum_{k \in \mathbf{\Omega}_{K}} \|\mathbf{C}_{k}^{(t)}- \mathbf{C}_{k}^{(t+1)}\|_{\text{F}}^2] \\ \geq \frac{1}{{\rho^{(\ell)}}^2} (\|\mathbf{J}_{\text{F}}^{(t+1)}\|_{\text{F}}^2 + \|\mathbf{J}_{\text{X}}^{(t+1)}\|_{\text{F}}^2 + \|\mathbf{J}_{C}^{(t+1)}\|_{\text{F}}^2) = \frac{\|\boldsymbol{J}^{(t+1)}\|_{\text{F}}^{2}}{{\rho^{(\ell)}}^2}, \end{gather}\tag{51}\] where \(\boldsymbol{J}^{(t+1)}\) is defined in Lemma 3. The sufficient descent property in Lemma 3 guarantees that the left-hand side of 51 is always nonnegative. By summing 51 from \(1\) to \(T\), we have \[\begin{gather} \label{eq:subgrad95and95T} M {\rho^{(\ell)}}^2 [\mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(1)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(T)}; \boldsymbol{U}^{(\ell)})] \\ \geq \sum_{t =1}^{T} \|\boldsymbol{J}^{(t)}\|_{\text{F}}^{2} \geq T \min_{t \in [1, T]} \|\boldsymbol{J}^{(t)}\|_{\text{F}}^{2}. \end{gather}\tag{52}\] The value of ALF is lower bounded, i.e., \(\mathcal{L}_{\rho}(\boldsymbol{z}; \boldsymbol{U}) > -\infty\), since the summation of Frobenius norms in the ALF is no less than zero and the term \(\|\mathbf{U}_{k}\|_{\text{F}}^2/(2\rho^{(\ell)})\) is finite for all \(k \in \mathbf{\Omega}_{K}\) due to the boundness of multipliers. The \(\varepsilon^{(\ell)}\)-stationary point means \(\min_{t \in [1, T]} \|\boldsymbol{J}^{(t)}\|_{\text{F}} \leq \varepsilon^{(\ell)}\). Then, 52 yields that \[T \leq \frac{M u {\rho^{(\ell)}}^2}{{\varepsilon^{(\ell)}}^2}, \nonumber\] where \(u = \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(1)}; \boldsymbol{U}^{(\ell)}) - \mathcal{L}_{\rho^{(\ell)}}(\boldsymbol{z}^{(T)}; \boldsymbol{U}^{(\ell)})\) is a positive number. Therefore, the sequences generated by Algorithm 4 will achieve an \(\varepsilon^{(\ell)}\)-stationary point after \(\mathcal{O}\left(\frac{{\rho^{(\ell)}}^2}{{\varepsilon^{(\ell)}}^2}\right)\) iterations.

9 Proof of Theorem 2↩︎

The convergence proof consists of two main steps, where the first step shows the feasibility of the limit point of the sequence generated by Algorithm 1 and the second step shows that any limit point is a feasible stationary point of problem 8 . We present these two steps in two subsections separately.

9.1 Feasibility of the Limit Point↩︎

According to the updating rule of the penalty parameter in 17 , either the penalty parameter keeps fixed after a finite number of iterations or \(\lim_{\ell \rightarrow \infty} \rho^{(\ell+1)} = +\infty\). The violation of the constraints will be decreased by a factor of \(\delta < 1\) after each iteration and thus tends to zero as \(\ell \rightarrow \infty\) if it is the former case.

Now we consider the condition where \(\rho^{(\ell)} \rightarrow \infty\). Denote \((\boldsymbol{C}^{(\infty)}, \mathbf{F}^{(\infty)}, \mathbf{X}^{(\infty)})\) as any limit point of the sequence generated by Algorithm 3. Based on the choice of the initial point in 23 and the previous analysis, we can see that condition 13 holds for any \(\ell \geq 0\). That is \[\begin{gather} f^{(\ell+1)} + \sum_{k \in \mathbf{\Omega}_{K}} \Big[\frac{\rho^{(\ell)}}{2} \|\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\|_{\text{F}}^2 \\ + \mathcal{R}\Big\{tr\Big[{\mathbf{U}_{k}^{(\ell)}}^{\text{H}} \Big(\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\Big)\Big]\Big\}\Big] \leq \zeta,\nonumber \end{gather}\] where \(f^{(\ell+1)} = f(\mathbf{F}^{(\ell+1)}, \mathbf{X}^{(\ell+1)})\). Dividing both sides of the above inequality by \(\rho^{(\ell)}\) yields \[\begin{gather} \sum_{k \in \mathbf{\Omega}_{K}} \frac{1}{2} \|\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\|_{\text{F}}^2 \leq \frac{1}{\rho^{(\ell)}} \Big[\zeta - f^{(\ell+1)} \\ - \sum_{k \in \mathbf{\Omega}_{K}} \mathcal{R}\Big\{tr\Big[{\mathbf{U}_{k}^{(\ell)}}^{\text{H}} \Big(\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\Big)\Big]\Big\} \Big]. \nonumber \end{gather}\] Since \(\zeta\) is a finite constant given in 11 , \(f(\mathbf{F}, \mathbf{X})\) is lower bounded by zero, and \(\mathbf{U}_{k}\) is also bounded according to 14 , it follows that the right-hand side of the above inequality tends to zero when \(\rho^{(\ell)} \rightarrow \infty\). Therefore, \[\lim_{\ell \rightarrow \infty} \sum_{k \in \mathbf{\Omega}_{K}} \|\mathbf{X}^{(\ell+1)} \mathbf{J}_{k} {\mathbf{F}^{(\ell+1)}}^{\text{H}} - \mathbf{C}_{k}^{(\ell+1)}\|_{\text{F}}^2 \leq 0, \nonumber\] which implies \(\|\mathbf{X}^{(\infty)} \mathbf{J}_{k} {\mathbf{F}^{(\infty)}}^{\text{H}} - \mathbf{C}_{k}^{(\infty)}\|_{\text{F}} = 0\) for all \(k \in \mathbf{\Omega}_{K}\). Therefore, any limit point of the sequence generated by Algorithm 3 is feasible.

9.2 The Limit Point is a Stationary Feasible Solution of Problem 8↩︎

Based on Definition 1 and Theorem 1, we know that as \(\varepsilon^{(\ell)} \to 0\) and \(\ell \to \infty\), any limit point generated by Algorithm 3 is a stationary point. In the previous subsection, we have shown that every such limit point is also a feasible solution to problem 9 . Hence, the pair \((\mathbf{F}^{(\infty)}, \mathbf{X}^{(\infty)})\) constitutes a feasible stationary point of problem 9 . Furthermore, according to [38], any feasible stationary point of the augmented Lagrangian problem is also a feasible stationary point of the original problem. Therefore, we conclude that every limit point generated by Algorithm 3 is a feasible stationary solution of the original problem 8 .

References↩︎

---

1. This work is supported in part by the Swedish Research Council (VR) through the project 6G-PERCEF under Grant 2024-04390, in part by National Natural Science Foundation of China under Grant 62471208, in part by Guangdong Provincial Natural Science Foundation under Grant 2024A151510098, in part by ZTE Industry-University-Institute Cooperation Funds under Grant No. IA20240610003, and in part by Shenzhen Science and Technology Program under Grant JCYJ20240813094627037. The work of Y.-F. Liu was supported in part by the National Natural Science Foundation of China (NSFC) under Grant 12371314 and Grant 12021001. (Corresponding Author: Ya-Feng Liu.)↩︎

2. Kecheng Zhang and Weijie Yuan are with the School of Automation and Intelligent Manufacturing, Southern University of Science and Technology, Shenzhen 518055, China (e-mail: zhangkc2022@mail.sustech.edu.cn; yuanwj@sustech.edu.cn)↩︎

3. Ya-Feng Liu is with the Ministry of Education Key Laboratory of Mathematics and Information Networks, School of Mathematical Sciences, Beijing University of Posts and Telecommunications, Beijing 100876, China (email: yafengliu@bupt.edu.cn)↩︎

4. Musa Furkan Keskin and Henk Wymeersch are with Department of Electrical Engineering, Chalmers University of Technology, 41296 Gothenburg, Sweden (e-mail: {furkan; henkw}chalmers.se?)↩︎

5. Zhongbin Wang and Shuqiang Xia are with the ZTE Corporation and the State Key Laboratory of Mobile Network and Mobile Multimedia Technology, Shenzhen 518055, China (e-mail: {wang.zhongbin; xia.shuqiang}zte.com.cn?)↩︎

6. Minimizing \(P_{\text{MUI}}\) is not equivalent to maximizing sum rate. We use this metric to simplify the optimization problem.↩︎

7. We adopt the full-power radar receive filter to ensure that the mainlobe energy remains as high as possible while suppressing the sidelobe level. Although this can be achieved by maximizing the sidelobe-to-mainlobe ratio, it leads to a fractional programming problem, increasing the complexity of the solution. In future work, we will consider using the sidelobe-to-mainlobe ratio as an optimization objective or constraint.↩︎

8. We attempted to identify quantitative metrics for evaluating radar image quality, but to the best of our knowledge, no formal or widely accepted metric has been established. As in prior classical works [29], [41], visual inspection remains the standard evaluation approach, which we also adopt in this study.↩︎