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

Authors: Joon-Ho Lee; Sung-Woo Cho; So-Hee Jeong
Journal: International Journal of Antennas and Propagation
 (2013)

---

## Abstract

We note that the cost function of the MUSIC (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.


---

## Body

## 1. Introduction

There has been much study on the determination of the AOA (angle-of-arrival) of incident signals using array antenna structure [1–15]. Array signal processing algorithm can also be extended to radar target recognition [16, 17] and detection [18, 19].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 two-dimensional 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 areM incident signals. The cost function of the maximum likelihood (ML) algorithm [1, 3–7] is M-dimensional for estimation of azimuth angle only and 2M-dimensional for simultaneous estimation of azimuth angle and elevation angle. For large value of M, the optimization of M-dimensional or 2M-dimensional cost function by exhaustive grid search is not feasible. That is why Newton-type 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 M. 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 the optimization of the two-dimensional cost function of the MUSIC algorithm.There has been much study on applying the Newton-type method to the optimization of the cost function of the ML DOA algorithm [3–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.


## 2. Standard MUSIC Algorithm

Consider the case of a uniform circular array (UCA) withN antenna elements. M 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)P(θ,ϕ)=1aH(θ,ϕ)UUHa(θ,ϕ).The array vector can be written as(2)a(θ,ϕ)=[exp(jψ1)exp(jψ2)⋯exp(jψN)]T,
where ψn for the UCA is defined as follows:
(3)ψn(θ,ϕ)=2πrλcosϕ[cos(θ-2π(n-1)N)]n=1,…,N.In the standard MUSIC algorithm, to evaluate (1), two-dimensional exhaustive grid search with search step of [Δθ,Δϕ] has to be implemented.


## 3. LM Implementation of the MUSIC Algorithm

LetF(θ,ϕ) be defined from the MUSIC spectrum in (1) as
(4)F(θ,ϕ)=1P(θ,ϕ)=aH(θ,ϕ)UUHa(θ,ϕ).It is quite clear that the local maxima ofP(θ,ϕ) correspond to the local minima of F(θ,ϕ).Let a vector valued functionf(θ,ϕ) be in the form of
(5)f(θ,ϕ)=[f1(θ,ϕ)f2(θ,ϕ)⋯f2(N-M)(θ,ϕ)]T.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(θ,ϕ)∥22 corresponds to the cost function of F(θ,ϕ) in (4):
(6)F(θ,ϕ)=∥f(θ,ϕ)∥22=∑m=12(N-M)‍|fm(θ,ϕ)|2.The MUSIC algorithm can be formulated as(7)[θ^,ϕ^]=argminθ,ϕF(θ,ϕ)=argminθ,ϕ∥f(θ,ϕ)∥22=argminθ,ϕ∑m=12(N-M)‍|fm(θ,ϕ)|2.We will show how to express the cost function in (4) in least squares form. From (4), we have
(8)F(θ,ϕ)=aH(θ,ϕ)UUHa(θ,ϕ)=∥UHa(θ,ϕ)∥2.In (8), note that UHa(θ,ϕ) is a column vector with N-M entries. Let the ith entry of UHa(θ,ϕ) be denoted by ei(θ,ϕ), i=1,2,…,N-M:
(9)ei(θ,ϕ)≡(UHa(θ,ϕ))i=∑n=1N‍(Re(ui)n-jIm(ui)n)×(cos(ψn(θ,ϕ))+jsin(ψn(θ,ϕ))).Note that, sincea(θ,ϕ) is complex-valued, ei(θ,ϕ) is also complex valued. Let e(θ,ϕ) be defined as follows:
(10)e(θ,ϕ)=[e1(θ,ϕ)⋯eN-M(θ,ϕ)]T=[u1Ha(θ,ϕ)⋯uN-MHa(θ,ϕ)]T.Using (8) and (10), it can be easily shown that the cost function of the MUSIC algorithm can be written as
(11)F(θ,ϕ)=∥e(θ,ϕ)∥22=∑i=1N-M‍|ei(θ,ϕ)|2.Letf2i-1(θ,ϕ) and f2i(θ,ϕ) in (5) be defined by the real part and the imaginary part of ei(θ,ϕ), respectively. From (9), we have, for i=1,2,…,N-M,
(12)f2i-1(θ,ϕ)≡Re[ei(θ,ϕ)]=Re(uiHa(θ,ϕ))=∑n=1N‍(Re(ui)ncos(ψn(θ,ϕ))+Im(ui)nsin(ψn(θ,ϕ))),f2i(θ,ϕ)≡Im[ei(θ,ϕ)]=Im(uiHa(θ,ϕ))=∑n=1N‍(Re(ui)nsin(ψn(θ,ϕ))-Im(ui)ncos(ψn(θ,ϕ))).From the definition off2i-1(θ,ϕ) and f2i(θ,ϕ), it is clear that the following is true:
(13)|ei(θ,ϕ)|2=(f2i-1(θ,ϕ))2+(f2i(θ,ϕ))2i=1,…,N-M.From (11) and (13), the cost function can also be written in terms of fm(θ,ϕ)′s as follows:
(14)F(θ,ϕ)=∑i=1N-M‍|ei(θ,ϕ)|2=∑i=12(N-M)‍(fi(θ,ϕ))2.The Jacobian matrix corresponding to the matrixf(θ,ϕ) in (5) is defined as
(15)J(θ,ϕ)=[∂f1∂θ∂f1∂ϕ⋮⋮∂f2(N-M)∂θ∂f2(N-M)∂ϕ].Using (12) in (15), the explicit expressions for the entries of the matrix J(θ,ϕ) are, for i=1,2,…,N-M,
(16)∂f2i-1∂θ=∑n=1N‍ψn(θ,ϕ)∂θ(-Re(ui)nsin(ψn(θ,ϕ))+Im(ui)ncos(ψn(θ,ϕ))),∂f2i-1∂ϕ=∑n=1N‍ψn(θ,ϕ)∂ϕ(-Re(ui)nsin(ψn(θ,ϕ))+Im(ui)ncos(ψn(θ,ϕ))),∂f2i∂θ=∑n=1N‍ψn(θ,ϕ)∂θ(Re(ui)ncos(ψn(θ,ϕ))+Im(ui)nsin(ψn(θ,ϕ))),∂f2i∂ϕ=∑n=1N‍ψn(θ,ϕ)∂ϕ(Re(ui)ncos(ψn(θ,ϕ))+Im(ui)nsin(ψn(θ,ϕ))).In Figure1, we outline the Levenberg-Marquardt algorithm. In the problem at hand, x can be written as x=(θ,ϕ), and k is the number of the iterations [20], and F(θ,ϕ) is defined in (8) or (14).Figure 1
Levenberg-Marquardt algorithm [20].In the LM algorithm, the updates of the estimates are obtained from(17)[J(θ,ϕ)TJ(θ,ϕ)+μI][θLMϕLM]T=-J(θ,ϕ)Tf(θ,ϕ).Due toμ>0, the coefficient matrix is positive definite, which makes [θLMϕLM]T be a descent direction. When μ is very large, we get
(18)[θLMϕLM]T≃-1μg
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].L
(
θ
LM
,
ϕ
LM
) in Figure 1 is defined as [20, 21]
(19)L(θLM,ϕLM)=F(θ,ϕ)+[θLMϕLM]J(θ,ϕ)Tf(θ,ϕ)+12[θLMϕLM]J(θ,ϕ)TJ(θ,ϕ)[θLMϕLM]T.


