Spherical subspace tracking for efficient, high performance adaptive signal processing applications

SIGNAL PROCESSING ELSEVIER

Signal Pro,:essing 50 11996~ I01 121

Spherical subspace tracking for efficient, high performance adaptive signal processing applications 1 Ronald D. DeGroat*, Eric M. Dowling, Hao Ye, Darel A. Linebarger Erik Jonsson School of Engineering and Computer Sctem'e, The Unirer ,ity of Texas at Dallas. P.O Box 830688. MS EC33. Richardson. TX 750b;.5-06,'18, USA

Received 29 March 1995: revised 20 September 1995

Abstract in this paper we develop and analyze multilevel spherical subspace updating algorithms. These algorithms track eigenvectors and/or eigensubspaces associated with time-varying data that arises in adaptive signal processing, direction of arrival (DOA) or frequency tracking problems. Spherical subspaces can be created by setting groups of cigenvalues to their average val-e. Because of the sphericity, the subspace basis vectors can be reoriented to deflate the update problem and reduce computation. In this paper, we focus on a four-level update that efficiently tracks the signal subspace and generates just enough eigenlevel information to make strongly consistent signal subspace rank estimates. The four-level method is O(nr) per update when n is the data dimension and r is the signal subspace tank. In order to analyze the multilevel algorithms, we use Ljung's method to derive ordinary differential equations (ODEs) that characterize the asymptotic eigenlevel and subspace projection operator trajectories. Under certain weak assumptions, we prove that the multilevel algorithms asymptotically converge with probability one to the true subspace and eigen!evel values. We also present a four-level MDL-based rank detection scheme and show it to be strongly consistent. Simulation studies support the theoretical results and demonstrate the performance of our methods with frequency tracking data. Zesammenfasseag In dieser Arbeit entwickeln und analysieren wir mehrstufige Algorithmen zur Nachfiihrung sphiirischer Unter~ume. Diese Algorithmen verfolgen Eigenvektoren und/oder Eigenunterr~iume, die in Zusammenhang mit zeitvarianten Daten stehen, wie sic bei der adaptiven Signalverarbeitung, bei Problemen der Sch~itzung tier Einfallsrichtung (IX)A) oder bei Frequcnzveffolgungsproblemen auftreten. Sph/irisehe Unterr~iume k6nnen entstehen, indem man Gruppen yon Eigenwerten auf ihren Mittelwert setzt. Wegen der sph/irischen Symmetric k6nnen die Basisvektoren des Unterraums neu ausgerichtet werden, so dab das Nachf'tihrproblem vereinfacht und der Recbenaufwand reduziert wird. In dieser Arheit konzentrieren wit ans auf einen vierstufigen Nachfiihralgorithmus, der den Signalunterraum wirksam veffolgt and gerade ausreicbende Information fiber die Eigenniveaus erzeugt, um cine streng konsistente Sch~itzung des Rangs des Signalunterraums zu erm6glichen. Die vierstufig¢ Methode ben6tigt O(nr) Rechenoperationen pro Schritt, wobei n di¢ Dimension der Daten und r der Rang des Stgnalunterraums ist. Um die mehrstufigen Algorithmen zu analysieren,

* Corresponding author. Tel.: 214-883-2894;e-mail: degroat@ utdatlas.edu. *This research was supported in part by the National Science Foundation Grant MIP-9203296 and the Texas Advanced Research Program Grant 009741-022 0165-1684/96/$15.00 i" 1996 Elsevier Science B.V. All rights reserved Pll S01 6 5 - 1 6 8 4 ( 9 6 ) 0 0 0 0 4 - 7

102

R.D. DeGroat et al. / Signal Processing 50 (1996) I01-1]1

verwenden wir Ljungs Verfahren zur Ableitung gew6hnlicher Differentialgleichungen (ODEs), die die asymptotischen Tmjektorien der Eigenniveaus und des Unterraum-Projektionsoperators charakterisieren. Unter bestimmten schwachen Annahmen beweiscn wir, dab die mehrstufigen AIgorithmen mit Wahrscheinlichkeit ! gegen den echten Unterraum und di¢ echten Werte der Eigenniveaus asymptotisch konvergiercn. Wir stellen auch ein vicrstufiges MDL-basienes Verfahma zur Rangdetektion vor und zeigen, da6 es streng konsistent ist. $imulationsstudien unterstiitzen die theoretischen Ergcbnisse und stellen die Leistungsf'dhigkeit unse~er Methoden an Hand yon Frequenzverfolgl.mgsdaten untcr Bewcis. gmanl6 Nous developpons et analysons dans oct article des algorithmes de mise /t jour dc type sous-espaoc sph~rique muRi-niveaux. Ces algorithmes poursuivent les vecteurs propres et/ou les sous-espaocs propres associ6s ~vec les donn6es variant dans le temps provenant de traitements adaptatifs de signaux, de probl6mes de peursuite de direction d'arriv6e ou dc fr~qucncc. Los sous-espaces sph6riqucs peuvent etre cr66s en mettant/~ leer valeur moyenne des groupes de valeurs proprcs. Du fait du caract6re sph6rique, Its vccteurs de base des sous-espaocs peuvent ~tre r~-oricnt6s de mani6rc /Lsimplifier lc probleme de mise ~ jour et r6daire its calculs. Dans cet article, nous nous conccntrons sur une mise ~ijour quatre niveaux qui pout'suit ellicaocment le sous-.cspaocde signal et g6n6re juste asscz d'information sur Ic nivcau des valeurs propres pour obtenir des estim&'s de rang du sous-espaoc signal fortcment consistentes. La m6thodc A quatrc niveaux ¢st O(nrl par raise/~ jour o/l nest la dimension des donn6es et rest le rang du sous-espaoc signal. Dans le but d'analyser les algorithmes multi-niveaux, nous utilisons ia m6thode de Ljung pour d6river les 6quations diff6rentielles ordinaires qui caract6risent le nivcau asymptotique des valeur pmpres et les trajectoires des op6rateurs de projection sur Ic sous-cspacc. Sous certaines hypothc~ses faibles, nous pronvons qu¢ ies aigorithmes multi-niveaux convergent avec probabilit6 un vers le vrai sous-cspacc et les vrais nivcaux des vaicurs propres. Nous pr~scntons 6galemcnt une technique dc d6tcction de rang quatrc-niveaux bas6c sur le MDL et montrons qu'clle et fortement consistente. Des 6tudes de simulation 6tayent les r6sultats th6,oriques et mettcnt en lumi~.rc lcs performances de nos m6thodcs sur des donn6es de poursuite de fr6quenoc. Kcgavords: Array signal proo~ssing; Subspace methods; Recursive estimation

1. lmtrodactiom Most modern high-resolution direction of arrival (DOA) and frequency estimation techniques extract information from a subspace of the correlation or data matrix associated with an observed vector process. These subspace based methods include Pisarenko's harmonic retrieval algorithm [45], MUSIC [49], ESPRIT [48], SVD-based reduced rank linear prediction [53], the minimum norm method [3g], the TLS method [97], and others. Eigen and subspace tracking methods are also beJag appfied to a variety of adaptive signal processing problems, e.g., [18,19, 56]. In order to achieve practical, real-time adaptive implementations, computationally efficient, numerically stable subspace tracking methods are needed. Many researchers have studied SVD, EVD and other subspace tracking-oriented problems. A broad overview of SVD-based tracking algorithms is presented by Comon and Golub in [7]. Golub [333 introduced the first eigen-updating

