% ccentroid.m
%
% Jason R. Blevins <jrb11@duke.edu>
%
% An implementation of the classical centroid method, algorithm 5.1 of
% Chu and Funderlic, SIAM J. MATRIX ANAL. APPL., 23, 1025-1044, 2002.
%
% $Id: ccentroid.m,v 1.4 2004/08/06 20:34:54 jrblevin Exp $

function [B,V,Z] = ccentroid(Y)

[m,l] = size(Y);

max_iter = 100;

B = [];
V = [];
Z = [];

p = rank(A);

R = Y*Y'/l;

for i = 1:min(p, max_iter)
    [R_new, z, b, v] = ccentroid_step(R);

    B = [B,b];
    V = [V,v];
    Z = [Z,z];
    R = R_new;
end


function [R_new,z,b,v] = ccentroid_step(R)

[m,m] = size(R);

z = ones(n,1);

% machine zero threshold
tol = 10^(-13);

P = R - diag(diag(R));

w = P*z;

% Loop until the corresponding components of w and z have the same
% sign.
while ( sign(w)'*sign(z) < n )
  k = -1;
  
  % Step 1
  for j = 1:n
    if ( sign(w(j)) ~= sign(z(j)) )
      if ((k == -1) | ( abs(w(j)) > abs(w(k)) ))
	k = j;
	
	% Step 2
	z(k) = -z(k);
        
	% Step 3
	w = w + 2*sign(z(k))*P(:,k);
      end
    end
  end
end

% Equation 4.4
v = A'*z/norm(A'*z);
% The loading vector is given by equation 4.6
b = (R*z)/(sqrt(z'*R*z));
% The new product moment matrix of the orthogonally reduced loading
% matrix A can be found using equation 4.7.
R_new = R - (R*z*z'*R)/(z'*R*z);