Transcript Introduction to Subspace Tracking IEEE AES Dallas Chapter Meeting Richardson, TX Feb 28, 2006 Scott Baker L-3 Communications Greenville, TX [email protected].

Introduction to Subspace Tracking
IEEE AES Dallas Chapter Meeting
Richardson, TX
Feb 28, 2006
Scott Baker
L-3 Communications
Greenville, TX
[email protected]
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
Data Model
 The time-series model for a sum of complex
sinusoids in noise may be written
Q
j q n
x ( n )   sq e
 e(n)
q 1
 This time series may be placed in the
vector form
Type?
Singular?
Type?
 x (1) 


x(2)
  A (ω )s  e ( n )
x(n)  
  


 x ( n ) 
 j  (1)
1
e
A (ω )  

 e j1 ( n )

 And in a matrix form
X  x  N  M  1, x  N  M ,  , x  N 
e
e
j  2 (1)
j 2 ( n )
Size?
Looks like?
j  Q (1) 



j Q ( n )

e

e
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)
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
X  U ΣV
H
R  X
where
UU
VV
H
X  V ΛV
H
What is U?
H
H
 UU
H
 VV
H
 IN
 IM
What is V?
Data Model
 The singular value decomposition can be
seen in more detail as
X U ΣV
H
 u1
u2

 1


uN 



2

0





M 

 v H
 1
 v2H


v H
 M







NxN
NxM
ui are unitary left singular vectors
vi are unitary right singular vectors
i are singular values
MxM
Data Model
 The eigenvalue decomposition can be seen
in more detail as
R  V ΛV
H
 v 1
v2

 l1

v M 



l2






l M 
 v H
 1
 v2H


H
v
 M







MxM
vi are unitary right eigenvectors
li are eigenvalues
MxM
MxM
How is the SVD and EVD related?
How are the singular values and eigenvalues related?
How are the singular vectors and eigenvectors related?
How can we interpret an eigenvalue and an eigenvector?
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)
1
Figure 3 – V (last 38)
0.25
0.8
0.2
0.6
0.15
0.4
0.1
0.2
0.05
0.5
0
0
0
-0.2
-0.05
-0.4
-0.1
-0.6
-0.15
-0.8
-0.2
-1
-0.25
0
10
20
30
40
50
60
70
80
90
100
-0.5
0
5
10
15
20
25
30
35
40
-1
0
5
10
15
20
25
30
35
40
Subspaces & Tracking
 From this example, consider a partition of the EVD
R  V ΛV
H
 V S
Signal
Eigenvectors
Span?
VS and A ?
Λ S
VN 



ΛN 
V H
 S
 V N H
Signal
Eigenvalues
Noise
Eigenvectors “large”



Noise
Eigenvalues
“small”
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 2 – V (1st two)
Figure 1 - x
3
2
Figure 3 – Eigenvalues
0.3
700
0.2
600
0.1
500
0
400
-0.1
300
-0.2
200
-0.3
100
1
0
-1
-2
-3
0
10
20
30
40
50
60
70
80
90
100
-0.4
0
5
10
15
20
25
30
35
40
0
0
5
10
15
20
25
30
35
40
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 VS VSH + VNVNH = I
 MUSIC



The MUSIC algorithm uses VN to estimate the
frequencies
It is possible to show VNH A = 0
Which the same as VNH [a1,a2,…,aQ] = 0
Subspaces & Tracking
 Recall that the frequency matrix had the form
 j  (1)
1
e
A (ω )  

 e j1 ( n )

e
e
j  2 (1)
j 2 ( n )
j  Q (1) 



j Q ( n )

e

e
 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 VNH aq = 0 and we know VN from the
decomposition of the data
 To find q, what can we do?
Subspaces & Tracking
 Yup, good guess!
 We only need to build a sweep vector a()
for  from DC to ½ the sample frequency
 We can plot the inverse instead so that
peaks are at the Q frequencies

S() = 1 / || VNH a()||
 S() is then the frequency spectrum with
peaks at the Q frequencies
 Way cool!
Subspaces & Tracking
 So how do we do this again?
 Step 1: Build a data matrix
X  x  N  M  1, x  N  M ,  , x  N 
 Step 2: Build a correlation matrix and find its EVD
H
H
R  X X  V ΛV
 Step 3: Sort the eigenvalues largest to smallest and
partition the EVD into signal portion & noise portion
V  V S
VN

 Step 4: Compute the MUSIC spectrum & locate Q
peaks
S() = 1 / || VNH a()||
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');
(type ‘help pmusic’)
pmusic(R,3,'whole')
MUSIC Spectrum
Ps eudos pec trum Es timate via MUSIC
60
50
40
Power (dB)
30
S()
20
10
0
-10
-20
-30
0
0.2
0.4
0.6
0.8
1
Normaliz ed Frequenc y

1.2
1.4
( rad/s ample)
1.6
1.8
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
Subspaces & Tracking
 Let’s slightly modify our data model so that it is
recursive
What type of update
is this called?
~
H
R   R  1   xx
 The factor  is called a fading factor and is slightly
