Global and Local Optimization Algorithms for Optimal Signal Set Design

The problem of choosing an optimal signal set for non-Gaussian detection was reduced to a smooth inequality constrained mini-max nonlinear programming problem by Gockenbach and Kearsley. Here we consider the application of several optimization algorithms, both global and local, to this problem. The most promising results are obtained when special-purpose sequential quadratic programming (SQP) algorithms are embedded into stochastic global algorithms.


Introduction
Consider the simple communication system model shown in Fig. 1. The goal is to transmit one of M possible symbols, i.e., an M -ary signaling system, over a memoryless additive noise channel. We will assume all signals are discrete-time with T samples. The transmitter assigns a unique signal s m : {1, . . . , T } → ‫ޒ‬ to each symbol m ʦ {1, . . . , M }. It is this signal that is sent through the channel. At the other end, the received signal is where n : {1, . . . , T } → ‫ޒ‬ is a noise process, and the job of the receiver is to decide which symbol was transmitted. Our goal is to design a set of signals s m , m = 1, . . . , M , which maximize, subject to constraints on the signals, the probability of a correct decision by the receiver given a particular channel noise distribution.
Of course, in order to design an optimal signal set, the action of the channel and the receiver must be completely specified. For the channel, we assume the noise process is independent and identically distributed (iid) with distribution p N . Further, we assume that the noise process is independent of the symbol being transmitted. Our detection problem falls into the class of M -ary Bayesian hypothesis testing problems where, for m = 1, . . . , M , the hypotheses are defined as follows, Finally, it is assumed that the receiver was designed using the minimum average probability of error criterion (or the uniform cost criterion). It is well known that is known since, given a signal set, this distribution is completely determined by the distribution on the channel noise. Specifically, in view of our assumptions, If the prior distribution were known, the quantity to be maximized could be expanded as As p (H m ) is not assumed to be known, the worst-case prior distribution will be used to compute p ({correct decision}) for any given signal set. In particular, let The goal will be to find signal sets which maximize It is not difficult to show that this is equivalent to maximizing min mʦ{1,...,M} p ({correct decision} | H m ). (1) A standard assumption in transmitter design is that the signals are restricted to be of the form where k : {1, . . . , T } → ‫,ޒ‬ k = 1, . . . , K , are given basis functions and ␣ m,k ʦ ‫,ޒ‬ m = 1, . . . , M , k = 1, . . . , K , are the free parameters. Finally, due to power limitations in the transmitter, the signals are forced to satisfy some type of power constraint, either peak amplitude or average energy. In this paper, we will assume a peak amplitude constraint, i.e., where C > 0 is given. Note that we could just as easily have considered an average energy constraint in our formulation. Our design problem is thus reduced to choosing parameters ␣ m,k in order to maximize Eq. (1), subject to the constraints Eq. (3).

