RANK-DEFICIENT QUADRATIC-FORM MAXIMIZATION OVER M-PHASE ALPHABET: POLYNOMIAL-COMPLEXITY SOLVABILITY AND ALGORITHMIC DEVELOPMENTS∗ Anastasios T. Kyrillidis
George N. Karystinos
School of Computer & Communication Sciences EPFL CH-1015 Lausanne, Switzerland E-mail: anastasios.kyrillidis@epfl.ch
Electronic & Computer Engineering Department Technical University of Crete Chania, 73100, Greece E-mail:
[email protected]
ABSTRACT
In the present work, we modify the algorithm in [6] to serve complex-domain optimization problems. Specifically, The maximization of a positive (semi)definite complex quadratic we present an algorithm for the efficient computation of form over a finite alphabet is N P-hard and achieved through the M -phase vector of length N that maximizes a rank-D exhaustive search when the form has full rank. However, if quadratic form with polynomial complexity for any M , D, the form is rank-deficient, the optimal solution can be comand N if D is independent of N . The algorithm uses 2D − 1 puted with only polynomial complexity in the length N of auxiliary hyperspherical coordinates that partition the multhe maximizing vector. In this work, we consider the general tidimensional hypercube into distinct regions of polynomial case of a rank-D positive (semi)definite complex quadratic size, each of which corresponds to a different candidate vecform and develop a method that maximizes the form with retor. Therefore, the method reduces the size of the candidate spect to a M -phase vector with polynomial complexity. The vector set from exponential to polynomial and the proposed proposed method efficiently reduces the size of the feasible algorithm turns out to be time and memory efficient, fully set from exponential to polynomial. We also develop an alparallelizable and rank-scalable. gorithm that constructs the polynomial-size candidate set in polynomial time and observe that it is fully parallelizable and 2. PROBLEM STATEMENT rank-scalable. We consider the computation of the M -phase vector that maximizes the quadratic form
1. INTRODUCTION Unconstrained complex quadratic maximization over a finite alphabet captures many problems that are of interest to the communications and signal processing community. Interestingly, it has been recently proven that the maximization of a quadratic form with a binary vector argument is no longer N P-hard if the rank of the form is not a function of the problem size [1] and can be computed optimally and efficiently in polynomial time through a variety of computational geometry (CG) algorithms, such as the incremental algorithm for cell enumeration in arrangements [2] and the reverse search method [3],[4]. However, it should be noted that, although the incremental algorithm is applicable to M -phase vectorargument alphabets with M > 4, it is not known whether it has practical importance in higher dimensions due to lack of parallelizability and memory management inefficiency. On the other hand, the reverse search method is highly parallelizable and speed/memory efficient but has been applied, until today, only for M = 2, 4. From a communicationssystems-optimization perspective, developments have led to algorithms that optimally solve the problem of maximization of a rank-D quadratic form over the M -phase alphabet for D = 1 for any M [5]. ∗ This
work was supported in part by Alexandros S. Onassis Foundation.
978-1-4577-0539-7/11/$26.00 ©2011 IEEE
3856
sopt arg max sH VVH s = arg max VH s s∈AN M
s∈AN M
(1)
where V ∈ CN ×D is a rank-D complex matrix, s ∈ AN is j2πm M a M -phase N -tuple vector argument, AM = e M m = 0, 1, . . . , M − 1 is the M -phase alphabet, and M ∈ {2k | k = 1, 2 . . . }. In the next section, we use the framework presented in [6] and propose a more generalized algorithm for the maximization of a rank-deficient quadratic form over any M -phase alphabet where M ∈ {2k | k = 1, 2 . . . }. 3. EFFICIENT RANK-DEFICIENT QUADRATIC-FORM MAXIMIZATION WITH A M -PHASE VECTOR ARGUMENT 3.1. Problem Reformulation W.l.o.g., we assume that each row of V has at least one nonzero element, i.e. Vn,1:D = 01×D , ∀n ∈ {1, 2, . . . , N }. Let φi:j [φi , φi+1 , . . . , φj ]T . To develop an efficient method for the maximization in (1), we introduce 2D −1 auxiliary hyperspherical coordinates φ1:2D−1 ∈ (− π2 , π2 ]2D−2 ×
ICASSP 2011
(−π, π] and define the hyperspherical real vector with unit radial coordinate ⎤ ⎡ sin φ1 ⎥ ⎢ cos φ1 sin φ2 ⎥ ⎢ ⎥ ⎢ cos φ1 cos φ2 sin φ3 ⎥ ⎢ ⎥ ⎢ . .. (2) c˜(φ1:2D−1 ) ⎢ ⎥ ⎥ ⎢
⎥ ⎢ 2D−2 cos φi sin φ2D−1 ⎥ ⎢ ⎦ ⎣ i=1
2D−2 cos φ cos φ i 2D−1 i=1
relocates the φ2D−1 of the hyperspheri angular coordinate π π jα in the interval (− M cal vector c(φ1:2D−1 )e ,M ] and results in the same value in (5). Thus, without loss of optimality, we choose α ∈ arg{AM } π π ,M ]. Then, (5) becomes such that φ2D−1 ∈ (− M sopt = arg max max sH Vc(φ1:2D−1 ) . (6) φ1:2D−1 ∈ s∈AN M π 2D−2 π π (− π ×(− M ,M ] 2,2]
2D×1
as well as the D×1 hyperspherical complex vector c(φ1:2D−1 ) c˜2:2:2D (φ1:2D−1 ) + j˜c1:2:2D−1 (φ1:2D−1 ).1 From Cauchy-Schwarz Inequality, we observe that for any a ∈ CD , aH c(φ1:2D−1 ) ≤ ac(φ1:2D−1 ) = a with equality if and only if φ1:2D−1 ∈ (− π2 , π2 ]2D−2 ×(−π, π] 2 are the hyperspherical coordinates of vector a, i.e. if a a c(φ1:2D−1 ) = a since aH c(φ1:2D−1 ) = aH a = a. Using the above, a critical equality for our subsequent developments is H sopt = arg max max s Vc(φ1:2D−1 ). (3) s∈AN M
φ1:2D−1 ∈ π 2D−2 (− π , ×(−π,π] 2 2]
with equality if and only if θˆ = θ arg aH c(φ1:2D−1 ) . We observe that Cauchy-Schwarz Inequality and (4) are simultaneously satisfied with equality if and only if φ1:2D−1 ∈ (− π2 , π2 ]2D−2 × (−π, π] are the hyperspherical coordinates of vector a and θˆ = θ. Then, applying some computations, (3) can be further transformed into: max sH Vc(φ1:2D−1 ) . (5) sopt = arg max φ1:2D−1 ∈ π 2D−2 (− π ×(−π,π] 2,2]
We note that given a hyperspherical complex vector c(φ1:2D−1 ) and φ2D−1 ∈ (−π, π], there always exists an 2πm angle α ∈ arg{AM } = m = 0, 1, . . . , M − 1 that
max
max s∗n Vn,1:D c(φ1:2D−1 ) . (7)
φ1:2D−1 ∈ sn ∈AM π π n=1 Φ2D−2 ×(− M ,M ]
We observe that the original problem in (1) is decomposed into symbol-by-symbol maximizations for a given φ1:2D−1 ∈ π π , M ]. For such a set of angles, the maximizaΦ2D−2 × (− M tion argument of the sum in (7), e.g. symbol sn , depends only on the corresponding row of matrix V. As φ1:2D−1 vary, the decision in favor of sn is maintained as long as a decision boundary is not crossed. Due to the structure of AM and given the definitions above, the M 2 decision boundaries for the determination of 2k+1 sn are given by Vn,1:D c(φ1:2D−1 ) = Aejπ M , A ∈ R, k = 0, 1, . . . , M 2 − 1, or, equivalently, 2k+1 M − 1. e−jπ M Vn,1:D c(φ1:2D−1 ) = 0, k = 0, 1, . . . , 2 (8) For n = 1, . . . , N and k = 0, 1, . . . , M 2 − 1, (8) becomes ˜ l,1:2D c˜(φ1:2D−1 ) = 0, l = 1, 2, . . . , M N (9) V 2 ˜ :,1:2:2D−1 = V ˆ and V ˜ :,2:2:2D = V ˆ with where V (M −1)π π 3π ˆ = V ⊗ [ej M ej M · · · ej M ]T and ⊗ denotes the V Kronecker product. Motivated by the statements above and the inner maximization rule in (7), for each D × 1 complex vector v we define the decision function s that maps φ1:2D−1 to AM according to s(vT ; φ1:2D−1 ) arg max {s∗ vT c(φ1:2D−1 )}.
M
1 The definition of c ˜(φ1:2D−1 ) and c(φ1:2D−1 ) differentiates the developments of this present work from the work in [7] where the poynomial solvability of (1) was treated in a similar way with the current work but led to an algorithm with higher complexity. In this work, the new definition of ˜ c(φ1:2D−1 ) and c(φ1:2D−1 ) helps us minimize the number of candidate vectors and overall complexity. 2 We observe that equality is also achieved for any rotated version of a for any ω ∈ (−π, π] since c(φ1:2D−1 ), i.e. ejω c(φ1:2D−1 ) = ejω a H jω a jω H jω a e c(φ1:2D−1 ) = a e a = e a = a. But, for clarity reasons and w.l.o.g., we present the case for ω = 0 in the above statement.
N
3.2. Decision Functions and Candidate Vector Set S(VN ×D )
Furthermore, we observe that for any a ∈ CD and any θˆ ∈ (−π, π], ˆ aH c(φ1:2D−1 )e−j θ ≤ aH c(φ1:2D−1 ) (4)
s∈AN M
We drop the arg operator, interchange the maximizations in (6), and obtain the equivalent problem (for Φ (− π2 , π2 ])
s∈AM
(10)
Furthermore, for the given N ×D complex observation matrix V, we can construct the vector decision function s using (10) π π where each point φ1:2D−1 ∈ Φ2D−2 × (− M ,M ] is mapped to a candidate M -phase vector according to ⎡ ⎤ s(V1,1:D ; φ1:2D−1 ) ⎥ ⎢ .. (11) s(VN ×D ; φ1:2D−1 ) ⎣ ⎦. . s(VN,1:D ; φ1:2D−1 )
3857
Computing s(VN ×D ; φ1:2D−1 ) for ∀φ1:2D−1 ∈ Φ2D−2 × π π (− M ,M ], we collect all M -phase candidate vectors into set N S(VN ×D ) φ1:2D−1 ∈ s(VN ×D ; φ1:2D−1 ) ⊆ AM . π π Φ2D−2 ×(− M ,M ]
π π Since φ1:2D−1 takes values from the set Φ2D−2 × (−M ,M ], H our problem in (1) becomes sopt = arg maxs∈S(V) V s , i.e. the M -phase candidate vector sopt that maximizes the expression above belongs into the set S(VN ×D ). In D the following, we (i) show that |S(VN ×D )| = d=1 i d−1 N N −i M 2(d−i)−2 M − 1 and (ii) dei=0 i 2 2 2(d−i)−1 velop an algorithm for the construction of S(VN ×D ) with 2D complexity O(( MN ). 2 )
3.3. Hypersurfaces and Cardinality of S(VN ×D ) ˜ MN determine MN According to (9), the rows of V 2 hy2 ×2D ˜ 1,1:2D ), . . . , H(V ˜ MN persurfaces H {H(V )} that 2 ,1:2D 2D−2 partition the (2D − 1)-dimensional hypercube Φ × π π ,M ] into K noninterleaving cells C1 , . . . , CK such that (− M π π the union of all cells equals Φ2D−2 × (− M ,M ] and the intersection of any two distinct cells, say Ck , Cj for k = j, is empty. Each cell Ck corresponds to a distinct sk ∈ AN M in the sense that s(VN ×D ; φ1:2D−1 ) = sk for any φ1:2D−1 ∈ Ck and sk = sj if k = j, k, j ∈ {1, . . . , K}. Let I2D−1 {i1 , . . . , i2D−1 } ⊂ {1, . . . , MN 2 } denote the subset of 2D − 1 indices that correspond to hypersurfaces ˜ i2D−1 ,1:2D ). We identify the follow˜ i1 ,1:2D ), . . . , H(V H(V ing cases. (a) Intersections of 2D − 1 hypersurfaces where at most two surfaces originate from the same row of V. (b) Intersections of 2D − 1 hypersurfaces where at least three surfaces originate from the same row of V. Two basic properties of such intersections follow. Proposition 1 (i) Each subset of H that consists of 2D − 1 hypersurfaces has either a single or uncountably many interπ π sections in Φ2D−2 × (− M ,M ]. (ii) Each combination of 2D − 1 hypersurfaces from the set H has a unique intersection point that constitutes a vertex of a cell if and only if no more than two hypersurfaces originate from the same row of the matrix V. According to Proposition 1, Part (ii), combinations of the form (b) do not have a unique intersection point but infinitely many intersection points; thus no cell is created and these combinations can be ignored. On the other hand, combinations of the form (a) have a unique intersection π π ,M ] that leads point φ(VN ×D ; I2D−1 ) ∈ Φ2D−2 × (− M Q cells, say C1 (VN ×D ; I 2D−1 ), C2 (VN ×D ; I2D−1 ), . . . , M 0 M 1 CQ (VN ×D ; I2D−1 ), Q ∈ ( M 2 − 1) , ( 2 − 1) , . . . , ( 2 − 1)D−1 and each cell is associated with a distinct M -phase candidate vector sq (VN ×D ; I2D−1 ), q = 1, 2, . . . , Q, in the sense that sq (VN ×D ; φ2D−1 ) = sq (VN ×D ; I2D−1 ) for all φ1:2D−1 ∈ Cq (VN ×D ; I2D−1 ) and φ(VN ×D ; I2D−1 ) is
a single point of Cq (VN ×D ; I2D−1 ) where φ2D−1 is minimized. The number of cells Q “led” by an intersection point depends on the number p of participating pairs of hypersurfaces that originate from the same row of matrix V and equals p to ( M 2 − 1) . Since each cell is associated with a distinct M -phase candidate vector, we can collect all these vectorsinto J (VN ×D ) I2D−1 ⊂{1,2,..., M N }, s(VN ×D ; I2D−1 ) ⊆ AN M . Tak2
φ(VN ×D ;I2D−1 )∈ π π ,M ] Φ2D−2 ×(− M
ing into consideration only cells into the region of interπ , π ], we observe that |J (VN ×D )| = est Φ2D−2 × (− M i M 2(D−i)−2 M D−1 N N −i M i=0 2 − 1 , i.e. there are i 2(D−i)−1) 2 π π |J (VN ×D )| candidate vectors s in Φ2D−2 ×(− M ,M ], associated with cells each of which minimizes φ2D−1 component at a single point that constitutes the intersection of the corresponding 2D − 1 hypersurfaces. It can also be shown that if π π ,M ], we take into consideration all regions in Φ2D−2 × (− M all candidates form the candidate set S given by S(VN ×D ) =
D−1
J (VN ×(D−d) ).
(12)
d=0
with cardinality |S(VN ×D )| = |J (VN ×D )|+· · ·+|J (VN ×1 )| i D d−1 N N −i M 2(d−i)−2 M = = d=1 i=0 i 2 2 − 1 2(d−i)−1 2D−1 O(( MN ). 2 )
4. ALGORITHMIC DEVELOPMENTS AND COMPLEXITY In this section, we present the basic steps of the proposed algorithm for the construction of S(VN ×D ) for arbitrary N, D ∈ N, D < N and M ∈ {2k | k = 1, 2 . . . }. From (12), we observe that the initial problem of the determination of S(VN ×D ) can be divided into smaller parallel construction problems of J (VN ×d ) for d = 1, . . . , D. Moreover, the construction of J (VN ×d) ) can be fully parallelized since the candidate vector(s) s(VN ×d ; I2d−1 ) can be computed independently for each I2d−1 . In the following, we assume a certain value for d ∈ {1, . . . , D} and a certain set of indices I2d−1 = {i1 , . . . , i2d−1 }. According to the derivations in the previous section, the com˜i ˜ i1 ,1:2d ), . . . , H(V ) bination of hypersurfaces H(V 2d−1 ,1:2d intersects at a single point φ(VN ×d ; I2d−1 ) that “leads” Q cells associated with Q different M -phase candidate vectors sq (VN ×d ; I2d−1 ), q = 1, 2, . . . , Q. It can be shown that the evaluation of the decision function in (10) at the intersection of the 2D − 1 hypersurfaces under consideration determines definitely the corresponding symbol sn if and only if no hypersurface originates from Vn,1:d . For the hypersurfaces that pass through the intersection, the rule in (10) becomes ambiguous. In such a case, we have constructed disambiguation rules that solve the ambiguity in polynomial
3858
30
80
10
10 RankD Proposed algorithm Exhaustive search Reverse Search Method
RankD Proposed algorithm Exhaustive search 70
10
25
10
60
10 20
10
50
Cardinality
Cardinality
10
15
10
40
10
30
10 10
10
20
10 5
10
10
10
0
10
0
0
5
10
15
20
25
30
Sequence length N
35
40
45
10
50
0
5
10
15
20
25
30
Sequence length N
35
40
45
50
Fig. 1. Cardinality of candidate set S(VN×D ) of proposed method, reverse search, and exhaustive search (D = 2, M = 4).
Fig. 2. Cardinality of candidate set S(VN×D ) of proposed algorithm and exhaustive search (D = 3, M = 32).
time with respect to the length algorithm visits inde N . The 2D−1 pendently |S(VN ×D )| = O ( MN intersections and ) 2 computes the candidate M -phase vector(s) associated with each intersection. For each I2d−1 , the cost of the algorithm . Therefore the overall complexity of the algois O MN 2 rithm for the computation S(V fixed of MN N ×D) with D < N MN 2D 2D−1 O = O ( . ) ) becomes O ( MN 2 2 2 We observe that the computation of the candidate vectors of S(VN ×D ) is performed independently from cell to cell, which implies that there is no need to store the data that have been used for each candidate and we only have to store the “best” vector that has been met. Therefore, the proposed method is fully parallelizable and rank-scalable and its memory utilization is efficiently minimized, in contrast to the incremental algorithm in [2]. Compared to previous works on the maximization of a complex rank-deficient quadratic form over a finite field, we recall that N −1 the reverse search method [3],[4] computes D−1 (as many as our proposed algorithm) candii=0 i 2D−1 dates for M = 2 and i=0 2Ni−1 (twice as many as our proposed algorithm) candidates for M = 4 [see Fig. 1]. Additionally, the corresponding complexity of the algorithm proposed in [3],[4] is of the order O(N 2D LP( MN 2 , 2D)) and O (2N )2D LP( MN , 2D) for M = 2, 4, respectively, where 2 LP( MN , 2D) is the time to solve a linear programming (LP) 2 inequalities in 2D variables. optimization problem with MN 2 It turns out that the complexity of reverse search method MN is O ( 2 )2D+1 for M = 2, 4, i.e. one order of magnitude higher than the proposed algorithm. In addition, the reverse search method is restricted to M = 2, 4. On the other hand, the incremental algorithm proposed in [2],[8] is a timeefficient algorithm that solves the maximization problem of interest but becomes impractical even for moderate values of D since it follows an “incremental” strategy to construct the candidate set: it solves the problem inductively and, thus, it is
too complicated to be implemented. Furthermore, the critical disadvantage of this method is its memory inefficiency since it needs to store all the extreme points, all faces and their incidences in memory. Finally, the algorithm proposed in [5] deals optimally the problem of the maximization of a rankdeficient quadratic form for any M = {2k |k = 0, 1, . . . } but only for D = 1.
3859
5. REFERENCES [1] J.-A. Ferrez, K. Fukuda, and T. M. Liebling, “Solving the fixed rank convex quadratic maximization in binary variables by a parallel zonotope construction algorithm,” Eur. J. Oper. Res., vol. 166, pp. 35-50, 2005. [2] H. Edelsbrunner, J. O’Rourke, and R. Seidel, ‘Constructing arrangements of lines and hyperplanes with applications,” SIAM J. Comput., vol. 15, pp. 341-363, May 1986. [3] D. Avis and K. Fukuda, “Reverse Search for enumeration,” Discrete Appl. Math., vol. 65, pp. 21-46, Mar. 1996. [4] V. Pauli, L. Lampe, R. Schober, and K. Fukuda, “Multiplesymbol differential detection based on combinatorial geometry,” IEEE Trans. Commun., vol. 56, pp. 1596-1600, Oct. 2008. [5] K. M. Mackenthun Jr., “A Fast Algorithm for Multiple-Symbol Differential Detection of MPSK,” IEEE Trans. Commun., vol. 42, pp. 1471-1474, Apr. 1994. [6] G. N. Karystinos and A. P. Liavas, “Efficient computation of the binary vector that maximizes a rank-deficient quadratic form,” IEEE Trans. Inform. Theory, vol. 56, pp. 3581-3593, July 2010. [7] D. S. Papailiopoulos and G. N. Karystinos, “Polynomialcomplexity maximum-likelihood block noncoherent MPSK detection,” in Proc. IEEE ICASSP 2008, Las Vegas, NV, Apr. 2008, pp. 2681-2684. [8] I. Motedayen, A. Krishnamoorthy, and A. Anastasopoulos, “Optimal joint detection/estimation in fading channels with polynomial complexity,” IEEE Trans. Inform. Theory, vol. 53, pp. 209-223, Jan. 2007.
