MaRLI: Attack on the Discrete-Phase WISL Minimization Problem

To facilitate target localization, active radar signals or sequences are designed to have low auto-correlation. This goal is typically achieved by the minimization of the auto-correlation integrated sidelobe level (ISL) metric, or the more general, weighted version of the ISL, known as the WISL metric. In this paper, we introduce two problem formulations to address <italic>unimodular quartic program</italic>, which arises in numerous applications, with WISL minimization being one of them. The first approach is almost equivalent, while the second scheme utilizes a bi-quadratic formulation, both of which simplify the main problem to a quadratic program. We propose efficient approaches to WISL minimization for unimodular sequence design that take advantage of the low-cost and easily implementable <italic>power method-like iterations</italic> (PMLI). Moreover, to design the discrete-phase unimodular sequences with low WISL, the <inline-formula><tex-math notation="LaTeX">$M$</tex-math></inline-formula>-<italic>a</italic>ry unimodular sequence design via <italic>R</italic>elaxed PM<italic>LI</italic> (MaRLI) is proposed by changing and adopting the PMLI approach to tackle this challenging problem. Several numerical results and comparisons between MaRLI and other well-known sequence design algorithms are presented to illustrate the effectiveness of the proposed method. As we show, MaRLI appears to achieve a better suppression of the finite-alphabet autocorrelation levels compared to the state-of-the-art alternatives.


I. INTRODUCTION
D ESIGNING sequences/codes with good correlation prop- erties is considered to be one of the most fundamental problems in wireless communications, radar signal processing and other active sensing schemes [1], [2], [3], [4], [5], [6], [7], [8].Sequences with impulse-like periodic correlation (called perfect sequences) have found interest in pulse compression and multi-user wireless communications [5], [9].They are required in typicalcode division multiple-access (CDMA) systems to The authors are with the Department of Electrical and Computer Engineering, University of Illinois Chicago, Chicago 60607 USA (e-mail: aeamaz2@uic.edu).
Digital Object Identifier 10.1109/TSP.2023.3339398handle the multi-user interference and are also used in the synthesis of orthogonal matrices for source coding as well as complementary coding [10].Moreover, in multiple access channels (MACs), having more sequences with good correlation properties expands the capacity of the communication system [6].While sequences with good periodic and good aperiodic correlation have a considerable set of common applications, there are also some cases in which merely good aperiodic correlation properties are of interest.For instance, in synchronization applications, while sequences with good periodic correlation properties are employed when the sequence can be transmitted several times in succession, sequences with good aperiodic correlation are required when the sequence can be used only once [11].We note that finding and studying sequences with good aperiodic correlation properties is usually a harder task than those with good periodic correlation.In particular, unlike the case of periodic correlations, it is not possible to construct sequences with exact impulsive aperiodic autocorrelation [7].There have seen significant efforts dedicated to finding low autocorrelation sequences in the statistical signal processing literature.We face two main categories of existing sequence design methods which are either analytical or computational.In the realm of analytical methods, sequences are constructed based on closed-form expressions, such as the Frank sequences [12], the Chu sequences [13], and the Golomb sequences [14].Traditional computational methods are either based on search algorithms, particularly stochastic search or evolutionary algorithms.As a result, they become computationally expensive when the sequence length grows larger than a few hundreds.In addition, the convergence of evolutionary (heuristic) methods is not guaranteed [15].Consequently, these sequence construction methods cannot be effectively used to design extremely long sequences with very low autocorrelation, and that is why there is a considerable demand for modern optimization techniques in this area [1], [4], [9], [16].
Radar and communication signal optimization with unimodularity (also known as constant-modulus) constraint has attracted a significant deal of interest in the last two decades owing to the minimal peak-to-average-ratio (PAR) of the generated sequences, which facilitates a uniform temporal and spatial power allocation, as well as the fact that the signal generation can be done using low-cost amplifiers, without being exposed to the perils of gain non-linearity [2], [3], [4], [5], [17], [18], [19], [20], [21].Due to such favorable properties, we will focus our design effort on unimodular sequences.Assume Ω N is the set of complex unimodular sequences defined as (1) In addition to the unimodularity constraint, many applications require polyphase or finite-alphabet sequences to ease system implementation [5], [15], [22], [23].The feasible set of polyphase sequences Ω N M where M presents the level of quantization, is defined as Because the quantized phases are at M -levels, these sequences are called M -ary as well.The continuous unimodular set (1) is in fact a special case or a generalization of the M -ary feasible set (2) when M → ∞ [23].
The autocorrelation function of the signal s is defined as: where k ∈ {0, • • • , N − 1}, in the periodic and aperiodic cases respectively, with the signal power calculated at k = 0.The problem of sequence design for good correlation properties arises when small out-of-phase (k = 0) autocorrelation lags are desired.Due to importance of the sequence design in the aperiodic scenario in various real-world applications and its more challenging nature, we will focus on aperiodic scenario.However, the derivations for the periodic case are similar.
To design sequences with good correlation properties, several metrics have been defined to measure the goodness of such sequences, for example in the aperiodic autocorrelation scenario, the peak sidelobe level the integrated sidelobe level (ISL) the weighted integrated sidelobe level (WISL) where {ω k } are nonnegative weights, are typically considered.The WISL metric is a more general form of the ISL criterion which is particularly useful when the suppression of the autocorrelation lags are not equally important; e.g., if only some of them are to be suppressed.

