Logistic PCA

logpca.m

% Adapted from the original, available online: 
% http://www.andrewschein.com/software/lpca_code.tar
% 

function [U,V] = logpca(X,k,tol,imax)
  % ensure we have a binary matrix.
  X = spones(X);
  
  [N,D] = size(X);
  U     = rand(N,k);
  V     = rand(k,D);
  
  oldU  = U;
  oldV  = V;
  oldLL = 0;
  
  for iter = 1:imax

    % Begin U update
    b = 2*X*V' - repmat(sum(V,2)', size(X,1), 1);

    parfor n = 1:N
      Z = U(n,:)*V;
      T = tanh(Z/2) ./ Z;

      A      = (repmat(T,k,1).*V) * V';
      U(n,:) = (A \ b(n,:)')';
    end

    % Begin V update.
    b = 2*X'*U - repmat(sum(U,1), size(X,2), 1);
    
    % Update log likelihood.
    newLL = 0;

    parfor d = 1:D
      Z = U*V(:,d);
      T = tanh(Z/2) ./ Z;

      A      = (repmat(T,1,k) .* U)' * U;
      V(:,d) = (A \ b(d,:)');
      
      ix0 = X(:,d) == 0;
      ix1 = X(:,d) == 1;
      z   = U*V(:,d);

      % update LL to remove 
      newLL = newLL - (...
        sum(log(1+exp(-z(ix1)))) + ...
        sum(log(1+exp(z(ix0)))));
    end
    
    fprintf('Iteration %f\tLog-likelihood: %f\n', iter, newLL);
    
    if newLL < 0 && oldLL < 0 && abs(oldLL - newLL) < tol
      break;
    end
    
    % detect numerical errors.
    if isnan(newLL)
      U = oldU;
      V = oldV;
      break;
    end
    
    oldU    = U;
    oldV    = V;
    oldLL   = newLL;
  end
end