#
#  function [A,X,Sigma]=als2a(H,maxit,a00,epsilon)
#
#   Normalized scores algorithm (Unequal variances)
#
#   Input:
#
#   H       Data matrix.
#
#   maxit, a00, epsilon, are optional arguments, to control
#           iteration.
#
#   Output:
#
#   A       weights
#   X       scores
#   Sigma   loss function
#
#           these three matrices contain all the intermediate
#           values, till convergence. The i-th column corresponds
#           to the i-th iteration.
#
function [A,X,Sigma]=als2a(H,maxit,a00,epsilon)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_epsilon=1.0e-4;
   default_maxit=50;
   default_a0=ones(m,1);
# ------------------------------------------------------------------
#
#   Compute D (the diagonal of H'*H)
#
   D=diag(diag(H'*H));
#
#   Check if maxit (=maximum number of iterations) was entered.
#   Otherwise, set a default value.
#
#
   if nargin==1,
      maxit=default_maxit;
   endif
#
#   Check if a0 (=initial weight vector) was entered.
#   if so, check its dimensions.
#   Otherwise, set a default value.
#
   if nargin>=3
      [m1,m2]=size(a00);
      if m1==m && m2==1,
      a0=a00;
      else
         if m1==1 && m2==m,
            a0=a00';
         else
            disp('*** als2a: Wrong sized a0. Set default');
            a0=default_a0;
         endif
      endif
   else
      a0=default_a0;
   endif
#
#   Initialize matrices for storing all the iterations.
#
   A=zeros(m,maxit);
   X=zeros(n,maxit);
   Sigma=zeros(1,maxit);
#
#   Initialize a and x
#
   x=zeros(n,1);
   a=a0;
   sigma=norm(H,'fro')/m;
   A(:,1)=a;
   X(:,1)=x;
#
   if nargin<4
      epsilon=default_epsilon;
   endif
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=2:maxit
         x=H*a;
         lx=norm(x);
         x=x/lx;
         a=inv(D)*H'*x;
         sigma=norm(x*a'-H,'fro')/m;
         A(:,i)=a;
         X(:,i)=x;
         Sigma(i)=sigma;
         if norm(A(:,i-1)-a)<epsilon
            disp('*** als2a: Convergence for a attained');
            break
         endif
    endfor
    num_iter=i;
    A=A(:,1:num_iter);
    X=X(:,1:num_iter);
    Sigma=Sigma(:,1:num_iter);
endfunction