Application of the Levenberg-Marquardt Scheme to the MUSIC Algorithm for AOA Estimation

We note that the cost function of theMUSIC (multiple signal classification) algorithm is quadratic in an array vector and that it can be expressed in a least squares form. Based on this observation, we present a rigorous Levenberg-Marquardt (LM) formulation of the MUSIC algorithm for simultaneous estimation of an azimuth and an elevation. We show a convergence property and compare the performance of the LM-based MUSIC algorithm with that of the standard MUSIC algorithm via Monte-Carlo simulation. We also compare the performance of the MUSIC algorithm with that of the Capon algorithm both for the standard implementation and for the LM-based implementation.

In this paper, we consider a simultaneous estimation of the azimuth and the elevation using the MUSIC algorithm [2].In standard MUSIC algorithm, we have to perform twodimensional exhaustive grid search, which can be computationally intensive when the search step is small.The idea proposed in this paper is to apply the LM algorithm to the numerical optimization of the cost function defined from the MUSIC algorithm.From the numerical results, it is shown that a performance of the proposed scheme is better than that of the standard MUSIC algorithm in terms of the computational burden and the estimation accuracy.
In this paper, we give a rigorous formulation for the improvement of the azimuth/elevation AOA estimation for the UCA structure by applying the LM scheme [20,21] to the cost function of the MUSIC algorithm.
Assume there are  incident signals.The cost function of the maximum likelihood (ML) algorithm [1,[3][4][5][6][7] is dimensional for estimation of azimuth angle only and 2dimensional for simultaneous estimation of azimuth angle and elevation angle.For large value of , the optimization of -dimensional or 2-dimensional cost function by exhaustive grid search is not feasible.That is why Newtontype algorithms including the LM algorithm have been widely applied to the optimization of the cost function of the ML algorithm.
Note that, in the MUSIC algorithm [1,2], the cost function is one-dimensional for estimation of azimuth angle and two-dimensional for simultaneous estimation of azimuth angle and elevation angle, which is independent of the actual number of the incident signals .When we are only concerned with the estimation of the azimuth angle, the optimization of one-dimensional cost function using exhaustive grid search even with small search step is quite feasible.But, for simultaneous estimation of the azimuth and the elevation, the computational complexity of two-dimensional optimization by exhaustive grid search is much larger than that of one-dimensional optimization by exhaustive grid search, which justifies why we apply the LM algorithm to International Journal of Antennas and Propagation the optimization of the two-dimensional cost function of the MUSIC algorithm.
There has been much study on applying the Newtontype method to the optimization of the cost function of the ML DOA algorithm [3][4][5][6][7].To the best of our knowledge, there has been no study on applying the LM algorithm to the optimization of the cost function of the MUSIC algorithm for simultaneous estimation of the azimuth angle and the elevation angle.

Standard MUSIC Algorithm
Consider the case of a uniform circular array (UCA) with  antenna elements. signals are incident on the array, and we are to estimate the azimuth angle and the elevation angle of each incident signal.Note that, for simultaneous estimation of the azimuth and elevation, a uniform linear array (ULA) cannot be adopted since ULA-based algorithm cannot uniquely estimate the azimuth and the elevation simultaneously due to the ambiguity pertinent to the ULA structure.
Let  and  represent the azimuth and the elevation, respectively.We denote the noise eigenvector matrix of the array covariance matrix U. Then the MUSIC spectrum can be usually defined as [4]  (, ) = 1 a H (, ) UU H a (, ) .
The array vector can be written as where   for the UCA is defined as follows: In the standard MUSIC algorithm, to evaluate (1), twodimensional exhaustive grid search with search step of [Δ, Δ] has to be implemented.

LM Implementation of the MUSIC Algorithm
Let (, ) be defined from the MUSIC spectrum in (1) as It is quite clear that the local maxima of (, ) correspond to the local minima of (, ).
Let a vector valued function f(, ) be in the form of To apply the LM algorithm to the optimization of the cost function of the MUSIC algorithm in ( 4), what we have to do first is to define f(, ) so that ‖f(, )‖ 2  2 corresponds to the cost function of (, ) in (4): The MUSIC algorithm can be formulated as We will show how to express the cost function in (4) in least squares form.From (4), we have In ( 8), note that U H a(, ) is a column vector with  −  entries.Let the th entry of U H a(, ) be denoted by   (, ),  = 1, 2, . . .,  − : Note that, since a(, ) is complex-valued,   (, ) is also complex valued.Let e(, ) be defined as follows: Using ( 8) and ( 10), it can be easily shown that the cost function of the MUSIC algorithm can be written as Let  2−1 (, ) and  2 (, ) in ( 5) be defined by the real part and the imaginary part of   (, ), respectively.From (9), we have, for  = 1, 2, . . .,  − , Figure 1: Levenberg-Marquardt algorithm [20].
From the definition of  2−1 (, ) and  2 (, ), it is clear that the following is true: From ( 11) and ( 13), the cost function can also be written in terms of   (, )  s as follows: The Jacobian matrix corresponding to the matrix f(, ) in ( 5) is defined as Using ( 12) in (15), the explicit expressions for the entries of the matrix J(, ) are, for  = 1, 2, . . .,  − ,  In Figure 1, we outline the Levenberg-Marquardt algorithm.In the problem at hand, x can be written as x = (, ), and  is the number of the iterations [20], and (, ) is defined in (8) or (14).
In the LM algorithm, the updates of the estimates are obtained from Due to  > 0, the coefficient matrix is positive definite, which makes [ LM  LM ]  be a descent direction.When  is very large, we get which is a short step in the steepest descent direction.Therefore, we choose large  when the current estimates are assumed to be very far from the solution.This is good if the current iterate is far from the solution.If  is very small, then the updates of the LM method are approximately equal to the updates of the Newton method, which is a good step in the final stages of the iteration [20].
can be chosen to be small value of  = 10 −6 if the initial estimate, [ 0 ,  0 ], is assumed to be close to the solution.Otherwise,  is chosen to be larger value of  = 10 −3 [20].

Comment on How to Get Initial Estimates
The alternating projection (AP) algorithm [22,23] is very popular algorithm for getting the initial estimates of the AOA  estimation since it is computationally efficient.The shortcoming of the AP scheme is that the initial estimates obtained from the AP scheme are less accurate in comparison with those obtained from the exhaustive search algorithm.
Recently, there is a study [15] on the extension of the original AP scheme provided in [22,23].
The parameters for the LM algorithm are defined as follows: both  1 and  2 are chosen to be 10 −5 . is set to be 1.0.The maximum number of the iteration which corresponds to  max is 100.The RMSE (root mean square error) is obtained from ten repetitions.
In the first example, the LM algorithm is initialized on a grid of azimuth-elevation to check which initial angles finally achieve the global convergence.Therefore, initial estimates used for the LM algorithm are chosen arbitrarily to see whether each of the initial estimates actually leads to the global minimum or not.The cost function is evaluated at each two-dimensional grid point with azimuth increment of   Δ = 4 ∘ and the elevation increment of Δ = 4 ∘ followed by the LM optimization.Since the scheme is based on the LM formulation, it is shown that once the initial estimates are within the main lobe, the initial estimates are supposed to converge to the local minimum corresponding to the main lobe which is actually a global minimum.
The initial azimuth angles and the initial elevation angles which result in the global convergence are marked in Figure 2. The criterion for the global convergence is that the RMSE is arbitrarily chosen to be less than two degrees.
For the fixed number of the antenna elements, as the radius of the antenna array becomes larger, the number of the spurious peaks of (, ) increases.The initial estimate close to the spurious peak of (, ) will converge to the spurious peak, which is not the global minimum of (, ).Therefore, the larger radius of the antenna array results in smaller domain of initial angles leading to global convergence as shown in Figure 2, which reflects that the proposed scheme is very efficient in the case of the smaller radius of the antenna array.
In Figure 3, for two incident signals, we compared the performance of the LM optimization following the MUSIC algorithm of coarse exhaustive grid search with that of the standard MUSIC algorithm of fine exhaustive grid search.
The purpose is to illustrate the fact that the proposed scheme still works for more than one incident signal.
[Δ coarse , Δ coarse ] and [Δ fine , Δ fine ] denote the search step of the coarse exhaustive grid search and the fine exhaustive grid search, respectively.We arbitrarily choose Δ coarse = Δ coarse and Δ fine = Δ fine .We can see that the LM implementation of the MUSIC algorithm is superior to the standard MUSIC algorithm for all three cases in terms of the RMSE and the computational complexity.
In Figure 4, we compare the performances of the MUSIC algorithm and the Capon algorithm [11].Each algorithm is implemented in two ways: the first way is the standard MUSIC algorithm and the standard Capon algorithm of fine exhaustive grid search.The second way is the LM optimization following the MUSIC and Capon algorithms of coarse exhaustive grid search.
In Figure 4(a), it is shown that the performances of the MUSIC algorithm and the Capon algorithm are nearly equal for the case that the incident signals are separated enough.In Figure 4(b), for the case that two incident signals are close to each other, it is shown that the performance of the MUSIC algorithm is better than that of the Capon algorithm, both for the standard implementation and for the LM-based implementation.The results in Figure 4(b) can be justified due to the fact that the MUSIC algorithm is a superresolution algorithm, while the Capon algorithm is not a superresolution algorithm.
In the standard implementation of the Capon algorithm, when two incident signals are close to each other in Figure 4(b), the Capon algorithm fails to resolves two closely incident signals at low SNRs of 4 dB and 8 dB, which is why the RMSEs of the first signal at SNRs of 4 dB and 8 dB are very large.On the other hand, the corresponding RMSEs for the MUSIC algorithm are much smaller than those for the Capon algorithm, which reflects that the MUSIC algorithm successfully resolve two closely incident signals at low SNRs, while the Capon algorithm cannot resolve.For high SNRs of 12 dB to 40 dB in increment of 4 dB, the RMSEs of the first signal for the Capon algorithm and for the MUSIC algorithm are small, which reflects that both algorithms can resolve two closely incident signals at high SNRs.

Conclusion
In this paper, we propose to apply the LM formulation to the optimization of the cost function of the MUSIC algorithm.Using the LM algorithm, we implemented the nonlinear multivariable optimization for azimuth/elevation AOA (angle-ofarrival) estimation based on the MUSIC algorithm.
We have compared the performance of the LM implementation of the MUSIC algorithm with that of the standard MUSIC algorithm in terms of the estimation accuracy and the computational complexity.
We illustrated that the performance of the MUSIC algorithm is better than that of the LM-based implementation as well as the standard implementation when there are two incident signals that are very close to each other.
We consider the array structure of the UCA to be able to uniquely estimate the azimuth and elevation simultaneously, but it is quite straightforward to extend the proposed scheme to an arbitrary array structure by simply modifying the array vector consistently with the specific array structure as long as the adopted array structure is able to uniquely estimate both the azimuth angle and the elevation angle.In addition, the proposed scheme can also be applied to the other AOA estimation algorithm by simply modifying the cost function consistently with the adopted AOA estimation algorithm as long as the corresponding cost function can be reformulated in the least squares form.

Figure 2 :
Figure 2: Initial azimuth angles and the initial elevation angles resulting in the global convergence for various radii at SNR = 4 dB.