scheme, and his ideas were developed and expanded by Bunch and his co-workers in [5, 41. The rank-one eigen update was simplified in [50], when Schreiber introduced a transformation that makes the core eigenproblem real. Based on an additive whi~¢ noise model, Karasalo suggested that the noise cigenvalues be replaced by their average value so that deflation could be used to significantly reduce computation [37]. DeGroat reduced computation further by extending this colicept to the signal subspace [9,10,14,13]. To make eigen updating more practical, DeGroat and Roberts developed stabilization schemes to control the loss of orthogonality due to the buildup of roundoff error [ 15]. Further work related to eige',vector stabilization is reported in [43,8,24,41J. Owsley pioneered orthogonal iteration and stochastic-based subspace trackers in [44]. Stochastic gradient methods were advanced, analyzed and made systolic in [61,59]. Subspace trackers based on the Lanczos and conjugate gradient methods appear in [7, 27, 54, 57, 58, 29, 32]. The

R,D. DeGroat et al. / Signal Proce,ssing .50 (1996) 101-121

Lanczos-based algorithm in [58] has a strongly consistent rank estimator used to track the number of sources. Efficient RLS and perturbation-based subspace trackers appear in [60, 6], respectively. Fast systolic SVD-based subspace *.rackers are pre~ented in [42, 26,17]. Subspace estimation based on URV or QR decompositions appears in [!, 2, 51, 25]. Spherical subspace (SS) tracking algorithms average sets of eig:nvalues and apply deflation to reduce computation [37, 9,10, 20, I 1, 19,16, 46]. In [19] it was shown that the two-level SS update [9, 10] is asymptotically a stochastic gradient (LMStype) algorithm. We call this update the two-level signal averaged (SA2) subspace tracker. SA2 averages all of the signal eigenvalues into a signal eigenlevel and all of the noise eigenvalues into a noise eigenlevel. SA2 has been analyzed using algorithm ordinary differential equations (ODEs) [39~. Since every subspace has a unique associated projection operator, the subspace tracking algorithm will have an associated asymptotic trajectory over a projection operator manifold in C n~n. In [20] w~ derived the SA2 algorithm ODE and used a Lyapunov function to show convergence with probability one (w.p.l). In fact, the SA2 subspace trajectories can be described with Lie bracket notation and follow an isospectral flow as described by Brockett's ODE [11, 3, 36]. This paper centers around three main results. First, we generalize the SA2 algorithm to multiple levels. In particular, we derive the four-level algorithm (SA4) to efficiently track the signal or noise subspace as well as the signal rank. This is important, as any practical subspace tracking system must be able to estimate the number of sc~rces in order to properly track them. The second main result is to generalize the SA2 algorithm ODE in [20] to the multilevel case. An important corollary result is that we can use the ODE approach to prove asymptotic convergence w.p.l for the SA4 algorithm. Also, asymptotic convergence w.p.l can now be shown for Karasalo'~ noise deflated update [37] which is widely used [7, 50,17, 30, 31, 46]. The third main result is that we present a strongly consistent SA4-based minimum description length (MDL) [55] detection scheme (SA4-MDL) to track the number of sources. The rest ol this paper is organized as follows, in Section 2 we present the signal model and lay the

mathematical framework for multilevel subspace tracking, in Section 3 we present the SA4 subspace tracking algorithm. In Section 4 we derive the multilevel ODE and define a Lyap~mov function over a projection operator manifold to prove asymptotic convergence w.p.I in a stationary environment. In Section 5 we present a customized SA4-MDL detection scheme used to track the number of sources. In Section 6 we show the SA4MDL rank tracking scheme to be strongly consistent. In Section 7 we present simulation results to veri~ our analytical results as well as to demonstrate the usefulness and performance of the SA4 subspace tracker using SA4-MDL in a nonstationary setting, in Section 8 we offer conclusions.

2. Signal model aMI n m l t i l ~ d framework

Consider the vector time series {x(k)} whose associated slowly time-varying correlation matrix may be estimated recursively according to R(k) = ~(k)R(k - 1) + (1 - ~t(k))x(k)x(k) n,

(1)

where 0 < ~(k) ~< 1 is a fading factor and ~(k)is the correlation matrix estimate at time index k. Suppose that the vector time series is generated by the standard narrow-band array processing model (for

further details see [52, 17, 21]), x(k) = ~ a(Os(k))bs(k) + R(k), j=!

(2)

where a(Oj(k))¢ C m is the j t h signal vector whose direction angle at time k is given by Os(k), hj(k) is the jth modulating process, and n(k)¢C ~ is additive white noise. Here the set of signal vectors, {a(Os(k))!Oj(k)¢ R , j - - 1. . . . . p} defines the signal subspace. If ~.e define the matrix A (O(k))¢ C* x• to have a(Oj(k)) as its flh column, then we can write

x(k) = A(O(k))b(k) + R(k),

(3)

where now thejth element orb(k) c C • is bs(k). With this matrix-vector formulation, the true instantaneous correlation matrix at time k may he written

a(k) = E Ex(k)x(k) H] = A (O(k))B(k)A "(O(k)) + ~ !, (4)

104

R.D. DeGroat et al. / Signal Processing 50 (1996) 101-121

where B(k) & E[b(k)b(k)a,l. Let the eigenvalue decomposition (EVD) of R(k) be given by

R(k) = ~ ~./(k)wdk)w~(k)= W(k)A(k)W"(k), (5) i=l

where w/are unit norm eigcnvectors and ~.~are the corresponding ¢igenvalues with W = [wt[ ... Iw~,l and A = diag(:q . . . . . ,~n). To simplify notation, we often drop the time index, k, but we remind the reader that the correlation matrix and its EVD may be slowly varying functions of time. in order to track the essential eigcninformation of(5) in a computationally efficient manner, spherical subspace updating can be used. An m-level sphericalized version of (5) can be defined as

• =S(R;s~,s, . . . . . s , _ , ) = ~ d~U~UH i=l

= ~ U/D/U~ = ~ diP/, r=l

(6)

i=l

where Ur is the ith spbericalized subspac¢, d/is the corresponding eigenlevel and Pi = U/U~ is the projection operator associated with U~. There are m spherical subspaces in (6). Each subspace, U~, is composed of m/ basis vectors with the s/ values defining the boundary rang~ of the basis vectors. More specifically, if we let Tr be an appropriately sized arbitrary unitary matrix, we can relate the sphericalized subspaces to the corresponding subspaces in (5) as U~ = W/T/= [w~. . . . . . . . w,,,l T/=

[,as.... ~ ... u J f o r i = l , . . . , m , s o = O < s t < ... < s , = n , d r = a v e (A. . . . . i . . . . . ":-.O -~ (l/m/)x ~ ' = ..... l ;tk, Dr = dJ,,, is an m~/mr diagonal matrix corresponding to the eigenlevel d~,m/= s ~ - s/_ ~ and ~;.7=tm/= n. Notice that the basis vectors within a given subspaee can be arbitrarily reoriented with any unitary transformation, T/. This arbitrariness can be exploited to deflate an associated updating problem. The sphericalizing function defined in (6) is most easily understood by studying the three examples that follow. We should also note that averaging the eigenlevels to achieve sphericalization is not necessary for subspaee convergence. As we will see in Section 4, subspace convergence is achieved even when the eigenlevels are set to distinct, arbitrary constant values. Another point to notice is that the U/

eigenspaces of (6), like the w/eigenvectors of (5), are ill-defined if the corresponding eigenlevels/eigenvalues are repeated. Since EVD and SS EVD updating is normally applied to finite amou,ts of data, we can genersIly assem¢ that the estimated eigenleveis/eigenvalues will be distinct. For the SA2 update [9,l, a two level sphericalized correlation matrix with an r dimensional dominant subspace can be obtained from an unsphericalized matrix (5) as • (sAz) ----S(R;r) = ~ d/U/U~,

(7)

/=1

where Ul = [ul ... u,] = Wl Tl, U2 = [u,+, ... u,] = WzT2, dl =avv().., . . . . . A,) and d2 = ave(~+ t . . . . . A,). The main advantage of SA2 uvdating is that deflation allows the computation to be reduced to O(nr) as compared to O(n 3) for full rank-one eigenupdating. SA2 updating tracks the dominant (or subdominant) subspace more slowly than full eigenupdating, but the SA2 subspace is more robust to noise perturbations. The main disadvantage of SA2 is that it does not provide enough eigenlevel information to easily detect the signal subspace rank. The signal eigenstructure (SE) update [37, 50, 9,l has r distinct eigenvalues in the signal subspace and one sphericalized noise level in the noise suhspace, i.e., v4"!

@~SE~= $(R; 1,2 . . . . . r) ffi Y. d/UrUH,

(8)

where Ui -- u / - - wrt/for i = 1, 2 . . . . . r with [tr[ -- 1,

U,+~=W,+rr,+l=[u,+l

... u,,l, dr=~./ i =

1,2 . . . . . r, and d,+l = ave(~,+t . . . . . A,). The SE update provides signal subspace tracking performance that is very similar to full rank-one eigenupdating, but the computational cost is O(nrZ). If r is greater than or equal to the true signal rank, the SE update also provides enough eigenvalue information to detect the signal subspace rank. The SA4 update, a major focus of this paper, involves a four-level sphericalized correlation matrix, @~SA4~= S ( R ; r - l,r,r + I ) = ~ diU/U~, r=1

(9)

R.D DeGroat et aL / Signal Processing 50 (1996) 101-121

where U , = [ u , ... u,_,], U2=u,, /-/3=u,~,. 154 = [u,+2 ... u,], dl = ave(iL, . . . . . 2,_,), d2 = L, d3 = 2,+1 and d, = ave(L+2 . . . . . ;.,). The SA4 update tracks the signal subspace much like the SE update, but at a computational cost that is more like SA2, i.e. O(nr). Also, as we will show later, the SA4 update provides just enough eigenvalue information for M D L based signal detection.

3. SA4 algorithm development In this section we develop the SA4 algorithm. Suppose that at time k - 1 we have a four-level sphericalized EVD estimate of the correlation matrix. As described in (9), a four-level EVD estimate can be obtained from an unsphericalized EVD by replacing the first r - I eigenvalues by their average value, and by replacing the last n - r - 1 eigenvalues by their average value. At time k - - 0 we can initialize the estimate as such, and from then on each rank one update can be resphericalized to produce a new four-level sphericalized EVD estimate. Hence, at time k - 1 we have the four-level sphericalized EVD estimate,

boundary eigenvectors. Since the signal and the noise tracking algorithms are nearly identical in structure, we will focus our development on the sig,'al subspace version, i.e., we will explicitly track

[0, O~0~].

We update the sphericalized correlation matrix, ~(k - i), in the same way that it(k - I) is updated in (1). The only differences are that the EVD update of ~ ( k - 1 ) can be deflated because of the sphericalized subspaces and after each update, the modified matrix must he rcsphericalized. To simplify the notation, we will drop the k time index and absorb the (I - ~) scalar term in (I) to write the sphericalized version of the correlation matrix update as = ~ 0 + (1 - ~)xx" =

~ODO" +

(l -- "t).,~ H

= 0 ( z 6 + pp")O", = UG(2/)

+

lx/T~-~O"x (HI)

where G = diag(pl/IP, I ..... P,/I/~.1) is a diagonal unitary transformation that makes the matrix in-

[

a , ; , _ ,

~(k - l) = U(k -- l)/)(k - I)U(k - !)" = [U, U,: 03 04]

~2

a3

L from which we wish to obtain the four-level sphericalized EVD estimate, ~(k). Recall from (9) that /-)2 -- ~, and 03 -- d,+ ~ are, respectively, spanned by a single vector. In (10), a~ and a4 represent the average eigenvalues and the columns of 01 and 04 are the associated spherical subspace basis vectors. To avoid expensive and redundant computations, the algorithm will actually only update either [ 0 1 0 2 0 3 ] or [ 0 2 0 3 0 4 ] . Later in this paper, we show that ifr = p, [02 0 , ] is an estimate of the signal subspace span, a2 and a3 are estimates of the smallest signal and largest noise eigenvalue with U2 = d, and 03 = d,+ 1 being the corresponding

p=

7~.'T)GHOH, 7 = G"~,

vWl o. I •

i

(to)

o-j

side the parenthesis real [50] by factoring out the complex phase information contained in ~. (Note:. In this paper, the tilde will generally be used to indicate updated or modified quantities.) We can now concentrate on finding the eigendecomposition of the completely real matrix ( a / ) + ~r/T). Based on the structure displayed in (10) and (! 1), we can deflate the spherical subspaccs 01 and 04 with a block Householder transformation [9,34] 01 H =

0

0 1,

I2

0

ott, j

(12)

106

R.D. DeGroat et aL / Signal Processing $0 (1996) 101-121

whose subblocks Ht ~ C ('- t,~(,- l~ and N , e C ( ' - ' - U~(m-r- t) are selected to annihilate the top (r - 2) and bottom n - r - 2 elements of 3'. Next we apply this Householder transformation as

--- (OGH)HX(~ + 77T)H(HXGa0a),

(13)

where the inner transforms are multiplied into the middle term, and the outer transforms properly align the sphericalized portions of the signal and noise subspaces contained in OGH= [ 0 t 6 1 H l 0262 0363 O,G,H,] SO that a deflated diagonalization can be accomplished. The middle term in (13) now takes the form HH(ff./~ + 3'7TIH = |

F

, (14)

main advantage of this notation is that we can now specify individual vectors from the intermediate quantities, i.e. ~l°H) denotes the jth column vector of

0 ~°') -- D a Y .

As previously mentioned, the above algorithm can be made more efficient if the larger subspace, 0t or O., is not explicitly tracked. In the following we explicitly track Or, but 04 could be just as easily tracked. It can be shown that the effect of applying the Householder reflections to Dt and 04 is to align the boundary vectors, ~.~0.) - , - t and ~.(o.) - , + 2 , with the projection of the current data vector, x(k), into Dt and 04, respectively. Hence, if we want to track [01 (/'2 0s'l without carrying along 04, we need only compute the projection ofx(k) into 04 according to ~(GH) -r+2 : ~/~4 where 4, =(,,-[~,

which contains the core 4 x 4 eigenproblem, r = a

a2

+

a,j

~ [ft

f2 i3

and

i,],(15)

L~,j

where [~t 12 i3 i , 3 ~ [~,-t ~, ~,+t ~,+2] are the nonzero elements of ~ --HH3'. Note that the subblocks of H cancel on the scaled identity subMocks of~0 in (10). Now the problem reduces to diagonalizing the all real 4 x 4 matrix, F. Let {~e R 4~* be the orthogonal matrix that diagnnalizes F. Then if we embed it as [z,_, Q = / o o

L

o

o

~

o

o

~,_,_~

, (16)

we can finish the diagnnalization according to

=(OGH~)QnHH(od~ + ~,'r) H~(QaHnGHOn) = (OGHQ)J~(QnHaGaDU),

(17)

so that J~ is a diagonal matrix and the newly updated subspaces are given by the appropriate column partitions of 0 = OGHQ. To denote more easily intermediate quantities in the updating process, we will use the following notation: ~ o ) = O~G~ and ~ o . , = 0~6~H~. The

O, 0~][0, ~,

(ts)

O,]")x

i ~ = II011~ = I]xll ~ - Y~ - iff - i ~ .

To

re-

sphericalize the EVD estimate, we must reaverage the eigenvalues at = ( ( r - 2)at + ;lt)/(r- !) and a , = ((n - r - 2)a4 + ~la)/(n - r - 1) to restore the form of (10). The overall algorithm has complexity O(nr) and is summarized below. 1. Compute a truncated version of p, i.e.,

If= [0, 02 031"x-2. Compute *.runcated versions of: G = diag(flt/IPtl ..... P.+ diP,+,l) and ~ = G " 8 - [~t . . . . . ~,+t]. 3. Absorb the complex phase of~: [~G, ~ a , i~3~,] = E0, 0~ 0~3 a . 4. Compute

5. Using [Yl . . . . . 7,-t ], find the Householder reflection Ht which annihilates all but the last element, i.e., [3't . . . . . 3',- t 1Ht = I'0. . . . . 0,i1 ]. Compute ~ou)__ t~,~O)Ht and 371 = 11[3'. . . . . . 3',-,31l. 6. Compute % ffi 114'11.Note: ~2 = 3', and i3 = 3',+t. 7. C o r n. -. .-.t. -. . ~(os, 2 -_ 4 , 1 i , . 8. Construct the core matrix F according to (15) and compute the 4 x 4 diagnnalizing rotation, ~, so that b = diag(~t,a2,t/3,a~,) = 0 F 0 , x.

R.D. DeGroat e~ al. / Signal Processing 50 (1996) I01-121

9. Update the affected eigenvectors, [ u , - t , u , , ~ kmr- | ,~r , r+l, r+ 2 J ~l~. 10. Resphericalize the eigenlevels, ~t = ((r - 2)at + ~[t)/(r - 1) and ~,t = ((n - r - 2)~.~ + ,~4)/ (n - r - 1). We should mention that almost any decomposition or rank one updating method can be applied to the core problem in (15). For example, perturbation-based eigenupdating, URV, or one-sided SVD approaches could be used. In particular, SVD versions of the SA2 and SA4 updates are described in [11, 23], respectively. A systolic array implementation of the one sided SA4 update is given in [22].

stochastic algorithm ODE method of Ljung; for details see [39, 40]. This basic technique has been used to prove convergence of SA2 in [20], and the SVD version of SA4 in [23]. In this paper, we use this analysis technique to show that, under certain weak assumptions, the subspaces and the eigenvalues of an m-level SS update asymptotically converge with probability one to the desired quantities. We should remind the reader that an m-level SS update is really a rank-oue EVD update in which m subsets of the eigenvalues have been pteprocessed to create m spherical subspaccs. Thus, we are really using the ODE analysis to prove that recursively sphericalizing the subsl~ces o f a rank-one updated correlation matrix asymptotically produces the same subspace spans as an unsphericalized rankone updated correlation matr/x. Heuristically, this

4. Multilevel spherical sUbSlmce convergence Conceptually, it is easy to generalize the fourlevel SS update to an m-level update. The key idea is that for a rank-one update [5], only one vector, ~l = ~I ° m , out of each subspace, ~Gu~ = OIG~H~ & [Oi, ~l], is modified according to

[t~t . . . . . UM] = [~t . . . . . um]O,

107

(19)

where ~ is the orthogonal matrix from the deflated m xm core eigenproblem, / ~ = O F O T, and F = D + ~7~7 s is an m x m generalization of (15). Recall that Pt = i ~ ' ~ is the projection operator associated with the ith subspace. For analysis purposes, it is also helpful to realize that ~ = Pix/llP, xH and ~t -IIPlxl[ [9]. To understand the asymptotic behaviour of the m-level SS EVD update, we will focus on the projection operators associated with the subspaces of interest. Based on (6), there are an infinite number of orthogonai bases for a particular sphericalized subspace, Ui = JFiTi, but the projection operator associated with the ith subspace is unique, i.e., P~ = W t W ~ = UjU~ (see [34, p. 75]). Thus, we can characterize adaptive subspaec evolution in terms of projection operator trajectories over projection operator manifolds. The trajectory itself is characterized by an ordinary differential equation (ODE) derived by taking limits and expectations on the updating direction associated with the updating rule of the adaptive subspace algorithm. ODEs for the eigenlevels can also be simultaneously derived. To perform the analysis, we follow the basic

should not comc as a surprise. Sphericalizing once at *he end, or recursively sphericalizing the subspaces after each EVD update should asymptotically produce the same subspaec spans because an EVD by definition decomposes that data into maximal power ranked orthogonal components. An SS update and an unsphericalized EVD update both decompose the data into successively maximal power ranked orthogonal components, but the SS EVD update does this with orthogonal spans (or subspaces), whereas the unsphericalized EVD date decomposes the data into individual orthogonai vectors. Both types of EYD updates are moving towards the same goal (a maximal power ranked ortho$onal decomposition), but the SS EVD update partitions the decomposition into m u l t i d i n g n ~ n a l subspaces rather than the usual one-dimensional eigenvectors. Ljung [39, 40] states that a general recursive algorithm can be written y(t) = yit -- 1) t- g(t)Q(t;y(t - 1~ ~p(t)),

(20)

where y(- ) is a sequence of n,-dimensional column vectors, which are the objects of interest. In our case, we have to generalize y(.) so that it can represent a sequence of matrices and scalar quantities. The y(- ) are referred to as "the estimates'. T h e gain sequence, g(t), is assumed to be a sequence of positive scalars. The n,-dimeusional vector cp(t) is the observation obtained at time t. The function Q(.;.,.) is a deterministic function that satisfies

108

R.D. DeGroat et al. / Signal Processing 50 (1996) 101-121

some regularity conditions [39, 40]. This function, together with the gain sequence define the algorithm. Instead of a completely general dependence of ~o(t) on y(.), the following structure for 8enerating ~o(.) is used in [39,40]: ~(t) = A ( y ( t -- l))¢p(t -- l) + B ( y ( t - 1)) e(t)), (21)

where A(') and B(" ) are n¢ × n# and n~ x ne matrix functions. The input e(-) is a sequence of ne:. 1 random vectors. There are also nine assumptions, AI-A9, that must he satisfied (see [39] for details on these and other assumptions that can be made). Verification of these assumptions will he addressed as we 80 along. For an ODE analysis of an m-level SS update, we satisfy At if e(k) ----x(k) is a sequence of stationary, independent vectors with a positive-definite correlation matrix, R--E[x(k)x(k)n]. Stationarity is needed to satisfy A5 which will be discussed later. To satisfy A2, we must also assume that le(t)l < C w.p.l for all t. In practice, any measurement is bounded, so this assumption is not restrictive. For our algorithm, we have ~0(t)= e(t) so A ffi 0 and B ffi I which satisfy the continuity requirements of A4. The regularity conditions A3 and A5 involve the function Q(.;., .) which will be dealt with later. Assumptions A6-A9 involve the gain sequence which will also he addressed later. For an ODE analysis of an m-level SS update, the estimates of interest, y(t), consist of two sets of parameters: (1) the eigenlevels, a(k) --[at(k) . . . . . a~,(k)], and (2) the subsp~ce spans, ~/~(k), i = 1. . . . . m, or equivalently (and less ambiguously), the projection operators, P(k) = [P,(k) . . . . . P,,(k)] where P,(k) -- O,(k)l~,(k) u. Thus, we let y(t) in (20) he defined as ¥(t) = [~(t),P(t)].

Since the m-level SS update is really a rank-one EVD update (with spherici*.y enforced on certain subspaces), wc do not have an explicit expression for Q(-;.,.) in the algorithm itself. However, Q(.,., .) defines the update direction of the algorithm at each step. Since the ODE approach is used to describe the asymptotic trajectory of the algorithm, we can use asymptotic perturbation analysis to find an explicit expression for Q(.; .,.). Using the fact that an m-level SS update contains an m × m symmetric eigenproblem at its core, we can easily

derive perturbation results for the eigenlevels and the associated projection operators. See Appendix A. Now, following Ljung's ODE approach to analyzing a recursive algorithm (20) we write the subspace update in a generalized form, that is ¥(k) = Y(k - 1) + g ( k ) Q ( Y ( k - 1),x(k)),

(22)

where Y(k) ~ l i t (k). . . . . Y2m(k)] ffi [al (k) . . . . . am(k~ P l ( k ) . . . . . P~(k)] contains the m eigenlevels and their associated projection operators. For our algorithm it is convenient to partition (22) and look at each of the eigenlevel and projection operator updates individually as aj(k) = a~(k - 1) + g(k)Q~(k; Y ( k - 1),x(k)), i ffi 1 . . . . . m,

(23)

and P~(k) ffi P~(k - 1) + g(k)Q,(k; r ( k - 1),x(k)),

i = 1 . . . . . m,

(24)

where Ql(';', ") is the partition of Q(';-, ") corresponding to Yl(') and o(k) is a gain sequence that satisfies assumptions A6-A9. One choice of g(k) that satisfies A6-A9 is g ( k ) = 1/k [39]. By letting u ( k ) - - 1 - g ( k ) - - 1 I l k , so that 45-• (k)~ + (1 - ~t(k))xx s -- (1 - 1 / k ) 4 + (l/k)xx n, we ultimately can obtain (22), as we will show later. Ljung [39] uses the gain sequence to prove convergence by introducing fictitious time, z = Y:k= 1 g(k). The multilevel SS algorithm ODE uses • as its independent 'time' variable. The next step in Ljung's approach is to find the (appropriately partitioned) function f(Y(z))= [ft (Y(~)). . . . . fm(Y(z)), fro+, (Y(~)). . . . . f2m(Y(~))], where fi(¥(x)) ffi lim E x [Q,(V(k - l),x(k))],

k--Qv i = 1. . . . . 2m,

(25)

is the crux of the partitioned ODE, d

Y~(z) =fi(Y(O),

i = 1. . . . . 2m,

(26)

which describe, the asymptotic trajectories of the eigenlevels and the associated projection operators.

R.D. DeGroat et al. / Si;;nal Processing 50 fl996) 101-121 The problem is really to find an expression for

f(Y(z)). Assumption A5 is satisfied by the existence of f(Y(z)) together with Lemmas 4. I and 4.2 (which appear later in this section). Fortunately, (25) involves expectation and limiting operations, so ~t is not too diffi~,ult to evaluate. Let us first consider the eigenlevels. If we: let u= 1-s where the perturbation parameter = g(k) = l/k, then as k --, o0, we can look at Lhe eigenlevel perturbation when the sphericalized correlation matrix is updated as 4 = ( 1 - e) ~ + txx a. To maintain sphericity, the eigenlevels must be reaveraged after the rank-one update. The new ith sphericalized eigenlevel is then given by

3,(k) = (1 - e)3,(k - 1) + erd,(k) + O(e 2),

(27)

where the function fdY(O) is the right-hand side of (30). To apply Ljung's O D E analysis to the subspaces, we make a slight (but equivalent) reformulation of the algorithm in terms of projection operators, and then determine the corresponding Q~(.;.,.). F r o m (~,~ we see that a simple relationship exists between ~he sphericalized subslraces and the associated proiection operators, in fact, the projection operators

are ltzlly defined by the EVD of the spaericalized correlation matrix at times k and k - 1, i.e. ~(k)

1) + (1 - ~t)xxn; st . . . . . sin) =

= S(~tC~(k -

)"~= 1 aiP,(k) and ~ ( k - I ) = ~ = t a, P l ( k - !). We can write the projection operator associated with the tth subspacc at time k - 1 as P,(k -

1)=/.~°m(k-

where = £I,(/"

6d,(k) = 3,(k) - g ( k - 1) _ - ~ ( l , l,~.dl2 a,(k - 1)) + O(e2 )

is the first-order perturbation of the ith cigenicvcl, 11~7J112= t r ( H n / . ~ i ( k - l ) " x x n ( f i ( k l)H~)= tr(t/'i(k -- 1) "xxn/.~/(k - 1)) = tr(Pi(k - l ) " x x n P ~ ( k - 1)) is the amount of energy that projects into the ith subspace and mi is the multiplicity of the ith eigenlevel. Since for each eigenlevel, only one eigenvalue out of the m~ repeated values can change under a rank-one modification, dividing the change by m~(as we have done in (28)) is equivalent to reaveraging the eigenlevel. According to Ljung's definition of z. e = g(k) = l/k -~ dz as k ~ o~. Based on (28). (25) and the definition of z, we can obtain an expression for an infinitesmal change in the ith eigenlevel, d~,(¢) = ,-~limE = ( ( I =

i

(29)

da,(~)=(ltr(P,(~)"RP,(z))-~,(~)),

i ffi 1 . . . . . m,

i =

ID a

1 .....

m,

(31)

= t f i G i H i , i = 1 . . . . . m, ~ is

the portion of ~Ga~ that is unchanged by a rank one EVD update and ~i is the one vector in ~ r m that is going to be altered by the m x m orthogonal transformation ~ which cont~-ins the eigenvectors from the core eigenproblem in the deflated update (19). After updating the appropriate components of each of the m subspnces as described by (19) the updated projection operators are given by

P,(k) = U~m(k)(U~Gn)(k)) n = ~l,~ly + ~ , i = 1. . . . . m.

(32)

For analysis purposes, we are interested in the projection operator update; hence we write P,(k) = P , ( k - 1) + A,(k),

i = l . . . . . m,

so that d d k ) ~- P~(k) - P d k - 1),

(33)

i = I . . . . . m, o r i = I . . . . . m.

(34)

which yields the eigenlevel ODEs, ODE(i):

~H ,

d,(k) = ~ , ( k ) ~ ( k ) - ~,(k)~n(k),

H,~,(k),[z - a i ( k ) ) ] d z

tr(P~lz)HRp~(~)) - ai(~) dz

I)(O}°m(k-

+ uiz,

where [17i ui) ~ ~ m

(28)

109

(30)

Comparing (24) with (33) shows that g(k)Qi(g(k), x(k)) ~-ddk) is the data-defined stochastic update direction for the ith orthogonal projection operator. Since g(k) Qi(g(k),x(k)) is the update direction as defined in (34), from (34) and the definition of we h a w dPdz) = lim E x [ ~ d k ) ~ ( k ) - ~dk)~"(k)]dz. tt~l-

(35)

R.D. DeGroat et al. / Signal Processing 50 (1996) 101-121

I10

We will use perturbation analysis to obtain the desired result (see [35, pp, 308-314]). In (19), the m × m core rotation, ~ (which diagonalizes the core ¢igenproblem), will update ~dk) to ~dk) according to ~i(k) = qle~l(k) + q2eu2(k) + "'" "4" qme~m(k),

(36)

where qo are the individual elements of ~.. For an infinitesimal (~<<1) update, i.e., 4 = ( 1 - ~)~ + ~ x a, it can be shown that to a first-order approximat[an [35], qo m ~

-k O(g2) for i~j---- l . . . . . m

(37)

with qee = 1. Substituting (37) into (36) yields

-

(38) Substituting (38) into (34), (34) into (35), replacing by dT, and ignoring higher-order terms that vanish in the limit of (35), we obtain dPdz) = lira E~ L #~._ ?~j(k)~:(k! ( ~ d k ) ~ ( k ) ",. . . . 1 de(k) -- dj(k) + ~,(k)~P(~))] d~.

(39)

for j :~ i ffi 1. . . . . m,

(40)

Letting 1

a,(~.)

-

a,(k)

and using the fact that ~e= ;'ex/llP, xll and Ye= IIP~xll [9]. we can write

f~÷.(r(O) = lim

~.

2vo(k)E.[~,(k)f,~(k)O~dk)~(k)

+ ,~Ak)l,s(k))] = lira

~.

2vo(k)F.,,[Pdk)x(k)x(k)aP.dk)

k~m j # i = l

+ Pj(k)x(k)x(k)nPdk)] =

~ j#e=!

2vo(¢){Pdz)RP~(z) +PjRPdz)},

ODE(i + m):

+P~(z)RP,(~)),

i = 1. . . . . m.

(42)

The ODEs in (42) and (30) could be stacked together into one massive first-order ODE, but we

#~, e'Te(k)'7~(k) ~jtk) + O(e2). ~e(k) = li(k) + J e= t a,(k) aAk)

vo(k ) A_ ~

where Pd¢) = Ud~)U~(¢) for i = 1. . . . . m. (Note: Eqs. (41) and (30) satisfy condition A5.) It can also be seen that the asymptotic expressions for Q~( ) are continuously differentiable w.r.t. ¥ and x which satisfies A3. A non-asymptotic proof of A3 can be easily derived from the partial derivatives given in [12]. Using (26), (35) and (41) the ODE for the ith proj~ction operator becomes

(41)

will find it convenient to look at each projection operator individually. Note that if m = 2, and we assume that the eigenlevels are fixed, this set of ODEs simplifies to the SA2 ODE described in [20]. Now that we have a set of ODEs for the m-level SS update that dcscrit~s the trajectories of each eigenlevel and its associated subspaee projection operator, we need to show that each projection asymptotically converges w.p.1 to the true subspaee projection operator and, similarly, the eigenlevels converge to the true eigenlevels. As we will see, the eigenlevel convergence is dependent on the subspace convergence, but not vice versa. In fact, the eigenlevels can be set to fixed values and the subspaces will still converge. However, the eigenlevel values affect the speed of subspace convergence. To prove subspac¢ convergence, we first show that the dominant subspace of R is the global stable attractor point of ODE(m + I). Then we show that the ith true subspace becomes a global attractor for ODE(0, i = m + 2 . . . . . 2m, when the more dominant subspaces approach the true subspaces closely enough. In other words, the more dominant subspace must approach convergence before the less dominant subspaces can converge. To begin, note that ODE(i + m) in (42) goes to zero if Pale) -- Pe -- Ul U a, i -- 1. . . . . m, where Pj is the true projection operator defined by the columns of /-]1¢ C ' x ' ' which form an orthonormal basis for the span of the true ith subspace of R. To see that this is true, write the EVD of the true

R.D. DeGroat et aL / Signal Processing .TO (19~,f) 101-121

correlation matrix as R -- 5"~'=t W s A s ~ . Now since span(U,)= span(We), we can substitute this expr,'ssion for R into (42) with Pdz) = U~U" and note how all terms go to zero. Now that we know the ODE has stable points at Pc(Z)--U~U~", i = 1. . . . . ,'n, we must analyze these points further to determine if they are global attractors. To do this, define the Lyapunov function V (PdO) = tr( Pd~)RP~(O ).

(43)

Define the raak r, n-dimensional projection.~ operator manifold, p~., A_ { p l p ~ C . . . . p2 = p = p a , Rank(Pj = r}.

III

Proof. For P~(~) to be in P,,,., a necessary condition is that Pd~) satisfy Pd¢) : (Pt(¢)) z. Differentiating both sides of this static constraint gives the dynamic constraint,

