function [y,w,e]=crls(x,M,epoch)
% This function extracts the first m principal components from x
% using the cascade least square neural network
% x is with dimension RxN, each row is an input vector.
% [y,w,e]=crls(x,M)
% M is the number of principal components desired, 
% y is the principal components, w is the principal vectors, e is the
% resudal error. y=wx
% Copyright by Ze Wang, all rights reserved by Shanghai Jiao Tong University. 2002
% If you use this code, please cite:
%     Ze Wang, Yin Lee, C.S. Leung, Y.S. Zhu, An Improved Sequential Method for Principal Component Analysis. Pattern Recognition Letters, vol. 24/9-10, pp. 1409 – 1415, 2003.
% 
if nargin<1
    disp('[y,w,e]=crls(x,M)');
    return;
end
try epoch;
    ;
catch
    epoch=100;
end
[R,N]=size(x);
epsilon=1e-5;
e=x;
y=zeros(M,N);
for j=1:M
    w(j,:)=eye(1,R);
    %w(j,:)=rand(1,R);
    %w(j,:)=w(j,:)./norm(w(j,:));
    if j>1
        w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
    end
    w(j,:)=w(j,:)./norm(w(j,:));
    eta=sum(sum(e.^2))/N;
    delta=200;
    for s=1:epoch
        
        for k=1:N
            
            y(j,k)=w(j,:)*e(:,k);
            eta=y(j,k).^2+eta;
            prew=w(j,:);
            w(j,:)=w(j,:)+(y(j,k)/eta)*(e(:,k)'-w(j,:)*y(j,k));
            delta=norm(w(j,:)-prew);
            if delta<epsilon
                sprintf('delta(w) is now less than epsilon for the %d-th component at epoch %d, k=%d',j,s,k)
                %
                %
                %break;
            end
%             if j>1
%                 w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
%             end
%  without this, the convergence is faster
        end
        if delta<epsilon
             %   sprintf('delta(w) is now less than epsilon,k=%d',k);
                %
                %
                %break;
        end
        if j>1
             w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
        end

    end
    if j>1
        w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
    end
    w(j,:)=w(j,:)./norm(w(j,:));
    y(j,:)=w(j,:)*e;
    e=e-w(j,:)'*y(j,:);
end
% another loop to reduce the error
% e=x;
% for j=1:M
%     eta=sum(sum(e.^2))/N;
%     delta=200;
%     for s=1:100
%         
%         for k=1:N
%             
%             y(j,k)=w(j,:)*e(:,k);
%             eta=y(j,k).^2+eta;
%             prew=w(j,:);
%             w(j,:)=w(j,:)+(y(j,k)/eta)*(e(:,k)'-w(j,:)*y(j,k));
%             %w(j,:)=w(j,:)+(y(j,k)/eta)*(e(:,k)');
%             delta=norm(w(j,:)-prew);
%             if delta<epsilon
%                 sprintf('delta(w) is now less than epsilon for the %d-th component at epoch %d, k=%d',j,s,k)
%                 %
%                 %
%                 break;
%             end
%         end
%         if delta<epsilon
%              %   sprintf('delta(w) is now less than epsilon,k=%d',k);
%                 %
%                 %
%                 break;
%         end
%         if j>1
%              w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
%         end
% 
%     end
%     if j>1
%         w(j,:)=w(j,:)-w(j,:)*w(1:j-1,:)'*w(1:j-1,:);
%     end
%     w(j,:)=w(j,:)./norm(w(j,:));
%     y(j,:)=w(j,:)*e;
%     e=e-w(j,:)'*y(j,:);
% end
%
%