less than unity
 The effective window size is B = 1 / (1-)
 The update assumes that we know the old EVD
R  VS
Λ S
VN  



ΛN 
V H
 S
 V N H



 And we wish to know the new EVD

~
~
R  VS
~
Λ
~
VN  S


~
V H
S
 ~
~
Λ N   V N H




Do we need
the whole
EVD?
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
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
~
~
~
l1  l1  l 2  l 2    l M  l M
original eigenvalue, largest
updated eigenvalue, largest
Gene Golub
(Stanford)
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
~
H


R   R  1   xx
Stick in the EVD
~
H
H
R  V Λ V
 1   xx
~
H
H
H
R  V [ Λ  1   V xx V ] V
~
H
H
R  V [ Λ  1   zz ] V
Because V is?
And z is?
What is left in the brackets then?
Recall the monotonicity relation
The Secular Equation Subspace Tracker
Note we have the form
D   Λ  1   zz
~
H where
R  VDV
Note that we can form the EVD
D QΒQ
H
Is D a diagonal matrix?
H
And so the resultant EVD of the updated correlation is
~
R  VQ
Β Q
~ ~
V Λ
H
V
H
Remember the EVD…
~
~ ~ ~H
R V ΛV
~
So how are the eigenvalue s of D and R related?
The Secular Equation Subspace Tracker
Answer: They are equal
Hence, to find the updated eigenvalues, we can find the eigenvalue of D.
D   Λ  1   zz
H
Recall from basic linear algebra, we can write
Db   b
These two expressions can be used to find a formula for the
updated eigenvalues
det( Db   b )  0
Combining these three equations leads to the secular equation
The Secular Equation Subspace Tracker
The secular equation has zeros & poles at the eigenvalue locations
1
VD i z
~
vi 
1
Di z
2
~
D i  D  li I
2
M
zi
 (  )  1  (1   ) 
i  1 (l i   )
Eigenvector
Computation
(i.i. not needed)
Original eigenvalue (p)
Efficient to curve fit
here between poles
Updated eigenvalue (z)
5
4
3
2
 ( )
1
0
-1
-2
-3
-4
-5
1
1.5
2
2.5
3
3.5

4
4.5
5
5.5
6
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
A Plane Rot. Subspace Tracker
 How might we proceed to update an SVD?
X U ΣV
H
Adding a row
X

~ 
X  
H 
 (1   ) x 
 U Σ V H
 
H
 (1   ) x



Similar trick
to what?
What structure does this
factor have?
U
 
0
0 
Σ
 H

H V
I   (1   ) x V 
U
 
0
0
Σ

H 
H H
QQ


H  ZZ V
I
 (1   ) x V 
U
 
0
0 ~ H H
Q Z V
I
What is Q? Z?
Is this in SVD form?
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
 c

 s
s *  a   r 
    
c  b  0 
| r | | a |  | b |
2
2
Real Given’s Rotation
A p p ly
2
G iv en's
b
 sa  cb  0
a
c 
a
2
c  s
2
2
1
2
 b
2
s 
b
a
a
c
B efor e
R ota tio n
r
A fter
R ota tio n
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
Q
Given’s
Matrix
a
b
0
*
* *
Original
Matrix
=
r
0
x
x
x
x

Transformed
Matrix
A Plane Rot. Subspace Tracker
*
Q1
Left
(Row)
Given’s
Rotations
First
*
*
*
*
*** *
*
*
Q3
0 0
=
0
*
*
*
**
=
x
*
*
*
** *
What is the
resultant
form?
x
*
x
*
x
*
0 0 0 0
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
A Plane Rot. Subspace Tracker
 So, going back to our SVD update
U
~
X  
0
0 ~ H H
Q Z V
I
 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
~
~ ~ ~H
X  U ΣV
Example Simulation
 Single signal
8
1.62
 SNR is 15 dB
1.6
1.58
 Sample Rate = 640MHz
1.56
 Linear FM signal
1.54
frequency
140 to 160 MHz
 20MHz / 1000 samples
 12.8 GHz / second
  = 0.975 (80 samples)
 TQR + root-MUSIC
 Spline Interpolation
 Polynomial = feature

Freq Track w/ MTH=5, SNR=15, EBS = 40
x 10
1.52
1.5
1.48
1.46
1.44
1.42
500
1000 1500 2000 2500 3000 3500 4000 4500 5000
time sample
8
x 10
1.6
1.58
1.56
This is a spline
fit to the first
ramp.
1.54
1.52
1.5
1.48
1.46
1.44
1.42
100
200
300
400
500
600
700
800
900
Other Examples
x
8
Freq Trac k
10
w/
MTH=5,
SNR=15,
EBS
=
40
1.6
frequency
1.55
1.5
1.45
1.4
0
x
100
8
200
300
Freq Trac k
10
400
500
600
time s ample
w/
MTH=5,
SNR=15,
700
EBS
800
=
900
1000
40
1.6
1.5
frequency
1.4
1.3
1.2
1.1
1
0
200
400
600
800
1000
1200
time s ample
1400
1600
1800
2000