Why does the first principal component show least difference in my PCA?

Question

I am using PCA on a set of facial images in MatLab. Creating an average face and randomizing others are working fine.

In my function vectorComparison I want to see the difference on each principal component vector when using the standard deviation. But when I use eig_face_index = 1 I see less of a difference than when I use 2, or 3 etc.

The higher indexes also seem to add more color, which could be due to noise in the eigenfaces, as I am using RGB space.

Why does my initial vector show the least difference. Shouldn't it be the other way around?

Here is all of the code I am using:

main.m

clear;clc;close all;
[imvecs,img] = loadImages();
meanval = meanValue(imvecs);
[T, D] = covarianceMatrix(imvecs, meanval);
[eigvecs, eigvals] = findEigVecs(imvecs, T, D);
eigenfaces = createEigenFaces(eigvecs, imvecs, img); 

%%
[mean_image] = createAverageFace(meanval, img);
%%
[stdev_vec] = createRandomFace(eigvals, eigvecs, imvecs, meanval, img);
%%
vectorComparison(meanval, eigvecs, stdev_vec, img, mean_image);

loadImages.m

function [imvecs,img] = loadImages()
images = dir('D:My\Path\*.png');  
imgPath = 'D:My\Path\';
img=imread([imgPath images(1).name]);
n=length(images);

for i = 1:n
    img = imread([imgPath images(i).name]);
    imvecs{i} = double(img(:));
end
return

meanValue.m

function meanval = meanValue(imvecs, imageNr)
    %Creates the mean value from our images.
    sumvec=imvecs{1};

    for i = 2:(size(imvecs,2))
        sumvec = sumvec + imvecs{i};
    end
    meanval = sumvec ./(size(imvecs,2));
return

covarianceMatrix.m

function [T, D] = covarianceMatrix(imvecs, meanval)
    D = [];
    for i = 1:size(imvecs,2),
       diff = imvecs{i} - meanval;
       D = [D, diff];
    end
    %Dimensionality reduction
    T = (D' * D) ./ (size(imvecs,2));
return

findEigVecs.m

function [eigvecs, eigvals] = findEigVecs(imvecs, T, D)
    [U,eigvals,V] = svd( T );
    eigvecs = [];
    for i = 1:size(imvecs,2),
        eigvec = D * U(:,i);
        eigvec = eigvec ./ sum(eigvec);
        eigvecs = [eigvecs, eigvec];
    end
return

createEigenFaces.m

function [eigenfaces] = createEigenFaces(eigvecs, imvecs, img)
    for i = 1:size(imvecs,2),
        eigface =  reshape(eigvecs( : , i), size(img));
        eigface = eigface - min(min(min((eigface))));
        eigface = eigface ./ max(max(max((eigface))));
        eigenfaces{i}=eigface;
        %figure;imagesc(eigface);
    end
return

createAverageFace.m

function [mean_image] = createAverageFace(meanval, img)
    mean_image = reshape(meanval, size(img));
    figure;imagesc(mean_image./255);
    title('Average Face')
return

createRandomFace.m

function [stdev_vec] = createRandomFace(eigvals, eigvecs, imvecs, meanval, img)
    stdev_vec = sqrt(diag(eigvals));
    t = (100 * rand(size(imvecs,2),1) - 50) .* stdev_vec;
    new_face1 = meanval + (eigvecs * t);
    new_face1 = reshape(new_face1, size(img));

    figure;imagesc(new_face1./255);
    title('Random Face')
return

vectorComparison.m

function [] = vectorComparison(meanval, eigvecs, stdev_vec, img, mean_image)
    t = zeros(17,1);
    eig_face_index = 1;
    t(eig_face_index) = 1000;
    t = t.*stdev_vec;

    new_face1 = meanval + (eigvecs * t);
    new_face1 = reshape(new_face1, size(img));

    new_face2 = meanval - (eigvecs * t);
    new_face2 = reshape(new_face2, size(img));

    figure;
    title('PCA Comparison')
    subplot(3,1,1), subimage(new_face1./255)
    subplot(3,1,2), subimage(mean_image./255)
    subplot(3,1,3), subimage(new_face2./255)
return

Answer 1

I found what was wrong here. In my function findEigVals I make each vector a unit vector (length 1) by deviding by the sum of the vector itself. This is possible as long as the content of the vector does range positive entries. However, as we cannot know if we have a vector in pos or neg direction (both equally valid), we cannot use sum here.

Instead we need to replace this with norm, matlab's way of stating normalization.

function [eigvecs, eigvals] = findEigVecs(imvecs, T, D)
    [U,eigvals,V] = svd( T );
    eigvecs = [];
    for i = 1:size(imvecs,2),
        eigvec = D * U(:,i);
        eigvec = eigvec ./ norm(eigvec);
        eigvecs = [eigvecs, eigvec];
    end
return

If anyone is using a version of the code above, note that the random face function will give too strong values in rgb-space. Replace the createRandomFace with the following:

stdev_vec = sqrt(diag(eigvals));
min_range = 0;
max_range = 2;
t = ((max_range - min_range)*rand(n,1)) .* stdev_vec + min_range;
new_face1 = meanval + (eigvecs * t);
new_face1 = reshape(new_face1, size(img));

figure;imagesc(new_face1./255);