A. Prior Art
A well-studied approach to the design of sequences with low (or good) correlation properties is through the minimization of the WISL criterion [2], [3], [4].The minimization of the WISL produces an optimization problem with a quartic objective in terms of the constrained transmit sequence which is deemed to be a difficult problem.Moreover, it deals with a considerable number of local optima, while many such local optima may in fact be of a good quality for deployment [17], [23].Moreover, since this problem is highly non-convex, it cannot be tackled by classical large-scale optimization techniques developed for convex problems with relatively simple constraints [24], [25].
In [4], the minimization of WISL has been accomplished by the weighted Cyclic Algorithm New (WeCAN).Although large sequence lengths have been accessible through WeCAN, the cost of design in terms of computational time can reach several hours (and even days) when the code length grows large [2], [3].In [26], the authors proposed a unimodular sequence design based on the alternating direction method of multiplier (ADMM) for WISL minimization.
To design polyphase sequences, in [27], the authors proposed an algorithm applying stochastic methods to find sequences that are good local minima for a base energy criterion.One of the well-known examples of polyphase sequences with optimal correlation properties is the Barker code, whose length in the binary case is limited to 13 [1].To design generalized polyphase Barker sequences, many Heuristic techniques such as simulated annealing (SA), threshold accepting (TA), and great deluge algorithm (GDA) have been proposed [27], [28], [29], [30], resulting in sequences lengths that are not greater than 77 [30].Note that this restriction is because of the nature of the target sequence.The iterative twisted approximation (ITROX) method was proposed in [6] to achieve discrete phase suboptimal sequence design with limited alphabet sizes.
The polyphase constraint is considered to be ill-posed, where the WISL minimization becomes more challenging with a plethora of poor local optima.To tackle the WISL minimization problem under the unimodularity and finite-alphabet constraints, majorization-minimization (MM)-based algorithms have been investigated in [2], [15].In particular, the authors of [15] provide an algorithm, based on the MM paradigm, capable of minimizing the p -norm of the autocorrelation sidelobes.Also, to minimize the combination of the PSL and WISL metrics with the same constraints, an iterative coordinate-descent (CD) framework was proposed in [23], [31].Another interesting approach to the design of finite-alphabet unimodular sequences with low WISL is the inexact alternating direction penalty method (IADPM) which was inspired by the method of multipliers, and in particular, the alternating direction penalty method (ADPM) framework proposed in [32].The proposed approach was used to tackle the quadratic optimizations under unimodularity and similarity constraint (SC) for both continuous and discrete scenarios.In each iteration, the IADPM method converts the WISL optimization problem into two subproblems while locally increasing the considered penalty factor considered within the framework [33], [34].The authors showed that, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
compared to the developed CD and MM framework of [15], the IADPM has a better performance in autocorrelation suppression but consumes more computational resources.Moreover, it was demonstrated that IADPM finds a feasible point with higher probability than ADMM.

B. Contributions of the Paper
Our main contributions in this paper are: 1) PMLI-based unimodular sequence design.This paper tackles the unimodular quartic programs (UQ 2 P) appearing in the many signal processing applications [2], [3], [19], [34].We address this problem using two distinct approaches.The first involves identifying a almost equivalent relaxed quadratic counterpart that transforms the problem into a quadratic program while introducing extra variables tied to the main variables.The second approach employs a bi-quadratic transformation, converting the initial objective into a bi-quadratic form where the main optimization variables become two parameters, allowing cyclic problem handling where each subproblem becomes a quadratic one.To explore efficient algorithms for the unimodular quartic program, we transform the well-known WISL minimization problem into a quartic program.We design sequences with good correlation properties as measured by the WISL criterion.The goal is accomplished by deploying one of the most efficient solvers currently available for unimodular quadratic programs (UQPs) and quadratic programs with polyphase variables, known as the power method-like iterations (PMLI).Inspired by the power method, the PMLI algorithm takes advantage of simple matrix vector multiplications, leading to a low computational cost algorithm [17], [35], [36].We transform the quartic formulation of the WISL metric to a quadratic alternative which can be readily solved via the PMLI.Particularly, the WISL metric will be overparametrized in such a way to enable its reformulation as a UQP.The resulting algorithm will be referred to as the Overparametrized WISL minimizer via the PMLI (OWL).Another novel algorithm, the Cyclic PMLI (CyPMLI), will be proposed which converts the WISL minimization into two UQP subproblems which are successively tackled via the PMLI.

2) PMLI-based polyphase unimodular sequence design.
To design polyphase unimodular sequences with low correlation properties, a problem that arises in a wide variety of communication and active sensing applications ranging from signal design for transmission to signal processing at receive side [37], [38], the conventional projection operator in the PMLI, namely the exponential function, will be replaced by a relaxation operator.Such an operator has helpful properties to ensure a converge to the desired discrete set-the set of polyphase unimodular sequences in this case.The effect of employing operators has been widely investigated in [6].By integrating the relaxation operator into the OWL and CyPMLI approaches, the M -ary (finite-alphabet) unimodular sequence design with the Relaxed PMLI (MaRLI) will be achieved.
3) Comparison with state-of-the-art algorithms.To evaluate the performance of the proposed algorithms, we will present several numerical results investigating the autocorrelation function and the required computational time in comparison with the existing algorithms such as We-CAN, MM and IADPM.By analysing the order complexity of proposed algorithms, we will show that OWL has a lower computational cost compared to that of the WeCAN and CyPMLI while the obtained correlation levels by the CyPMLI algorithm appear to be even smaller than WeCAN and alternatives.In the finite-alphabet scenario, we demonstrate that the MaRLI algorithm has a better performance in terms of the achieved correlation results compared to that of the MM and IADPM.Moreover, the required CPU time results for MaRLI appear to be lower than that of WeCAN and IADPM, but higher than that of the MM approach.

