Regularized Alternating Least Squares

als.m

function [U,V] = als(X,k,tol,imax,lambda)

  [d,n] = size(X);
  U     = rand(k,d);
  V     = rand(k,n);
  I     = eye(k);

  % convergence monitors.
  olderr = 0;
  
  for iter = 1:imax
    
    % update V.
    parfor j = 1:n
      ix  = find(X(:,j) > 0);
      num = numel(ix);

      if num == 0
        continue;
      end

      x = X(ix,j)';
      Q = U(:,ix);

      V(:,j) = (x*Q') / (Q*Q' + lambda*num*I);
    end
    
    % update cost
    newerr = 0;
    
    % update U.
    parfor i = 1:d
      ix  = find(X(i,:) > 0);
      num = numel(ix);

      if num == 0
        continue;
      end

      x = X(i,ix);
      Q = V(:,ix);

      U(:,i) = (x*Q') / (Q*Q' + lambda*num*I);
      
      err     = U(:,i)'*V(:,ix);
      newerr  = newerr + sum((X(i,ix) - err).^2);
    end
    
    % regularization terms.
    sums = sum(X > 0, 2);
    for i = 1:d
      newerr = newerr + ...
        sums(i,1) * lambda * U(:,i)'*U(:,i);
    end

    sums = sum(X > 0, 1);
    for j = 1:n
      newerr = newerr + ...
        sums(1,j) * lambda * V(:,j)'*V(:,j);
    end
    
    fprintf('Iteration %f\tCost function: %f\n', iter, newerr);
    
    if newerr > 0 && olderr > 0 && olderr - newerr < tol
      break;
    end
    
    olderr = newerr;

  end

  U = U';
  
end