## 4. 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].


## 5. Numerical Results

In numerical results, the UCA is used. Therefore,f(θ,ϕ) is obtained from (12), and J(θ,ϕ) is given by (16) by substituting (3) for ψn(θ,ϕ).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 kmax 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 Figure2. The criterion for the global convergence is that the RMSE is arbitrarily chosen to be less than two degrees.Initial azimuth angles and the initial elevation angles resulting in the global convergence for various radii at SNR = 4 dB.
(a)
r
=
0.2
λ
(b)
r
=
0.4
λ
(c)
r
=
0.6
λ
(d)
r
=
0.8
λFor the fixed number of the antenna elements, as the radius of the antenna array becomes larger, the number of the spurious peaks ofP(θ,ϕ) increases. The initial estimate close to the spurious peak of P(θ,ϕ) will converge to the spurious peak, which is not the global minimum of F(θ,ϕ). 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 Figure3, 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.Comparison of the standard MUSIC algorithm with fine exhaustive grid search and the MUSIC algorithm with coarse exhaustive grid search followed by LM optimization (r=0.679λ, five antennas, [θ1(true)ϕ1(true)]=[-21∘30∘],  [θ2(true)ϕ2(true)]=[120.5∘72.5∘]).
(a)
[
Δ
θ
fine
,
Δ
ϕ
fine
]
=
[
2
∘
,
2
∘
], [Δθcoarse,Δϕcoarse]=[4∘,4∘]
(b)
[
Δ
θ
fine
,
Δ
ϕ
fine
]
=
[
2
∘
,
2
∘
], [Δθcoarse,Δϕcoarse]=[7∘,7∘]
(c)
[
Δ
θ
fine
,
Δ
ϕ
fine
]
=
[
4
∘
,
4
∘
], [Δθcoarse,Δϕcoarse]=[7∘,7∘][
Δ
θ
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 Figure4, 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.Comparison of the standard MUSIC algorithm and the standard Capon algorithm with fine exhaustive grid search and the MUSIC algorithm and the Capon algorithm with coarse exhaustive grid search followed by LM optimization (r=0.679λ, five antennas, [Δθfine,Δϕfine]=[4∘,4∘], [Δθcoarse,Δϕcoarse]=[7∘,7∘]).
(a)
[
θ
1
(
true
)
ϕ
1
(
true
)
]
=
[
-
2
1
∘
3
0
∘
], [θ2(true)ϕ2(true)]=[120.0∘72.0∘]
(b)
[
θ
1
(
true
)
ϕ
1
(
true
)
]
=
[
85
.
5
∘
63
.
5
∘
], [θ2(true)ϕ2(true)]=[120.5∘72.5∘]In Figure4(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 Figure4(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.


## 6. 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-of-arrival) 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.



---

*Source: 402842-2013-08-27.xml*