function [eigvect,eigval,Kt]=kpca_n(K,n)
% [eigvect,eigval,Kt] = kpca_n( K, n )
%
%   Computes the first n principal components of the
%   kernel matrix K.  Uses some computational tricks to speed
%   things up.
%
%   Also returns Kt, which is the kernel matrix which has
%   been centered in feature space.
%
%   Eigenvalues are sorted in descending order, and their
%   corresponding eigenvectors are also sorted.  Eigenvectors
%   are returned as column vectors.
%
%   D. Wingate
%

nx = size( K, 1 );

%
% Center the kernel matrix.
%
% This is the fast way to do it -- we 
% leverage the fact that the necessary matrix multiplies
% are actually rank-one!
%

o = ones(nx,1);
oon = 1/nx;
ko = oon*o*(o'*K);
Kt = K - ko - ko' + oon*oon * o*(o'*(K*o))*o';
Kt = min( Kt, Kt' ); % ensure symmetry

%
%  Find the eigenvectors
%

%[ eigvect, eigval ] = jdqr( Kt, n, 'LM' ); % diagonalizing Kt;
[ eigvect, eigval ] = eigs( Kt, n, 'LM' ); % diagonalizing Kt;

% just in case we had some slightly negative eigenvalues...
eigvect = real(eigvect);
eigval = real(eigval);

eigval = diag( eigval );

%
% eigenvectors are already normalized to unit length...
% we want it so that 1 = lambda_i * e_i' * e_i.
%

nbeigval = length( eigval );
for i=1:nbeigval % normalizing eigenvector
    eigvect(:,i)=eigvect(:,i)/sqrt(eigval(i));
end;

%
% sort the results
%

[ aux, ind ] = sort( -eigval );
eigval = eigval( ind );
eigvect = eigvect( :, ind );