Pl(z) = PI(z)Pi(Q + P,(T)P,(T).

(47)

It iS easy to substitute O D E ( / + m) of (42) for Pe(z) on the right-hand side of (47) to verify that (42) satisfies (47~. This shows that the algorithm ODE will propagate a projection operator quantity, it remains to be shown that the ODE will not admit a trajectory that leads from a rank m, projection operator to a rank q projection operator where me =/=q. To do this note that after d~ 'time' units, the projection operator has moved

(44) Note that for the positive-definite matrix R (see (6) for the definition of si and m~) max

V(Pi(z))=tr(UiU~RUeU~I ) =

P,(r)ePm,..

~ )=sl- 1 +

~.s, I

(45) with Pe(r) subject to the constraint

Pe(~)PS(O = PS(~)P,(~) = 0, j=l

..... (i-l)

and

i=2 ..... m

(46)

and Pz (T) is unconstrained. The maximal points are obtained when Pe(~) -- UeU~ which shows that the stable points are unique maximum points provided A,, > ).,,+ t, i.e., the boundary eigenvalues of R are distinct. A potential problem with the above analysis is that the ODE might take P~(T) to a point off the rank m~ projection operator manifold. It is clear that points off the manifold could increase the value of V(Pd¢)) well beyond the constrained maximum of (45). The following lemma shows that ODE trajectories that start on a rank me projection operator manifold will always stay on the manifold. The lemma is needed to support the subsequent Lyapunov function based convergence theorem. Lemma 4.1. I f Pi(z) is a solution to (42) with initial condition PdO) ~ Pmj,n then Pe(~)¢ P~,.. for every ~ E~o, oo).