C. Organization and Notations
The rest of the paper will be organized as follows: In Section II, we formulate the WISL metric in its quartic formulation.The associated WISL minimization problem will thus be called the UQ 2 P. Section III is devoted to the overparametrization of the WISL metric in such a way to enable the formulation of its almost equivalent quadratic objective.The resulting UQP is efficiently tackled by taking advantage of the PMLI approach.Next, by using a bi-quadratic transformation, the WISL minimization problem will be separated to two sub-UQP problems that will be tackled by the CyPMLI algorithm in Section IV.Section V is dedicated to the MaRLI algorithm that benefits from projections on discrete signal sequences via relaxation operators, which can rather efficiently tackle the discrete-phase WISL minimization problem in the finite-alphabet unimodular signal design scenario.In Section VI, we perform a comprehensive investigation of the performance of the proposed algorithms for the task of the unimodular sequence design.Finally, Section VII concludes the paper.
Notation: Throughout this paper, we use boldface lowercase, boldface uppercase, and calligraphic letters for vectors, matrices, and sets, respectively.The notations C, R, and Z represent the set of complex, real, and integer numbers, respectively.We may represent a vector x = [x i ] in terms of its elements {x i } or (x) i .We use (•) and (•) H to denote the vector/matrix transpose and the Hermitian transpose, respectively.The identity matrix of size N is , where b rs is the (r, s)-th entry of B. The function diag(•) outputs a diagonal matrix with the input vector as its main diagonal.The p -norm of a vector b is b p p = i b p i .For any real number x, the function [x] yields the closest integer to x (the largest is chosen when this integer is not unique) and {x} = x − [x].The s-dimensional all-ones vector, all-zeros vector, and the identity matrix of size s × s are 1 s , 0 N , and I s , respectively.The minimum and maximum eigenvalues of B are denoted by λ min (B) and λ max (B), respectively.The real, imaginary, and angle/phase components of a complex number are Re(•), Im(•), Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

II. WISL METRIC: A QUARTIC FORMULATION
The autocorrelation formulations (3) can be written in their compact form as: where J P k and J AP k are the periodic and aperiodic shift matrix operators that shift the vector argument by k indices.For the sake of brevity, in this paper we are going to focus on the aperiodic case that is deemed to be more challenging, and denote the associated shift matrix as J k .Here, J k = J −k is given by [5], It is straightforward to verify that the derivations for the periodic case would be similar.
Recall the weighted integrated sidelobe level (WISL) of the autocorrelation defined in (6): with weights ω k ≥ 0 that are predetermined based on the emphasis given to minimizing different autocorrelation lags.It is easy to see that the WISL formulation subsumes the ISL metric as a special case by simply considering identical {ω k } [2].Note that f (s) can be written as Moreover, let Therefore, the WISL formulation can be rewritten as The problem of interest is to design unimodular sequences s that minimize the WISL metric: It is clear that f (s) is quartic in the unimodular sequence s.Therefore, ( 13) is a unimodular quartic programming (UQ 2 P).
In particular, (13) may be recast as maximize s H G(s)s with G(s) = λ m I − G(s), where a diagonal loading with λ m I has no effect on the solution of ( 13) due to the fact that s H G(s)s = λ m N − s H G(s)s.Such a transformation will be helpful in our later use of the PMLI approach.
In the following sections, we will propose novel algorithms to efficiently tackle the obtained UQ 2 P.

III. WISL METRIC: A QUARTIC TO QUADRATIC TRANSFORMATION
In this section, at first, we overparametrize the WISL metric in such a way to enable a quadratic reformulation.Next, we propose a cyclic algorithm to efficiently solve the resulting optimization problem with easily implementable PMLI.We will refer to the emerging algorithm as the Overparametrized WISL minimizer via the PMLI (OWL).

A. Problem Formulation
The WISL metric (10) can be equivalently written as where are Hermitian and skew-Hermitian matrices, respectively, defined as It is readily concluded that s H J (r) k s are real and imaginary values, respectively.It follows from (15) that where j = √ −1 and jJ k is a Hermitian matrix.Note that using the diagonal loading process in the form k + λI, and k + λI, where λ ∈ R is given by [17], [36], guarantees the positive-definiteness of J(r) k and J(i) k .It follows from ( 16) and ( 17) that the WISL metric can be reformulated based on the proposed diagonal loading as As discussed in Section II, f (s) is a quartic function with respect to unimodular sequence s, which is difficult to minimize.Nevertheless, similar to the proposed transition in [4] and [36], we can devise a quadratic alternative: For any matrix Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where the overparametrized WISL metric is given by where Ũ(r) k and Ũ(i) k are the Hermitian square roots of where . As a result, the WISL minimization problem boils down to the following quadratic program: By using the diagonal loading process, we can recast ( 22) as a maximization problem, i.e., where R = λ I − R, with λ being larger than the maximum eigenvalue of R. It is easy to verify that the objective of P I can be simplified as