The Optimization Problem
The details of casting the design problem introduced in the previous section in such a way that it may be solved using standard nonlinear programming algorithms are presented in this section. The design of optimal signal sets under the assumption of Gaussian noise has been well studied (see, e.g., Ref. [22]). In fact, a gradient-based first-order algorithm was developed and analyzed in Ref. [5] for the case of Gaussian noise, K = 2 basis functions, and an average energy constraint on the signals. The performance of optimal detectors in the presence of non-Gaussian noise as a function of the signal set was first studied by Johnson and Orsak in Ref. [10]. It was shown in Ref. [10] that the dependence of detector performance on the signal set is related to the Kullback-Leibler (KL) distance between distributions for the various hypotheses. Based on this work, Gockenbach and Kearsley Ref. [7] proposed the nonlinear programming (NLP) formulation of the signal set design problem which is considered here.
Given our assumptions on the noise process, the loglikelihood ratio may be written Note that, since randomness only enters the received signal through the additive noise process, when hypothesis H m is true, the receiver computes and, for m' m , Thus, upon substitution, the statistic of interest to us is i.e., the KL distance between the noise distribution and the noise distribution shifted by Ϫ␦ . Note that if we assume a symmetric distribution for the noise (this is not a restrictive assumption), then K N (и) will be an even function. It is not difficult to show that Substituting the expansion Eq. (2), we see that, under our assumptions, the signal set design problem is equivalent to solving the optimization problem It is only necessary to consider m' > m since K N (и) is an even function.
The problem (SS ) is already in a form which may be handled by algorithms which tackle inequality constrained mini-max problems. Such algorithms have been developed by, e.g., Kiwiel in Ref. [13], Panier and Tits in Ref. [16], and Zhou in Ref. [25], all in the context of feasible iterates. In Ref. [25], Zhou extends the nonmonotone line search-based algorithm of Ref. [26] to handle the constrained case. The algorithm of Ref. [16] extends the feasible SQP algorithm of Ref. [17] to handle mini-max objective functions. A recent algorithm for the constrained mini-max problem which does not generate feasible iterates is the augmented Lagrangian approach of Rustem and Nguyen Ref. [23]. The problem may be converted into an equivalent single objective smooth nonlinear programming problem by adding an additional variable. This is the approach taken by Gockenbach and Kearsley in Ref. [7]. Additionally, they add a "regularization" term (possibly zero) to the converted objective. Specifically, consider where ⑀ r is small (possibly zero). It turns out that problems (SS ) and (SS' ) are difficult to solve using standard off-the-shelf nonlinear programming codes. To begin with, it was observed in Ref. [7] that outside of the feasible region, the linearized constraints for problem (SS ) are often inconsistent, i.e. no feasible solution exists. This can be a big problem for sequential quadratic programming (SQP) based algorithms, generally considered the most successful algorithms available for NLP problems with a reasonable number of variables. Of course, with feasible iterates, the linearized constraints are always consistent and the solutions of the QP sub-problems are always well-defined. As an alternative to feasible iterates, there are infeasible SQP-based algorithms which use special techniques to deal with inconsistent QPs (see, e.g., Ref. [24,12,4]). Another difficulty in applying a local NLP algorithm is that problem (SS ) has many local solutions which may prevent convergence to a global solution.

Local Algorithms
Sequential Quadratic Programming (SQP) has evolved into a broad classification encompassing a variety of algorithms. When the number of variables is not too large, SQP algorithms are widely acknowledged to be the most successful algorithms available for solving smooth nonlinear programming problems. For an excellent recent survey of SQP algorithms, and the theory behind them, see Ref. [4].
In general, an SQP algorithm is characterized as one in which a quadratic model of the problem is formed at the current estimate of the solution and is solved in order to construct the next estimate of the solution. Typically, in order to ensure global convergence, a suitable merit function is used to perform a line search in the direction provided by the solution of the quadratic model. While such algorithms are potentially very fast, the local rate of convergence is critically dependent upon the type of second order information utilized in the quadratic model as well as the method by which this information is updated.

Infeasible SQP
Numerous SQP algorithms that do not require iterates to remain feasible have been suggested by researchers (e.g., Refs. [3,6,24] among others). Because of the nonlinear nature of the constaints appearing in this specific class of problems, the linearizations employed by SQP are frequently inconsistent. As a result, constraint perturbation techniques must be employed to ensure that the algorithm can generate a descent direction at every iteration. This, at least in part, explains the superior performance of hybrid SQP algorithms reported on here (see Tables 2 and 3). These algorithms are particularly desirable because they do not require a feasible starting point.

Feasible SQP
In Ref. [19,17,18,9], variations on the standard SQP iteration are proposed which generate iterates lying within the feasible set. Such methods are sometimes referred to as "Feasible SQP" (or FSQP) algorithms. It was observed that requiring feasible iterates has both algorithmic and application-oriented advantages. Algorithmically, feasible iterates are desirable because • The Quadratic Programming (QP) subproblems are always consistent, i.e., a feasible solution always exists, and • The objective function may be used directly as a merit function in the line search.
State of the art SQP algorithms typically include complex schemes to deal with inconsistent QPs. Further, the choice of an appropriate merit function (to enforce global convergence) is not always clear. Requiring feasible iterates eliminates these issues. In order to overcome the problem of inconsistent QPs in this work, we use the CFSQP implementation Ref. [14] of the FSQP algorithm proposed and analyzed in Ref. [18] as our local algorithm.

Stochastic Approach
In one attempt to overcome the problem of local solutions, we use a stochastic two-phase method (see, e.g., Ref. [1]) where random initial points are generated in the feasible region and a local method is repeatedly applied to a subset of these points. Such an approach may be thought of as simply a "smart" way of generating many initial points for our fast local algorithm with the hopes of eventually identifying a global solution. Specifically, we will use the Multi-Level Single Linkage (MLSL) approach Ref. [1,11], which is known to find, with probability one, all local minima (hence the global minima) in a finite number of iterations.
Let ᏹ denote the cumulative set of local minimizers identified by the MLSL algorithm. At iteration ᐉ , for some integer N > 0 fixed, we generate N points ␣ (ᐉϪ1)N+1 , . . . , ␣ ᐉN distributed uniformly over the feasible set X for (SS ). For each of the points ␣ i ʦ {␣ 1 , . . . , ␣ ᐉN } we check to see if there is another point ␣ j within a "critical distance" r ᐉ of ␣ i which also has a smaller objective value. If not, then the local algorithm, call it ᏸ is applied with initial point ␣ i and the computed local minimizer is added to the set ᏹ. After all points are checked, r ᐉ is updated, ᐉ is set to ᐉ + 1 and the process is repeated. At any given iteration, the local maximizer with the smallest objective value is our current estimate of the global solution. For ease of notation, define the (mini-max) objective Further, let ᏸ(␣ ) denote the local minimizer obtained when a local algorithm is applied to problem (SS ) with initial point ␣ . The following algorithm statement is adapted from Ref. [1].
It remains to specify how we select the critical distance r ᐉ , the definition of the metric ||и|| s we use for signal sets (as parameterized by ␣ ), and how we generate the sample points. Following Ref. [1], we use where n is the number of variables (MK for our problem), m (X ) is the volume of the feasible region, and > 2. To compute m (X ), note that in view of symmetry with respect to the signals, where A is the volume of the feasible region for the parameters corresponding to one signal (recall, M is the number of signals). The quantity A is easily estimated using a Monte Carlo technique.
Note that, for our problem, as far as the transmitter is concerned, a given signal set is unchanged if we were to swap the coefficients ␣ m 1 ,k , k = 1, . . . , K , with ␣ m 2 ,1 , k = 1, . . . , K , where m 1 m 2 . The distance "metric" we use in Algorithm MLSL should take this symmetry into account. Consider the following procedure for computing the distance between signal sets parameterized by ␣ 1 and ␣ 2 .
This is not a metric in the strict sense of the definition, though it suffices for our purposes and we will set To aid the generation of sample points, before starting the MLSL loop we compute the smallest box which contains the feasible set X . By symmetry with respect to the signals, we can do this by solving 2K linear programs. Specifically, let ␣ k ʦ ‫,ޒ‬ k = 1, . . . , K be defined as the optimal value of the linear program (LP) and let ␣ k ʦ ‫,ޒ‬ k = 1, . . . , K be defined as the optimal value of the LP Then, it should be clear that Using standard random number generators, it is a simple matter to choose samples from the uniform distribution on the box ᑜ. Thus, for Step 1 of Algorithm MLSL, we repeatedly generate samples from the uniform distribution on ᑜ, discarding those which do not lie in X , until we find N which do lie in X . It should be clear that such a procedure is equivalent to drawing N samples from the uniform distribution on X .

Expanded Space Approach
In this section we describe a reformulation of the problem along the lines of that considered in Ref. [8] in the context of molecular conformation. The essence of the approach is to add variables in such a way that the global solution of the new equivalent problem has an enlarged basin of attraction. While the method does not guarantee convergence to a global solution, it does increase the likelihood.
To use such an approach in our context, we assign to each signal M Ϫ 1 sets of "auxiliary" coefficients. Each set will be used exclusively for computing the interaction with one particular other signal. For a given signal, the method will asymptotically force these sets of auxiliary coefficients to coalesce, as they must in order for the solution to be feasible for the original problem. Thus, the problem has grown from KM to KM (M Ϫ 1) variables. Let ␣ denote the vector of all such coefficients. As mentioned above, the coefficients ␣ m,ᐉ,k and ␣ ᐉ,m,k , k = 1, . . . , K will be used only to compute the interaction between signals m and ᐉ . Correspondingly, define the objective function which encapsulates this interaction as A schematic representation of the auxiliary coefficients and their interactions is given in Fig. 2 for the case M = 3.
The constraint that the sets of auxiliary coefficients for each signal must all be equal at a feasible solution is easily expressed as a set of linear equality constraints. In particular, it is not difficult to show that there exists a linear operator W , a KM (M Ϫ 2) ϫ KM (M Ϫ 1) matrix, whose kernel is precisely the set of all vectors satisfying this constraint. Finally, in view of this constraint, it is only necessary to enforce the power constraint on one set of auxiliary coefficients for each signal. Thus, our "equivalent" expanded space problem is Following Ref. [8], instead of attempting to solve this problem directly, we incorporate the equality constraints into the objective functions as a quadratic penalty term. This allows us to approach solutions from "infeasible" points by carefully increasing the penalty parameter and asymptotically forcing the auxiliary coefficients to coalesce for each signal. Specifically, define the penalized objective For a fixed value of > 0, we can solve (using, e.g., a mini-max SQP-type algorithm) the problem Let ᑦ (␣ 0 , ) denote the solution of PES ( ) obtained using a local algorithm (such as those discussed in Sec. 3) starting with the initial point ␣ 0 . Using the solution just obtained as the next initial point, the process may be repeated after appropriately increasing the penalty parameter . At any time, a candidate global solution for the original problem (SS ) may be computed by "projecting" a solution ␣ = ᑦ (␣ ' , ) of PES ( ) into ‫ޒ‬ MK and solving (SS ) using a local algorithm with the projected vector as the initial point. We will denote the projection operation Proj(␣ ) and, using the notation from Sec. 4.1, the computed local solution for (SS ) starting from the initial point Proj(␣ ) will be denoted ᑦ(Proj(␣ )). One possible method to compute the required projection is to simply take the arithmetic average of the corresponding components of the auxiliary coefficients, i.e., if ␣ = Proj(␣ ), then componentwise It remains to specify how we update the penalty parameter . At "major" iteration 1 i , the penalty parameter will be increased by a multiplicative factor, call it ␦ i , i.e., In addition, we will decrease the factor ␦ i when the projection for the current major iteration produces an estimate which does not improve upon the previous estimate. If i is increased too fast for a given problem we could get trapped in a local solution. The precise algorithm statement follows.

Homotopy Approach
Suppose that for a given set of basis functions (and fixed problem size) we know the globally optimal signal set for one particular noise distribution, say p 1 (и). For example, in many cases it is possible to analytically compute the globally optimal signal set for Gaussian noise. In this section, we discuss a method for computing candidate globally optimal signal sets for a different noise distribution, say p 2 (и), based on this knowledge.
The so-called homotopy approach relies upon the observation that, for ʦ [0,1 ], p N ( ; ) = (1 Ϫ ) p 1 ( ) + p 2 ( ) is a valid noise distribution. Gradually increasing and repeatedly applying a local algorithm, the computed optimal signal set should trace out a continuous "path" from that for p 1 to that for p 2 . Let ␣ ( , ␣ 0 ) denote the computed optimal solution of (SS ) obtained via a local algorithm starting with the initial point ␣ 0 and using the noise distribution p N (и ; ). At iteration i , we compute and the homotopy parameter i is updated. The proper updating of i is critical. Clearly, we should have Further, the rate of convergence of the method is directly related to how quickly i is increased. On the other hand, if i is increased too quickly then ␣ i+1 may "jump" to a new path of minimizers. For preliminary tests we simply increment i by a small, fixed, amount.
Numerically, one of the biggest challenges associated with this approach is the computation of the KL distance K N and the derivative of the KL distance K' N . For the distributions considered in this work these functions are readily obtained analytically. Unfortunately, this is no longer possible for convex combinations of these distributions as considered in this section. Consequently, we are forced to turn to numerical quadrature. For the preliminary numerical implementations we use a simple infinite trapezoid rule, i.e. we use the approximation For those integrands with a slowly decaying tail we use the change of variables (t ) = e t Ϫ e Ϫt .

Numerical Results
Following Ref. [7], we consider the noise distributions p N listed in Table 1  In order to demonstrate the need for special-purpose SQP codes, we attempted to solve the problem using VF02AD 2 from the Harwell subroutine library Ref. [15], a standard SQP code based on Powell's algorithm Ref. [21] and a hybrid SQP code recently developed Ref. [3] and analyzed Ref. [2]. These codes do not directly solve mini-max problems, we used the formulation (SS' ) suggested in Ref. [7]. In Tables 2 and 3, we list the number of times VF02AD and the BKT SQP algorithm Ref. [3] successfully converged to a local minimizer out of 20 trials for a given noise distribution and basis (and regularization parameter). In parenthesis we report on the number of times each algorithm succeeded in converging to the global minimizer. For each trial the initial point was drawn from the uniform distribution over the feasible set.
It is clear from the table that the standard SQP algorithm had little success converging to a local solution. The failures were essentially always due to inconsistent constraints in the QP sub-problem. In our trials, neither the Feasible SQP algorithm nor the hybrid BKT SQP algorithm failed to converge to, at least, a local solution.
We ran Algorithm MLSL for 20 different problem instances. The algorithm was stopped after it appeared that no better local minimizers would be found (i.e., the estimate of the global minimum remained constant for 2 Certain commercial equipment, instruments, or materials are identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose.  several MLSL iterations). In Tables 4 and 5 we list our computed minimum values for instances of (SS ) with M = 8 and M = 16, respectively. Note that our solutions respectively. Note that our solutions agree with those reported in Ref. [7]. In all cases, execution was terminated after no more than 10 to 15 minutes. In Figs. 3 through 8 we show the optimal signal constellations for several of the instances of (SS ) corresponding to the optimal values listed in Tables 4 and 5.

Conclusions
In this paper we have presented an important and difficult optimization problem along with a broad arsenal of numerical optimization algorithms and modern enhancements that can be used to solve it. These problems are not "large-scale" by modern computing standards but they are, nonetheless, extremely difficult problems to solve efficiently.
Numerical solutions to these problems were located using SQP methods embedded into stochastic global algorithms. Numerous numerical tests suggest that this embedding procedure can significantly improve the performance of the SQP method.
Because there are so many different algorithms and implementations for the solution of nonlinear programming problems there is a need to create an accepted collection of test problems arising from the application of optimization. Because of the difficulties it poses, this family of problems is a natural candidate to use as a practical and significant test problem.