pca.rar_Principle Component_pca_pca-sift

function [new_data, eigvector, eigvalue] = pca(data, eigvector_n) %PCA Principal component analysis % Usage: % [NEW_DATA, EIGVECTOR, EIGVALUE] = PCA(DATA, EIG_N) % % DATA: Rows of vectors of (zero-mean) data points % EIG_N: No. of selected eigenvectors % NEW_DATA: The data after projection % EIGVECTOR: Each column of this matrix is a eigenvector of DATA'*DATA % EIGVALUE: Eigenvalues of DATA'*DATA % % Note that DATA must be zero-mean'ed before calling this function. % % Type "pca" for a self-demo. % Roger Jang, 970406, 990612, 991215 if nargin == 0, selfdemo; return; end if nargin < 2, eigvector_n = min(size(data)); end if size(data, 1) >= size(data, 2), [eigvector, d] = eig(data'*data); eigvalue = diag(d); % ====== Sort based on descending order [junk, index] = sort(-eigvalue); eigvalue = eigvalue(index); eigvector = eigvector(:, index); new_data = data*eigvector; if eigvector_n < size(data, 2), eigvalue = eigvalue(1:eigvector_n); eigvector = eigvector(:, 1:eigvector_n); new_data = new_data(:, 1:eigvector_n); end else % This is an efficient method which computes the eigvectors of % of A*A^T (instead of A^T*A) first, and then convert them back to % the eigenvectors of A^T*A. [eigvector1, d] = eig(data*data'); eigvalue = diag(d); % ====== Sort based on descending order [junk, index] = sort(-eigvalue); eigvalue = eigvalue(index); eigvector1 = eigvector1(:, index); % Eigenvectors of A*A^T eigvector = data'*eigvector1; % Eigenvectors of A^T*A eigvector = eigvector*diag(1./(sum(eigvector.^2).^0.5)); % Normalization new_data = data*eigvector; if eigvector_n < size(data, 1), eigvalue = eigvalue(1:eigvector_n); eigvector = eigvector(:, 1:eigvector_n); new_data = new_data(:, 1:eigvector_n); end end function selfdemo % ====== Demo for 2D data data_n = 1000; x = 5+randn(data_n, 1); y = 10+randn(data_n, 1)/2; x_mean = mean(x); y_mean = mean(y); theta = pi/6; z = ((x-x_mean) + i*(y-y_mean))*exp(i*theta)+x_mean+i*y_mean; x = real(z); y = imag(z); plot(x, y, '.'); axis image; t = [x-x_mean y-y_mean]; [new_data, v, d] = feval(mfilename, t); v1 = -v(:, 1); v2 = -v(:, 2); arrow_x = [-1 0 nan -0.1 0 -0.1]+1; arrow_y = [0 0 nan 0.1 0 -0.1]; arrow = arrow_x + j*arrow_y; arrow1 = 2*arrow*(v1(1)+i*v1(2))*d(1)/data_n+x_mean+i*y_mean; arrow2 = 2*arrow*(v2(1)+i*v2(2))*d(2)/data_n+x_mean+i*y_mean; line(real(arrow1), imag(arrow1), 'color', 'c', 'linewidth', 2); line(real(arrow2), imag(arrow2), 'color', 'g', 'linewidth', 2); title('Axes for PCA'); % ====== Demo for Iris data load iris.dat data = iris(:, 1:4); % Take only the input part data_n = size(data, 1); data = data - ones(data_n, 1)*mean(data); % Make data zero-mean [new_data, eigvector, eigvalue] = feval(mfilename, data); index1 = find(iris(:,5)==1); index2 = find(iris(:,5)==2); index3 = find(iris(:,5)==3); figure; plot( new_data(index1, 1), new_data(index1, 2), '*', ... new_data(index2, 1), new_data(index2, 2), 'o', ... new_data(index3, 1), new_data(index3, 2), 'x'); legend('Class 1', 'Class 2', 'Class 3'); title('PCA projection of IRIS data onto the first 2 eigenvectors'); axis equal; axis tight; figure; plot( new_data(index1, 3), new_data(index1, 4), '*', ... new_data(index2, 3), new_data(index2, 4), 'o', ... new_data(index3, 3), new_data(index3, 4), 'x'); legend('Class 1', 'Class 2', 'Class 3'); title('PCA projection of IRIS data onto the last 2 eigenvectors'); axis equal; axis tight;