Daniel@0: function [evals, evec] = eigdec(x, N)
Daniel@0: %EIGDEC	Sorted eigendecomposition
Daniel@0: %
Daniel@0: %	Description
Daniel@0: %	 EVALS = EIGDEC(X, N computes the largest N eigenvalues of the
Daniel@0: %	matrix X in descending order.  [EVALS, EVEC] = EIGDEC(X, N) also
Daniel@0: %	computes the corresponding eigenvectors.
Daniel@0: %
Daniel@0: %	See also
Daniel@0: %	PCA, PPCA
Daniel@0: %
Daniel@0: 
Daniel@0: %	Copyright (c) Ian T Nabney (1996-2001)
Daniel@0: 
Daniel@0: if nargout == 1
Daniel@0:    evals_only = logical(1);
Daniel@0: else
Daniel@0:    evals_only = logical(0);
Daniel@0: end
Daniel@0: 
Daniel@0: if N ~= round(N) | N < 1 | N > size(x, 2)
Daniel@0:    error('Number of PCs must be integer, >0, < dim');
Daniel@0: end
Daniel@0: 
Daniel@0: % Find the eigenvalues of the data covariance matrix
Daniel@0: if evals_only
Daniel@0:    % Use eig function as always more efficient than eigs here
Daniel@0:    temp_evals = eig(x);
Daniel@0: else
Daniel@0:    % Use eig function unless fraction of eigenvalues required is tiny
Daniel@0:    if (N/size(x, 2)) > 0.04
Daniel@0:       [temp_evec, temp_evals] = eig(x);
Daniel@0:    else
Daniel@0:       options.disp = 0;
Daniel@0:       [temp_evec, temp_evals] = eigs(x, N, 'LM', options);
Daniel@0:    end
Daniel@0:    temp_evals = diag(temp_evals);
Daniel@0: end
Daniel@0: 
Daniel@0: % Eigenvalues nearly always returned in descending order, but just
Daniel@0: % to make sure.....
Daniel@0: [evals perm] = sort(-temp_evals);
Daniel@0: evals = -evals(1:N);
Daniel@0: if ~evals_only
Daniel@0:    if evals == temp_evals(1:N)
Daniel@0:       % Originals were in order
Daniel@0:       evec = temp_evec(:, 1:N);
Daniel@0:       return
Daniel@0:    else
Daniel@0:       % Need to reorder the eigenvectors
Daniel@0:       for i=1:N
Daniel@0:          evec(:,i) = temp_evec(:,perm(i));
Daniel@0:       end
Daniel@0:    end
Daniel@0: end