use PCA to reduce the dimensionality of hand written digit data from 256 to 2, and achieve better classification rate with 1-Nearest_Neighbor classifier

Author: Robert

bonus.m

@@ -0,0 +1,24 @@
+clear all;
+load digits;
+
+
+m = [2,5,10,20];
+numIter = size(m,2);
+avgErrTest = zeros(numIter,1);
+
+for iter=1:numIter
+    [avgErrTest(iter,1)] = pca_classifier(m(1,iter));    
+end
+
+
+
+figure;
+hold on;
+
+plot(m,avgErrTest,'r-');
+
+title('avg classification error rate VS PC');
+xlabel('number of Principle Component');
+ylabel('avg classification error rate');
+
+

digits.mat

@@ -0,0 +1,17 @@
+function [label_test] = knn(k, data_train, label_train, data_test)
+
+error(nargchk(4,4,nargin));
+
+D = size(label_train, 1);
+
+dist = l2_distance(data_train, data_test);
+[sorted_dist, nearest] = sort(dist);
+nearest = nearest(1:k,:);
+
+label_test = zeros(D, size(data_test, 2), k);
+for i=1:k
+    label_test(:,:,i) = label_train(:, nearest(i, :));
+end
+
+label_test = mean(label_test,3);
+label_test = label_test == repmat(max(label_test, [], 1), D, 1);

knn.m

@@ -0,0 +1,72 @@
+function d = L2_distance(a,b,df)
+% L2_DISTANCE - computes Euclidean distance matrix
+%
+% E = L2_distance(A,B)
+%
+%    A - (DxM) matrix 
+%    B - (DxN) matrix
+%    df = 1, force diagonals to be zero; 0 (default), do not force
+% 
+% Returns:
+%    E - (MxN) Euclidean distances between vectors in A and B
+%
+%
+% Description : 
+%    This fully vectorized (VERY FAST!) m-file computes the 
+%    Euclidean distance between two vectors by:
+%
+%                 ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
+%
+% Example : 
+%    A = rand(400,100); B = rand(400,200);
+%    d = distance(A,B);
+
+% Author   : Roland Bunschoten
+%            University of Amsterdam
+%            Intelligent Autonomous Systems (IAS) group
+%            Kruislaan 403  1098 SJ Amsterdam
+%            tel.(+31)20-5257524
+%            [email protected]
+% Last Rev : Wed Oct 20 08:58:08 MET DST 1999
+% Tested   : PC Matlab v5.2 and Solaris Matlab v5.3
+
+% Copyright notice: You are free to modify, extend and distribute 
+%    this code granted that the author of the original code is 
+%    mentioned as the original author of the code.
+
+% Fixed by JBT (3/18/00) to work for 1-dimensional vectors
+% and to warn for imaginary numbers.  Also ensures that 
+% output is all real, and allows the option of forcing diagonals to
+% be zero.  
+
+if (nargin < 2)
+   error('Not enough input arguments');
+end
+
+if (nargin < 3)
+   df = 0;    % by default, do not force 0 on the diagonal
+end
+
+if (size(a,1) ~= size(b,1))
+   error('A and B should be of same dimensionality');
+end
+
+if ~(isreal(a)*isreal(b))
+   disp('Warning: running distance.m with imaginary numbers.  Results may be off.'); 
+end
+
+if (size(a,1) == 1)
+  a = [a; zeros(1,size(a,2))]; 
+  b = [b; zeros(1,size(b,2))]; 
+end
+
+aa=sum(a.*a); bb=sum(b.*b); ab=a'*b; 
+d = sqrt(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);
+
+% make sure result is all real
+d = real(d); 
+
+% force 0 on the diagonal? 
+if (df==1)
+  d = d.*(1-eye(size(d)));
+end

l2_distance.m

@@ -0,0 +1,15 @@
+function [ avg_err_rate ] = pca_classifier( k )
+
+load digits;
+
+[base,mean,data_train] = pcaimg([train2,train3],k);
+[base,mean,data_test] = pcaimg([test2,test3],k);
+
+label_train = [ones(1, 300) zeros(1,300); zeros(1,300) ones(1,300) ];
+label_test = [ones(1, 300) zeros(1,300); zeros(1,300) ones(1,300) ];
+
+result = knn(1, data_train, label_train, data_test);
+total = sum( abs( label_test(1,:) - result(1,:)));
+avg_err_rate = total/600;
+end
+

pca_classifier.m

@@ -0,0 +1,28 @@
+function [base,mean,projX] = pcaimg(X,k)
+%% pca analysis of image data
+%% input: X: vectorized image data
+%%        k: number of eigenvectors you want to keep
+%% output: base: eigenvectors
+%%         mean: mean of data
+%%         projX: the projected data in the low-dimensional 
+%%                space.
+
+disp('eigendecomposition...');
+[xdim,ndata] = size(X);
+mean = sum(X,2)/ndata;  % compute mean of data
+X = X-repmat(mean,1,ndata); % substract the mean
+cov = X*X'/ndata;  % form the covariance matrix
+
+[ev,ed]=eig(cov);  %eigendecomposition
+ed = diag(ed);
+
+% Sort eigenvectors by eigenvalues
+[foo,p]=sort(-ed);
+ed = ed(p);
+ev = ev(:,p);
+
+% Take the top k eigenvectors
+base = ev(:,1:k);
+
+% project the data into low-dim space
+projX = base'*X;