#
#  function [A,X,Sigma]=als2b(H,maxit,x00,epsilon)
#
#   Normalized weights algorithm (Unequal variances)
#
#   Input:
#
#   H       Data matrix.
#
#   maxit, x00, 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]=als2b(H,maxit,x00,epsilon)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_epsilon=1.0e-4;
   default_maxit=50;
   default_x0=ones(n,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 x0 (=initial scores vector) was entered.
#   if so, check its dimensions.
#   Otherwise, set a default value.
#
   if nargin>=3
      [n1,n2]=size(x00);
      if n1==n && n2==1,
      x0=x00;
      else
         if n1==1 && n2==n,
            x0=x00';
         else
            disp('*** als2b: Wrong sized x0. Set default');
            x0=default_x0;
         endif
      endif
   else
      x0=default_x0;
   endif
#
#   Initialize matrices for storing all the iterations.
#
   A=zeros(m,maxit);
   X=zeros(n,maxit);
   Sigma=zeros(1,maxit);
#
#   Initialize a and x
#
   a=zeros(m,1);
   x=x0;
   sigma=norm(H,'fro')/m;
   A(:,1)=a;
   X(:,1)=x;
#
   if nargin<4
      epsilon=default_epsilon;
   endif
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=2:maxit
         a=inv(D)*H'*x;
         la=sqrt(a'*D*a);
         a=a/la;
         x=H*a/m;
         sigma=norm(x*a'-H,'fro')/m;
         A(:,i)=a;
         X(:,i)=x;
         Sigma(i)=sigma;
         if norm(A(:,i-1)-a)<epsilon
            disp('*** als2b: 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