PHASE RETRIEVAL USING COORDINATE DESCENT TECHNIQUES

Abstract

Coordinate descent is applied to recover a signal-of-interest from only magnitude information. In doing so, a single unknown value is solved at each iteration, while all other variables are held constant. As a result, only minimization of a univariate quartic polynomial is required, which is efficiently achieved by finding the closed-form roots of a cubic polynomial. Cyclic, randomized, and/or a greedy coordinate descent technique can be used. Each coordinate descent technique globally converges to a stationary point of the nonconvex problem, and specifically, the randomized coordinate descent technique locally converges to the global minimum and attains exact recovery of the signal-of-interest at a geometric rate with high probability when the sample size is sufficiently large. The cyclic and randomized coordinate descent techniques can also be modified via minimization of the l.sub.1-regularized quartic polynomial for phase retrieval of sparse signals-of-interest, i.e., those signals with only a few nonzero elements.

Claims

1. A method for obtaining phase information for a signal-of-interest from only magnitude information for the signal-of-interest, the method comprising: expressing the magnitude information as an objective function, f(x) the objective function being in the form of a multivariate quartic polynomial; selecting a coordinate index, i.sub.k, according to one or more predefined coordinate descent rules; and iteratively minimizing, according to the selected coordinate index, only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value.

2. The method of claim 1 where selecting the coordinate index, i.sub.k, comprises: cyclically assigning i.sub.k sequential values from 1 to 2N for each of the variable components of the multivariate quartic polynomial.

3. The method of claim 1 where selecting the coordinate index, i.sub.k, comprises: randomly assigning i.sub.k values from 1 to 2N, where each value is assigned with equal probability, for each of the variable components of the multivariate quartic polynomial.

4. The method of claim 1 where selecting the coordinate index, i.sub.k, comprises: assigning i.sub.k as: $i_k=\arg .Math..Math.\underset{i}{\min }.Math..Math.f_i\left({\stackrel{\_}{x}}^k\right).Math.,.Math.f_i\left({\stackrel{\_}{x}}^k\right)=\frac{f\left({\stackrel{\_}{x}}^k\right)}{{\stackrel{\_}{x}}_i^k}.$

5. The method of claim 1 where iteratively minimizing only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value comprises: for the kth iteration, minimizing f(x.sup.k) with respect to the i.sub.kth (i.sub.k {1, . . . , 2N}) variable while maintain remaining 2N1 {x.sub.i.sup.k}.sub.ii.sub.k, variable components at a constant value.

6. The method of claim 1 where iteratively minimizing only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value comprises: solving the only one variable component while maintain every other variable component at a constant value.

7. The method of claim 1 where the one or more predefined coordinate descent rules is selected from the group consisting of: a cyclic rule; a random rule; and a greedy rule.

8. The method of claim 1 where the signal-of-interest is a sparse signal and one or more of the predefined coordinate descent rules is obtained by minimizing a l.sub.1-regularized quartic polynomial.

9. A system for obtaining phase information for a signal-of-interest from only magnitude information for the signal-of-interest, the system comprising: a memory; one or more processors coupled to the memory, the one or more processors: expressing the magnitude information as an objective function, f(x), the objective function being in the form of a multivariate quartic polynomial; selecting a coordinate index, i.sub.k, according to one or more predefined coordinate descent rules; and iteratively minimizing, according to the selected coordinate index, only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value.

10. The system of claim 9 where selecting the coordinate index, i.sub.k, comprises: cyclically assigning i.sub.k sequential values from 1 to 2N for each of the variable components of the multivariate quartic polynomial.

11. The system of claim 9 where selecting the coordinate index, i.sub.k, comprises: randomly assigning i.sub.k values from 1 to 2N, were each value is assigned with equal probability, for each of the variable components of the multivariate quartic polynomial.

12. The system of claim 9 where selecting the coordinate index, i.sub.k, comprises: assigning i.sub.k as: $i_k=\arg .Math..Math.\underset{i}{\min }.Math..Math.f_i\left({\stackrel{\_}{x}}^k\right).Math.,.Math.f_i\left({\stackrel{\_}{x}}^k\right)=\frac{f\left({\stackrel{\_}{x}}^k\right)}{{\stackrel{\_}{x}}_i^k}.$

13. The system of claim 9 where iteratively minimizing only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value comprises: for the kth iteration, minimizing f(x.sup.k) with respect to the i.sub.kth (i.sub.k {1, . . . , 2N}) variable while maintain remaining 2N1 {x.sub.i.sup.k}.sub.ii.sub.k variables at a constant value.

14. The system of claim 9 where iteratively minimizing only one variable component of the multivariate quartic polynomial while maintaining every other variable component of the multivariate quartic polynomial at a constant value comprises: solving the only one variable component while maintain every other variable component at a constant value.

15. The system of claim 9 where the one or more predefined coordinate descent rules is selected from the group consisting of: a cyclic rule; a random rule; and a greedy rule.

16. The system of claim 9 where the signal-of-interest is a sparse signal and one or more of the predefined coordinate descent rules are obtained by minimizing a regularized quartic polynomial.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0009] For a more complete understanding of the present invention, reference is now made to the following descriptions taken in conjunction with the accompanying drawings, in which:

[0010] FIG. 1 shows the normalized reduction of objective function versus number of iterations/cycles at SNR=20 dB according to one or more embodiments;

[0011] FIG. 2 shows relative recovery error versus number of iterations/cycles at SNR=20 dB according to one or more embodiments;

[0012] FIG. 3 shows empirical probability of success versus number of noise-free measurements according to one or more embodiments;

[0013] FIG. 4 shows NMSE of recovered signal versus SNR according to one or more embodiments;

[0014] FIG. 5 shows the magnitudes of recovered signals according to one or more embodiments, where lines notated with a star (*) and lines notated with a hollow circle () denote the recovered and true signals, respectively;

[0015] FIG. 6 shows the probability of success versus number of noise-free measurements for sparse phase retrieval according to one or more embodiments;

[0016] FIG. 7 shows scatter plots of constellations of received signal and equalizer outputs according to one or more embodiments;

[0017] FIG. 8 shows ISI versus number of iterations/cycles according to one or more embodiments;

[0018] FIG. 9 shows a phase retrieval system according to one or more embodiments; and

[0019] FIG. 10 shows a phase retrieval method according to one or more embodiments.

DETAILED DESCRIPTION OF THE INVENTION

[0020] Described embodiments provide for recovering a complete signal-of-interest, i.e., the signal's phase information and magnitude information, from only intensity measurements of the signal. Doing so is generally referred to herein as phase retrieval. Conventionally, phase retrieval involves restoring a time-domain signal from power spectrum observations of the signal-of-interest, although the Fourier transform can be generalized to any linear mappings. Nevertheless, phase retrieval is a nonconvex optimization problem where minimizing a multivariate fourth-order polynomial is required. This problem is generally considered to be NP-hard, where NP-hardness (i.e., nondeterministic polynomial-time hardness), in computational complexity theory, is a class of problems that are, informally, at least as hard as the hardest problems in NP. More precisely, a problem, H, is NP-hard when every problem, L, in NP can be reduced in polynomial time to H. That is, a given solution for problem L can be verified to be a solution for problem H in polynomial time. As a consequence, finding a polynomial algorithm to solve any NP-hard problem would give polynomial algorithms for all the problems in NP, which is unlikely as many of them are considered hard.

[0021] According to the inventive concepts described herein, coordinate descent techniques are utilized to solve the foregoing problems more efficiently than previously known approaches. That is, a single unknown is solved at each iteration while all other variables are held constant. As a result, only minimization of a univariate quartic polynomial is required, which is efficiently achieved by determining the closed-form roots of a cubic polynomial. Embodiments can utilize one or more of three (3) coordinate descent techniques, referred to herein as (1) cyclic, (2) randomized, and (3) greedy coordinate descent techniques. As discussed, these three (3) coordinate descent techniques globally converge to a stationary point of the nonconvex problem, and specifically, the randomized coordinate descent technique locally converges to the global minimum and attains exact recovery of phase information at a geometric rate with high probability when the sample size is sufficiently large.

[0022] Further, the cyclic and randomized coordinate descent techniques can be modified by minimizing the l.sub.1-regularized quartic polynomial for phase retrieval of sparse signals, i.e., signals with only a few nonzero elements. Also, a specific application of the three (3) coordinate descent techniques to blind equalization in digital communications is described herein. Experimental evaluations reveal that described embodiments have been demonstrated to be superior to known approaches in terms of computational simplicity and/or recovery performance.

[0023] To better illustrate the inventive concepts described herein, general mathematical principles of phase retrieval are first discussed. Mathematically, the problem of phase retrieval is to find the complex-valued signal, x .sup.N, from M phaseless observations, b.sub.m . This is shown by equation (1):

b.sub.m=|a.sub.m.sup.Hx|.sup.2+v.sub.m, m=1, . . . , M equation (1)

where a.sub.m .sup.N are known sampling vectors, v.sub.m are additive zero-mean noise terms, and M>N. It should be appreciated that in optical imaging applications, a.sub.m is the discrete Fourier transform vector. If b.sub.m=|a.sub.m.sup.Hx|+v.sub.m, it can be converted to b.sub.m.sup.2. The x can be recovered up to a global phase, [0, 2), because e.sup.ix is also a solution. Adopting the least squares criterion, x is determined from equation (2):

