Exponential Family PCA (Poisson Distribution)

exppca.m

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

  [d,n] = size(X);
  
  % random initialization.
  U = rand(d,k);
  V = rand(k,n);
  I = eye(k);
  
  % for numerical errors.
  oldU  = U;
  oldV  = V;
  
  % monitor convergence.
  olderr = 0;
  
  for iter = 1:imax
    
    % update U.
    parfor i = 1:d
      x  = X(i,:);
      
      h      = U(i,:)*V;
      g      = exp(h);
      D      = diag(g);
      num    = (h + (x - g) / D) * D*V';
      dom    = V*D*V' + 10^-5*I;
      U(i,:) = num / dom;
    end
    
    % monitor cost.
    newerr = 0;
    
    % update V.
    parfor j = 1:n
      x  = X(:,j);
      
      h      = U*V(:,j);
      g      = exp(h);
      D      = diag(g);
      num    = U'*D * (h + D \ (x - g));
      dom    = U'*D*U + 10^-5*I;
      V(:,j) = dom \ num;
      
      % update cost.
      q      = U*V(:,j);
      newerr = newerr + sum(log(1 ./ factorial(x)) + (x.*q) - exp(q));
    end
    
    fprintf('Iteration %f\tCost function: %f\n', iter, newerr);
    
    if newerr < 0 && olderr < 0 && newerr - olderr < tol
      break;
    end
    
    % detect numerical errors.
    if isnan(newerr)
      U = oldU;
      V = oldV;
      break;
    end

    oldU    = U;
    oldV    = V;
    olderr  = newerr;
  end

end