B. OWL Algorithm
We resort to a task-specific alternating optimization (AO) or cyclic optimization algorithm [5], [35], [39], wherein we tackle P I for s, {u However, it was demonstrated in [17], [35], [36] that UQP solutions of interest can be efficiently approximated by deploying the PMLI approach.
Assume R is positive definite and s(t+1) ∞ t=0 be a sequence of unimodular codes where s(t+1) is the minimizer of the following criterion: minimize where s = s 1 .The desired vector s(t+1) of ( 25) is readily evaluated by PMLI in each iteration as Note that If s (t+1) = s (t) and R 0, which implies as Re s(t+1)H Rs (t) > s(t)H Rs (t) [17].Therefore, the UQP objective is increasing through the PMLI.Moreover, we have Due to the fact that the sequence s(t)H Rs (t) is convergent, (30) implies that s(t+1) − s(t) 2 is also converging to zero through the PMLI iterations.For fixed s, the minimization problem in (19) with respect to auxiliary variables {u (r) k }, and {u In order to see why the overparametrization in ( 19) is useful, note that where are the solutions at the t-th iteration.According to (32), the objective is a decreasing function with respect to auxiliary variables u (r)(t) k , and in each iteration.In addition, the objective of the UQP is monotonically decreasing with respect to s (t) with the PMLI algorithm which leads to conclude that the OWL algorithm can Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The steps of OWL is summarized in the Algorithm 1.
Note that the minimum eigenvalues of J share the same set of eigenvalues.The minimum eigenvalue of these matrices is given by Proof: The matrices J (r) k and J (i) k have the following structure, which makes them symmetrically sparse tridiagonal (SST) Toeplitz matrices [40]: The eigenvalues of a SST Toeplitz matrix with center diagonal a 0 and two off-diagonals a ω and a −ω where the number of zeroes between center diagonals and off-diagonals is ω − 1, are obtained as [40] where with β = mod (N, ω), N ω is the quotient of N/ω and θ( 1) In the case of , the eigenvalues are evaluated as To prove the theorem, it is solely needed to show that It follows from ( 37) that max{θ j1 } = Nωπ Nω+1 , and Nω+2 .Also, since the value of N is even in our experiments, the maximum value for N ω is N/2.Therefore, (39) is rewritten as which holds for N ≥ 2. It is worth noting that the offdiagonals {a ω , a −ω } of J (r) k and J (i) k are {0.5, 0.5}, and {j0.5, j0.5}, respectively, leading to same eigenvalues obtained from (35), which proves the theorem.
The guaranteed convergence of OWL depends on the PMLI approach and a monotonically increasing objective function with respect to s.To this end, R must be positive definite which is achieved by choosing a proper diagonal loading parameter λ .The following theorem is useful in selecting λ for our algorithm: Theorem 2: To guarantee the positive definiteness of R, the following lower bound for λ must hold: where R t and u t are defined in (21), and N is the sequence length.
Proof: According to the Schur complement, R is positive definite if and only if [41], [42]: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The first condition implies that λ ≥ λ max (R t ).To achieve a boundary for λ from the second condition, one can utilize the below relation: By choosing λ in a such way to guarantee the lower term in ( 42) is positive, the second condition of the aforementioned Schur complement is met as well.Thus, the proper value of λ will satisfy the following inequality as a result of ( 42): Considering the two bounds λ ≥ λ max (R t ), and ( 43), the lower bound in (41) will be established for λ .Note that k is a banded Toeplitz matrix whose eigenvalues are efficiently evaluated as discussed in [43].
According to Theorems 1 and 2, the parameters of the OWL algorithm, namely λ and λ may be chosen with a negligible computational cost.

IV. UNIMODULAR BI-QUADRATIC PROGRAMMING
In this section, we will use a bi-quadratic transformation to reformulate the quartic WISL metric as a bi-quadratic function.By taking advantage of the bi-quadratic transformation with respect to new variables s 1 and s 2 , whose convergence to the optimal sequence s will be investigated.

A. WISL-Minimized Sequence Design via UBQP
To find a quadratic alternative to f (s) defined in (10), we define It is also interesting to observe that if either s 1 or s 2 are fixed, minimizing g (s 1 , s 2 ) with respect to the other variable can be done via a UQP formulation [17], [35]: which is originally a unimodular bi-quadratic programming (UBQP) problem.Assume λ m be the maximum eigenvalue of G(s i ), where λ m I G(s i ).Thus, G(s i ) = λ m I − G(s i ) is positive definite.As a result, (45) can be reformulated as maximize Note that a diagonal loading with λ m I has no effect on the solution of (45) due to the fact that s In order to guarantee that a minimization of g (s 1 , s 2 ) leads to the minimization of f (s), a connection must be made to show that s 1 and s 2 obtained from (46) are convergent to the same sequence.We may consider: By adding the norm-2 error between s 1 and s 2 as a penalty with the Lagrangian multiplier to (46), we have the following regularized Lagrangian problem1 : where η is the Lagrangian multiplier.The penalty is a quadratic function as well.Consequently, the objective function of ( 48) is recast as Thus, the final UBQP formulation for ( 13) is given by where λ M is the maximum eigenvalue of G(s i ), and

B. CyPMLI
The desired vector s(t+1) j of P II is readily evaluated by PMLI in each iteration as or equivalently, It is worth pointing out that the difference of recursions for (51) and (52) compared to that of (46), i.e., s is merely in including more momentum or the effect of η at each iteration.In other words, the update process (52) requires only a simple matrix vector multiplication while leveraging information about previous update s (t) i through the momentum term.This update process resembles that of the gradient descent projection with a heavy ball momentum, where in each iteration, the solution information from the previous step is incorporated using a momentum term [45], [46].Such power method-like iterations are already shown to be convergent in terms of the signals, implying that s 1 and s 2 will be converging to each other as well.Moreover, according to (29), the error between two successive UBQP objectives is given by where Γ = λ M I − G.This difference can be further simplified as, Due to the symmetry of s with respect to s i and s j , based on (30), we can write where φ 1 = arg (s ) and φ 2 = arg (s Besides the minimum eigenvalue, the lower bound of the UBQP problem is dependent on the momentum compared to the conventional setting in (30).We call our algorithm Cyclic PMLI (CyPMLI), whose steps are summarized in Algorithm 2.
To guarantee the convergence of PMLI and a monotonically increasing objective function, G(s (t) i ) must be positive definite which is achieved by choosing a proper λ M through the following result: Theorem 3: To guarantee the positive definiteness of G(s (t) ) in each iteration, its maximum eigenvalue must meet the following condition: where η is the momentum parameter, N is the sequence length, and λ max is the maximum eigenvalue of G defined in (11).
Proof: According to the Schur complement, G(s i ) is positive definite if and only if [41], [42]: The first condition implies that λ M ≥ λ max (G), where λ max (G) is the maximum eigenvalue of G(s i ).To achieve a boundary for λ M from the second condition, one can utilize the below inequality: The objective function is g (s i , sj ) = sH j G(s i )s j .3: for t = 0 : E − 1 do 6: while until convergence do 7: 8: i ← i + 1.

16:
end while 1 .18: end for 19: return s ← s By choosing λ M in a such way to guarantee the lower term in (57) remains positive, the second condition of the aforementioned Schur complement is met as well.Thus, the optimal value of λ M may satisfy the following inequality, which is obtained from (57) considering s H s = N : It is straightforward to verify that to satisfy both λ M ≥ λ max (G), and (58), the lower bound (56) must be satisfied for λ M .Also, it is easy to verify that the second term of the max (•) in ( 56) is greater than λ max as follows: where the last inequality is trivial and completes the proof.According to Theorem 3, λ M should not only be larger than the maximum eigenvalue of G but should also satisfy the tighter inequality discussed in (56).Although the theorem proposes a condition on λ M , it in effect ensures there is enough momentum η to guarantee the solutions s 1 and s 2 are convergent.
Note that the objective of the UBQP problem P II will converge as follows in conjunction with the convergence of the optimization variables: where sj = s j 1 , and λ M is chosen based on Theorem 3. We can expand the objective of P II as follows: which may be simplified considering s i = s j = s as where the second part is the WISL metric which is being minimized here.In most practical cases, where our proposed algorithm suppresses the WISL to virtually zero or near zero values (as can be observed in the numerical examples), we can approximate the value of the objective at convergence as (60).In fact, this claim is numerically investigated in Section VI.

C. Convergence Analysis
Based on the cyclic structure inherent in the proposed algorithms and equations ( 27)-( 30), (32), and ( 53)-( 55), it becomes evident that the sequence of objective values, specifically the WISL, computed through these algorithms, consistently trends in a non-increasing manner.Additionally, it is readily apparent that these objectives possess a lower bound of 0. Consequently, the sequence of objective values is ensured to progressively converge towards a finite value.Despite substantial endeavors documented in existing literature [17], [47], the convergence rate of the PMLI algorithm remains an ongoing subject of investigation.However, the confirmation of its convergence to a stationary point has been rigorously established in [17].Therefore, we possess a fundamental level of assurance that each constituent subproblem (which constitutes a UQP) within the CyPMLI or OWL methodologies converges towards a stationary point.This fact has been both theoretically and numerically validated in [17].As a result, it follows that the ultimate subproblem that the algorithm addresses upon convergence inherently represents a stationary point of the overarching problem.It is noteworthy that the convergence analysis of cyclic algorithms continues to be an actively studied realm within the literature, and under reasonably moderate conditions, the propensity for convergence towards a stationary point can be demonstrated [48].In order to elucidate the concept of a stationary point within our context, we initially present a first-order optimality condition derived from [49].This condition pertains to the minimization of a smooth function across a diverse constraint set.A point s is said to be a stationary point of the problem (50) (with objective g(s) = g(s 1 , s 2 )), if it satisfies the first-order optimality condition as follows where T s (s) denotes the tangent cone at s. Since the sequences s (t) are bounded, we know that it has at least one limit point.Consider a limit point s (∞) , we have i.e., s (∞) is a global minimizer of the problem.It is easy to show that minimizing g s 1 , s (∞) 2 is equivalent to the following problem: Define u (s 1 ) as the objective function, we have ∇u s As mentioned in the paper, the matrix G is made positive semi definite by the loading process, consequently, we have It is straightforward to verify that ∇u s implying that s (∞) is a stationary point of the problem.Analogous outcomes can also be attained for the OWL algorithm and its problem.

V. WISL MINIMIZATION FOR POLYPHASE SEQUENCE DESIGN
Previously, we proposed two novel algorithms to design WISL-minimized unimodular sequences.In this section, in addition to the unimodularity constraint, a finite or discrete phase constraint will be considered within the WISL minimization problem as follows: This is an NP-hard problem in general whose global optimum could be a needle in the haystack of poor local optima [23], [37].In such cases, one generally needs to perform an exhaustive search to find the polyphase sequences that poses good correlation properties [6].To approximate the solution to this optimization problem, we will apply a relaxation operator to OWL and CyPMLI algorithms.The resulting algorithm is referred to as M-ary unimodular sequence design via Relaxed PMLI (MaRLI).

A. Relaxation Operators
To adapt the PMLI-based algorithms such that they can handle constrained alphabet, instead of employing exponential function e j arg(.) in each iteration, we utilize a relaxation operator defined as follows [6]: Definition 1: Consider a function h(x, t) : C × Z → C; as an extension, for every vector x let h(x, t) be an elementwisely monotonic operator iff for any x ∈ C, both |h(x, t)| and arg(h(x, t)) are monotonic in t.A set X is converging to a set X under a function h iff for every x ∈ X and for every x ∈ X , there exists an element x ∈ X such that The function h is called identity iff for any x ∈ X and x ∈ X satisfying (70), x is the closest element of X to x.The sequence of sets X (t) ∞ t=0 where Accordingly, we propose a relaxation operator for the design of M -ary unimodular sequences, corresponding to X = C − {0} and X = Ω N M .In particular, one can employ the following function which facilitates the converging of X to X : where s i are elements of s, ν 1 and ν 2 are real positive numbers.
To see how, verify that |h(s i , t)| and arg(h(s i , t)) are monotonic in t.Moreover, we can write The function h is identity since according to Definition 1, for any s i ∈ X and s i ∈ X satisfying lim t→∞ h(s i , t) = s i , s i is the closest element of X to s i , and the sequence of sets at iteration t, where φ is a projected operator of the PMLI on X .

B. Proposed Algorithm: MaRLI
The development of the MaRLI algorithm is centered around applying the relaxation operator defined in (71) at each iteration of the PMLI method, hence paving the way for the design of WISL-minimized M -ary unimodular sequences.This algorithm is inspired by [6], where relaxation operators were utilized along with the alternating projections approach.
The optimization problem (25) may be rewritten by considering the discrete phase constraint as follows: minimize The desired finite-alphabet sequence s(t+1) is therefore given at each iteration as These iterations will be referred to as the MaRLI-I algorithm.Moreover, to tackle the UBQP problem P II with the feasible set s i , s j ∈ Ω N M , the MaRLI algorithm will yield the following sequence at each iteration: are the solutions at the t-th iteration. 2: 3: Obtain R from ( 23). 6: j ← 0.

14:
end while 15: s (t+1) ← s (j) . 16: 18: end for where sj = λ M I − G(s i , and i = j ∈ {1, 2}.The MaRLI algorithm used for the UBQP problem P II is referred to as MaRLI-II algorithm. It is worth noting that both versions of the MaRLI algorithm can also be applied to the continuous phase unimodular sequence at which M → ∞.In this scenario the relaxation operator utilized in each iteration is formulated as h(s, t) = |s| e −νt e j arg(s) .The steps of MaRLI-I and MaRLI-II algorithms are summarized in Algorithms 3 and 4, respectively.We refer the interested reader to [6] for a convergence analysis of sequence design when relaxation operators are employed.

VI. NUMERICAL INVESTIGATIONS
In this section, we evaluate the performance of the proposed algorithms, i.e., CyPMLI and OWL in designing continuous unimodular sequences, and MaRLI (MaRLI-I and MaRLI-II) in designing polyphase (or finite-alphabet) unimodular sequences.In continuous unimodular sequence design, we compare CyPMLI and OWL with WeCAN [4], while for the finitealphabet sequence design, we compare MaRLI with IADPM Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
1 and s (t) 2 are the solutions at the t-th iteration.

18:
end while 19: and the MM-based approach [15], [33].All methods are initialized by randomly generated sequences.The metric for comparisons is the correlation level, defined as

A. Continuous Unimodular Sequence Design
Herein, we consider the design of a continuous unimodular sequence of length N , where N ∈ {8, 16, 32, 64, 128, 512}.Suppose that we are interested in suppressing the correlation terms {r 1 , • • • , r c }. Therefore, ω k in Eq. ( 6) will have the following form: We compare CyPMLI and OWL with WeCAN and PS algorithm proposed in [31].The correlation level of the designed sequences are shown in Fig. 1 for different values of N , whose values are averaged over 1000 experiments.The momentum parameter η is considered to be η = 2 in all experiments for CyPMLI.In this example, we afford all four algorithms the same design time for fair comparison.It appears that the obtained unimodular sequences by CyPMLI and OWL algorithms have lower correlation sidelobes at the required lags comparing to that of WeCAN and PS algorithm.Furthermore, comparing the performance of CyPMLI and OWL algorithms in Fig. 1 reveals the effectiveness of the penalty term s 1 − s 2 2 2 in CyPMLI.In other words, including more momentum or the effect of η in the updating process of CyPMLI described in (52) boosts the performance of CyPMLI compared to that of OWL.To compare the complexity cost of CyPMLI and OWL algorithms with WeCAN and PS algorithm, we designed a unimodular sequence of length N = 100.In this case, we imposed the stopping criterion such that the designed sequences by all four algorithms have the same correlation levels.The CPU time for WeCAN and PS algorithm were 23 s and 1.3 s, respectively, while CyPMLI and OWL algorithms required 0.3 s and 0.05 s, respectively, whose values are averaged over 1000 experiments.Therefore, one may conclude that the computational cost of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.To further numerically scrutinize the performance of the CyPMLI and OWL algorithms, we present the behaviour of the WISL metric (6), in Fig. 2(a) and 2(b), respectively for 1000 arbitrary unimodular sequence design experiments with length N = 100.Note that the momentum parameter η for the CyPMLI algorithm is considered to be η = 2 in all experiments.We show that both approaches can significantly suppress the WISL metric when the total number of iterations for each algorithm is fixed.It is important to mention that the number of iterations for each algorithm is set such that the WISL metric achieves the value around 1e − 07 or less which leads to 2500 and 2e + 04 for the CyPMLI and OWL algorithms, respectively.Moreover, as can be observed in Fig. 2(b), the WISL metric behaves monotonically decreasing during the iterations of the OWL algorithm.This is an interesting observation since one may not theoretically comment on the behaviour of the WISL during the iterations of the OWL algorithm (we only proved in (32) that the overparametrized WISL f r (s) decreases monotonically during the iterations of the OWL algorithm).For CyPMLI, we also present the decreasing behaviour of the regularizer s 1 − s 2 2 2 in Fig. 2(c).As can be observed, in the iterative process of the CyPMLI algorithm, the sequences s 1 and s 2 are getting closer and converge to the optimal sequence s .The behaviour of the CyPMLI objective presented in (50) during the iterations of the CyPMLI algorithm is shown in Fig. 2(d).
According to (60), the CyPMLI objective approximately converges to λ M (N + 1) = 4.0585e + 04, by proper minimization of WISL.The decreasing behaviour of f r (s) presented in (21) during the iterations of the OWL algorithm is illustrated in Fig. 2(e).This result validates the theoretical decreasing nature of f r (s) during the iterations of OWL discussed in (32).
To visually investigate the effect of tuning the momentum parameter η in CyPMLI, we consider the design of a unimodular sequence of length N = 64 with different values of the momentum parameter η ∈ {0.5, 1, 2, 4}.As shown in Fig. 3(a), momentum governs acceleration; however, choosing very large values for η places too much weight on previous information s (t) i , and the update process will be less effective.On the other Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.[33] and MM [15] in designing a polyphase sequence of length N = 100 with M = 32 and M = 2048, respectively.hand, choosing a small momentum leads to reduced incorporation of past information as well as emphasis on the convergence of the two sequences s 1 and s 2 .To further study the effect of the momentum parameter η, Fig. 3(b) illustrates the value of the WISL metric at the last iteration of the CyPMLI algorithm respect to the momentum parameter η when the total number of iterations is fixed.Note that the WISL of 1000 arbitrary experiments are plotted in Fig. 3(b).In Fig. 3(b), the CyPMLI algorithm appears to effectively minimize the WISL when the momentum parameter is η = 1.

B. Polyphase Unimodular Sequence Design
We compare the performance of MaRLI algorithms (MaRLI-I and MaRLI-II) with IADPM and MM in designing a finite-alphabet unimodular sequence of length N = 100 and M ∈ {2, 32, 2048}.Specifically, we consider ν 1 = 1.2 and ν 2 = 1e − 09 for MaRLI-I and ν 1 = 1.3, ν 2 = 1e − 04 and the momentum parameter η = 2 for MaRLI-II.As was discussed in [33], the computational complexity of IADPM is O I a N 2 + cN log 2 N , where c and I a are the maximum range cell of interest defined in (77) and the iteration number of IADPM, respectively.Also, in [15], the complexity of the MM-based algorithm is given as O (I m N log N ) where I m is the total number of iterations employed in the MMbased algorithm, which is due to the use of FFT (IFFT) for WISL metric.The computational complexity of MaRLI-II is O I b (N + 1) 2 , where I b is the total iteration number employed in the MaRLI-II due to using matrix vector multiplications.Because of the SST structure of matrices J Note that for clarity of presentation in the case of the binary sequence design (M = 2), we compare the performance of MaRLI-I, MaRLI-II and IADPM in a pair-wise manner as illustrated in Fig. 4(a)-4(c) and omit the results related to the MM algorithm since IADPM performs better in terms of the correlation level compared to that of MM.As can be observed in Fig. 4, it appears that the obtained unimodular sequences by MaRLI-I and MaRLI-II have lower correlation sidelobes at the required lags compared to that of MM and To further investigate the performance of MaRLI-I and MaRLI-II, we present the behaviour of WISL metric (6), in Fig. 5  that of the continuous case (M → ∞), thus leading to a larger minimized WISL.Interestingly, as can be observed in Fig. 5(a), it seems that the behaviour of the WISL metric during the iterations of MaRLI-I may be monotonically decreasing similar to the continuous phase scenario.For MaRLI-II, we further present the decreasing behaviour of the regularizer s 1 − s 2 2 2 in Fig. 5(c).As can be observed, in the iterative process of the MaRLI-II, the sequences s 1 and s 2 become closer until they converge to the optimal sequence s .The behaviour of the overparametrized WISL f r (s) presented in (21) during the iterations of the MaRLI-I is illustrated in Fig. 5(d).Unlike what we had in the OWL algorithm regarding the theoretical monotonically decreasing behaviour of f r (s), one may not necessarily comment on the behaviour of f r (s) during the iterations of MaRLI-I.Interestingly, however, as can be seen in Fig. 5(d), f r (s) has a monotonically decreasing behaviour in the iterations of the MaRLI-I.In Fig. 5(e), the behaviour of the objective function in (50) during the iterations of the MaRLI-II is illustrated.One can easily verify that the result in (60) also holds for the case of the polyphase sequences.Considering this fact, the MaRLI-II objective approximately converges to λ M (N + 1) = 4.0618e + 04.

VII. DISCUSSION
We embark on a two-step approach.Initially, we transformed the WISL metric into a quartic formulation, subsequently exploring diverse avenues to address it.Our focus was directed towards both continuous and discrete phase scenarios.For continuous phase situations, we introduced two distinct solutions.One is rooted in an almost equivalent quadratic problem, while the other hinges on a bi-quadratic problem formulation.Furthermore, our consideration extends to discrete phase scenarios.To navigate past the challenges posed by poor local minima in the discrete WISL minimization problem, we adopted the relaxation operators.Interestingly, these operators, initially proposed for various optimization problems, exhibit remarkable effectiveness even within the context of our specific design problem.
We introduced two original formulations for unimodular quartic programs.These novel formulations pave the way for the application of efficient and easily implementable approaches, such as the power method that solely involves matrixvector multiplication.
Our selection of the WISL metric derives from its representation as a weighted norm-2 of autocorrelation values.This quadratic nature enables us to transform the problem into a unimodular quartic programming challenge, which holds relevance across diverse applications.The primary motivation driving this paper is the proposition of innovative algorithms tailored to address the complexities inherent in such problems.These algorithms can extend beyond WISL minimization to any quartic programming context.
We made a comparison of two distinct perspectives for solving quartic programming problems.The first approach involves formulating a relaxed quadratic problem through the introduction of additional variables.The second approach entails transforming the main variable into two separate Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
variables, thereby splitting the problem into two sub quadratic problems.These subproblems can be tackled using quadratic programming solvers like the power method.Our numerical results demonstrated the effectiveness of the proposed approaches compared to that of MM and coordinate descent algorithms.

Manuscript received 5
January 2023; revised 31 August 2023 and 30 November 2023; accepted 1 December 2023.Date of publication 7 December 2023; date of current version 18 December 2023.This work was supported in part by the National Science Foundation under Grant ECCS-1809225.An earlier version of this paper was presented at the International Conference on Acoustics Speech and Signal Processing (ICASSP) 2023 [DOI: 10.1109/ICASSP49357.2023.10096466].The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Pascal Vallet.(Arian Eamaz and Farhang Yeganegi contributed equally to this work.)(Corresponding author: Arian Eamaz.) (r) k }, and {u (i) k } cyclically.The optimization problem P I is a unimodular quadratic program (UQP) with respect to the unimodular sequence s.Due to their NP-hard nature, UQPs have a highly multi-modal optimization objective.

Algorithm 1
The OWL algorithm for WISL minimization.Input: Initialization values for s (0) and u (r)(0) k and u (i)(0) k , parameters λ and λ , total number of iterations Γ. Output: Designed unimodular sequence s * , optimized parameters u (r) k , and u (i) k may be analytically obtained through the following result: Theorem 1: The matrices J (r) k and J (i) k

) Algorithm 3
The MaRLI-I algorithm for overparametrized WISL problem P I with discrete phase constraint.Input: Initialization values for s(0) and u (r)(0) k and u (i)(0) k , parameters λ and λ , parameters of the relaxation operator ν 1 , and ν 2 , total number of iterations Γ. Output: Designed M -ary unimodular sequence, optimized parameters u

Fig. 2 . 2 2
Fig. 2. Sequence design with OWL and CyPMLI.The behaviour of the WISL metric for 1000 arbitrary design experiments with sequence length N = 100 during the iterations of (a) the CyPMLI algorithm with η = 2, and (b) the OWL algorithm, when the number of iterations is considered to be fixed for each algorithm to achieve the WISL around 1e − 07 or less.As can be observed, by increasing the number of iterations, both CyPMLI and OWL algorithms properly minimize the WISL criterion.(c) The behaviour of the regularizer s 1 − s 2 2 2 in CyPMLI algorithm.As can be seen, by increasing the number of iterations, s 1 and s 2 converges to each other and to the optimal sequence s .(d) The behaviour of the CyPMLI objective presented in (50) during the iterations of the CyPMLI algorithm.Based on (60), the CyPMLI objective approximately converges to λ M (N + 1) = 4.0585e + 04.(e) The decreasing behaviour of f r (s) presented in (21) during the iterations of the OWL algorithm which has been theoretically investigated in (32).

Fig. 3 .
Fig. 3. (a) The effect of momentum parameter η in WISL minimization for the sequence design with length N = 64 using the CyPMLI algorithm when the number of iterations is fixed.(b) The plot of the WISL values in the last iteration of the CyPMLI algorithm for various η, with a sequence length of N = 64.It appears that η = 1 is leading to the lowest WISL values with CyPMLI.

k
in(15), one can easily verify that the complexity of MaRLI-I is O (((4c + 2)I c − 1) N ), where I c denotes the number of iterations for MaRLI-I.The correlation levels of the designed sequences are shown in Fig.4for different values of M , the values of which are averaged over 1000 experiments.
(a) and 5(b), respectively for 1000 arbitrary polyphase unimodular sequence design experiments with length N = 100 and M = 512.Note that for MaRLI-I, we set the parameters as ν 1 = 1.2 and ν 2 = 1e − 10, and for MaRLI-II we consider ν 1 = 1.2, ν 2 = 1e − 04, with the momentum parameter set to η = 2 for all experiments.As can be observed in Fig. 5(a) and 5(b), both algorithms can significantly suppress the WISL metric up to values near 1e − 04 when the total number of iterations for each algorithm is fixed.The minimized WISL values obtained by MaRLI-I and MaRLI-II are larger compared to that of the results reported in Fig. 2(a) and 2(b) for the continuous case.The reason behind this is the finite-alphabet restriction imposed within both MaRLI-I and MaRLI-II.In other words, the M -ary feasible set (2) is smaller compared to

Fig. 5 . 2 2
Fig. 5. Polyphase sequence design with MaRLI-I and MaRLI-II.The behaviour of the WISL metric for 1000 arbitrary polyphase sequence design experiments with sequence length N = 100 and M = 512 during the iterations of (a) the MaRLI-I algorithm with parameters ν 1 = 1.2 and ν 2 = 1e − 10, and (b) the MaRLI-II algorithm with parameters ν 1 = 1.2, ν 2 = 1e − 04 and η = 2, when the number of iterations is considered to be fixed for each algorithm to achieve the WISL around 1e − 04 or less.As can be observed, by increasing the number of iterations, both MaRLI-I and MaRLI-II properly minimize the WISL (6), (c) Decreasing behaviour of the regularizer s 1 − s 2 2 2 in MaRLI-II algorithm.As can be seen, by increasing the number of iterations, s 1 and s 2 converges to each other and to the optimal sequence s .(d) The behaviour of f r (s) presented in (21) during the iterations of the MaRLI-I described in Algorithm 3.Although one cannot theoretically comment on the decreasing nature of f r (s) in the case of the polyphase sequence design, the reported numerical results show the monotonically decreasing behaviour of f r (s).(e) The behaviour of the objective function in (50) during the iterations of the MaRLI-II.Based on (60), the MaRLI-II objective will also approximately converge to λ M (N + 1) = 4.0618e + 04.

Algorithm 2
The CyPMLI algorithm for WISL minimization. 2: IADPM methods.For all algorithms, the correlation level of the designed sequences decreases as well as the value of M grows large since a larger feasible set is acquired for each algorithm.The reason of the superior performance of MaRLI-II compared to that of MaRLI-I in terms of the correlation level is the same as discussed in the continuous design case.The corresponding CPU times are summarized in TableI.The MaRLI-I and MaRLI-II can achieve lower correlation levels in much less CPU time compared to that of IADPM.The CPU time of MM algorithm appears to be less than that of the MaRLI-I and MaRLI-II.However, increasing the number of iterations (CPU time) in MM algorithm is not effective in the design procedure since it cannot further improve the objective function (WISL), and subsequently the correlation levels.As can be observed in TableI, by increasing the level of phase quantization M , the CPU time required by the MaRLI-I and MaRLI-II algorithms are reduced owing to the fact that the M -ary feasible set (2) is expanding and thus, proposed algorithms may more easily reach good local optima.