Efficient AM Algorithms for Stochastic ML Estimation of DOA

The estimation of direction-of-arrival (DOA) of signals is a basic and important problem in sensor array signal processing. To solve this problem, many algorithms have been proposed, among which the Stochastic Maximum Likelihood (SML) is one of the most concerned algorithms because of its high accuracy of DOA. However, the estimation of SML generally involves the multidimensional nonlinear optimization problem. As a result, its computational complexity is rather high. This paper addresses the issue of reducing computational complexity of SML estimation ofDOAbased on theAlternatingMinimization (AM) algorithm. We have the following two contributions. First using transformation of matrix and properties of spatial projection, we propose an efficient AM (EAM) algorithm by dividing the SML criterion into two components. One depends on a single variable parameter while the other does not. Second when the array is a uniform linear array, we get the irreducible form of the EAM criterion (IAM) using polynomial forms. Simulation results show that both EAM and IAM can reduce the computational complexity of SML estimation greatly, while IAM is the best. Another advantage of IAM is that this algorithm can avoid the numerical instability problem which may happen in AM and EAM algorithms when more than one parameter converges to an identical value.


Introduction
The localization of multiple signal sources by a passive sensor array is of great importance in a wide variety of fields, such as radar, geophysics, radio-astronomy, biomedical engineering, communications, and underwater acoustics.The basic problem in this context is to estimate direction-of-arrival (DOA) of narrow-band signal sources located in the far field of the array.
Among these techniques, the ML technique is well known for its high accuracy of DOA.There are two famous ML criterions, that is, Deterministic or Conditional Maximum Likelihood (DML) [1,[3][4][5]7] and Stochastic or Unconditional ML (SML) [2,[5][6][7].The difference between them lies in their signal models.In particular, the SML shows much higher resolution than other techniques.Furthermore, the SML technique can solve coherent signals without any preprocessing, such as the spatial smoothing [9].However, the SML technique does not become popular in practice because the estimation of SML generally involves the multidimensional nonlinear optimization problem, and it requires high computational complexity.
To solve the multidimensional nonlinear optimization problem, the Alternating Minimization (AM) algorithm is one of the most classic techniques.It is an iterative technique and usually needs one-dimensional global search in the updating process.There are also some other efficient techniques, such as Alternating Projection (AP) [3], Expectation Maximization (EM) [13], Space Alternating Generalized EM (SAGE) [14], genetic algorithm [15], ant colony algorithm [16], and Particle Swarm Optimization (PSO) [17,18].Generally, they are all iterative techniques and usually are formulated for explicit criterions (e.g., these algorithms formulated for the criterions of DML and WSF).As for the SML 2 International Journal of Antennas and Propagation criterion, the AM algorithm is still the most commonly used algorithm although its computational complexity is a little high [6].
Therefore, this paper addresses the issue of reducing computational complexity of SML estimation of DOA based on the conventional AM algorithm.
Firstly we show a brief description of SML estimation of DOA and the conventional solving method, AM algorithm.Then using transformation of matrix and properties of spatial projection, we propose an efficient AM (EAM) algorithm by dividing the AM criterion of SML into two components.One depends on a single variable parameter while the other does not.The computational complexity of EAM can be greatly reduced compared to AM algorithm.However, numerical instability may happen in calculation of EAM criterion when more than one parameter converges to an identical value.As a result, oscillation may happen in the convergence phase and extra calculation is needed.To solve this problem and reduce computational complexity further, based on the EAM criterion we get the irreducible form of the EAM criterion (IAM) using a uniform linear array.In this way, the EAM criterion can be written into polynomial forms.The common zeros can be easily canceled in numerator and denominator of the EAM criterion.Furthermore, the computational complexity is also reduced.Finally, simulation results are also shown to demonstrate the validity of the proposed EAM and IAM algorithms.
The rest of this paper is organized as follows.In Section 2 we introduce the problem formulation of DOA.A brief introduction of exact definition of SML estimation for incoherent signals is shown in Section 3. In Section 4, we introduce the conventional AM algorithm and our proposed two efficient algorithms, that is, EAM and IAM algorithms.Simulation results are shown in Section 5 and conclusion is drawn in Section 6.

Problem Formulation
Without loss of generality, consider that there are  sensors and  narrow-band sources far from the array, centered around a known frequency, impinging on the sensor array from distinct directions  1 ,  2 , . . .,   , with respect to a reference point, respectively.
Note that the received signals may be coherent because of multipath propagation.In the case where there are signals coherent, the independent signal number is less than .The task of DOA estimation is to detect all  directions.Also note that here we assume that the signals are narrow-band.For wideband signals, the CSM algorithms [19] can be used as a preprocessing technique to change them into the narrowband.
Furthermore, the sensor configuration can be arbitrary and we assume that all the sensors are omnidirectional and not coupled [6,7].
Using complex envelope representation, the -dimensional vector received by the array can be expressed as where   () is the th signal received at a certain reference point.n() is a -dimensional noise vector.a() is the "steering vector" of the array towards direction , which is represented as where   () is the amplitude response of the th sensor to a wave-front impinging from the direction .  () is the propagation delay between the th sensor and the reference point.The superscript  denotes the transpose of a matrix.
In the matrix notation, (1) can be rewritten as Suppose that the received vector x() is sampled at  time instants  1 ,  2 , . . .,   and define the matrix of the sampled data as where The problem of DOA finding is to be stated as follows.Given the sampled data X, obtain a set of estimated directions of  1 ,  2 , . . .,   .

SML Estimation
In this section, brief descriptions of the exact SML criterion are shown [6,7].
To solve the problem of ML estimation of DOA, we make the following assumptions.
(A1) The array configuration is known and any  steering vectors for different  directions are linearly independent; that is, the matrix A(Θ) has full rank.(A2) , , and  satisfy the condition that a unique solution of DOA exists in the noise-free case.As for a uniform linear array,  ≤ 2/(2 + 1) and  ≥  are the necessary conditions of the uniqueness [20], where  is the rank of S. (A3)  and  are known.(A4) The noise samples n(  ) are statistically independent Gaussian random vectors with zero mean and the covariance matrix  2 I  , where I  is a  ×  identity matrix.
International Journal of Antennas and Propagation 3 (A5) s(  ) are independent samples from a complex Gaussian random vector which has zero mean and signal covariance matrix with rank  and is independent of the noise.
According to [6], the SML criterion can be derived based on the process of the array covariance matrix R defined as where S R = {SS  } is the signal covariance matrix.
The SML criterion is shown as follows [6,7]. where is the sample covariance matrix of the sampled data.V S (Θ) is a  ×  matrix composed of an orthonormal system of the signal subspace spanned by A(Θ).V N (Θ) is a ×(−) matrix composed of an orthonormal system of the noise subspace, which is an orthogonal complement of the signal subspace.
The  ×  matrix, R SS , corresponds to the covariance matrix of the components for x() in the signal subspace.R NN is the covariance matrix of the components for x() in the noise subspace.
From ( 8) and ( 9), we can see that the estimation of SML is to find a set of Θ which minimizes  SML in (9).This is a multidimensional nonlinear optimization problem.
Note that there are literatures [21,22] discussing the preciseness of SML criterion derived above.This paper focuses on the problem of reducing the computational complexity of SML estimation.Our efficient algorithms are derived based on (9).

Efficient AM Algorithms for SML Estimation
The AM algorithm is the most classic estimation algorithm for multidimensional nonlinear optimization problem in DOA estimation.In this section, we will introduce the conventional AM algorithm firstly, and then we will introduce our proposed efficient AM algorithms.

Conventional AM Algorithm.
The AM method is a popular iterative technique for solving a nonlinear multivariate minimization problem with a multimodal criterion [6].It can be applied to the SML criterion of DOA in the following manner.
Let   ( Θ() ) be a cost function of ( 9) for which the signal number is assumed to be  instead of , where Θ() = { θ1 , θ2 , . . ., θ }.Convergence Phase.Repeat the following updating process until all parameters are converged.At each updating process, let one parameter, say θ , be variable and let all other parameters be held fixed.Find θ minimizing the criterion   ( Θ() ) by one-dimensional global search with respect to θ .Change the index  of the parameter to be updated into

Initialization
Although a global minimum is not guaranteed in the AM algorithm, global solutions can be obtained in most cases because of one-dimensional global searches performed in each update process.

Efficient AM (EAM) Algorithm for SML Estimation.
In this section, we propose an efficient version of the AM algorithm.We call it the EAM algorithm.Using transformation of matrix and properties of spatial projection, the EAM algorithm divides the SML criterion into two components.One depends on a variable parameter and the other one is independent of the variable parameter.In order to simplify expressions, parameters to be estimated are represented without the accent hat, and the argument Θ or Θ is omitted.
In each updating process, let   be a variable parameter and define V S l is an orthonormal system of the subspace spanned by {A  }.
V N l is an orthonormal system of the orthogonal complement of the subspace spanned by {A  }.
Note that when the value of   changes, only v  (  ) varies and all others above are fixed.From these definitions, we have span where ⊕ represents the direct sum of subspaces.It follows from ( 19) that there exist a ( −  + 1) × ( −  + 1) unitary matrix T 1 and a  ×  unitary matrix T 2 which satisfy Substituting ( 20) into ( 17), we have Taking the trace of both sides in (22), we have Substituting the definition of P ⊥ A l and v  (  ) into (23), we have Substituting ( 21) into (11), we have Taking the determinant of both sides in (25) and using the definition of ( 16), we have Substitute the definition of P ⊥   and v  (  ) into (26), and define and then we have and then the final form of the proposed algorithm, EAM algorithm, can be derived from ( 24) and (28) as follows: In ( 29), (30), and (31), all quantities except for u(  ) are fixed and can be computed before starting the onedimensional search with respect to   .Therefore main computations of each step in the one-dimensional search are a product of the matrix V  N l and the vector a(  ) in (29) and evaluation of two Hermitian forms in (30) and (31).Evaluation of u(  ) requires ( × ( −  + 1)) multiplications.Then evaluation of ‖u(  )‖ 2 requires (( −  + 1) 2 ) multiplications.Therefore, at each step of one-dimensional search with respect to the variable parameter,   , ( × ( − +1)) multiplications are required.Therefore, computational complexity can be greatly reduced since, in the conventional AM algorithm, we need to calculate the whole SML criterion in each step of one-dimensional search.
(1) Oscillation Problem in EAM Algorithm.Define the subspace The linear combination a (  ) − a (  ) always belongs to U(Θ).Even when   →   , we have dim{U(Θ)|   →  } = , where dim{} represents the dimension of the subspace.On the other hand, when   =   , we have dim{U(Θ)|   =  } < .This implies that the criterion has discontinuous points in the parameter space.
Next, we consider calculation of the EAM criterion v  (  ) in (18).When the variable parameter   approaches the value of a fixed parameter   , v  (  ) vanishes and both the numerator and denominator in (18) become zero.Then v  (  ) becomes indefinite.Therefore the calculation of v  (  ) becomes numerically instable.This can be verified in Section 5.
In the case that the DOA can be solved, the numerical instability does not occur, since each parameter in the convergence phase of the EAM criterion comes apart from others.However, at the threshold region, when more than one signal approaches an identical value, the numerical instability becomes significant.In practice, when this case happens, the sequence of DOA obtained in the convergence phase of the EAM criterion shows oscillation that is because the estimated directions can not converge well due to the numerical instability and would oscillate around that identical value.Let us give an example.In the simulation there are two sources located in 0 and 8 degrees (the true DOAs).When SNR = 0 dB and with specific noise samples, the solutions of SML criterion are 3.999 and 4.001 degrees (the solutions of SML).However, in this case, when calculating  1 =  2 = 4 degrees, the oscillation happens.As a result, the estimated DOAs may be 4 and 4 degrees (the estimated DOAs).In this case, because of oscillation, extra computation is required.Furthermore, the estimated DOAs are wrong (because they are not the solutions of SML criterion) although they are very close to the solutions of SML criterion.As a result, the estimation accuracy is also affected.These explanations will be shown in Section 5.
To solve these problems, the key point is to cancel the common zeros in the numerator and denominator of the EAM criterion.Next we try to establish the irreducible form of the EAM criterion using a uniform linear array.

Irreducible Form of Efficient AM Criterion (IAM).
In this section, we derive the irreducible form of the EAM criterion using a uniform linear array.With a uniform linear array, the EAM criterion can be written into polynomial forms.Then we can easily cancel the common zero in both numerator and denominator.Thus numerical instability never happens in IAM criterion.Furthermore, we can find that the IAM algorithm can reduce the computational order of EAM criterion from square to one in each updating process.
The array configuration is uniform linear array composed of omnidirectional sensors, of which steering vector is represented as where  is the wavelength of signals impinging on the array and Δ is the sensor spacing between two adjacent sensors.As a necessary condition that a unique direction  is determined by the phase parameter , Δ ≤ /2 is imposed on the array configuration.In this paper, Δ = /2.Using the uniform linear array, we derive the irreducible form of (30) and (31).Define which are the varying parts in (30) and (31).First we derive the irreducible form of (30).Substituting ( 15) and (29) into   (), we get Using the uniform linear array defined above, the steering vector a() in (34) can be represented as follows: Then   () can be rewritten into the form of a rational function where   * () is the paraconjugate of   () defined as and the superscript * is the complex conjugate of a complex number.
Let   be variable and all other parameters are held fixed as well as in last subsection.As we have discussed in the problem of EAM, when   becomes equal to   , both the polynomials () and () have double zeros at  =  −(  ) in the complex z-plane, since it holds   ( −(  ) )P ⊥   = 0. Without canceling these common zeros,   () is indefinite at  =   .
Since the column vectors in W  are all orthogonal to a(  ), the projection matrix P ⊥ A l can be written as Here, note that   () represents a polynomial, while W  represents a matrix.Using the expressions where   () is the common zero factor at  =  −() ,  ≤ ,  ̸ = , and As for the irreducible form of f (), it can be derived like this form similarly, where Define the following polynomials where and similarly define ñ () to calculate the matrix Ñ which realizes the same function as   ().
Then we have the final form of   () and f () shown as follows: where Re{} represents the real part of the complex value.
Therefore the irreducible form of efficient AM criterion (IAM) is shown like this.
Furthermore, we can find that R N l , R S l , N  , Ñ , and D  can be calculated before each updating process.Therefore, at each updating process we only need to evaluate Re{  (  )}, Re{ñ  (  )}, and Re{  (  )}.Evaluation of them only needs ( − ) multiplications.Compared to (( − )) multiplications in each updating process of EAM algorithm, the computational complexity can be reduced further.
Here we should note that the IAM algorithm is only applicable to the ULA that is because the IAM algorithm is formulated based on (34) and (35).The steering vector can only be written into this form when the array is uniform linear array (ULA).Based on (34) and ( 35), the EAM criterion can be written into polynomial form.And then the irreducible form of EAM can be derived by canceling the common zero in both numerator and denominator.Therefore, the IAM is only applicable to the ULA.

Simulations
In this section, we show some simulation results to demonstrate the validity of the EAM and IAM algorithms.In simulation, the array configuration is a uniform linear array as discussed above.
The SNR is defined as The Root-Mean-Square-Error (RMSE) is defined as where θ, is the estimation of   at the th trial.The number of trails in our simulation is 100.
In Figure 1, the scenario is  = 3,  = 2,  = 2, SNR = 10,  = 100.Two true sources are independently located at 0 and 8 degrees.The estimated bearing  2 is fixed at 10 degrees, while  1 varies from 9.99999 to 10.00001.The dashed line represents the value of SML with EAM criterion, while the solid line represents the IAM criterion.It shows clearly that numerical instability occurs when the EAM criterion is used.In particular, the value changes violently around the point  2 =  1 = 10 degrees because it is indefinite.As for the IAM criterion, we can find that the value becomes monotonic and it is numerically stable.Due to the numerical instability, oscillation may happen as Figure 2 shows.The scenario is the same as Figure 1.In Figures 2(a) and 2(b), the estimated two bearings  1 and  2 converge to an identical value represented by the solid line and the dashed line.The convergence condition is that when the variation of each bearing is less than 10 −5 or when the iteration reaches the maximum number, 800.When the EAM criterion is used, the iteration does not stop until it reaches the maximum number because of oscillation.As for the IAM criterion, it converges well for less than 100 iterations.Figure 2(c) shows the oscillation rate of the two criterions in 100 independent trials.We can find that there is no oscillation when the IAM criterion is used.Note that the oscillation rate of EAM decreases when SNR increases as shown in Figure 2(c) that is because the oscillation happens when more than one parameter converges to an identical value.In simulations, two sources are independently located in 0 and 8 degrees.When SNR is relatively low, the case that the estimated two bearings are extremely close (e.g., 3.990 and 4.010) or even equal exists.In this case, oscillation happens.When SNR increases, the estimated two bearings will significantly separate each other.As a result, the oscillation rate will decrease.
Figure 3 shows the comparison of RMSE between SML, MUSIC, and ESPRIT.The scenario is the same as Figure 2 except for SNR.In simulation we use the original AM and our proposed EAM and IAM algorithms for SML estimation of DOA.We find that the RMSE of all the three algorithms are completely the same.The reason is that the proposed IAM and EAM algorithms are just transformations for SML criterion and there is no any modification.All these algorithms can find the solution of SML successfully (of course we have to delete the oscillation cases for AM and EAM criterions.Obviously, oscillation affects the estimation accuracy because the bearings do not converge to their optimal value).Also Figure 3 shows that the solution of SML is much better than that of MUSIC and ESPRIT.Figures 4 and 5 and Table 1 show the comparison of computational complexity of AM, EAM, and IAM.Figures 4 and 5 show the average amount of operations of each algorithm.In Figure 4,  = 3,  = 2,  = 2,  = 100.The true DOA is 0 and 8 degrees.In Figure 5  = 8,  = 5,  = 3,  = 100.The true DOA is 0, 8, 16, 24, and 32 degrees.Note that in Figure 5  < , which means that there are correlated signals (coherent case).
In Figures 4 and 5, "amount of arithmetic operations" represents the summation of all the complex addition, subtraction, multiplication, and division.These two figures show clearly that both EAM and IAM can reduce the computational complexity of SML estimation greatly, while IAM is the best.
Table 1 shows the detailed comparison of computational complexity including average iteration times and main calculation process of Figure 5 when SNR = 5 dB for 30 independent trials.For all these three algorithms, in the one-dimensional global search of each updating process, we use different step-sizes for searching.At first, we have a relatively large step-size (0.5 degrees) for rough search (note that this step-size should not be too large; otherwise the onedimensional global minimum may be missed) and then much smaller step-sizes (0.001 and 0.00001 degrees) for fine search.For example, for the AM algorithm in the one-dimensional global search of an updating process, the searching range for the variable parameter (  ) is from −90 ∘ to 90 ∘ .As a result, the rough search needs 180/(0.5)times of calculation of  SML .And then the fine search needs (0.5 * 2/(0.001)+ 0.001 * 2/(0.00001))times of calculation of  SML .Since the number of iterations is 127, the main computational complexity of AM algorithm is 83,820 times of calculation of  SML .For the EAM and IAM algorithms, the main computational complexity is also shown in Table 1 in the same manner.As we have discussed above, for the EAM and IAM algorithm, before each updating process, we can calculate many components in advance (the detailed components are shown in Table 1).We do not need to calculate the whole cost function  SML every time.As a result, the computational complexity of EAM and IAM is much smaller than that of AM algorithm.Furthermore, from Figures 4 and 5, we can find two phenomena.First, the efficiency of EAM and IAM is more obvious than that of AM when the number of parameters increases.Here the number of parameters represents the number of incident signals, that is, .The task of DOA is to find the directions of all the incident signals.When there are more parameters, the AM algorithm needs more iterations to solve the multidimensional nonlinear optimization problem of SML.For our proposed algorithms (EAM and IAM), in each updating process the computational complexity is greatly reduced.As a result, the efficiency is more obvious when the number of parameters increases (in Figure 5 there are 5 parameters to be estimated, and in Figure 4 there are 2).Second, the computational complexity of all these three methods changes with SNR that is because the computational complexity of all these methods mainly depends on the iteration times for the initial value converges to the estimated value.As we show in Section 4.1, the initial value is determined by the [Initialization Phase].Obviously, the initial value will change according to different SNRs.Therefore, the total iteration times will change with SNR.As a result, computational complexity also changes with SNR.

International Journal of Antennas and Propagation
In a huge number of simulations, we have confirmed the efficiency of our proposed EAM and IAM algorithms, while IAM is the best.

Conclusions
In this paper, to reduce the computational complexity of SML estimation, based on the AM algorithm, we propose two more efficient algorithms, that is, EAM and IAM algorithms.The EAM algorithm mainly uses transformation of matrix and properties of spatial projection to divide the SML criterion into two components.One is variable, while the other is fixed.Computational complexity can be greatly reduced since the fixed part can be calculated once in advance and only the varying part should be calculated in each onedimensional updating process.To avoid numerical instability of EAM (note that because of the numerical instability, wrong estimation of DOA may be got) and to reduce computational complexity further, we derive the irreducible form of EAM, that is, IAM algorithm.The main idea of IAM is to use a uniform linear array and rewrite the EAM criterion into polynomial forms.Then the irreducible form can be got by canceling the common zero factor of the polynomial form.Simulation results show that the IAM algorithm can avoid the numerical instable problem of EAM and reduce computational complexity further.
Phase.For a certain value of SNR, first assuming a single signal,  = 1, find θ1 minimizing  1 ( Θ(1) ) by onedimensional global search with respect to θ1 .Next, assuming two signals,  = 2, and fixing θ1 at the value obtained for a single signal, find θ2 minimizing  2 ( Θ(2) ) by onedimensional global search with respect to θ2 .Continue in this fashion until all the initial values for θ ,  = 1, 2, . . .,  are computed.

Figure 1 :
Figure 1: Numerical instability in EAM and stability in IAM criterion.

Figure 2 :
Figure 2: Oscillation in EAM algorithm while not in IAM.