[00001]$\begin{array}{cc}\underset{x{}^N}{\min }.Math.f\left(x\right):=\underset{m=1}{\overset{M}{.Math.}}.Math.{\left({.Math.a_m^H.Math.x.Math.}^2-b_m\right)}^2.& \mathrm{equation}.Math..Math.\left(2\right)\end{array}$

[0024] In terms of real-valued representations, equation (2) is equivalent to the expression shown at equation (3):

[00002]$\begin{array}{cc}\underset{\stackrel{\_}{x}{}^{2.Math.N}}{\min }.Math.f\left(\stackrel{\_}{x}\right):=\underset{m=1}{\overset{M}{.Math.}}.Math.{\left(q_m\left(\stackrel{\_}{x}\right)-b_m\right)}^2..Math.\mathrm{where}& \mathrm{equation}.Math..Math.\left(3\right)\\ \stackrel{\_}{x}=\left[\begin{array}{c}\mathrm{Re}\left(x\right)\\ \mathrm{Im}\left(x\right)\end{array}\right].Math..Math.\mathrm{and}& \mathrm{equation}.Math..Math.\left(4\right)\\ q_m\left(\stackrel{\_}{x}\right)={\stackrel{\_}{x}}^T.Math.{\stackrel{\_}{A}}_m.Math.\stackrel{\_}{x},.Math.{\stackrel{\_}{A}}_m=\left[\begin{array}{cc}\mathrm{Re}\left(x\right)& -\mathrm{Im}\left(A_m\right)\\ \mathrm{Im}\left(x\right)& \mathrm{Re}\left(A_m\right)\end{array}\right]{}^{2.Math.N2.Math.N},.Math.A_m=a_m.Math.a_m^H{}^{NN}.& \left(5\right)\end{array}$

[0025] The objective function f (x) is a multivariate quartic polynomial of x=[x.sub.1 . . . x.sub.2N].sup.T, implying that its minimization is generally NP-hard. Again, NP-hardness (i.e., non-deterministic polynomial-time hardness) is a class of problems that are, informally, at least as hard as the hardest problems in NP. More precisely, a problem, H, is NP-hard when every problem, L, in NP can be reduced in polynomial time to H. That is, a given solution for problem L can be verified to be a solution for problem H in polynomial time. As a consequence, finding a polynomial algorithm to solve any NP-hard problem would give polynomial algorithms for all the problems in NP, which is unlikely as many of them are considered hard.

[0026] According to described embodiments, the coordinate update strategy is exploited to minimize f(x). Coordinate descent is an iterative procedure that successively minimizes the objective function along coordinate directions. The result of the kth iteration can be denoted as x.sup.k=[x.sub.1.sup.k . . . x.sub.2N.sup.k].sup.T. In the kth iteration, f(x.sup.k) is minimized with respect to the i.sub.k th. (i.sub.k {1, . . . , 2N}) variable while keeping the remaining 2N1 {x.sub.i.sup.k}.sub.ii.sub.k constant. This is equivalent to performing a one-dimensional search along the i.sub.kth coordinate, which can be expressed as:

[00003]$\begin{array}{cc}{}_k=\arg .Math..Math.\underset{}{\min }.Math.f\left({\stackrel{\_}{x}}^k+.Math..Math.e_{i_k}\right)& \mathrm{equation}.Math..Math.\left(6\right)\end{array}$

where e.sub.i.sub.k is the unit vector with the i.sub.kth entry being one and all other entries being zero. From the foregoing, x is updated by:

x.sup.k+1=x.sup.k+.sub.ke.sub.i.sub.k equation (7)

which implies that only the k.sub.kth component is adjusted:

x.sub.i.sub.k.sup.k+1x.sub.i.sub.k.sup.k+.sub.k equation (8)

while all other components remain constant. Since x.sup.k is known, f(x.sup.k+e.sub.i.sub.k) is a univariate function of . Therefore, equation (6) is a one-dimensional minimization problem.

[0027] A coordinate descent solution is shown in the following:

TABLE-US-00001 Algorithm 1 CD for Phase Retrieval Initialization: Determine x.sup.0 .sup.2N using the spectral method, referenced herein: for k = 0, 1 . . . , do Choose index i.sub.k {1, . . . , 2N}: [00004]$\begin{array}{c}{}_k.Math.=.Math.a.Math..Math.\mathrm{rg}.Math..Math.\underset{\mathrm{aR}}{.Math.\min }.Math..Math.f.Math..Math.\left({\stackrel{\_}{x}}^k+e_{i_k}\right);\\ {\stackrel{\_}{x}}_{i_k}^{k.Math.+.Math.1}{x}_{i_k}^k+a_k;\end{array}$ Stop if f (x.sup.k) f (x.sup.k + 1) < is satisfied, where > 0 is a small tolerance parameter. end for

[0028] According to the inventive concepts, prior to presenting a solution to equation (6), one or more of three (3) selection rules for the coordinate index, i.sub.k, are followed by the phase retrieval system described herein. According to one embodiment, a cyclic rule is followed where i.sub.k is first assigned a value of 1, then 2, and so forth through 2N. These steps are then repeated, starting with i.sub.k=1 again. That is, i.sub.k is cyclically assigned a value from {1, . . . , 2N}. Accordingly, one cycle corresponds to 2N iterations. According to another embodiment, a random rule is followed where i.sub.k is randomly selected from {1, . . . , 2N} with equal probability. According to yet another embodiment, a greedy rule is followed where i.sub.k is chosen as

[00005]$\begin{array}{cc}{i}_k=\arg .Math..Math.\underset{i}{\max }.Math..Math.f_i\left({\stackrel{\_}{x}}^k\right).Math.,.Math.f_i\left({\stackrel{\_}{x}}^k\right)=\frac{f\left({\stackrel{\_}{x}}^k\right)}{{\stackrel{\_}{x}}_i^k}.& \mathrm{equation}.Math..Math.\left(9\right)\end{array}$

[0029] That is, in following the greedy rule the coordinate with the largest absolute value of the partial derivative is chosen. Therefore, in following the greedy rule, the phase retrieval system must compute the full gradient at each iteration, which is not required in following the cyclic and random rules. The three (3) coordinate descent techniques: cyclic, random, and greedy rules are referred to as CCD, RCD, and GCD, respectively. Also, experimental evaluations show that the GCD converges faster than CCD and RCD at the expense of the extra full gradient calculation.

[0030] Next, the closed-form solution for equation (6) can be derived. For notational simplicity, the superscript and subscript k in equation (6) are omitted. That is, given x for any current iteration and the search direction e.sub.i, equation (10) can be followed:

[00006]$\begin{array}{cc}\underset{}{\min }.Math.\left(\right):=f\left(\stackrel{\_}{x}+.Math..Math.e_i\right).& \mathrm{equation}.Math..Math.\left(10\right)\end{array}$

[0031] Employing equation (3), () is expressed as

[00007]$\begin{array}{cc}\left(\right)=\underset{m=1}{\overset{M}{.Math.}}.Math.{\left(q_m\left(\stackrel{\_}{x}+.Math..Math.e_i\right)-b_m\right)}^2& \mathrm{equation}.Math..Math.\left(11\right)\end{array}$

where the mth term is

.sub.m()=(q.sub.m(x+e.sub.1)b.sub.m).sup.2. equation (12)

[0032] Expanding the quadratic function in equation (12) yields:

q.sub.m(x+e.sub.i)=.sup.2e.sub.i.sup.TA.sub.me.sub.i+2e.sub.i.sup.TA.sub.mx+x.sup.TA.sub.mxc.sub.2,i.sup.m.sup.2+c.sub.1,i.sup.m+c.sub.0.sup.m equation (13)

where c.sub.2,i.sup.m, c.sub.1,i.sup.m and c.sub.0.sup.m are the coefficients of the univariate quadratic polynomial. Note that c.sub.0.sup.m has no relation with i. According to equation (5), the coefficient of the quadratic term can be simplified to:

