Standard PCA code, finidng the eigenvalues of a non square matrix

In a PCA algorithm, you need to find the eigenvalues of a covariance matrix in order to derive the 1st algebraic solution to PCA using linear algebra. Quick synopsis:

X is a data set (mxn).

mxn are not always different, m is # of measurement types, n is the # of sampless.

we find some orthonormal matrix P where Y = PX such that the covariance matrix C (matrix of variances between two data sets) = (normalization constant)*Y*Y transpose is diagonalized.

Now I need to find the eigenvalues of the covariance matrix but in this code, the matrix appears to not be square, which can not be if this code needs to do fulfill its purpose. Can someone tell me, who may or may not be familiar with PCA algorithm what am I misunderstanding or have incorrect? The problem is in the very last line of the code, previous is included for background.

PCA Code:

I = imread('image');

      [irow, icol] = size(I);
     data = reshape(I',irow*icol,1); 
  % [sR ,sC ,eR ,eC] = deal(1, 3, 2, 4);

% Compute the sum over the region using the integral image.

% PCA1: Perform PCA using covariance.

% data - MxN matrix of input data

% (M dimensions, N trials)

% signals - MxN matrix of projected data

% PC - each column is a PC

% V - Mx1 matrix of variances

[M,N] = size(data);

% subtract off the mean for each dimension

mn = mean(data,2);

data = double(data)

data = data - repmat(mn,1,N);

% calculate the covariance matrix

covariance = 1 / (N-1) * (data) .* (data);

% find the eigenvectors and eigenvalues

[PC, V] = eig(covariance);

[jgg]

I posted this reply on your other problem, which you should probably delete:

I'm not sure if you're doing exactly what you think. Your data matrix is is Mx1 since you're only loading in a single image. You can't do PCA on a sample of exactly one image, since this will (at most) give you a single component which is pretty useless.

For instance, at the end of the day, your covariance matrix should be square. You don't need to find the eigenvalues of non-square matrix. Here:

covariance = 1 / (N-1) * (data) .* (data);

should probably be:

 covariance = (1 / (N-1)) * (data'*data);

instead.

[Neo]

I am confused by what you mean by single image means Mx1. This describes the dimensions of a matrix, so wouldn't a single matrix be composed of some data matrix MxN or MxM? And shouldn't that line be:

covariance = (1 / (N-1)) * (data) .* (data.');

[Neo]

I understand your last statement, but I am not sure why you say "For, instance" because I am not sure where the beginning of that logic began, might you mind explaining a little more? I think I am missing something small because I thought I had a pretty good grasp of PCA.

[jgg]

You read in an image which is like 200px by 200 px. This is read in as a 200 by 200 matrix.

This code data = reshape(I',irow*icol,1); converts it into a 40,000 by 1 vector. This is the first step in doing PCA for images; take the image, then treat each pixel as a variable. The 4000 items are the variables you're going to do PCA on; each image creates one vector.

In your code, you only read in a single image, which creates one vector. You cannot do PCA on a single image for exactly this reason; the maximum number of principal components is one less than the number of images used to calibrate the PCA, which in this case is zero.

For instance, at this point in your code:

 [M,N] = size(data);

you should find that M is a large number and N is 1. This step:

 data = data - repmat(mn,1,N);

should return a vector of zeros, since you are computing mn by the average of a single observation ( data ).

I think if you think carefully about why you cannot perform PCA with a single image, you will understand what is wrong with your code.

[Walter Roberson]

Is this the same discussion as http://uk.mathworks.com/matlabcentral/answers/261394-standard-pca-code-finidng-the-eigenvalues-of-a-non-square-matrix ?

[Neo]

Yes, deleted redundant discussion on other thread.

[Neo]

jgg, I see now what you are saying. I did not realize that the vector is literally given by each "type of measurement" given by an image. Here is my revised code, I am about to run it but I am not sure if this will allow the matrix to be square so that I can now allow the calculation of the eigenvectors/eigenvalues:

 I = imread('image');
I2 = imread('image2');

    [irow, icol] = size(I);
    [irow2, icol2] = size(I2);

   data = reshape(I',irow*icol,1); 
   data2 = reshape(I2',irow*icol,1); 
% [sR ,sC ,eR ,eC] = deal(1, 3, 2, 4);
% Compute the sum over the region using the integral image.
% PCA1: Perform PCA using covariance.
% data - MxN matrix of input data
% (M dimensions, N trials)
% signals - MxN matrix of projected data
% PC - each column is a PC
% V - Mx1 matrix of variances
[M,N] = size(data);
[M2,N2] = size(data2);
% subtract off the mean for each dimension
mn = mean(data,2);
mn2 = mean(data2,2);
data = double(data);
data2 = double(data2);
data = data - repmat(mn,1,N);
data2 = data2 - repmat(mn2,1,N);
% calculate the covariance matrix
covariance = 1 / (N-1) * (data) .* (data.');
covariance2 = 1 / (N-1) * (data2) .* (data2.');
% find the eigenvectors and eigenvalues
[PC, V] = eig(covariance);
[PC2, V2] = eig(covariance2);

[Neo]

jgg?

[Christine Tobler]

I think you want to do matrix multiplication of data and data transposed. In this case, you should use

data * data'

instead of

 data .* data'

the .* stands for element-by-element multiplication.

[jgg]

Again, I still don't think you understand what you're trying to do properly. I'm pretty sure at this point: data = data - repmat(mn,1,N); you will have a matrix of zeroes.

Basically, you're trying to calculate the covaraince of a single observation; this doesn't make sense.

What you want to do is something like this:

 I = imread('image');
I2 = imread('image2');
    [irow, icol] = size(I);
    [irow2, icol2] = size(I2); %this needs to be the same as I
   d = reshape(I',irow*icol,1); 
   d2 = reshape(I2',irow*icol,1);
   data = [d, d2];
   mn = mean(data,2);
   [M,N] = size(data);
   data = data - repmat(mn,1,N);
   covariance = 1 / (M-1) * (data')*(data);

You can then calculate the eigenvectors from a transform of this, or from 1 / (N-1) * (data)*(data'); (which is usually not a good idea for images)

You probably want way more than two images, though. It's pretty impossible to isolate the principle components of a dataset with two observations.

I think that this lecture is a very clear illustration of the procedure if you want a reference.

[Star Strider]

See also: Matrix dimension must agree to use .*

This seems to be a related problem.