Some Experiments with Eigenanalysis of Data Matrices

June 28, 2007

This post gives some examples of the use of singular value decomposition (SVD) in some data separation problems. Could be useful in radar, sonar, computer tomorgaphy, etc…

Let’s start with having a look where to find some info on this… The problem at hand is also called the principal components analysis (PCA). So, the book by J V Stone is an excellent start. Lately, I have found the following reference. A lot of ideas (and inspiring Matlab code) came from the HOSA (higher-order spectral analysis) Matlab package and its accompanying manual. There will be probably much more (possibly better) references to study. If you know any good ones, please let me know.

The examples are accompanied with a Scilab code. Scilab is a free equivalent of Matlab. I would encourage to try the code, simply by pasting it into a workspace window. You could either use Matlab (then you need to make some changes) or you can try Scilab (it doesn’t take too much effort and disk space to install it).

Example 1

Consider the following example: A two-dimensional harmonic function

 

f(x,y)=A\exp{\left[-\jmath (x k_x + y k_y)\right]}

can be also expressed as

 

f(x,y)=A\exp{(-\jmath x k_x)} \exp{(-\jmath y k_y)},

 

where x,y are variables and k_x,k_y are some constants. A sampled version would be:

 

f(x_i,y_j)=A\exp{(-\jmath x_i k_x)} \exp{(-\jmath y_j k_y)}

 

One can then write that

 

f(x_i,y_j)=\{f\}_{i,j}=A\exp{(-\jmath \mathbf{x} k_x)} \exp{(-\jmath \mathbf{y}' k_y)},

 

where \mathbf{x},\mathbf{y} are column vectors.

 

A Scilab code that would generate such a function is shown in this table:

 


clear
cj=sqrt(-1);m=16;n=16;
x=linspace(-1,1,m)'*%pi;
y=linspace(-1,1,n)'*%pi;

s=zeros(m,n);
kx=3;ky=1;
g=exp(-cj*x*kx);
h=exp(-cj*y'*ky);
s=g*h;

 

 

The vectors \mathbf{x},\mathbf{y} form an orthogonal basis for the function f. Singular value decomposition factors a matrix into a product of three matrices: \mathbf{U},\mathbf{V} and \mathbf{D}. Matrices \mathbf{U},\mathbf{V} are orthogonal or unitary square matrices, \mathbf{D} is a diagonal matrix. Hence, it should not come as a surprise that the SVD performed on matrix \{f\}_{i,j} will give the following result:

 

\{f\}_{i,j}=\mathbf{U\Sigma V}',

 

where

\mathbf{U} = \left[\,\exp{(-\jmath \mathbf{x} k_x)} \quad \mathbf{0} \quad \dots\right]

\mathbf{\Sigma} = \left[\,\mathbf{a} \quad \mathbf{0} \quad \dots\right]

 

\mathbf{V} = \left[\,\exp{(-\jmath \mathbf{y} k_y)} \quad\mathbf{0}\quad \dots\right]

\mathbf{a} is a coulmn vector containing a single eigenvalue a in the first entry, the rest are zeros. \mathbf{o} is a column vector with zero entries.

 

 

A code that could demonstrate this is in the table below. We also plot the Fourier transform of the signal s. This will allow us to observe separation properties of SVD, later on.

 


[U,D,V]=svd(s);
pad=96;
S=fftshift(fft2(...
fftshift([s,...
zeros(m,pad-n);...
zeros(pad-m,pad)])));

fig1=scf(1);
map=jetcolormap(256);
xset("colormap",map)
Matplot(log10(...
abs(S)+1e-12)*50)
fig2=scf(2);
plot(diag(D),'.')


fig3=scf(3);
plot(x,U(:,1),...
x,g/sqrt(D(1,1)))


fig4=scf(4);
plot(y,V(:,1),...
y,h/sqrt(D(1,1)))

 

Example 2

In the second example, we add one more signal of much smaller amplitude than signal s. For instance, like so:


s1=zeros(m,n);kx=-2;ky=5;
w=exp(-cj*x*kx);
z=exp(-cj*y'*ky);
s1=.1*w*z;
s=s+s1

After running the code from example 1 again, one can observe that now there will be two singular values corresponding to two singular vectors in matrices \mathbf{U} and \mathbf{V}. It can be further observed that again, the first two columns of \mathbf{U} will be very close to functions g and w, respectively. Just as well as columns of matrix \mathbf{V} will be close to functions h and z. Clearly, now the signals s and s can be separated simply by nulling their corresponding singular values. The singular values are in descending order. Thus, the signals will be ordered depending on their amplitude.

A different situation occurs, if signals s and s1 have comparable amplitudes. Matrix \mathbf{D} will still contain two singular values, but the singular vectors will no longer correspond to basis functions of the original signals. In such cases, it is said that the original signals are mixed. Canceling one of the singular values will not produce one of the two original signals, as can be easily verified. The way out of this leads onto so-called independent component analysis, or ICA for short. Fortunately, we will see later that there are some cases of interest where PCA (that is, the method described above) will be sufficient.

 

Example 3

In the last example, let us consider the following special case: Suppose that both signals in our simulation have the same constant k_y, say equal to 1. What will the SVD look like? Even though two signals are present, only one singular value is calculated. The reason is quite plain:

\sum_i A_i\exp{\left[-\jmath (x k_{x,i} + y k_y)\right]} = g(y)\cdot f(x)

 

where

 

f(x) = \sum_i A_i\exp{(-\jmath x k_{x,i})}

g(y) = \exp{(-\jmath y k_y)}

 

and

 

k_{x,i} \in [-K_x,K_x], \qquad K_x=\mbox{const.}

 

So, now we have a signal that can be decomposed into two singular vectors f(\mathbf{x}) and g(\mathbf{y}).

 

Our simulation is performed for i=2, but the general case can be easily verified:

 


s=zeros(m,n);
f=0;
Kx=7;ky=1;
kx=linspace(-Kx,Kx,10);
h=exp(-cj*y'*ky);
for i=1:10,
g=exp(-cj*x*kx(i));
f=f+rand(1,1)*g;
end
s=f*h;

It becomes obvious that adding one more signal such as s1 will produce the same behavior as above.

Conclusions

One may ask why all this exposition when all the operations described can be performed in the Fourier domain as the figures clearly show. Well, the reason is that the most interesting cases are when the number of samples is so small that the Fourier transform using the FFT algorithm fails to resolve the signals in question. It can be verified by setting smaller numbers m and n, that while capability of FFT analysis to resolve the signals decreases, SVD analysis remains unaffected by the number of samples.

 

Posted by doctuh4iah
Filed in Scientific Stuff
Leave a Comment »

Leave a Reply