Introduction to Subspace Tracking IEEE AES Dallas Chapter Meeting Richardson, TX Feb 28, 2006

1. Introduction to Subspace TrackingIEEE AES Dallas Chapter MeetingRichardson, TXFeb 28, 2006 Scott Baker L-3 Communications Greenville, TX scott.baker@L-3com.com

2. Outline • Introduction • Time Series Data Model • Subspaces & the Subspace Tracking Idea • Subspace Tracking Categories • A Secular Equation Subspace Tracker • A Plane Rotation Subspace Tracker • Simulations

3. Data Model • The time-series model for a sum of complex sinusoids in noise may be written • This time series may be placed in the vector form Type? Singular? Type? • And in a matrix form Size? Looks like?

4. Data Model • Two Typical Goals in Time Series Estimation • 1)Parameter estimation • Number of signals Q • The frequencies themselves • The signal amplitudes • 2)Noise reduction • Find a cleaner version of x(n) • Typical assumptions • Signal model is stationary over a block • Noise is pre-whitened as needed • N >= M (more snapshots than sensors)

5. Data Model • The data matrix is often place in form of a correlation matrix • Data matrix forms lead to an SVD • Correlation matrix forms lead to an EVD where What is U? What is V?

6. Data Model • The singular value decomposition can be seen in more detail as N x N N x M M x M ui are unitary left singular vectors vi are unitary right singular vectors si are singular values

7. Data Model • The eigenvalue decomposition can be seen in more detail as M x M M x M vi are unitary right eigenvectors li are eigenvalues M x M How is the SVD and EVD related? How are the singularvalues and eigenvalues related? How are the singular vectors and eigenvectors related? How can we interpret an eigenvalue and an eigenvector?

8. Data Model Box 1 M = 40; x = sin(2*pi*1 * [1:100]/22); % 1 Hz sinewave, figure; plot(x) X = datamtrx(x, M); R = X'*X; [V, L] = svd(R); figure; plot(V(:,1));hold on; plot(V(:,2)) Period = 22 samples N=61 Box 2 >> disp(L(1:3,1:3)) 693.4222 0 0 0 531.3492 0 0 0 0.0000 Why not one? Figure 1 - x Figure 2 – V (1st two) Figure 3 – V (last 38)

9. Subspaces & Tracking • From this example, consider a partition of the EVD Signal Eigenvectors Signal Eigenvalues “large” Noise Eigenvalues “small” Noise Eigenvectors Span? VS and A ?

10. Subspaces & Tracking • Extend our previous example, by adding noise Box 1 M = 40; x = sin(2*pi*1 * [1:100]/22); % 1 Hz sinewave, x = x + randn(1,100); figure; plot(x) X = datamtrx(x, M); R = X'*X; [V, L] = svd(R); figure; plot(V(:,1));hold on; plot(V(:,2)) Only one change Effect? Meaning? Figure 1 - x Figure 2 – V (1st two) Figure 3 – Eigenvalues

11. Subspaces & Tracking • Knowledge of the signal eigenvectors VS gives us knowledge of the frequency matrix A • Knowing the span of VS is equivalent to knowing the span of VN • One can show that VSVSH + VNVNH = I • MUSIC • The MUSIC algorithm uses VN to estimate the frequencies • It is possible to show VNHA = 0 • Which the same as VNH [a1,a2,…,aQ] = 0

12. Subspaces & Tracking • Recall that the frequency matrix had the form • So, we know the structure of a column and if we had the qth frequency, we could construct the qth column • We do know that VNHaq = 0 and we know VN from the decomposition of the data • To find wq, what can we do?

13. Subspaces & Tracking • Yup, good guess! • We only need to build a sweep vector a(w) for w from DC to ½ the sample frequency • We can plot the inverse instead so that peaks are at the Q frequencies • S(w) = 1 / || VNHa(w)|| • S(w) is then the frequency spectrum with peaks at the Q frequencies • Way cool!

14. Subspaces & Tracking • So how do we do this again? • Step 1: Build a data matrix • Step 2: Build a correlation matrix and find its EVD • Step 3: Sort the eigenvalues largest to smallest and partition the EVD into signal portion & noise portion • Step 4: Compute the MUSIC spectrum & locate Q peaks S(w) = 1 / || VNHa(w)||

15. Subspaces & Tracking randn('state',1); n=[0:99]; x=exp(i*pi/2*n)+2*exp(i*pi/4*n)+exp(i*pi/3*n)+randn(1,100); R=corrmtx(x,12,'mod'); pmusic(R,3,'whole') (type ‘help pmusic’) MUSIC Spectrum S(w) w