dPi(z ) =

~.

2Vu(T)(PdT)JibwJ(T)+ Pj(T)JI[Pe('c)dT).

j~i=l

(48) Using the Frobenius norm we have IlPdOIIv = for each Pdz)e P ..... Now suppose the trajectory took Pe(T) to another matrix, P e pc.,, q ~ me. Then IlPlle = ~ ~ ~ implying that IIPe(~) + dPiIIF ----IIPIIF-Butsince dPi is vanishingly small, IIPdr) + dPilJF = IIPd¢)IIF -- v / ~ l , Thus, it is impossible to jump from a rank m~ projection to a rank q =~ me projection while moving along a continuous trajectory that is constrained to lie on a projection operator manifold. ]n terms of the Frobenius norm, the P.,.. manifold lies on an mr radius hypersphere while P,.. ]ks on a q-radius hypersphere. The idea is that you cannot take a 'quantum ]cap' from one projection operator hypersphere to another by adding an infinitesimal perturbation. I-1 Lemma 4 . 7 . I f PdO)rj(o)=o and Pd~), i = l . . . . . m is a solution to (42), then Pd~)Ps(z) = Ofor ~ >~0. Proof. Differentiating both sides of the static constraint, Pi(~)Ps(~) = e, gives the dynamic constraint P,(r)Ps(r) + Pd~)PS(~) = O.

