[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: how can I do Principal Components Analysis with octave?

From: Vic Norton
Subject: Re: how can I do Principal Components Analysis with octave?
Date: Wed, 12 May 2004 09:31:01 -0400

Hi Mario,

I'll give you my slant on Principal Components Analysis. I prefer
singular value decomposition to the covariance-eigen-value route. I
understand the former better, and it is more efficient.

Start with an m x n data matrix X and an m x 1 vector of weights w. We
require w(i) > 0 for i = 1,..,m and sum(w) = 1. In the normal
application w(i) = 1/m for i = 1,..,m.

   E = sum(diag(w) * X)
is the 1 x n matrix of expected values of the columns of X,
   V = (X - ones(m, 1) * E)' * diag(w) * (X - ones(m, 1) * E)
is the n x n covariance matrix corresponding of X
(V = ((m - 1)/m) * cov(X) in the normal case),
   S = sqrt(diag(V))'
is the 1 x n matrix of standard deviations of the columns of X, and
   C = diag(1 ./ S) * V * diag(1 ./ S)
is the correlation matrix. (C = corrcoef(X) in the normal case.)

Now let's geometrize these computations. First reassign
   w = sqrt(w)
and put
   X = diag(w) * X.
Then the expected values of the columns of the original X are
   E = w' * X.
One can think of E as the perpendicular projection of the new X
onto w-axis (more or less) with
   F = X - w * E
the complementary projection. Now
   S = sqrt(sum(F .^ 2))
is the 1 x n matrix of standard deviations of the columns of the
original X, and
   V = F' * F
is the covariance matrix above. To normalize the situation set
   F1 = F * diag(1 ./ S).
   C = F1' * F1
is the correlation matrix above.

OK, let's get to principal component analysis. We work with F or
F1. It is not necessary to compute the covariance or correlation

If the units of measurement are important, do
   [U, Sig, W] = svd(F, 1)
otherwise do
   [U1, Sig1, W1] = svd(F1, 1).
The singular values, the coefficients of the diagonal matrices Sig and
Sig1, measure the significance of the various factors:
   Sig(1) >= Sig(2) >= ... >= 0
      with  sum(diag(Sig) .^ 2) = trace(V) = total variance
   Sig1(1) >= Sig1(2) >= ... >= 0
      with  sum(diag(Sig1) .^ 2) = trace(C) = n.
These singular values are, in effect, principal standard deviations.
Their squares are the eigenvalues of the covariance and correlation

(Note in passing. Put  Y = Sig * W' and Y1 = Sig1 * W1'.
Then Y and Y1 are n x n and  V = Y' * Y, C = Y1' * Y1.)

For two dimensional plots, plot the n columns of
   Y([1, 2], :) = diag( [Sig(1), Sig(2)] ) * W(:, [1, 2])'
   Y1([1, 2], :) = diag( [Sig1(1), Sig1(2)] ) * W1(:, [1, 2])',
point by point. These points represent your data. The horizontal, 1-axis
is the most significant here. It is the axis of maximum variance
(correlation?) in the data.



At 9:58 AM +0100 5/11/04, Ted Harding wrote:
 On 11-May-04 rino mailing wrote:
 > I'd like to do Principal Components Analysis with octave
 > What are the command I ave to write?
 > How to plot the result?

 You find the principal components by first finding the eigenvalues
 and eigenvectors of the covariance or correlation matrix (depending
 on whether you want the results to be independent of the units in
 which the variables are measured).

 So let C be (say) the covariance matrix. Then

   [V,D] = eig(C)

 gives you the eigenvectors as columns of V, and the eigenvalues
 as the diagonal elements of D, in increasing order of magnitude.

 So, therefore, ifr X is the N by k matrix of data, C=cov(X),
 [V,D]=eig(C), then X*V(:,k) gives you the first principal component,
 X*V(:,k-1) the second, and so on.

 I'm not sure what you mean by "plot the result". Some people want
 to see a graph which shows the proportions of variance accounted
 for by the successive principal components; for this you can plot
 the diagonal elements diag(D) of D in reverse order, perhaps in the
 form 100*diag(D)/sum(diag(D)). Or the amount of variance accounted
 for by the first, first+second, first+second+third, etc; for this
 you can plot the (reverse-order) cumsum of diag(D).

 Or you may want to say plot the first principal component against
 the second, so: PC1=X*V[:,k];PC2=X*V[:,k-1];plot(PC1,PC2)

 This may be useful in exploring possible clustering; or you may
 already have an index say Grp which identifies each row of X
 as belonging to a group by Grp==1, Grp==2, Grp==3, ... , in which
 case you can plot these subsets of (PC1,PC2) in different colours,
 to see where the groups fall in the plot.

 And so on ...

 Hoping this helps,

 E-Mail: (Ted Harding) <address@hidden>
 Fax-to-email: +44 (0)870 167 1972
 Date: 11-May-04                                       Time: 09:58:49

Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:
How to fund new projects:
Subscription information:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]