Maximum Likelihood DOA Estimation Based on
the Cross-Entropy Method
Yen-Chih Chen and Yu T. Su
Department of Communications EngineeringNational Chiao Tung University
1001 Ta Hsueh Road, Hsinchu 30056, TAIWAN joechen.cm88g@nctu.edu.tw, ytsu@mail.nctu.edu.tw
Abstract— In this paper, we propose two simulation based maximum likelihood (ML) methods to estimate the direction of arrival (DOA) by a novel combination of the Cross-Entropy (CE) method and the polynomial parameterization scheme. The CE method is an efficient stochastic approximation method for solving both discrete and continuous optimization problems. We bridge the ML approach and the stochastic search algorithm by properly randomizing the desired parameters. Numerical results show that the proposed CE-based algorithms yield highly accurate DOA estimation with fast convergence rate while requiring only linear processing complexity. Compared with the conventional iterative quadratic maximization likelihood (IQML) method, the proposed algorithms can alleviate the error propagation effect in low signal to noise ratio (SNR) region and asymptotically approach the Cram´er-Rao bound in high SNR region.
I. INTRODUCTION
The problem of DOA estimation has been extensively investigated in the last few decades [1]. Among various DOA estimation approaches, the class of maximum likelihood (ML) methods is both unbiased and asymptotically efficient. Moreover, the ML methods can alleviate the mutual can-cellation effect caused by the coherent sources and provide high estimation accuracy. However, the multi-dimensional optimization nature of the ML methods also require high computing complexity. Bresler and Macovski [2] suggested a polynomial parameterization approach for conditional ML estimation and applied iterative projection to reduce the computation burden. The resulting ML algorithm, abbreviated IQML, is an efficient and near-optimal solution. However, IQML suffers from slow convergence rate in low SNR region and it does not give satisfactory performance when the sample size is small.
In this paper, we suggest two novel simulated-based opti-mization approach to solve the ML DOA estimation problem. The Cross-Entropy (CE) method developed by Rubinstein [3] was originally aimed to solve the rare event simulation problem. It was later proposed as an unified approach to perform global search for combinatorial and continuous multi-extremal optimization problems. Under the CE method frame-work, deterministic signals or parameters, e.g. DOA in array signal processing, are randomized according to some
para-metric distributions related with a set of dynamic parameters. Through iteratively minimizing the Kullback-Leibler distance (cross entropy) between the the parametric distribution and the target distribution, we can update the set of parameters and improved the estimate quality. Thanks to the accurate and robust nature of the CE method [3], we are able to obtain highly accurate ML DOA estimation and eliminate the error propagation effect incurred by the iterative projection in the conventional IQML algorithm.
To solve the complexity issue associated with the simulated searching approach, we approximate source information by sample covariance matrix and adopt the polynomial parame-terization procedure used in the IQML algorithm to translate the original problem into one that is a computational much more efficient. In short, this paper presents a new estimate method that not only links the conventional conditional ML and the simulated approximation approaches but also suggest a novel bridging procedure so that the advantages of both approaches are preserved.
The organization of this paper is as follows. In Section II, we describe the system model of array signal processing and review the conventional IQML method. In Section III, we describe the generic concept of the CE method. In Section IV, we present the proposed maximum likelihood estimation based on the CE method. In Section V, a polynomial param-eterizing CE method is proposed to lower the computation complexity. In Section VI, we provide some simulation ex-amples and discuss the numerical behavior. Finally, in the last Section, we give summary remarks on the proposed new methods and suggest some future works.
II. SYSTEMMODEL ANDIQMLMETHOD A. System Model
Assume thatD narrowband plane waves imping on an uni-form linear array (ULA) equipped withM antenna elements. At time instantn, the multi-dimensional received vector y(n) of size M × 1, i.e., y(n) = [y1(n), y2(n), . . . , yM(n)]T, can
be modelled as
wherex(n) = [x1(n), x2(n), . . . , xD(n)]T is aD × 1 source
vector,T being the transpose operator, and w(n) is an M ×1 noise vector with independent zero-mean Gaussian random components and covariance matrix Kw = σw2I. Matrix A
denotes the array manifold matrix with D array manifold vectorsa(θi), i = 1, 2, · · · , D,
A = [a(θ1), a(θ2), · · · , a(θD)] , (2)
where
a(θi) =
1, ejπsin(θi), · · · , ejπ(M −1)sin(θi)T (3) and θi is the DOA associated with the ith impinging plane
wave.
Stacking upNs consecutive independent received vectors
and defining Y def= [y(1), · · · , y(Ns)], X def= [x(1), · · · ,
x(Ns)], and W def
= [w(1), · · · , w(Ns)], we obtain
Y = AX + W, (4)
As W is Gaussian distributed, the conditional ML function
can be written as f(Y|X, A) = N s n=1 f(y(n)|x(n), A) = 1 πM Ns|K w|Nse −Nsn=1y(n)−Ax(n)2K−1 w , (5) where x2 K−1def= xHK−1x (6) denotes the generalized 2-norm of x, and H denotes the
conjugate transpose. From (5), the conditional ML estimate can be written as max {θi,x} − Ns n=1 y(n) − Ax(n)2 2 . (7)
The solution, though can be found by exhaustic search, is computational intractable, especially when the array size M and data block Ns are large. An alternative scheme is to
substitute the source vectorx(n) by its least squared solution
for any fixed A, and transform (7) into
min {θi} Ns n=1 y(n) − AA†y(n)2 2 = min {θi} Ns n=1 yH(n)P⊥ Ay(n) (8) = min {θi}Ns trP⊥ ARˆy , (9)
where† denotes the matrix pseudoinverse operation, ⊥ the or-thogonal projection,P⊥
Adef= I − AA†, and ˆRydef= N1sYYH
is the sample correlation matrix.
B. IQML method
Using the coefficients of the polynomial b(z)def= D
i=1
(z − zi) = b∗0+ b∗1z + . . . + b∗DzD (10)
where {zi def= ejπsin(θi)} are its roots and ∗ represents
complex conjugate, Bresler and Macowski [2] construct the M × (M − D) Toeplitz matrix B B = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ b0 .. . b0 bD ... ... bD ... b0 ... ... bD ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (11)
They further proved thatP⊥
A= PB, and (8) is equivalent to min b b H Ns n=1 YT(n)(BTB∗)−1Y∗(n) b def= min b b HQb (12) subject to: b0= 1, (13)
where b = [b0, b1, . . . , bD]T and Y(n) is the accumulated
data matrix Y(n)def= ⎡ ⎢ ⎣ yD+1(n) · · · y2(n) y1(n) .. . ... ... ... yM(n) · · · yM −D+1(n) yM −D(n) ⎤ ⎥ ⎦ . The constraint (13) is imposed to guarantee a non-trivial b.
Instead of optimizing (12) directly, [2] searches the solution by iteratively projecting on the temporary solution space.
III. CROSS-ENTROPYMETHOD
For convenience of reference, we give a brief review of Rubinstein’s CE method. Let Θ be the domain of variable θ, and S be the score function of θ defined on Θ. The CE method attempts to solve the following optimization function
arg max
θ∈ΘS(θ) (14)
by relating the above problem to a rare event simulation problem. A popular approach to reduce the required sample size for estimating the rare event probability is the so-called importance sampling technique. To find the optimal impor-tance sampling density within a class of densitiesf(θ; v), one iteratively adapts the parameterv so that the Kullback-Leibler
distance (cross entropy) between the associated density and the optimal importance sampling density is minimized. In short, a generic CE method can be described by the following two steps [3].
1. Generate samples from the importance density specified by the parameters from the previous iteration.
2. Update the parameters for next iteration according to the order of the score values associated with new samples and the CE criterion.
IV. ML DOAESTIMATION BASED ONCEMETHOD The score function of (7) is the total sum of the negative 2-norm of the error vectors. To apply the CE method, we have to include a set of important density functions over the deterministic A and X. At the kth iteration, for BPSK
signals, a natural choice of density function for each element {{x(k)i (n)}Di=1}Nn=1s inX(k)is the Bernoulli distribution, i.e.,
x(k)i (n) ∼ Ber(p(k)i,n), such that P (x(k)i (n) = 1) = 1 − P (x(k)i (n) = −1) = p(k)i,n.
Besides, we have to model the D-dimensional DOA in-formation in the continuous space [−π, π) for the manifold matrix A(k). To be compatible with the results for spatially
distributed source [4][6], each element{sin(θi(k))}D
i=1is
mod-elled as truncated Gaussian distributed, i.e., sin(θ(k)i ) ∼ N (µ
(k) i , σ
(k)
i ) (15)
given the constraint
−1 ≤ sin(θi(k)) ≤ 1, for i = 1, 2, . . . , D. (16)
In order to sample {sin(θ(k)i )}Di=1 in the constraint region,
we adopt the acceptance-rejection method described in [7]. By defining the approximated important density, the CE-ML algorithm for DOA estimation is described in Table I in which the following three recursive relations are used.
p(k+1)i,n = α · 1 U u=1I{S(X(k) u ,A(k)u )≥γ(k)} · U u=1 I{S(X(k) u ,A(k)u )≥γ(k)}I{x(k)i,u(n)=1} + (1 − α)p(k)i,n (17) µ(k+1)i = α · 1 U u=1I{S(X(k) u ,A(k)u )≥γ(k)} · U u=1 I{S(X(k) u ,A(k)u )≥γ(k)}sin(θ (k) i,u) + (1 − α)µ(k)i (18) (σi(k+1))2= α · 1 U u=1I{S(X(k) u ,A(k)u )≥γ(k)}− 1 · U u=1 I{S(X(k) u ,A(k)u )≥γ(k)} sin(θi,u(k)) − µ (k+1) i 2 + (1 − α)(σ(k)i )2 (19)
While the algorithm iterates, the generated standard deviation σ(k)i decreases to zero. When the algorithm converges, the
approximated important density degenerates to the impulse function, i.e.N (θ∗
i, 0), matching the important density of the
deterministic case in (7).
By sampling from the exact signal space, we obtain perfor-mance superior to that of the conventional covariance matrix based methods, which are corrupted by unknown signal phase. However, when the data block is large, sampling over such a high dimensional signal space ofX(n) becomes impractical.
V. CE-PPMLMETHOD
To ease the computation burden of the pure CE method given in section IV, we propose an alternative CE method that does not generate samples from high dimensional X(n). We
use the information generated from least square projection of the sample covariance matrix ˆRy in (9). This approximation
greatly reduces the sampling complexity especially when the block size,Ns, is much larger than the number of impinging
plane waves, D. Furthermore, we invoke the technique of polynomial parameterization described in Section II-B to simplify the calculation ofPB = P⊥Ain (9). Since the matrix
(BHB) is a banded Hermitian matrix, its inverse can be
efficiently computed via the inverse of aD × D matrix rather than an(M −D)×(M −D) matrix [8]. Although the required inverse matrix size D × D is the same as P⊥
A, we can avoid
calculating the cumbersome(M − 1)D exponential functions in A after applying polynomial parameterization, especially
when the array sizeM and the sample size U are large. The resulting CE-PPML algorithm is given in Table II.
As the DOA information of the polynomial parameteri-zation scheme is represented by b, the proposed CE-PPML
algorithm first generates a list of samples {bu}Uu=1 from an
auxiliary phase distribution. Following the basic framework of the CE method, we then replace the iterative projecting steps of the IQML algorithm by the iterative important sampling procedure in Table I. As will be proved by simulation, the divergence phenomenon of the conventional IQML method is greatly reduced by this method especially in low SNR region. The main ideas of the proposed CE algorithm can be sum-marized by the two design criteria: (i) minimizing the number of target parameters by using deterministic signal processing techniques with performance and complexity constraints; and (ii) randomizing the desired parameters and applying the CE method. By following the criteria we are able to take the advantage of the global searching ability of the CE method while maintaining linear computing complexity.
VI. SIMULATIONRESULTS
We consider a standard ULA equipped with M = 10 antenna elements. Two adjacent elements of the ULA are separated by half wavelength. A single signal source with impinging arrival angle equal to sin(θ) = −0.3 is assumed. The block number Ns is 10. The sample size used in
CE-based methods is 500 for all SNR regions and, except for Fig. 3, the smooth parameter α is 0.6. Maximum iteration
number is set to be 200. When the iteration number exceeds 200, the algorithm is regarded as being fail to converge.
Fig. 1 shows the root mean squared error (RMSE) per-formance of the IQML and the proposed (ML, and CE-PPML) algorithms using a short block size. The stopping cri-terion used in Table I and Table II is based on the convergence of the score value in each iteration. The algorithms terminate when the difference of the score values of five successive iterations are within a predefined threshold, sayη = 10−8. It
is clearly that the two proposed CE-based methods outperform the IQML algorithm and the CE-ML algorithm yields the best overall average performance. The Cram´er-Rao bounds for array systems with unknown complex and real signal source using the conditional ML approach [5], [9]
CRBcomplex(ψ) = 6 · SNR −1 Ns(M2− 1)Mπ2, ψ = sin θ (20) CRBreal(ψ) = 3 · SNR −1 Ns(M − 1)M(2M − 1)π2 (21) are also included in the figure for comparison purpose.
Since signal sources are estimated by the sample correlation matrix ˆRy, the IQML and the proposed CE-PPML methods
can only achieve the bound for unknown phase (complex) signal sources. On the other hand, the CE-ML method, having drawn samples from the signal domain thus avoiding possible phase error, can achieve the bound for known phase (real) signal sources. This performance advantage is obtained at the expense of increasing the required sample size for large data blocks. In other words, when Ns is small, since the
sample covariance matrix ˆRy is not good approximation of
its ensemble counterpart, the proposed CE-ML method gives better performance than the covariance matrix based methods. The large block size (Ns= 100) RMSE performance of the
CE-PPML and IQML algorithms are depicted in Fig. 2. Other algorithmic parameters used are the same as the previous example. By using a larger data block size, the sample covariance matrix ˆRybecomes a closer approximation of the
ensemble covariance matrixRyand the performance of both
algorithms improves accordingly. The CE-PPML algorithm still outperforms the IQML method for −10 ≤ SNR ≤ 10 dB. The convergence behavior comparison is given in Fig. 3. Fig. 3 indicates that the average iteration numbers for the two proposed CE methods remain relatively low while the IQML solution tends to diverge in the low SNR region.
The computational complexity of the CE-based algorithm is proportional to the product of iteration number and sample size while that of the IQML algorithm is proportional to the iteration number. Nevertheless, the required iteration number of the latter is much larger than that of the former algorithms in low SNR region. Furthermore, samples used in CE-based algorithms can be generated in parallel and then the process-ing delay will mainly depend on the iteration number.
VII. CONCLUSION
We present two CE-based methods to perform the maxi-mum likelihood DOA estimation. Our first CE-based approach searches the ML solution in the full signal and angle space. It gives accurate DOA estimation even if the block sizeNsis
small. However, when Ns is large, the algorithm is
compu-tational infeasible. To make the CE-based solution practical, we introduce the covariance matrix approximation and the polynomial parameterization scheme in our algorithm. The resulting CE-PPML algorithm alleviates the error propagation arose from the iterative projection process in the IQML approach and is much more robust in the low SNR region.
Although we consider only the case when the impinging plane waves have deterministic DOA, the proposed method can be easily extended to solve DOA estimation problems for spatially distributed source with arbitrary constellation by including additional system parameters, e.g., angular spread.
REFERENCES
[1] H. L. Van Trees, Optimum array processing, Wiley, New York, 2002. [2] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter
estimation of superimposed exponential signals in noise,” IEEE Trans.
ASSP., pp. 1081-1089, Oct. 1986.
[3] R. Y. Rubinstein and D. P. Kroese, The Cross-Entropy Method, Springer, 2004.
[4] O. Besson and P. Stoica, “Decoupled Estimation of DOA and Angular Spread for a spatially Distributed Source,” IEEE Trans. Signal Process., July 2000.
[5] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cram´er-Rao bound,” IEEE Trans. ASSP., pp. 720-741, May 1989.
[6] J. B. Anderson, “Distributions of phase derivative in mobile communi-cations,” IEE Proc.-H, pp. 197-201, Aug. 1990.
[7] G. S. Fishman, Monte Carlo : Concepts, Algorithms, and Applications, Springer 1995.
[8] R. Kumaresan, et al., “An algorithm for pole-zero modelling and spectral analysis,” IEEE Trans. ASSP., pp. 637-640, June 1986.
[9] J. Li and R. T. Compton, “Maximum Likelihood Angle Estimation for Signals with Known Waveforms,” IEEE Trans. Signal Process., pp. 2850-2862, July 1993. −15 −10 −5 0 5 10 15 20 −30 −25 −20 −15 −10 −5 0 SNR (dB)
Root Mean Squared Error (dB)
IQML CE−PPML CE−ML CRB (complex) CRB (real)
Fig. 1. Root Mean squared error (MSE) versus SNR:M = 10, Ns = 10, D = 1, α = 0.6, Trials = 500. The DOA is set at sin(θ) = −0.3.
−15 −10 −5 0 5 10 15 20 −35 −30 −25 −20 −15 −10 −5 0 SNR (dB)
Root Mean Squared Error (dB)
IQML CE−PPML CRB (complex)
Fig. 2. Root Mean squared error (MSE) performance versus SNR:M = 10, Ns = 100, D = 1, α = 0.6, Trials = 1000. The true DOA is at sin(θ) = −0.3. −150 −10 −5 0 5 10 15 20 20 40 60 80 100 120 140 SNR (dB)
Average Iteration Ratio
IQML CE−PPML CE−ML
Fig. 3. Average Iteration Ratio (AIR) behavior as a function SNR:M = 10, Ns = 10, D = 1, α = 0.9, Trial = 1000. AIR is defined as
total number of iterations
total number of trials . Hence, a higher AIR implies a lower convergence speed.
TABLE I
CEMETHOD FOR MAXIMUM LIKELIHOODDOAESTIMATION 1. Initializep(0)i,n= 0.5, µ(0)i = 0, and σ(0)i = 1. Set the
quantile parameterρ, and define a smooth parameter 0 < α < 1.
2. Draw U samples x(k)i,u(n) ∼ Ber(p (k)
i,n). Draw U
samples sin(θ(k)i,u) ∼ N (µ(k)i,u, σi(k)) by rejecting the samples that are outside the range[−1, 1]. Construct {X(k)u }Uu=1, and{A
(k)
u }Uu=1from the generated
sam-ples.
3. Calculate the score functions {S(X(k)u , A(k)u )}
ac-cording to (7).
4. Set a quantile parameter ρ, such that there is a γ(k)
satisfying γ(k) = arg max
γP (S(Z) ≥ γ) ≥ ρ for
Z ∈ {X(k)u , A(k)u }Uu=1.
5. Update the estimated Bernoulli parameterp(k+1)i,n for x(k+1)i (n) based on (17).
6. Stop at iteration k = K if the pre-defined stopping criterion is met. Output{µ(K+1)i }D
i=1as the estimate
of {sin(θi)}Di=1. Otherwise, let k = k + 1, update
µ(k+1)i ,(σi(k+1))2 using (18), (19) and repeat Steps
2 − 4.
TABLE II CE-PPMLALGORITHM
1. Initialize µ(0)i = 0, and σ(0)i = 1. Set the quantile parameter ρ, and a smooth parameter 0 < α < 1. 2. Draw U samples sin(θi,u(k)) ∼ N (µ(k)i,u, σi(k)). Reject
samples outside the range [−1, 1]. Find b(k)u and
construct the corresponding B(k)u described in (11)
from the generated samples. 3. Calculate the score function
S(b(k) u ) def = Nstr P(k)BuRˆy U u=1 where P(k)Bu def = B(k)u (B(k)u H B(k)u )−1B(k)u H . 4. Set a quantile parameter ρ, such that there is a γ(k)
satisfying
γ(k)= arg max
γ P (S(Z) ≥ γ) ≥ ρ
for Z ∈ {b(k)u }Uu=1.
5. Generate the estimated Gaussian parameter µ(k+1)i and σ(k+1)i for next iteration according to (18) and (19) with{X(k)u , A(k)u } being replaced by {b(k)u }.
6. Stop at iteration k = K if the pre-defined stopping criterion is met. Output{µ(K+1)i }Di=1as the estimate
of{sin(θi)}Di=1. Otherwise, setk = k+1, and repeat