(49)

To see that (42) propasates orthogonal projection operators, we substitute (42) into (49) and note that

112

R.D DeGroatet al. / Signal Processing50 (1996) 101-121

by I,emma 4.1, certain terms cancel and the remaining terms become zero, i.e, (49) is satisfied, when Ps(T)Pj(~) = 0 for • >I 0. Thus, if the projection operators are initially orthogonal, (42) will propagate orthogonal projections. I'-1

Theorem 4.1. The estimated etoenlevels and subspaces of the m-level rank-one SS update asymptotically converge with probability one to their true values as defined in (6),i.e.

Note: The above analysis assumes infinite precision and that the subspaces are orthogonal. With finite precision, the subspaces can become nonorthogonai due to accumulated roundoff error. Consequently, if large numbers of updates are to be computed, a computatlonally efficient, partial orthogonalization scheme should be used to maintain working precision orthogonality [15,41]. If machine precision orthogonality is maintained, i.e., IlfP(@~(k)")fp(Ds(k)) - l,,,ll < O(~) where fp(A) denotes a finite precision representation of the infinite precision quantity, A, then it is not difficult to show that liPs(k) - fP(Pdk))ll < O(~). For Bunch's algorithm (which is probably less stable than most EVD updates), analysis shows that the orthogonality error (and the projection operator error) after a rank-one update, but before reorthogonalization, is upper bounded by O(v/~) [15]. However, the orthogonality error can be brought back down to O(~) with an efficient, partial orthogonalization scheme ['15, 41]. Thus, the projection operators of the finite precision algorithm will remain within O(e) of the desired projection operator manifold as long as working precision orthogonality is maintained. Moreover, it is not possible for the finite precision 'projection operator' to jump from one rank to another since at no point in the finite precision update (either before or after reorthogonalization) is the projection operator error O(1) (which would indicate a possible change in rank). In other words, the finite precision algorithm will satisfy Lemmas 4.1 and 4.2 to within O(e), if working precision orthogonality is maintained. This is the best that we can hope for. Technically speaking, the proofs in this paper only apply to the infinite precision SS algorithm. A more complete finite precision analysis of the algorithm is a topic of future research. With the above results established, we are in position to prove the multilevel SS convergence theorem.

and

a~(k)--,ds

a s k - - , oo, w.p.l, i = l . . . . . m,

(50)

[Is(k) - ' Us = Wi Ti as k --, ~ , w.p. 1,

i = 1, ... ,m,

(51)

provided Ljuno's conditions AI-A9 are satisfied [39]. These conditions can be satisfied by assumino that input vectors are independent and stationary with a positive-definite correlation matrix, R = EEx(k)x(k)"], IIx(k)ll is bounded and the fadino factor is oe(k) = 1 - Ilk.

ProoL We start by showing that PI = UI U H, the projection operator defined by the true dominant subspace, is a global attractor point for PI(T) = UI(r)U~I(~) in ODE(m + 1) of (42). This is accomplished by showing that the Lyapunov function increases in the direction of any trajectory leading to PI --- Ut U~ . This, together with Lemmas 4.1 and 4.2 show that the algorithm dynamics are such that the expected trajectory Pt(z) converges to Pt -- U~U~. To show this, we take a derivative of the Lyapunov function along an arbitrary trajectory defined by (42). This turns out to be a straightforward calculation, d"~ v(Pl(z)) =

(trlPllT)RPl(z)))

=tr[dps(t)RPs(z)+P,(z)Rdp,(z)]

=

~.

2V,l(z)tr[Pil~)RP~(z)RP,(~)

j~i=!

+ Pj(x)RP~(z)RP,(x) + P,(~)RPi(z)RP~(~) + PI(z)RPj(T)RP~(z)].

(52)

For i = I, (d/dr) V(Pl(z)) I> 0 with equality only when Pt (3) = Pt, the true dominant projection operator. This is true because each of vu-terms is

R.D. DeGroatet al. / SignalProcessing50 (1996) 101-121 positive and the trace is positive when Pl (~) ~ P, because the true correlation matrix R = Y.~'=t djP~is positive definite. More specifically, the first and last terms inside the trace are nonnegative definite, e.g., zuPd~)RPj(~)RPdT)z=(zaPd~)RUdz)) (U~(z)RPdt)Z) ~ 0 for any nonzero z with equality only when Pt (z) -- Pt. The trace of these terms is therefore nonnegative. The trace of the middle two terms is zero because tr(Pl(z)RPi(~)RPdz))= tr(Pd~)Pj(~)RPdz)R) = 0 since Pd~)Pj{~) = 0. Thus, PI = UI U~ is a global attractor point for Pl(Z) = Ut(T)U~(z) in ODE(m + 1) of (42). For i = m + 2 , one of the v2f terms is negative, i.e, v21(T) -- l/(~2(z) - at(r)) < o. However, (d/dz)V(Pz(z)) > 0 for some ~ such that IIPx(~)- PI II < ~ because the negative term in the sum goes to zero. In other words, P2(r) = UzUa2 becomes a globally attractive point for ODE(m + 2) when PI (~) gets close enough to the true operator Pj = UaU~. This line of reasoning can be continued for the progressively less dominant subspaces. Thus, the subspaces of the m-level SS update asymptotically converge w.p.l to the true subspaces. To prove that the eigenlevel ODEs in (30) converge to the true eigenleveis, d~, i = 1. . . . . m, as defined in (6), we must show that d, i = 1. . . . . m, are global attractor points for the ODEs in (30). First, it is rather trivial to see that if the subspaces converge to the true subspaces, the true eigenleveis are stable points for (30). To show that these points are globally attractive, we define the Lyapanov function, V(~i(z)) = ((l/mi)tr(Pi(r)aRpi(~)) -- ~i(z)) 2. Clearly, as • - , oo, the projection operators converge to their true values and (1/mi)tr(Pd¢)nx RPdz)) --, d,, which implies V (~i(~)) -* (d, - ai(~) 2. Thus (d/dz) V(ad~)) = (c3V/O~) (~/~) = - (d~ - ~(z)) 2 ~< 0 as • ~ o~ with equality only at the true eigenlevel values. Consequently, the true eigenlevel values become globally attractive when the subspaces have sufficiently converged. Thus, the eigenlevels of the SS update asymptotically converge to the true values w.p.l. []

Theorem 4.2. The estimated eigenlevels and subspaces of the m-level rank-one SS update asymptotically converoe with probability one to their true values as defined in (6), provided Ljuno's conditions C 1-C7 are

satisfied [39]. The.w conditions can be satisfied provided x(k) is ergodic and random enough that x(k)n Pt(k)x(k) ~ O. w.p.l, and the fadino factor is ~t(k) = I - I/k. Proof. Due to space limitations, the proof is given in [12]. (Note: This theorem allows for Gaussian data which is needed in the next section.)

Corollary 4,1. The SE update converges to w.p.I to the true eigenvectors and eioenvalues of the siflnal subspace. Corollary 4.2. The SA4 update converoes w.p.l to the true subspaces ond eioenlevels.

5.

SA4-MDL malt I r a c k i ~ system

In order for the SA4 method to track time° varying subspaces for high-resolution algorithms such as MUSIC, minimum norm, and total least squares, it must also track the number of signals° In this section we present an M D L detection scheme [55] that is based on the SA4 update (SA4-MDL). In the next section we show it to be strongly consistent. SA4-MDL uses the four ¢iget~levels retained by the SA4 algorithm to deride if the subspace dimension should be increased, decreased or remain unchanged, or if there is no signal at all. At least four levels are needed to provide a single, well-defined eigenvector on both sides of the signal and noise subspace boundary. A good estimate of the boundary eigenvectors allows the signal subspace to be increased or decreased in a well defined manner. Let ,~i, i = 1. . . . . n, he the eigenvalues of the sample covariance matrix k(N) = (I/N) y.s= lx(k ) x(k) a and let x(q) he the number of degrees of freedom associated with the correlation matrix in a rank q estimation problem. Then the classical MDL cost function is given by [55] MDL(q) = - Nlog(L(q)) + ½;t(q)!og(N~, q = 0, ... ,n -- 1,

153)