[00008]$\begin{array}{cc}{c}_{2,i}^m=e_i^T.Math.{\stackrel{\_}{A}}_m.Math.e_i={\left[{\stackrel{\_}{A}}_m\right]}_{i,i}=\{\begin{array}{c}{.Math.{\left[a_m\right]}_i.Math.}^2,i=1,\mathrm{.Math.}.Math.,N\\ {.Math.{\left[a_m\right]}_{i-N}.Math.}^2,i=N+1,\mathrm{.Math.}.Math.,2.Math.N\end{array}& \mathrm{equation}.Math..Math.\left(14\right)\end{array}$

where [A].sub.i,i and [].sub.i represent the (i, i) entry of A and ith element of , respectively. It is also revealed from equation (5) that

[00009]$\begin{array}{cc}{\stackrel{\_}{A}}_m.Math.\stackrel{\_}{x}=\left[\begin{array}{c}\mathrm{Re}\left(A_m.Math.x\right)\\ \mathrm{Im}\left(A_m.Math.x\right)\end{array}\right]=\left[\begin{array}{c}\mathrm{Re}\left(\left(a_m^H.Math.x\right).Math.a_m\right)\\ \mathrm{Im}\left(\left(a_m^H.Math.x\right).Math.a_m\right)\end{array}\right].& \mathrm{equation}.Math..Math.\left(15\right)\end{array}$

[0033] Therefore, the coefficient of the linear term in equation (13) is computed as

[00010]$\begin{array}{cc}{c}_{1,i}^m=2.Math.e_i^T.Math.{\stackrel{\_}{A}}_m.Math.\stackrel{\_}{x}=\{\begin{array}{c}\mathrm{Re}\left({\left(a_m^H.Math.x\right)\left[a_m\right]}_i\right),i=1,\mathrm{.Math.}.Math.,N\\ \mathrm{Im}\left({\left(a_m^H.Math.x\right)\left[a_m\right]}_{i-N}\right),i=N+1,\mathrm{.Math.}.Math.,2.Math.N\end{array}.& \mathrm{equation}.Math..Math.\left(16\right)\end{array}$

[0034] The constant term is

c.sub.0.sup.m=x.sup.TA.sub.mx=x.sup.HA.sub.mx=|a.sub.m.sup.Hx|.sup.2. equation (17)

[0035] Since q.sub.m(x+e.sub.i) is quadratic, .sub.m() of equation (12) is a univariate quartic polynomial of . Inserting equation (13) into equation (12) results in:

.sub.m()=d.sub.4,i.sup.m.sup.4+d.sub.3,i.sup.m.sup.3+d.sub.2,i.sup.m.sup.2+d.sub.1,i.sup.m+d.sub.0.sup.m equation (18)

where

d.sub.4,i.sup.m=(c.sub.2,i.sup.m).sup.2

d.sub.3,i.sup.m=2c.sub.2,i.sup.mc.sub.1,i.sup.m

d.sub.2,i.sup.m=(c.sub.1,i.sup.m).sup.2+2c.sub.2,i.sup.m(c.sub.0.sup.mb.sub.m)

d.sub.1,i.sup.m=2c.sub.1,i.sup.m(c.sub.0.sup.mb.sub.m)

d.sub.0.sup.m=(c.sub.0.sup.mb.sub.m).sup.2. equation (19)

Note that d.sub.0.sup.m is not related to i.

[0036] Recalling ()=.sub.m=1.sup.M.sub.m(), it is clear that the coefficients of the quartic polynomial

()=d.sub.4,i.sup.4+d.sub.3,i.sup.3+d.sub.2,i.sup.2+d.sub.1,i+d.sub.0 equation (20)

are the sum of those of {.sub.m()}.sub.m=1.sup.M:

[00011]$\begin{array}{cc}{d}_0=\underset{m=1}{\overset{M}{.Math.}}.Math.d_0^m,d_{j,i}=\underset{m=1}{\overset{M}{.Math.}}.Math.d_{j,i}^m,j=1,\mathrm{.Math.}.Math.,4.& \mathrm{equation}.Math..Math.\left(21\right)\end{array}$

[0037] The minimum point of () is obtained by setting its derivative to zero:

()=4d.sub.4,i.sup.3+3d.sub.3,i.sup.2+2d.sub.2,i+d.sub.1,i=0. equation (22)

[0038] Equation (22) refers to finding the roots of a univariate cubic polynomial, which is easy and fast because there is a closed form solution. Due to the real-valued coefficients of the cubic polynomial, there are only two possible cases on the roots. The first case is that equation (22) has a real root and a pair of complex conjugate roots. In this case, the minimizer is the unique real root because the optimal solution of a real-valued problem must be real. The second case is that equation (22) has three real roots. Then the optimal a is the real root associated with the minimum objective. Once the coefficients of equation (22) are determined, the complexity of calculating the roots of a cubic polynomial is merely O(1).

[0039] In view of the foregoing, the complexity requirement of the coordinate descent techniques is now examined according to embodiments described herein. The most significant computational cost at each iteration is calculating the coefficients {d.sub.j,i}.sub.j=1.sup.4, or equivalently, computing and c.sub.2,i.sup.m, m.sub.1,i.sup.m, and c.sub.0.sup.m, m=1, . . . M, because {d.sub.j,i}.sub.j=1.sup.4 can be easily calculated from c.sub.2,i.sup.m, c.sub.1,i.sup.m, and c.sub.0.sup.m according to equation (19) and equation (21). From equation (14), c.sub.2,i.sup.m is the squared modulus of [a.sub.m].sub.i and can be pre-computed in advance before iteration, which requires O(M) multiplications for calculating all M coefficients {c.sub.2,i.sup.m}.sub.m=1.sup.M. According to equation (16) and equation (17), {a.sub.m.sup.Hx}.sub.m=1.sup.M must be computed to obtain {c.sub.1,i.sup.m}.sub.m=1.sup.M and {c.sub.0.sup.m}.sub.m=1.sup.M. This involves a matrix-vector multiplication Ax, where the mth row of the matrix A is a.sub.m.sup.H, that is:

[00012]$\begin{array}{cc}A=\left[\begin{array}{c}{a}_1^H\\ \mathrm{.Math.}\\ a_M^H\end{array}\right]{}^{MN}.& \mathrm{equation}.Math..Math.\left(23\right)\end{array}$

[0040] At first glance, the matrix-vector multiplication requires a complexity of O(MN). However, this complexity can be reduced to O(M) per iteration for the coordinate descent technique. By observing equation (7), it is determined that only a single element changes between two consecutive iterations. Specifically, this is seen from:

[00013]$\begin{array}{cc}{x}^{k+1}-x^k=\left(x_j^{k+1}-x_j^k\right).Math.e_j,j=\{\begin{array}{c}{i}_k,\mathrm{if}.Math..Math.jN\\ i_k-N,\mathrm{otherwise}\end{array}& \mathrm{equation}.Math..Math.\left(24\right)\end{array}$

which yields
where A.sub.:,j represents the jth column of A. It is only required to compute Ax.sup.0 before an iteration. Afterward, the matrix-vector product can be efficiently updated from that of the previous iteration by an efficient computation of a scalar-vector multiplication, (x.sub.j.sup.k+1x.sub.j.sup.k)A.sub.:,j, which merely requires O(M) operations.

[0041] In sum, the complexity of CCD and RCD is on the order of O(M) per iteration. Therefore, the complexity in 2N iterations, for a cycle in which CCD is followed, is the same as that of the WF using full gradient descent. While for GCD, an extra cost for computing the full gradient is needed, which results in a complexity of O(MN).

[0042] Further, as previously mentioned, the three (3) coordinate descent techniques converge to a stationary point from any initial value, and the RCD converges to the global minimum and achieves exact retrieval with high probability provided that the sample size is large enough.

[0043] According to embodiments described herein, the foregoing concepts can also be utilized to more efficiently retrieve phase values of sparse signals. If x is sparse, then the real-valued x is also sparse. The following l.sub.1-regularization, which is commonly used in compressed sensing, for sparse phase retrieval is adopted:

[00014]$\begin{array}{cc}\underset{\stackrel{\_}{x}{}^{2.Math.N}}{\min }.Math.g\left(\stackrel{\_}{x}\right):=\underset{m=1}{\overset{M}{.Math.}}.Math.{\left({\stackrel{\_}{x}}^T.Math.{\stackrel{\_}{A}}_m.Math.\stackrel{\_}{x}-b_m\right)}^2+.Math.{.Math.\stackrel{\_}{x}.Math.}_1& \mathrm{equation}.Math..Math.\left(26\right)\end{array}$

where x.sub.1=.sub.i|x.sub.i| is the l.sub.1-norm, r>0 is the regularization factor, and g(x)=f(x)+rx.sub.1. As there is no gradient for equation (26), the GCD cannot be implemented because it requires a gradient for index selection. Therefore, only CCD and RCD are discussed for the l.sub.1-regularization, and they are referred to as l.sub.1-CCD and l.sub.1-RCD, respectively.

[0044] The steps of the coordinate descent technique for solving equation (26) are similar to those in shown above for recovery general signals. The only difference is that an l.sub.1-norm term is added to the scaler minimization problem of equation (10), which is shown as

[00015]$\begin{array}{cc}\underset{}{\min }.Math.\left\{\left(\right)+.Math.{.Math.\stackrel{\_}{x}+.Math..Math.e_i.Math.}_1\right\}.& \mathrm{equation}.Math..Math.\left(27\right)\end{array}$

[0045] By ignoring the terms independent to , equation (27) is equivalent to

[00016]$\begin{array}{cc}\underset{}{\min }.Math.\left\{\left(\right)+.Math..Math.+{\stackrel{\_}{x}}_i.Math.\right\}.& \mathrm{equation}.Math..Math.\left(28\right)\end{array}$

[0046] Making a change of variable =+x.sub.i, substituting =x.sub.i into equation (20), and ignoring the irrelevant components, we obtain an equivalent scalar minimization problem

[00017]$\begin{array}{cc}\underset{}{\min }.Math.\left(\right):=u_4.Math.{}^4+u_3.Math.{}^3+u_2.Math.{}^2+u_1.Math.+.Math..Math..Math.& \mathrm{equation}.Math..Math.\left(29\right)\end{array}$

where the coefficients of the quartic polynomial {u.sub.j}.sub.j=1.sup.4 are calculated as

u.sub.4=d.sub.4,i

u.sub.3=d.sub.3,i4x.sub.id.sub.4,i

u.sub.2=d.sub.2,i3x.sub.id.sub.3,i+6x.sub.i.sup.2d.sub.4,i

u.sub.1=d.sub.1,i2x.sub.id.sub.2,i+3x.sub.i.sup.2d.sub.3,i4x.sub.i.sup.3d.sub.4,i. equation (30)

[0047] Although () is non-smooth due to the absolute term, the closed-form solution of its minimum can still be derived. The minimizer of () is evaluated in two intervals, namely, [0, ) and (, 0). Next, the set S.sup.+ containing the stationary points of () in the interval [0, ) is defined. That is, S.sup.+ is the set of real positive roots of the cubic equation:

4u.sub.4.sup.3+3u.sub.3.sup.22+u.sub.2+(u.sub.i+r)=0, 0. equation (31)

[0048] The S.sup.+ can include zero (0) elements, or can include one (1) or three (3) elements since the cubic equation of (31) may have zero (0), one (1), or three (3) real positive roots. Similarly, S.sup. is the set that contains the stationary points of () in (, 0), that is, the real negative roots of

4u.sub.4.sup.3+3u.sub.3.sup.2+2u.sub.2+(u.sub.1r)=0, <0. equation (32)

[0049] Again, S.sup. can include zero (0) elements, or can include one (1) or three (3) elements. The minimizer of () in [0, ) must be the boundary, i.e., 0, or one (1) element of S.sup.+. The minimizer in (, 0) must be an element of S.sup.. In sum, the minimizer of equation (29) is limited to the set {0 S.sup.+ S.sup.}, which includes a maximum of seven (7) elements, i.e.,

[00018]$\begin{array}{cc}{}^{*}=\arg .Math..Math.\underset{}{\min }.Math..Math.\left(\right),\left\{0.Math.{}^{+}.Math.{}^{-}\right\}.& \mathrm{equation}.Math..Math.\left(33\right)\end{array}$

[0050] Therefore, () must be evaluated over a set of, at most, seven (7) elements only, whose computation is straightforward and not computationally intensive. The coordinate of the l.sub.1-regularized coordinate descent technique is updated as x.sub.i.sub.k.sup.k+1*. If S.sup.+ S.sup.=, then x.sub.i.sub.k.sup.k+1=*=0, which makes the solution sparse.

[0051] As discussed, certain embodiments can be applied to several different fields, including blind equalization in a wireless communication system. To further illustrate the inventive concepts, these embodiments, as applied blind equalization in a wireless communication system, as described herein: Consider a communication system with discrete-time complex baseband signal model:

r(n)=s(n)*h(n)+v(n) equation (34)

where r(n) is the received signal, s(n) is the transmitted data symbol, h(n) is the channel impulse response, v(n) is the additive white noise, and * denotes convolution. The received signal is significantly distorted due to the inter-symbol interference (ISI) induced by the propagation channel

[0052] Blind equalization aims at recovering the transmitted symbols without knowing the channel response. If the equalizer is defined with P coefficients w=[w.sub.0 . . . w.sub.P1].sup.T and r.sub.n=[r(n) . . . r(nP+1)].sup.T, the equalizer output is

[00019]$\begin{array}{cc}y\left(n\right)=\underset{i=0}{\overset{P-1}{.Math.}}.Math.w_i^{*}.Math.r\left(n-i\right)=w^H.Math.r_n.& \mathrm{equation}.Math..Math.\left(35\right)\end{array}$

[0053] As many modulated signals in communication, such as phase shift keying (PSK), frequency modulation (FM), and phase modulation (PM), are of constant modulus, the constant modulus criterion is applied to obtain the equalizer:

[00020]$\begin{array}{cc}\underset{w}{\min }.Math.\underset{n}{.Math.}.Math.{\left({.Math.w^H.Math.r_n.Math.}^2-\right)}^2,\mathrm{where}& \mathrm{equation}.Math..Math.\left(36\right)\\ =\frac{\left[{.Math.s\left(n\right).Math.}^4\right]}{\left[{.Math.s\left(n\right).Math.}^2\right]}& \mathrm{equation}.Math..Math.\left(37\right)\end{array}$

is called the dispersion constant. It is obvious that the problem of equation (36) has the same form as the phase retrieval of equation (2). Both of them involve multivariate quartic polynomials. The only difference between phase retrieval and blind equalization is that the decision variable of the former is the unknown signal x while that of the latter is the equalizer w. Therefore, the WF and coordinate descent methods can be applied to solve the blind equalization problem of equation (36). By defining the composite channel-equalizer response as v=(n)=h(n)*w(n), the quantified ISI, which is expressed as

[00021]$\begin{array}{cc}\mathrm{ISI}=\frac{\underset{n}{.Math.}.Math.{.Math.v\left(n\right).Math.}^2-\underset{n}{\max }.Math.{.Math.v\left(n\right).Math.}^2}{\underset{n}{\max }.Math.{.Math.v\left(n\right).Math.}^2}& \mathrm{equation}.Math..Math.\left(38\right)\end{array}$

reflects the equalization quality. From the foregoing, smaller values of ISI imply better equalization. Further, all methods use the same initial value obtained from the spectral method. The measurement vectors satisfy a complex standard i.i.d Gaussian distribution.

[0054] The convergence behavior of the three (3) coordinate descent techniques is now discussed. The signal x and noise v.sub.m in equation (1) are i.i.d. Gaussian. In experimental evaluations, the following parameters were set: N=64 and M=6N. The results of WF are included for comparison. Note that it is fair to compare 2N iterations (one cycle) for the coordinate descent technique with one WF iteration because the computational complexity of the CCD and RCD per cycle is the same as the WF per iteration. The GCD has a higher complexity for every 2N iterations than WF, CCD, and RCD.

[0055] Two quantities are plotted to evaluate the convergence rate. The first quantity is the reduction of the objective function normalized with respect to b.sup.2:

[00022]$\begin{array}{cc}\frac{f\left({\stackrel{\_}{x}}^k\right)-f\left({\stackrel{\_}{x}}^{*}\right)}{{.Math.b.Math.}^2},b={\left[b_1.Math..Math.\mathrm{.Math.}.Math..Math.b_M\right]}^T& \mathrm{equation}.Math..Math.\left(39\right)\end{array}$

where x* denotes the optimal solution. The second quantity is the relative recovery error:

[00023]$\begin{array}{cc}\frac{{.Math.\left({\stackrel{\_}{x}}^k\right)-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right).Math.}^2}{{.Math.{\stackrel{\_}{x}}^{*}.Math.}^2}& \mathrm{equation}.Math..Math.\left(40\right)\end{array}$

where T.sub..sub.k(x*) is extracted from x.sup.ke.sup.j.sup.kx* has the minimum value, which reflects the convergence speed to the original signal. The signal-to-noise ratio (SNR) in equation (1) is defined as

[00024]$\begin{array}{cc}S.Math..Math.N.Math..Math.R=\frac{E\left[{.Math.b.Math.}^2\right]}{M.Math..Math.{}_v^2}& \mathrm{equation}.Math..Math.\left(41\right)\end{array}$

where .sub.v.sup.2 is the variance of FIGS. 1 and 2 show the normalized objective reduction and recovery error, respectively, at SNR=20 dB. Fifty (50) independent trials were performed while the averaged results are also plotted with thick lines. It can be clearly seen that the three coordinate descent techniques converge faster than the previously known WF technique. As seen, among the observed techniques, the convergence speed of the GCD is the fastest.

Statistical Performance

[0056] For purposes of discussing statistical performance, the experimental settings are the same as those discussed above, except M and SNR vary. The performance of GS method is also examined. The empirical probability of success and normalized mean square error (NMSE), which is the mean of the relative recovery error in equation (40), are used to measure the statistical performance. In this case, all results are averaged over two hundred (200) independent trials. In the absence of noise, if the relative recovery error of a phase retrieval scheme is smaller than 10.sup.5, it is considered a successful exact recovery.

[0057] FIG. 3 shows the empirical probability of success versus the number of measurements M. It is observed that the GCD is slightly better than WF, while CCD and RCD are slightly inferior to the WF.

[0058] FIG. 4 plots the NMSEs versus SNR and it is seen that the three (3) coordinate descent techniques and WF have comparable NMSEs, while each is superior to the GS method.

Phase Retrieval of Sparse Signal

[0059] Phase retrieval of a sparse signal with K nonzero elements is also discussed with respect to described embodiments. In addition to WF, the two convex relaxation based methods, namely, PhaseLift and PhaseCut, and sparse GS methods using hard-thresholding are examined for comparison. The regularization factor was set to r=2.35M for the l.sub.1-CCD and l.sub.1-RCD. In this case, the support of the sparse signal is randomly selected from [1, N] where N=64. The real and imaginary parts of the nonzero coefficients of x are drawn as random uniform variables in the range [2/{square root over (2)}, 1/{square root over (2)}] [1/{square root over (2)}, 2/{square root over (2)}].

[0060] FIG. 5 shows the recovered signal with K=5 and M=2N in the noise-free case. As seen, the WF, PhaseLift, PhaseCut, and sparse GS methods cannot recover the signal when M is relatively small while the l.sub.1-CCD and l.sub.1-RCD work well.

[0061] FIG. 6 plots the probability of success versus number of measurements, which demonstrates the superiority of l.sub.1-CCD and l.sub.1-RCD.

Blind Equalization

[0062] The application of the coordinate descent techniques to blind equalization is discussed with respect to described embodiments and compared to WF method. The results of the conventional super-exponential (SE) algorithm are also included. The transmitted signal adopts quadrature PSK (QPSK) modulation, namely, s(n) {1, 1, j, j}. A typical finite impulse response (FIR) communication channel with impulse response {0.4, 1, 0.7, 0.6, 0.3, 0.4, 0.1} is adopted.

[0063] FIG. 7 shows the constellations of the received signal and equalizer outputs of 1000 samples with 20 dB zero-mean Gaussian noise. The received signal is severely distorted due to the channel propagation. The SE, WF, and coordinate descent techniques succeed in recovering the transmitted signal up to a global phase rotation. Clearly, the CCD and GCD have a higher recovery accuracy.

[0064] FIG. 8 plots the ISI versus number of iterations/cycles with 25 dB zero-mean Gaussian noise when there are 2000 samples. The ISI is averaged over 100 independent trials. In addition to faster convergence than the WF and SE methods, the coordinate descent techniques, especially CCD, arrive at a lower ISI. This means that the coordinate descent techniques achieve a more accurate recovery.

[0065] With regard to the convergence analysis of the coordinate descent techniques, the inventive concepts described herein prove that the three (3) coordinate descent techniques globally converge to a stationary point from any initial value. Further, the sequence of the iterations generated by the RCD locally converges to the global minimum point in expectation at a geometric rate under a mild assumption. This implies that in the absence of noise, the RCD achieves exact phase retrieval under a moderate condition.

Global Convergence to Stationary Point

[0066] Lemma 1: Given any finite initial value x.sup.0 .sup.2N and f(x.sup.0)=f.sub.0, the sublevel set of f(x)

S.sub.f.sub.0={x|f(x)f.sub.0}equation (42)

is compact, i.e., bounded and closed. The iterates of the three coordinate descent techniques, namely, x.sup.k, k=0, 1, . . . , are in the compact set S.sub.f.sub.0.

[0067] Proof: If x.fwdarw., then f(x).fwdarw. on since f(x) is quartic. The converse-negative proposition implies that f(x)f.sub.0< guarantees x for all x S.sub.f.sub.0. Hence, S.sub.f.sub.0 is bounded. Now it is clear that all the points in the sublevel set satisfy f(x) [0, f.sub.0] since f(x)0.

[0068] Since the mapping f(x) is continuous and the image [0, f.sub.0] is a closed set, the inverse image {x|0f(x)f.sub.0} is also closed. This completes the proof that S.sub.f.sub.0 is compact. The coordinate descent technique monotonically decreases f(x), meaning that f(x.sup.k)f(x.sup.k1) . . . f(x.sup.0)=f.sub.0. Therefore, all the iterates {x.sup.k} must be in the compact set S.sub.f.sub.0.

[0069] Lemma 1 guarantees the analysis can be limited in the compact set S.sub.f.sub.0 rather than the whole domain of f(x). The following lemma states that the partial derivatives of f(x) are locally Lipschitz continuous on S.sub.f.sub.0, although they are not globally Lipschitz continuous over the whole domain .sup.2N.

[0070] Lemma 2 On the compact set S.sub.f.sub.0, the gradient f(x) is component-wise Lipschitz continuous. That is, for each i=1, . . . , 2N, it follows

|.sub.if(x+te.sub.i).sub.if(x)|L.sub.i|i|, t equation (43)

for all x,x+te.sub.i S.sub.f.sub.0, where L.sub.i>0 is called the component-wise Lipschitz constant on S.sub.f.sub.0. Further, we have

[00025]$\begin{array}{cc}f\left(\stackrel{\_}{x}+{\mathrm{te}}_i\right)f\left(\stackrel{\_}{x}\right)+t.Math..Math.{}_i.Math.f\left(\stackrel{\_}{x}\right)+\frac{L_i}{2}.Math.t^2.& \mathrm{equation}.Math..Math.\left(44\right)\end{array}$

[0071] Proof For t=0, both sides of (43) are equal to 0 and (43) holds. For t0, it means that x+te.sub.ix=t0. Since the function

[00026]$\begin{array}{cc}\frac{.Math.{}_i.Math.f\left(\stackrel{\_}{x}-{\mathrm{te}}_i\right)-{}_i.Math.f\left(\stackrel{\_}{x}\right).Math.}{.Math.\stackrel{\_}{x}+{\mathrm{te}}_i-\stackrel{\_}{x}.Math.}& \mathrm{equation}.Math..Math.\left(45\right)\end{array}$

is continuous on the compact set S.sub.f.sub.0, its minimum over S.sub.f.sub.0 is L.sub.i, which is attained by Weierstrass' theorem. Then:

[00027]$\begin{array}{cc}\max .Math.\frac{.Math.{}_i.Math.f\left(\stackrel{\_}{x}-{\mathrm{te}}_i\right)-{}_i.Math.f\left(\stackrel{\_}{x}\right).Math.}{.Math.\stackrel{\_}{x}+{\mathrm{te}}_i-\stackrel{\_}{x}.Math.}=L_i& \mathrm{equation}.Math..Math.\left(46\right)\end{array}$

immediately elicits equation (43). The following second-order partial derivative

[00028]$\begin{array}{cc}{}_{i,i}^2.Math.f\left(\stackrel{\_}{x}\right)={\left[{}^2.Math.f\left(\stackrel{\_}{x}\right)\right]}_{i,i}={\left[\frac{{}^2.Math.f\left(\stackrel{\_}{x}\right)}{\stackrel{\_}{x}.Math.{\stackrel{\_}{x}}^T}\right]}_{i,i}=\frac{{}^2.Math.f\left(\stackrel{\_}{x}\right)}{{\stackrel{\_}{x}}_i^2}& \mathrm{equation}.Math..Math.\left(47\right)\end{array}$

is well defined since f(x) is twice continuously differentiable, Noting that .sub.i,i.sup.2f(x) is the partial derivative of .sub.if (x) with respect to x, and shown by equation (43),

[00029]$\begin{array}{cc}{}_{i,i}^2.Math.f\left(\stackrel{\_}{x}\right)=\underset{t.fwdarw.0}{\lim }.Math.\frac{{}_i.Math.f\left(\stackrel{\_}{x}-{\mathrm{te}}_i\right)-{}_i.Math.f\left(\stackrel{\_}{x}\right)}{t}{L}_i& \mathrm{equation}.Math..Math.\left(48\right)\end{array}$

is obtained, which holds for all x S.sub.f.sub.0. Applying Taylor's theorem and equation (48), there exists a [0, 1] with x+te.sub.i S.sub.f.sub.0 such that

[00030]$\begin{array}{cc}\begin{array}{c}f\left(\stackrel{\_}{x}-{\mathrm{te}}_i\right)=.Math.f\left(\stackrel{\_}{x}\right)+t.Math.{f\left(\stackrel{\_}{x}\right)}^T.Math.e_i+\frac{t^2}{2}.Math.e_i^T.Math.{}^2.Math.f\left(\stackrel{\_}{x}-.Math..Math.{\mathrm{te}}_i\right).Math.e_i\\ =.Math.f\left(\stackrel{\_}{x}\right)+t.Math.f_i\left(\stackrel{\_}{x}\right)+\frac{t^2}{2}.Math.{}_{i,i}^2.Math.f\left(\stackrel{\_}{x}-.Math..Math.{\mathrm{te}}_i\right)\\ .Math.f\left(\stackrel{\_}{x}\right)+t.Math.f_i.Math.\left(\stackrel{\_}{x}\right)+\frac{L_i}{2}.Math.t^2.\end{array}& \mathrm{equation}.Math..Math.\left(49\right)\end{array}$

[0072] The component-wise Lipschitz constant L.sub.i is not easy to compute or estimate because the partial derivatives are complicated multivariate polynomials. However, the coordinate descent techniques do not require L.sub.i. This quantity is just used for theoretical convergence analysis. The minimum and maximum of all the component-wise Lipschitz constants are L.sub.min=min.sub.1i2N L.sub.i and L.sub.max=max.sub.1i2N L.sub.i, respectively. Employing similar steps of Lemma 2, we can prove that the full gradient f(x) is Lipschitz continuous:

f(x)f(z)Lxzequation (50)

with the full Lipschitz constant L. It is not difficult to show L.sub.i L.sub.i and thus L2N L.sub.max is obtained.

[0073] Theorem 1: The CCD, RCD, and GCD globally converge to a stationary point of the multivariate quartic polynomial from an arbitrary initialization.

[0074] Proof: Based on the component-wise Lipschitz continuous property of equation (44), it is derived that:

[00031]$\begin{array}{cc}\begin{array}{c}f\left({\stackrel{\_}{x}}^{k+1}\right)=.Math.\min .Math..Math.f\left({\stackrel{\_}{x}}^k+.Math..Math.e_{i_k}\right)\\ .Math.f\left({\stackrel{\_}{x}}^k+.Math..Math.e_{i_k}\right).Math.{|}_{=-{}_{i_k}.Math.f\left({\stackrel{\_}{x}}^k\right)/L_{i_k}}\\ =.Math.f\left({\stackrel{\_}{x}}^k-\frac{{}_{i_k}.Math.f\left({\stackrel{\_}{x}}^k\right)}{L_{i_k}}.Math.e_{i_k}\right)\\ .Math.f\left({\stackrel{\_}{x}}^k\right)-\frac{{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2}{L_{i_k}}+\frac{L_{i_k}}{2}.Math.\frac{{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2}{L_{i_k}^2}\\ =.Math.f\left({\stackrel{\_}{x}}^k\right)-\frac{1}{2.Math.L_{i_k}}.Math.{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2\\ .Math.f\left({\stackrel{\_}{x}}^k\right)-\frac{1}{2.Math.L_{\max }}.Math.{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2\end{array}& \mathrm{equation}.Math..Math.\left(51\right)\end{array}$

from which we obtain a lower bound on the progress made by each coordinate descent technique iteration

[00032]$\begin{array}{cc}f\left({\stackrel{\_}{x}}^k\right)-f.Math.\left({\stackrel{\_}{x}}^{k+1}\right)\frac{1}{2.Math.L_{\max }}.Math.{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2.& \mathrm{equation}.Math..Math.\left(52\right)\end{array}$

[0075] For different rules of index selection, the right-hand side of equation (52) will differ. For GCD, it chooses the index with the largest partial derivative magnitude. Hence, by equation (9), it follows

[00033]$\begin{array}{cc}{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2=.Math.{(f\left({\stackrel{\_}{x}}^k\right).Math.}_{}^2\frac{1}{2.Math.N}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2.& \mathrm{equation}.Math..Math.\left(53\right)\end{array}$

[0076] Plugging equation (53) into equation (52) leads to the following lower bound of the progress of one GCD iteration

[00034]$\begin{array}{cc}f\left({\stackrel{\_}{x}}^k\right)-f.Math.\left({\stackrel{\_}{x}}^{k+1}\right)\frac{1}{4.Math.{\mathrm{NL}}_{\max }}.Math.{.Math.f_{i_k}\left({\stackrel{\_}{x}}^k\right).Math.}^2.& \mathrm{equation}.Math..Math.\left(54\right)\end{array}$

[0077] This means that one GCD iteration decreases the objective function with an amount of at least f(x.sup.k).sup.2/(4N L.sub.max). Setting is k=0, . . . , j, in equation (54) and summing over all inequalities yields

[00035]$\begin{array}{cc}\underset{k=0}{\overset{j}{.Math.}}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^{2}4.Math.{\mathrm{NL}}_{\max }\left(f\left({\stackrel{\_}{x}}^0\right)-f\left({\stackrel{\_}{x}}^{j+1}\right)\right)4.Math.{\mathrm{NL}}_{\max }.Math.f_0& \mathrm{equation}.Math..Math.\left(55\right)\end{array}$

where f(x.sup.j+1)0 is used. Taking the limit as j.fwdarw. on equation (55), a convergent series is obtained

[00036]$\begin{array}{cc}\underset{k=0}{\overset{}{.Math.}}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^{2}4.Math.{\mathrm{NL}}_{\max }.Math.f_0.& \mathrm{equation}.Math..Math.\left(56\right)\end{array}$

[0078] If a series converges, then its terms approach to zero, which indicates f(x.sup.k) is a zero vector as k.fwdarw., that is, the GCD converges to a stationary point.

[0079] For RCD, since i.sub.k is a random variable, f(x.sup.k+1) is also random:

[00037]$\begin{array}{cc}E\left[f\left({\stackrel{\_}{x}}^{k+1}\right)\right]E\left[f\left({\stackrel{\_}{x}}^k\right)-\frac{1}{2.Math.L_{\max }}.Math.{\left(f_{i_k}\left({\stackrel{\_}{x}}^k\right)\right)}^2\right]=f\left({\stackrel{\_}{x}}^k\right)-\frac{1}{2.Math.L_{\max }}.Math.\underset{i=1}{\overset{2.Math.N}{.Math.}}.Math.\frac{1}{2.Math.N}.Math.{\left(f_i\left({\stackrel{\_}{x}}^k\right)\right)}^2=f\left({\stackrel{\_}{x}}^k\right)-\frac{1}{4.Math.{\mathrm{NL}}_{\max }}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2& \mathrm{equation}.Math..Math.\left(57\right)\end{array}$

where the fact that i.sub.k is uniformly sampled from {1, . . . , 2N} with equal probability of 1/(2N) is employed. Then the RCD at least obtains a decrease on the objective function in expectation

[00038]$\begin{array}{cc}f.Math.\left({\stackrel{\_}{x}}^k\right)-E\left[f\left({\stackrel{\_}{x}}^{k+1}\right)\right]\frac{1}{4.Math.{\mathrm{NL}}_{\max }}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2.& \mathrm{equation}.Math..Math.\left(58\right)\end{array}$

[0080] Following similar steps in the GCD, it is sufficient to prove that the expected gradient of the RCD approaches to the zero (0) vector and thus converges to a stationary point in expectation.

[0081] For CCD, i.sub.k takes value cyclically from { 1, . . . , 2N}. Applying the property of CCD, a lower bound of the decrease of the objective function after 2N CCD iterations can be derived:

[00039]$\begin{array}{cc}f\left({\stackrel{\_}{x}}^k\right)-f\left({\stackrel{\_}{x}}^{k+2.Math.N}\right)\frac{{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2}{4.Math.L_{\max }\left(1+2.Math.{\mathrm{NL}}^2/L_{\min }^2\right)}.& \mathrm{equation}.Math..Math.\left(59\right)\end{array}$

[0082] Setting k=0, 2N, . . . , 2jN, in equation (59) and summing over all the inequalities yields

[00040]$\begin{array}{cc}\begin{array}{c}\underset{k=0}{\overset{j}{.Math.}}.Math.{.Math.f\left({\stackrel{\_}{x}}^{2.Math.\mathrm{kN}}\right).Math.}^2.Math.\frac{\left(f\left({\stackrel{\_}{x}}^0\right)-f\left({\stackrel{\_}{x}}^{2.Math.\left(j+1\right).Math.N}\right)\right)}{4.Math.L_{\max }\left(1+2.Math.{\mathrm{NL}}^2/L_{\min }^2\right)}\\ .Math.\frac{f\left({\stackrel{\_}{x}}^0\right)}{4.Math.L_{\max }\left(1+2.Math.{\mathrm{NL}}^2/L_{\min }^2\right)}.\end{array}& \mathrm{equation}.Math..Math.\left(60\right)\end{array}$

[0083] Taking the limit as j.fwdarw. of (60) yields a convergent series. Thus, the gradient approaches to zero, indicating that the CCD converges to a stationary point.

Local Convergence to Global Minimum

[0084] If x* is an optimal solution of equation (2), then all the elements of the following set

.sub.c:={e.sup.jx*, [0, 2)}equation (61)

are also optimal solutions of equation (2). The distance of a vector z .sup.N to .sub.c is

[00041]$\begin{array}{cc}\mathrm{dist}\left(z,{}_c\right)=\underset{}{\min }.Math..Math.z-e^{j.Math..Math.}.Math.x^{*}.Math.& \mathrm{equation}.Math..Math.\left(62\right)\end{array}$

and the minimum of equation (62) attains at =(z). Similarly, the set of all optimal solutions of the real-valued problem (3) is defined as

[00042]$\begin{array}{cc}:=\left\{\left[\begin{array}{c}\mathrm{Re}\left(e^{j.Math..Math.}.Math.x^{*}\right)\\ \mathrm{Im}\left(e^{j.Math..Math.}.Math.x^{*}\right)\end{array}\right].Math.\stackrel{}{=}.Math.T_{}\left({\stackrel{\_}{x}}^{*}\right),\left[0,2.Math.\right)\right\}& \mathrm{equation}.Math..Math.\left(63\right)\end{array}$

where x*=[Re(x*).sup.T Im(x*).sup.T].sup.T is a global minimizer of (3). That is, T.sub.(x*) denotes the effect of a phase rotation to x*. The projection of onto is the point in closest to x.sup.k, which is denoted as T.sub..sub.k(x*) where

[00043]$\begin{array}{cc}{}_k=\arg .Math..Math.\underset{}{\min }.Math..Math.{\stackrel{\_}{x}}^k-T_{}\left({\stackrel{\_}{x}}^{*}\right).Math..& \mathrm{equation}.Math..Math.\left(64\right)\end{array}$

[0085] Then the distance of x.sup.k to is

[00044]$\begin{array}{cc}\begin{array}{c}\mathrm{dist}\left({\stackrel{\_}{x}}^k,\right)=.Math.\underset{}{\min }.Math..Math.{\stackrel{\_}{x}}^k-T_{}\left({\stackrel{\_}{x}}^{*}\right).Math.\\ =.Math..Math.{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right).Math..\end{array}& \mathrm{equation}.Math..Math.\left(65\right)\end{array}$

[0086] Proving that dist(x.sup.k,).fwdarw.0 further explains the benefits of embodiments described herein. The following lemma, which essentially states that the gradient of the objective function is well behaved, is crucial to the proof.

[0087] Lemma 3: For any z .sup.N with dist (z, .sub.c), the regularity condition

Re(f(z), ze.sup.j(z)x*) dist(z, .sub.c)+f(z).sup.2, >0, >0 equation (66)

holds with high probability if the number of measurements satisfies MC.sub.0N log N with C.sub.0>0 being a sufficiently large constant.

[0088] Although the regularity condition of Lemma 3 corresponds to the complex-valued case, the real-valued version according to equation (63) and equation (64) is obtained at once. For x.sup.k satisfying dist(x.sup.k, ), we have:

f(x.sup.k), x*) dist(x.sup.k, )+f(x.sup.k).sup.2. equation (67)

[0089] Theorem 2: Assume that the sample size satisfies MC.sub.0N log N with a sufficiently large C.sub.0. The RCD with a slight modification, in which the one-dimensional search is limited to a line segment, that is:

[00045]$\begin{array}{cc}{}_k=\arg .Math..Math.\underset{}{\min }.Math.f\left({\stackrel{\_}{x}}^k+.Math..Math.e_{i_k}\right),s.t..Math..Math..Math.2.Math..Math..Math.f_{i_k}\left({\stackrel{\_}{x}}^k\right).Math.& \mathrm{equation}.Math..Math.\left(68\right)\end{array}$

converges to in expectation with high probability at a geometric rate

[00046]$\begin{array}{cc}\left[{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^{k+1},\right)\right]{\left(1-\frac{.Math..Math.{}_{\max }}{N}\right)}^k.Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^0,\right),.Math.{}_{\max }>0.& \mathrm{equation}.Math..Math.\left(69\right)\end{array}$

[0090] Proof: The updating equation of the CD, i.e., x.sup.k+1=x.sup.k+.sub.ke.sub.i.sub.k, is equivalently expressed as

x.sup.k+1=x.sup.k.sub.kf.sub.i.sub.k(x.sup.k)e.sub.i.sub.k equation (70)

where .sub.k=.sub.k/f.sub.i.sub.k(x.sup.k). It requires >0 to ensure f(x.sup.k+1)<f(x.sup.k). Hence, |.sub.k|2|f.sub.i(x.sup.k) means 0<.sub.k2. Beginning from equation (65), it follows

[00047]$\begin{array}{cc}\begin{array}{c}{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^{k+1},\right)=.Math.{.Math.{\stackrel{\_}{x}}^{k+1}-T_{{}_{k+1}}\left({\stackrel{\_}{x}}^{*}\right).Math.}^2\\ .Math.{.Math.{\stackrel{\_}{x}}^{k+1}-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right).Math.}^2\\ =.Math.{.Math.{\stackrel{\_}{x}}^k-{}_k.Math.f_{i.Math.}\left({\stackrel{\_}{x}}^k\right).Math.e_{i_k}-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right).Math.}^2\\ =.Math.{.Math.{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right).Math.}^2+{{}_k^2\left(f_{i_k.Math.}\left({\stackrel{\_}{x}}^k\right)\right)}^2-\\ .Math.2.Math.{}_k.Math.f_{i_k.Math.}\left({\stackrel{\_}{x}}^k\right).Math.e_{i_k}^T\left({\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right)\right)\\ =.Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^k,\right)+{{}_k^2\left(f_{i_k.Math.}\left({\stackrel{\_}{x}}^k\right)\right)}^2-\\ .Math.2.Math.{}_k.Math.{f_{i_k.Math.}\left({\stackrel{\_}{x}}^k\right)\left[{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right)\right]}_{i_k}.\end{array}& \mathrm{equation}.Math..Math.\left(71\right)\end{array}$

[0091] It is already known that [(f.sub.i.sub.k(x.sup.k)).sup.2]=f(x.sup.k).sup.2/(2N) by (57). This also yields

[00048]$\begin{array}{cc}\left[{f_{i_k}\left({\stackrel{\_}{x}}^k\right)\left[{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right)\right]}_{i_k}\right]=\frac{1}{2.Math.N}.Math.\underset{i=1}{\overset{2.Math.N}{.Math.}}.Math.{f_{i_k}\left({\stackrel{\_}{x}}^k\right)\left[{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right)\right]}_i=\frac{1}{2.Math.N}.Math.f\left({\stackrel{\_}{x}}^k\right),{\stackrel{\_}{x}}^k-T_{{}_k}\left({\stackrel{\_}{x}}^{*}\right)\frac{}{2.Math.N}.Math.\mathrm{dist}\left({\stackrel{\_}{x}}^k,\right)+\frac{}{2.Math.N}.Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2& \mathrm{equation}.Math..Math.\left(72\right)\end{array}$

where the last line follows from equation (67). Combining equation (71) and equation (72) yields

[00049]$\begin{array}{cc}\left[{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^{k+1},\right)\right]\left(1-\frac{.Math..Math.{}_k}{N}\right).Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^k,\right)+\frac{{}_k}{2.Math.N}.Math.\left({}_k-2.Math.\right).Math.{.Math.f\left({\stackrel{\_}{x}}^k\right).Math.}^2\left(1-\frac{.Math..Math.{}_k}{N}\right).Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^k,\right)& \mathrm{equation}.Math..Math.\left(73\right)\end{array}$

where the last inequality follows from 0<.sub.k2. Successively applying equation (58), the following is obtained:

[00050]$\begin{array}{cc}\left[{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^{k+1},\right)\right]\underset{j=1}{\overset{k}{.Math.}}.Math..Math.\left(1-\frac{.Math..Math.{}_j}{N}\right).Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^0,\right){\left(1-\frac{.Math..Math.{}_{\min }}{N}\right)}^k.Math.{\mathrm{dist}}^2\left({\stackrel{\_}{x}}^0,\right),.Math..Math.{}_{\min .Math.}=\underset{1jk}{\min }.Math.{}_j.& \mathrm{equation}.Math..Math.\left(74\right)\end{array}$

[0092] As seen from the foregoing, embodiments of the phase retrieval system enable recovery of a signal using intensity measurements and, in doing so, recover the phase of the signal given only the magnitude-square component. This is done in a fast manner even when the sample number is small. Certain embodiments can also be applied for the recovery of sparse signals. These features are enabled by a novel framework for employing coordinate descent techniques. These techniques include three variants, namely, CCD, RCD, and GCD, which are devised for phase retrieval of general signals-of-interest. At each iteration, the coordinate descent technique only needs to find the minimum of a univariate quartic polynomial, which is achieved by finding the closed-form roots of a cubic polynomial. Thus, the coordinate descent technique is computationally simple and easy to implement.

[0093] Although the phase retrieval problem can be expressed as minimization of a multivariate quartic polynomial, it is still NP-hard. However, by applying coordinate descent, only the minimum of a univariate quartic polynomial at each iteration is determined. This yields a closed-form solution. The described coordinate descent techniques for phase retrieval, namely, CCD, RCD and GCD, are computationally simple and converge faster than WF, which is the state-of-the-art approach.

[0094] To further explain these features, it is theoretically proven that the three (3) coordinate descent techniques converge to a stationary point from any initial value. It is further proved that, according to certain embodiments, the RCD converges to the global minimum and achieves exact retrieval with high probability provided that the sample size is sufficiently large. Also, with the use of l.sub.1-regularization, CCD and RCD are extended for phase retrieval of sparse signals, referred to a l.sub.1-CCD and l.sub.1-RCD, respectively. The l.sub.1-CCD and l.sub.1-RCD outperform previously known techniques, namely, WF, PhaseLift, and PhaseCut, particularly when the number of intensity measurements is small.

[0095] According to other embodiments, CCD, RCD and GCD can be applied to blind equalization of constant-modulus communication signals. These techniques are superior to the benchmarks of WF and SE methods in terms of faster convergence and smaller ISI.

[0096] FIG. 9 shows a block diagram of phase retrieval system 900 for recovering a full signal-of-interest. Phase retrieval system 900 can be used with, or be part of, a system(s) that include, e.g., a telescope, a microscope, or other optical system using phase retrieval techniques and utilized in various applications, including astronomy, computational biology, crystallography, digital communications, electron microscopy, neutron radiography and optical imaging.

[0097] System 900 comprises signal capture block 901, which captures signals, e.g., transmitted or reflected by object 902. Signal capture block 901 receives or otherwise obtains signals or signal data relating to object 902, which can be signals-of-interest which are recovered by system 900.

[0098] Object 902, which may be a celestial star, a transmission antenna, or other signal source, emits or reflects signals to form incoming wave front 903. Wave front 903 may be a spherical wave front, a plane wave front, or the like. Wave front 903 propagates from object 902 to phase retrieval system 900. System 900 iterates a series of operations and has an input and an output. A data set having a plurality of elements, each element containing magnitude and phase information, is received at the input. After completing an iteration, the method outputs a new approximation of the received data set, and this approximation forms the basis for the input to the next iteration. It is intended that each iteration is a better approximation than the last iteration.

[0099] Within system 900, wave front 903 is applied to optics including lens 904, which can be a Fourier transform lens having its focus at screen 905. For instance, according to embodiments in which lens 904 is a Fourier transform lens, lens 904 performs a frequency-space transformation to produce a reconstruction of wave front 903 at the screen 905 in the spatial domain.

[0100] In operation of system 900, phase retrieval engine (PRE) 906 receives a full signal-of-interest from intensity or magnitude information utilizing the iterative coordinate descent techniques described herein. PRE 906 is in commutation with processors of system 900 and other components to enable the functions described herein.

[0101] Display unit 907 displays one or more signals-of-interest based on the coordinate descent techniques utilized by system 900 in a given embodiment. Display unit 907 can be electronically accessible and communicatively connected to other components of phase retrieval system 900. According to some embodiments, each pixel of a display device comprising display unit 907 is electronically accessible to other components of phase retrieval system 900. Receiving data from components of phase retrieval system 900, display unit 907 can display the signals-of-interest at or near real time.

[0102] Display unit 907 and components of phase retrieval system 900 operate in conjunction with one another to present the signals of interest to a user. Display unit 907 can be, or can comprise, a display device(s), such as an SLM display device or an LCoS display device. In other embodiments, display unit 907 can comprise one or more of high-resolution LCDs, 3D television (TV) displays, high-resolution LCoS display devices, high-resolution SLM display devices, or other desired display devices suitable for displaying the signals-of-interest.

[0103] It should be appreciated that signal-of-interest can be communicated over wired or wireless communication channels to display unit 907 or other display components (e.g., remote display components, such as a 3D TV display) to facilitate generation and display of the 3D signals-of-interest.

[0104] Phase retrieval system 900 can comprise additional components that enable system 900 to execute functions described herein. For example, phase retrieval system 900 further comprises communication block 908 that communicates information between components of phase retrieval system 900 (e.g., display component(s), capture device(s), processor component(s), user interface(s), data store(s), etc.). The information includes, for example, a source object, signal-of-interests, information relating to defined signal-of-interest generation criterion(s), information relating to an algorithm(s) that facilitate generating signals-of-interest, etc.

[0105] Phase retrieval system 900 also comprises aggregator 909 that aggregates data received from various entities (e.g., capture device(s), display component(s), processor component(s), user interface(s), data store(s), etc.). Aggregator 909 correlates respective items of data based at least in part on type of data, source of the data, time or date the data was generated or received, point with which data is associated, image with which data is associated, pixel with which a transparency level is associated, etc., to facilitate processing of the data.

[0106] Phase retrieval system 900 also can comprise a processor 910 that can operate in conjunction with the other components to perform the various functions of phase retrieval system 900. Processor 910 can employ one or more processors (e.g., central processing units (CPUs), GPUs, FPGAs, etc.), microprocessors, or controllers that can process data, such as data relating to parameters associated with phase retrieval system 900 and associated components, etc., to perform operations relating to generating signal-of-interests.

[0107] Phase retrieval system 900 can also comprise data store 911 that stores data structures (e.g., user data, metadata); code structure(s) (e.g., modules, classes, procedures), commands, or instructions; information relating to generating a signal-of-interest, parameter data; algorithms; criterion(s) relating to signal-of-interest generation; and so on.

[0108] Processor 910 can be functionally coupled to data store 911 to store and retrieve information desired to operate and/or confer functionality, at least in part, to other system components and/or substantially any other operational aspects of phase retrieval system 900. It should be appreciated that the various components of phase retrieval system 900 can communicate information between each other and/or between other components associated with phase retrieval system 900 as desired to carry out operations of phase retrieval system 900. It should be further appreciated that respective components of phase retrieval system 900 each can be a stand-alone unit, can be included within phase retrieval system 900 (as depicted), can be incorporated within another component of phase retrieval system 900 or a component separate from phase retrieval system 900, and/or a combination thereof.

[0109] Phase retrieval system 900 can further comprise intelligence block 912 that can be associated with (e.g., communicatively connected to) processor 910 and/or other components associated with system 900 to facilitate analyzing data, such as current and/or historical information, and, based at least in part on such information, can make an inference(s) and/or a determination(s) and the like. For example, based in part on current and/or historical evidence, intelligence block 912 can infer or determine a process or algorithm to use.

[0110] Intelligence block 912 communicates information related to the inferences and/or determinations to processor 910. Based on the inference(s) or determination(s) made by intelligence block 912, processor 910 can take one or more actions to facilitate generating a signal-of-interest.

[0111] Intelligence block 912 reasons about or infers states of system 900, its environment, and/or users from a set of observations as captured via events and/or data. Inferences can be employed to identify a specific context or action, or can generate a probability distribution over states, for example. The inference can be probabilisticthat is, the computation of a probability distribution over states of interest based on a consideration of data and events. Inferences can also refer to techniques employed for composing higher-level events from a set of events and/or data. Such inference results in the construction of new events or actions from a set of observed events and/or stored event data (e.g., historical data), whether or not the events are correlated in close temporal proximity, and whether the events and data come from one or several event and data sources. Various classification (explicitly and/or implicitly trained) schemes and/or systems (e.g., support vector machines, neural networks, expert systems, Bayesian belief networks, fuzzy logic, data fusion engines . . . ) can be employed in connection with performing automatic and/or inferred action in connection with the disclosed subject matter.

[0112] System 900 also can comprise presenter 913, which is connected to processor 910. Presenter 913 provides various types of user interfaces to facilitate interaction between a user and any component coupled to processor 910. Presenter 913 can be a separate entity that can be utilized with processor 910 and associated components. However, it is to be appreciated that presenter 913 and/or similar view components can be incorporated into processor 910 and/or can be a stand-alone unit. Presenter 913 provides one or more graphical user interfaces (GUIs) (e.g., touchscreen GUI), command line interfaces, and the like. For example, a GUI can be rendered that provides a user with a region or means to load, import, read, etc., data, and can include a region to present the results of such. These regions can comprise known text and/or graphic regions comprising dialogue boxes, static controls, drop-down-menus, list boxes, pop-up menus, as edit controls, combo boxes, radio buttons, check boxes, push buttons, and graphic boxes. The user can interact with one or more of the components coupled to and/or incorporated into the processor 910.

[0113] The user can also interact with the regions to select and provide information via various devices such as a mouse, a roller ball, a keypad, a keyboard, a touchscreen, a pen and/or voice activation, for example. Typically, a mechanism such as a push button or the enter key on the keyboard can be employed subsequent entering the information in order to initiate the search. However, it is to be appreciated that the claimed subject matter is not so limited. For example, merely highlighting a check box can initiate information conveyance. In another example, a command line interface can be employed. For example, the command line interface can prompt (e.g., via a text message on a display and an audio tone) the user for information via providing a text message. The user can than provide suitable information, such as alpha-numeric input corresponding to an option provided in the interface prompt or an answer to a question posed in the prompt. It is to be appreciated that the command line interface can be employed in connection with a GUI and/or API, in addition, the command line interface can be employed in connection with hardware (e.g., video cards) and/or displays (e.g., black and white, and EGA) with limited graphic support, and/or low bandwidth communication channels.

[0114] In accordance with one embodiment of the disclosed subject matter, the processor 910 and/or other components, can be situated or implemented on a single integrated-circuit chip. In accordance with another embodiment, processor 910, and/or other components, can be implemented on an application-specific integrated-circuit (ASIC) chip. In yet another embodiment, processor 910 and/or other components, can be situated or implemented on multiple dies or chips.

[0115] FIG. 10 shows a flow diagram of an exemplary method performed according to the concepts described herein. Method 1000 is performed to efficiently retrieve phase information for a signal-of-interest and can be executed by a phase retrieval system, such as system 900 shown at FIG. 9.

[0116] At step 1001, magnitude information for a signal-of-interest is expressed as an objective function, f(x). As previously discussed, the objective function, f(x), can be expressed in the form of a multivariate quartic polynomial.

[0117] At step 1002, a coordinate index, i.sub.k, is selected according to one or more predefined coordinate descent rules. The coordinate descent rules can be one or more of a cyclic rule, a random rule, or a greedy rule. In one embodiment, the coordinate index, i.sub.k, is cyclically assigned sequential values from 1 to 2N for each of the variable components of the multivariate quartic polynomial. In another embodiment, the coordinate index, i.sub.k, is randomly assigned values from 1 to 2N, where each value is assigned with equal probability, for each of the variable components of the multivariate quartic polynomial. In yet another embodiment, the coordinate index, i.sub.k, is assigned a value of:

[00051]$\begin{array}{cc}{i}_k=\arg .Math..Math.\underset{i}{\min }.Math..Math.f_i\left({\stackrel{\_}{x}}^k\right).Math.,.Math.f_i\left({\stackrel{\_}{x}}^k\right)=\frac{f\left({\stackrel{\_}{x}}^k\right)}{{\stackrel{\_}{x}}_i^k}.& \mathrm{equation}.Math..Math.\left(75\right)\end{array}$

[0118] At step 1003, one variable component of the multivariate quartic polynomial is minimized while every other variable component of the multivariate quartic polynomial is held at a constant value. Preferably, this is done on an iterative basis where, at each iteration, a different component is minimized while the other components are held at constant values. According to one embodiment, for the kth iteration, this is done by minimizing f(x.sup.k) with respect to the i.sub.kth (i.sub.k {1, . . . , 2N}) variable while maintaining the remaining 2N1 {x.sub.i.sup.k}.sub.ii.sub.k variables at a constant value. Also, step 1003 can involve solving only one variable component while maintain every other variable component at a constant value.

[0119] The aforementioned systems and/or devices have been described with respect to interaction between several components. It should be appreciated that such systems and components can include those components or sub-components specified therein, some of the specified components or sub-components, and/or additional components. Sub-components could also be implemented as components communicatively coupled to other components rather than included within parent components. Further yet, one or more components and/or sub-components may be combined into a single component providing aggregate functionality. The components may also interact with one or more other components not specifically described herein for the sake of brevity, but known by those of skill in the art.

[0120] Although the present invention and its advantages have been described in detail, it should be understood that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims. Moreover, the scope of the present application is not intended to be limited to the particular embodiments of the process, machine, manufacture, composition of matter, means, methods and steps described in the specification. As one of ordinary skill in the art will readily appreciate from the disclosure of the present invention, processes, machines, manufacture, compositions of matter, means, methods, or steps, presently existing or later to be developed that perform substantially the same function or achieve substantially the same result as the corresponding embodiments described herein may be utilized according to the present invention. Accordingly, the appended claims are intended to include within their scope such processes, machines, manufacture, compositions of matter, means, methods, or steps.