16. Subspaces & Tracking • So are we done yet? • In stationary cases, yes • In non-stationary cases, no • The Q frequencies are often time-varying • Hopping sigals • Chirping signals • Any kind of frequency modulated signal • It is too expensive to compute an EVD at every instant of time • Instead, use a subspace tracker to track the changes in VS (or equivalently VN) as new data comes in

17. Subspaces & Tracking • Let’s slightly modify our data model so that it is recursive What type of update is this called? • The factor a is called a fading factor and is slightly less than unity • The effective window size is B = 1 / (1-a) • The update assumes that we know the old EVD • And we wish to know the new EVD Do we need the whole EVD?

18. Subspace Tracking Categories • There are about 250 published papers on subspace tracking • It sky-rocketed after the invention of the MUSIC algorithm • There are several broad categories • Secular equation (e.g., Bunch) • Plane rotation (e.g., Businger) • Optimization • Lanczos • Perturbation • Power Method & Subspace Iteration • Eigenvector computation with a priori eigenvalues • Analog or Continuous Time • Hybrid Techniques

19. The Secular Equation Subspace Tracker • Golub discovered the secular equation in 1971 • Can be considered the first eigen- based subspace tracker (Businger wrote the first overall ST in 1970) • Reveals a nice relation between original and updated eigenvalues • Seven years later, Bunch, Nielsen and Sorensen introduced the classic secular equation subspace tracker • A keen observation is that the eigenvalues of the new correlation are larger assuming the data has non-zero energy Gene Golub (Stanford) original eigenvalue, largest updated eigenvalue, largest

20. The Secular Equation Subspace Tracker • This observation is important because it was noted that the rank-1 modification of the correlation can be written as the rank-1 modification of a diagonal matrix Stick in the EVD Because V is? And z is? What is left in the brackets then? Recall the monotonicity relation

21. The Secular Equation Subspace Tracker Note we have the form where Is D a diagonal matrix? Note that we can form the EVD And so the resultant EVD of the updated correlation is Remember the EVD…

22. The Secular Equation Subspace Tracker Answer: They are equal Hence, to find the updated eigenvalues, we can find the eigenvalue of D. Recall from basic linear algebra, we can write These two expressions can be used to find a formula for the updated eigenvalues Combining these three equations leads to the secular equation

23. The Secular Equation Subspace Tracker The secular equation has zeros & poles at the eigenvalue locations Eigenvector Computation (i.i. not needed) Original eigenvalue (p) Efficient to curve fit here between poles Updated eigenvalue (z)

24. A Plane Rot. Subspace Tracker • The first subspace tracker developed in 1970 • P. Businger of Bell Labs • Used plane rotations • Impetus was updating a least squares problem • SVD based solution • Businger acknowledged the problem was posed by G. Dahlquist and related by G. Golub

25. How might we proceed to update an SVD? A Plane Rot. Subspace Tracker Adding a row What structure does this factor have? Similar trick to what? What is Q? Z? Is this in SVD form?

26. A Plane Rot. Subspace Tracker • The Q and Z matrices are called plane rotation matrices • The specific type of plane rotation used here is called a Given’s rotation Real Given’s Rotation

27. A Plane Rot. Subspace Tracker • So how do we use Given’s rotations to diagonalize the matrix? • We have to be careful • Eliminating one element in a column / row may perturb another column / row • Need to apply Given’s rotations to both rows and columns a * x 0 r x Q = * * x b 0 x Given’s Matrix Original Matrix Transformed Matrix

28. A Plane Rot. Subspace Tracker * * x * * * = * Q1 * * Left (Row) Given’s Rotations First * * * * * * * 0 What is the resultant form? * * x * * x * * x = * * Q3 * * 0 0 0 0 0 0

29. A Plane Rot. Subspace Tracker • The bidiagonal form is reduced to a diagonal form • The process is iterative • Start with column rotations (right) • Then a set of row rotations (left) • Repeat until the energy is minimal off the main diagonal

30. A Plane Rot. Subspace Tracker • So, going back to our SVD update • Q and Z are series of Given’s rotation matrices • Q1Q2...QL • Z1Z2….ZP • And so our final result then is in the proper SVD form

31. Example Simulation • Single signal • SNR is 15 dB • Sample Rate = 640MHz • Linear FM signal • 140 to 160 MHz • 20MHz / 1000 samples • 12.8 GHz / second • a = 0.975 (80 samples) • TQR + root-MUSIC • Spline Interpolation • Polynomial = feature This is a spline fit to the first ramp.

32. Other Examples