R.D DeGroat el al. / Sign¢l Processing 50 (1996) 101-121

where

e(q)

I-I:=¢+, ~,

L(q)=A(q) ((l/(n-q)Y.'~=,+,,~,)'-" q : 0 . . . . . n - 1,

(54)

and N is the number of observation vectors. Ideally, the cost function should be minimized when its integer argument, q, equals the signal rank, or p in (2). If the eigenvalues are distinct and q is the assumed signal rank, then the number of degrees of freedom are given by z(q) = q(2n - q) + 1. MDLbased rank estimators simply search for the integervalued q which minimizes MDL(q) and then use the minimizing q as an estimate for the signal rank, p,

l0 can still take on any legitimate M D L value, i.e., 0 ~
i= 1,...,r-

h, = a2,

h,+ 1 = a3,

hi=a4,

i=r+2

1, (55)

. . . . . n,

so that the SA4-MDL cost function is given by

i . e . 1~ = q . i . .

M D L must be modified to work with SA4. Since all of the cigenvaluc information is not retained by SA4, (53) cannot be directly applkxi. Recall that for SA4, at is the average cigcnlevel associated with the r - 1 dominant basis vectors contained in ~s, a, and as are eigenvalue estimates for the rth and (r + l)st eisenv~tors contained in Dz and D3, respectively, and a,, is the eigcnlevel associated with the n -- r - I subdominant basis vectors contained in 19,. Since we assume the dominant ['i91C/=] and subdominant [D3/.Y,,] subspaces have exactly two dgenlevels associated with each of them, r is restricted as 2 < r ~ n - 2. The SA4-MDL criterion may not directly generate a signal rank estimate, 1~= qms,, but rather it can generate one of the following outcomes based on the four ¢igcnlevel estimates: 1. Do not change dominant subspaee size (al I> a~ > a3 ~ a,). Let ~ : r and r = r. 2. Increase dominant subspace size (a!/> a , / > as >a4). Let 10= r + 1 and r = r + 1. 3. Decrease dominant subspace size (a I > az a3 w, a4). Let,0 = r - 1 and r = r 1. 4. No signal (as ~ a2 ~ a3 ~ a4). Let /~= 0 if • = 2, otherwise let ~ = r - 1 and r = r - 1. Thus, when • > 2, it is only necessary to test for three possible outcomes, i.e., q = r - 1, r,r + 1 and when • = 2, we test for four possible outcomes, i.e., q =0,1,2,3. Although r is restricted as 2 < • ~ n - 2, the SA4-MDL method can step its way towards the global minimum, ~ = qmin and

MDLsA4(q) -- - N Iog(LsA4(q)) + ½ZsA4(q)log(N), q = 0, r -- I r,r + 1,

(56)

where BSA4(q)

LsA4(q)

l-l~=q + 1 hi

ASA4(q) /(1/(n -- q) Y-,=q+ " ~h ' x,,,-¢, ) q ----0,r -- l , r , r + 1,

(57)

and XsA4(q) represents the number of degrees of freedom in the SA4-MDL cost function. Since there is only one degree of freedom in the r - 1 repeated eigenvalues of the dominant subspace, we must subtract r - 2 from the total number of degrees of freedom in the SA4 EVD, i.e., !, ZSA4(q)=

q = 0,

q(2n -- q) + I -- (r -- 2). q = r -- 1, r,r + I.

(58) For the exponentially faded correlation matrix estimate in (1), the effective correlation window length can be approximated by 1

N ~ 1~"

(59)

Thus, we can use (58) and (59) in (53) to ob~,ain an M D L criterion fez SA4.

!15

R.D. DeGroatet al. / SignalProcessing~0 (1996) 101-i21

6. gA4-MDL strong consistency

Defining

In this section we demonstrate the SA4-MDL detection scheme to be strongly consiatent in a stationary signal environment. The proof is based on the result in Section 4 that the eigenlevels estimated by SA4 are strongly consistent estimates of the true (averaged) eigenlevels. Theoretically, we show that the strong consistency of SA4-M DL holds under an assumption of stationarity. However, simulation studies show that SA4-MDL is practical for slowly

( (l/(n- r-

•

I))~=,+2~1

)....

C

,

(62)

[!./=,+ z ,~l we can verify that fi (~...~._r_.~ =~, z ,~t) ' - ' - 1 BSA4(r + 1) = hi = i=r+2

= C ~ II i=r+2

time-varying environments too.

i

£~,

(63)

or

Theorem 5.1. The M D L criterion based on the SA4 estimated eigenlevels is strongly consistent provided the assumptions of Theorem 4.2 are satisfied and, additionally, the data is Gausslan.

-BSA4(r - ~ !} = C. B(r + 1)

Proof. We will show that MDLsA4(q) differs from

BsA,,(r)= hr+ 1BsA4(r + 1) = C B(r) ~,+~B(r + 1)

MDL(q)onlybyaconstantatqfr-l,r,r+l.

From Theorem 4.2, the SA4 estimated eigenlevels are the strongly consistent estimates of the true eigenlevels. Therefore, asymptotically, we have al = ~

1

(6o)

/-,~ ~j,

,_q AsA4(q) = ( ~ - ~ _q =~+ hi) t l

=

I ( '+l~ h , + = ~ + 2 h , ) ) ' - ' n----'q i=q+ 1 i

A(q).

h, BsAa(r)

= c.

(~)

q=r-l,r,r+l,

(67)

~

LsA4(q) = C , L(q)

MDLsA4(q) = MDL(q) - N log C -- ½(r -- 2) log N,

where ~i, i = 1 . . . . . n, are the eigenvalues of the unsphericalized correlation matrix estimate, k. For q = r - 1, r, r + 1, we have

1 =\n-~q\i=q+l

i~

or, equivalently,

~4 ~ n - - r - - l i ~ r + 2

(

Bs^4(r -- i ) Hence,

a~ = ,~,+,,

a~ = ~,,

(65)

and ~7;:

r-I ~,~ ~l,

l

(~)

Furthermore,

q = r - 1, r,r + 1.

space t o w a r d s r f 2 ,

~i+

i (61)

(68)

We must consider the q = 0 case separately. When r = 2, LsA4(O)= CL(O) which is consistent with (67) and (68). When r > 2, we may assume that qmi, > 0 and tbe q = 0 case should not be evalu° ated. If by chance, q,.l. = 0 when r > 2, MDI.,sa,t will step the dimension of the SA4 dominant suband then evaluate the q f O

case when r = 2 is reached. Since SA4-MDL only uses the M D L ~ , function a t q = r - - l,r, r4. 1 a n d q - - 0when r - - 2, and for these values, MDLsA, differs from the classical MDL ft:action only by a constant, the strong consistency of the SA4-MDL method is derived from

R.D. DeGroat et al. / Signal Processing .$0 (1996) 101-121

the strong consistency of the classical MDL method [62"1. If the initial estimate of the signal subspace rank,/~, is incorrect, MDLsA4 will (asymptotically) correctly tell us that we need to increase or decrease the size of the signal subspace. When the subspace size is stepped up or down by one, the asymptotic assumption for strong consistency is no longer met, and we must wait for SA4 update to reconverge enough so that an asymptotic state is achieved. Once the asymptotic assumption is met again, MDLsA4 will correctly tell us if we need to take another step. In this way, MDLsA4 will eventually step its way towards the true global minimum of the MDL function. []

misadjustment thereafter. The SNR is 20dB for all signals in Figs. 1-3. In this simulation, the data dimension (or number of .sensors) is set to n = 10 (similar results are obtained with n = 50), and the rank of the signal subspace is assumed to be known a priori as/; -- r = 4. Thus, the SA4 signal subspaee consists of a three-dimensional sphericalized subspace, ~)t, and a separately tracked fourth

oh

.

.

.

.

.

.

.

.

0.s

Hence, we see that the SA4-MDL asymptotically converges to the correct number of sources under stationary conditions with a diminishing gain sequence. And, in the next section, simulations show that SA4-MDL also works well m;der slowly timevarying conditions and a small but fixed gain.

7. Simdsfion results

We now present simulation results that verify the analysis in this paper. In fact, the simulations demonstrate that some relaxation of the assumptions (A1-A9 or CI-C7) used i , the ODE analysis is possible. Specifically, we use time series Toeplitz data matrices to generate the input vectors, so assumption AI is not satisfied because the input vectors are not independent. Nevertheless, t~e simulations work as the ODE analysis predicts. The simulations also demonstrate the ability of the SA4 tracker to track the signal subspace and the signal rank in a nonstationary environment. In the first set of simulations the SA4 algorithm is initialized with a spheric~lized correlation estimate 4~o = 0.011. The tracker is then turned on with an input consisting of four complex sinusoidal signals. To obtain a diminishing gain sequence, g(k) ffi ( 1 - ~t(k)), that asymptotically approaches zero, we use a fading factor that starts out small and gradually builds up to one, i.e., ~t(k)= 1 - l/(0.2k + 1). This particular choice of lading provides a good balance between rapid convergence at the beginning and a gradually diminishing

Fig. 1. ConvergenceofSA4frequencyestimatestotrueftequencies, with Oo = O.OIL a lading factor ore(k) = I - l/(0.2k + I), SNR = 2 0 dB for all sig~.als, n = I0, ~ = r = 4 , ]'1 = - 0 . 1 , f2 = 0.05,f3 = 0.17 and f4 = 0.2.

- -

uJ

--

U_2

.....

IU_S U_q

nm~.k Fig. 2. Convergence of SA4 subspace estima:es to true subspaces, with same signal scenario as Fig. I.

R.D. DeGroat et al. / Signal Proccss#,g 50 ll996) I01-121 0s¢-

117

,

oa

.

..

: i

:

;

o,

i i

i

" !

.., i. .i:.

l°.;t-it~........ ...~--.:L ... ~:-:~!~i=-i---:i~:: :..~ ' ~i-:-t ]

/ il il o,

Fig. 3. Convergence of SA4 eigenlevelestimates to true eigenlevels, with the same signal scenario as Fig. 1. The true eigenlevels are indicated with a straight horizontal line. component, /~'2, that should asymptotically approximate the fourth d o m i n a n t eigenvector of the true correlation matrix. The frequencies of the signals are normalized to lie in a range of - 0.5 to 0.5, and are given by ,I"1 = - 0.1, f2 = 0.05, f3 = 0.17 a n d f4 = 0.2. Fig. 1 illustrates the convergence of the frequency estimates. N o t e that within approximately 50 samples, SA4 has locked o n t o all four frequencies. Fig. 2 illustrates the convergence of the subspaces. We measure the distance between subspace U and V as dist(U,V)= IIUU"-VVHIIF. Three curves are plotted, but only two are visible. The solid curve displays IIUtUP -utU~llr, the dashed curve displays I I t ; , ~ - U+U~IIF= Ila,+," - u , u , " l I v and the dotted line displays convergence of the noise subspace, I'/Js/J4]. This curve falls primarily on top of the U2 subspace convergence curve. These two curves are similar because the bulk of the signal and noise subspace error occurs between these two subspaces, and as one corrects, so does the other. We d o not plot the subspace distance convergence curves for the d o m i n a n t and s u b d o m i n a n t c o m p o n e n t s of the noise subspace since the true noise subspace is spherical, the individual eigenvectors, u++t and Ur+,, are not uniquely defined. Fig. 3 shows the convergence of the four SA4 eigenlevei estimates, ~+-*di for i= 1 . . . . . 4. The straight horizontal lines represent the true eigenlevel values.

o,

illllllllll i oo

°o,.

,2.,,,

illl ,..

,.

.,¢

Fig. 4. SA4 tracking scenario with nonstationary sources and cha,ging signal rank. Tl~ starling and stopping frequenciesare indicated v~itha large "X', SNR = 10dB for all signals, n = 10 and :¢ = 0.95. The frequencies are estimated wit SA4 baled root-MUSIC and the signal rank is estimated with a low pass filtered SA4-MDL signal rank estimator.

Note that since the noise is white, a3 and a4 both converge to the same level, hence we only see three distinct convergence levels. Figs. 2 and 3 graphically illustrate the O D E - b a s e d convergence theorems. Fig. 4 illustrates a nonstationary tracking scenario involving SA4 updating and S A 4 - M D L signal source detection. To avoid startup transients, we initialize this simulation with 1000 d a t a vectors and a correct signal rank estimate o f ~ = • = 2. The S N R for all signals in this simulation is 10dB, the data dimension is set at n = 10 and the fading factor is fixed at ~, = 0 ~3. The number of signals starts at 2, builds to 5, ramps d o w n to 0, then back up to 2. As many as five signals are present at one time, two of which are nonstationary, in this simulation, the tracker must adaptively estimate ~, the estimated number of signals. Simultaneously, the /~-dimensional signal subspace is adaptively estimated and used with root M U S I C to estimate the frequencies. In this figure, frequency estimates and signal rank estimates are plotted on the same graph. The left vertical axis labels the normalized frequency range ( - 0.5 to 0.5.), the right vertical axis labels the signal rank estimate, and the horizontal axis is the time index. At time zero, two

118

R.D. DeGroat et al. / Signal Processing 50 (1996) 101-121

tracks are present starting at normalized frequencies off~ = 0.15 and f2 = 0.4. These tracks have frequency slopes of A f t / A k = O and A f 2 / A k - - -0.000004. The dark lines represent the frequency tracks and the integer-valued stairstep function indicates the SA4-MDL source number estimate. We see that at time k -- 0, SA4-MDL correctly estimates the source number to be two. At time k -- 1000, signal three turns on with an initial frequency off3 -- - 0.2 and a frequency slope of A f s / A k = 0.000004. A large 'X' is used to indicate the actual starting and ending points of the signals. Between the time a signal is turned on or off and when the change is detected by the SA4-MDL method, there is a small delay (on average, of 21 updates for signals turned on and 59 updates for signals turned off). It takes less time for SA4-MDL to detect an eigenvalue that is rising out of the noise level than it does for one that is fading its way towards the noise level. Signals four and five turn on at times k = 3000 and k = 5000 with start fiequencies f4 ffi 0.5 and fs ffi - 0.3, respectively. These signals are stationary, that is A f , , / A k = A f s / A k ffi 0. Note that when new signals are added, the S A ~ M D L stair-step function increments, and when a source disappears, it decrements. Near k = 3000 and k ffi 5000, SA4-MDL makes temporary errors which appear as 'glitches' (or spikes) in the stair-step curve. In order to smooth out the rank estimates, we used a single pole lowpass filter, H(z) ffi 0.05/(1 - 0.95z-t), with the output rounded to the nearest integer. The SA4-MDL method without rank smoothing would have yielded about nine "glitches', but with smoothing, only two 'glitches' remain as shown in Fig. 4. Of course, the price of smoothing involves an additional delay in detection. For this particular ~mulation, the additional delay amounted to 13 updates. Our experience generally indicates that SA4 performs more like a full EVD update. This is probably due to the fact that the subspace evolution is primarily governed by the closely spaced eigenvalues on the border of the signal and noise subspaces. Since these two cigenvalues in SA4 asymptotically converge to the true eigenvalue~, the SA4 subspaces evolve in a manner that is similar to the full EVD update.

8. Conclusion In this paper we developed and analyzed multilevel spherical subspace (S$) tracking algorithms. These algorithms are based on sphericalized eigendecompositions. Given a correlation matrix and its eigendecomposition whose eigenvahies are ordered from largest to smallest, the sphericalized estimate of the correlation matrix is obtained by setting groups of adjacent eigenvalues to a constant (usually average) value. Once done, the hypereilipsoid associated with the subspace becomes a hypersphere, and the subspace is said to be spherical. An interesting point is that while the hyperellipsoid has changed, the subspace spans have not. So DOA and frequency estimates obtained from high-resolution algorithms such as MUSIC, minimum norm, and TLS are invariant under eigendecomposition sphericalization. The advantage to using spherical subspaces is that the subspace bases can be reoriented to deflate the update computation. The basic idea of spherical subspace updating is to start with a sphericalized eigendecomposition estimate, make a rank one eigenupdate of the sphericalized estimate, and then resphericalize the updated estimate. This seemingly ad hoc approach to reducing computations stands up to rigorous stocha~:ic analysis. Ljung's ODE approach [39] shows that in a stationary environment the subspaces of an SS update asymptotically converge with probability one to the subspaces obtained from the true correlation matrix cigendecomposition. For many DOA estimation methods, the subspaces are all that is needed to obtain high-resolution DOA estimates. Yet the eigenvalues are needed for source detection. Hence, if the eigenvalues are sphericalized to just two levels, as in SA2, the algorithm will not produce enough information to easily estimate the number of sources. This was the primary motivation for deriving the SA4 update and other multilevel updates. In particular, it was important to show that SA4-MDL produces strongly consistent source rank estimates. Possibly, the most practical contribution of this paper was to develop and analyze the SA4 algorithm. This algorithm requires O(nr) computations and uses internally generated eigenlevel

R.D. DeGroat et aL / Signal Processing .50 (1996) 101-121 m f o r m a t i o n to track the signal rank w i t h o u t a d d i n g significantly to the algorithm's complexity, T h e analysis showed that the a l g o r i t h m is a stochastic a p p r o x i m a t i o n that converges to the true s u b s p a c e s with probability one. It also showed that the rank estimates are strongly consistent. T h e O D E m e t h o d s previously used to analyze adaptive filters are extended to projection o p e r a t o r manifolds to characterize s u b s p a c e convergence. T h e s t r o n g consistency results were based on the O D E based convergence theorem. While the analysis holds only for stationary envir o n m e n t s , we see via simulation that the results hold quite well for slowly time-varying environm e n t s as well.

Acknowledgements W e w o u l d like to t h a n k the reviewers for helping to i m p r o v e this p a p e r significantly.

Appendix A W e state the following l e m m a which deals with the effect of p e r t u r b a t i o n s o n the eigenlevels a n d s u b s p a e e s o f a sphericalized H e r m i t i a n matrix (i.e., a m a t r i x that h a s multiple eigenvalues). It is a direct extension o f t h e case for distinct eigenvalues [28, 35"1. L e m n m A.I. I f a sphericalized Hermitian matrix • i s perturbed as • = • + 6@ = @ + e,F, where E is a constant Hermitian matrix and [ ~l << I, then the eigenspaces and eigenlevels o f • are given by Oi = Ui + 6U~ + O{r, 2) and Di = Di + ODi + 0(~2), where aD,

=

H ~ U ~ 6q~U~Hi

(Al)

and ~" U H H]~U~6@UiHI 6U,--/., s s ff._--~. j~ti ~i - - t~j with Hi chosen to diagonalize ,SDi.

(A2)

119

References

[I] G. Adams, M.F. Griffinand G.W. Stewart,"Direction of arrival estimation using the rank rcvcalin$ URV decomposition",proc. Im~'~mt. Conf.. Acoust. Speech Signal Process.-9/. Toronto, Canada, 1991, pp. 1385-1388. [2] C. Bischofand G. Shroff, "On updating sigmtl subspac~", IEEE Trans. Signal Process.,Vol. SP-40, January 1992,pp. 96-105. [3] R.W. Brockctt, "Dynamical systems that sort lists, diagonalizc matrices and solve linear prolp'ammins problems". 27th Conf..on Decision and Contro/, 1988. [4] J.R. Bunch and C.P. Nielsen,"Updating the singular value decomposition", Numce. Mm&, VUI.31,1978, pp. I 11-129. [5] J.R. Bunch, C.P. Nilscn and D.C. Sorcuscn, -Rank-one modification of the s]mlmetri¢ eigenp¢oblem', Numer. Math., Vol. 31, 1978, pp. 31-48. [6] B. Champagne, "Adaptive eigendecompmition of data covariance matrices based on first-ordar perturbations', IEEE Trans. Signal. Process., Vol. SP-42, October 1994, pp. 2758--2770. [7] P. Comon and G.H. Golub. "Tracking a few extreme singular values and vectors in signal processing", Proc IEEE, Vol. 78, August 1990, pp. 1327--1343. [8] R.D. DeGroat, ~Authofs r©rAy",IEE£ Trans. Sional Procezs.. Vol. SP-39, August 1991, pp. 1913--1914. [9] R.D. ~ r o a t , "Non-iterative subslmce tracking", IEEE Trans. Sional Process.,Vol. 40, March 1992,pp. 571--577. [101 R.D. DeGroat and E.M. Dowling, "Non-iterative subspace updatiag," in: SPIE Adv. Sio. Proc. AIo&,Arehs., and Impis. II. San Diego, California, July 1991, pp. 376-387. Ill1 R.D. DeGroat and E.M. Dowling. ~Spherical sabslmce tracking: Analysis, convergence and detection schemes', 26th Asilomar Conf, October 1992, pp. 561-565. [12] R.D. DeGroat, E. Dowling and D. Linebarser, PtoofofSS convergence, UTD Tedmicul Relmrt, 1995. [13] R.D. IX'Groat and R.A. Roberts, "A family of rank~me eigenstructure updating methods", Proc. 21st Astlorr,~ Conf.., 1987, PiP. 564-568. [14] R.D. [kf;roat and B.A. Robert& "A family of rank-ott¢ subspace updating methods", in: E.F. Deptettere, ed., SVD and Signal Processi~y:.Agl~, Apld~ and Arcks, N o r t h - H o l . land, Amsterdam, 1988, pp. 277-300. [151 R.D. DcGroat and B.A. Roberts, " E ~ . ~ ; , nuraeri~lly stabilized rank-one eigmlstru~un~updating", IEEE TFans. Acoust. Seeech $/qna/Process~ VUI. ASSP-38, February 1990, pp. 301-316. [16] R.D. DeGroat, H. Ye and E.M. Dowfin& "An asympto~ analysis of spherical subspace tracking', 27th As//omar Conf., November 1993. [17]EM. Dowling, LP. Ammann and ILD. 13eGtrmL"A TQRiteration based SVD for real time angle and fgqum~ Iradiins', IEEE Trims. Sigma/Process.,April 1994,pp. 914-92..5. [18] E.M. Dowlins aral R.D. DeGroat, "Recursiv¢ total least squares adaptive filtering', in: S. Haykin, ed.. SPIE Adaptive Signal Processir,q, San Diego, California, July 1991, pp. 35-46.

120

R.D. DeGroat et al. / Signal Processing 50 (1996) I01-121

[19] E.M. Dowlin8 and R.D. DeGroat, "A spherical subspace based adaptive filler", Prec. imernat. Conf. Ac'ou~c, Speech Signal Process., April 1993, pp. II1: 504-507. r20] E.M Bowling and R.D. DeGroat, "Adaptation dynamics ~f the spherical subspace tracker", IEEE Trans. Signal Proce:s., October 1992, pp. 2599-2602. [21] E.M. Dowiing, R.D. DeGroat and D. Linebarger, "Efficient, high performance subspace based tracking problems", in: F.T. Luk, ed., Adv. Sig. Prec. Alos., Archs. and Appls. VL Vol. 2563, SPIE, 1995, pp. 253-264. [22] E.M. Dowling, R.D. DeGrnat, D.A. Linebarger and Z. Fu, "Real time architectures for spbericalized SVD updating", in: SVD and Signal Processing: Algorithms, Applications and Architectures III, North-Holland, Amsterdam, 1995. [23] EM. Dowiing, R.D. DeGroat, D.A. Linebarger and H. Ye, "Spbericalized SVD updating for subspace tracking". in: SVD and Signal Processing: Algorithms, Applications and Architectures IlL North-Holland, Amsterdam, 1995. [24] A. Ecklman and G. Stewart. "Scaling for orthogonality ", IEEE Trans. Signal Process., Vol ASSP-41, April 1993, pp. 1676-1677. [25] M.P. Fargnes, "Tracking moving sources using the rank revealing QR factorization', 25th Asilomar Conf on Sigs., Sys., and Comps., November 1991, pp. 292-295. [26] W. Ferzali and J.G. Proakis, "Adaptive SVD algorithm and applications", in: S VD and Signal Processing !1, Elsevier, Amsterdam, 1992. [27] D.R. Fuhrmann,'~'An algorithm for subspacecomputation, with applications in signal processing", SlAM J. Matrix Anal Appl., Vol. 9, April 1988, pp. 213-220. [28] K. Fukunagn, Introduction to Statistical Pattern Recognitim, 2nd Edition, Academic Press, New York, 1990. [29] Z. Fu and E.M. Dowlin8, "Conjugate gradient projection subspace tracking', 281h Ann. Asilomar Conf. on Signals, Systems, ant Computers, 1994, pp. 612-616. [30] Z. Fu and E.M. Dowling, "Conjugate gradient projection subspace tracking and systolic implementation", IEEE Trans. Signal Process., Submitted. [31] Z. Fu and E.M. Dowling, "Conjugate gradient eigenstructurn tracking adaptive spectral estimation", IEEE Trans. Signal Process., Vol. 43. May 1995,1151-1160. [32] Z. Fu and E.M. Dowling, ~Conjugate gradient eigenstructare tracking for adaptive spectral estimation", IEEE Trans. Signal Processing, Vol. 43, May 1995, pp. 1151-1160. [33] G.H. Golub, "Some modified matrix eigenvalue problems', SIAM Roy., Vol. 15, 1973, pp. 318-334. [34] G.H. Golub and C.F. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, Baltimore, Maryland. 1989. [35] W.W. Hager, Applied Numerical Linear Algebra. Prentice Hail, Engk-wood Cliffs, N J, 1988. [36] U. Helmke, "tsospectral flows on symmetric matrices and the Riceati equation", Systems Control Lett.. Vol. 16,1991, pp. 159-165.

[37] I. Karasalo, "Estimating the covariance matrix by signal • ub~pa¢i aiei'aliiif, lEE[ Trans. Acouat, Speech Bignal Process., Vol. ASSP-34, February 1986, pp. 8-12. ['38] R. Kumaresan and D.W. Tufts, "Estimating the angles of arrival of multiple plane waves", lEE[ Trans. Aerospace Electronic Syster~, V01. AES-19, January 1982, pp, 134-139, [39] L. Ljung, ~Analysis of recursive stochastic algorithms", IEEE Trans. Automat. Control, Vol. AC-22, 1977. [40] L. Ljung and T. Soderstrnm, Theory and Practice of Recursire Identification, MIT Press, New York, 1983. [41] R. Mathias, "Analysis of algorithms for orthogonalizing products of unitary matrices", J. Numer. Linear Algebra Appls., To appear. [42] M. Moonen, P. VanDooren and J. Vanderwalle, "Updating singular value decompositions. A parallel implementation", in: SPIE Advanced Algorithms and Architectures for Signal Processing, San Diego, California, 1989, pp. 80-91. [43] M. Moonen, P. VanDooren and J. Vanderwalle, "A note on efficient stabilized rank-one eigenstructure updating", lEE[ Trans. Signal Process., V01. SP-39, August 1991, pp. 1911-1913. [44] N.L. Owsley, "Adaptive data orthogonalization", Prec. lnternat. Conf. Acoust. Speech Signal Process., 1978. [45] V.F. Pisarenko, "The retrieval of harmonics from a covariance function", Geophys. J. Roy. Astron. Sac., Vol. 33, 1973, pp, 347-366. [463 D.J. Rabideau and A.O. Steinhardt, "Fast subspace tracking", Prec. 7th IEEE Spectrum Estimation Workshop, June 1994, pp. 353-356. [47] M.A. Rahman and K. Yu, "Total least squares approach for fsequen~ estimation ~:sing linear prediction", IEEE Acoust. Speech Signal Process., Vol. ASSP-34, October 1987, pp. 1440-1454. [48] R. Roy, A. Paulraj and T. Kailath, "ESPRIT - A subspace rotation approach to estimation of parameters of cisoids in noise", IEEE Trans. Acoust. Speech Signal Process., October 1986, pp. 1340-1342. [49] R. Schmidt, "Multiple emitter location and signal parameter estimation", lEE[ Trans. Antennas Propag.. Vol. AP-34, March 1986. [50] R. Schreiber, "Implementation of adaptive array algorithms", lEE[ Trans. Acoust. Speech Signal Process., Vol. ASSP-34, October 1986, pp. 1038-1045. [5 I'l G.W. Stewart, "An updating algorithm for subspace tracking", IEEE Trans. Signal Process., Vol. SP-40, June 1992, pp. 1535--1541. [52] p. Stolen and A. Nehorai, "MUSIC, maximum likelihood, and Cramer-Rao bound", IEEE Trans. Acoust. Speech Signal Process., Vol. ASSP-34, June 1986, pp. 585-602. [53] D.W. Tufts and R. Kumaresan, "Singular value decomposition and improved frequency estimation using linear prediction", IEEE Trans. Acous~. Speech Signal Process., Vol. 30, August 1982, pp. 671-675. [54] S. VanHuffel and H. Zha. "The restricted total least squares problem: formulation, algorithm, and properties", SIAM J. MAA., Vol. 12, April 1991, pp. 292-309.

R,D DeGroat el al. / Signal Processing 50 fl996) 101-121 [55] M. Wax and T. Kailath, "Detection of signals by informa:Ioil theoretic criteria", I EEE Trans. Acoust. Speech Sianal Process., ¥ol. ASSP-33, April 1985, pp. 387-392. [56.] Y.G. Xu and T. Kailath, "Application of fast subspace decomposition to signal processing and communication problems", IEEE Trans. Signal Process., Vol. 42. June i994, pp. 1453-146i. [57] G. Xu and T. Kailath, "Fast subspac,¢ dccomposir.on", IEEE Trans. Signal Process., Vol. 42, March 1994, pp. 539-551. [58] G. Xu, H. Zha, G. Golub and T. Kailath, "Fast and robust algorithms for updating signal subspaccs', IEEE Trans Circuits and Systems, Vol. 41, June 1994,

121

[59] B. Yang. "Gradient based subspace tracking algorithms and systolic implementation", Int~nat. J. High Speed Electron. System.s, Vol. 4, r~o. 2, 1993, pp. 203-218. [60] B. Yang, "Projection approximation subspece tracking'. IEEE Trans. Signal Process., VoL $P-43. January 1995, pp. 95 107 [_61] J.F. Yang and M. Kay©h, "Adal~iVc ¢ilb'nsubsp~ccs algorithms for direction or frequency estimation and tracking", IEEE Trans. Acoust. Speech Signal Process., Vo|. ASSP36, February 1988, pp. 241-251. [62] Y. Yin and P. Krishnaiah, "On some nonparamettic methodsfor detectionof the nurr,l~r of signals",IEEE Trans, ,Signal Process., Vol. 35, November 1987, pp. 1533-1538.