fttuck03.m

Tucker, PARAFAC, GRAM, RAFA and misc. 2nd order models with a test data set (old version now covered

function [A,B,C,G]=fttuck03(X,r1,r2,r3,W1,W2,W3,startmode,convlim,Ai,Bi,Ci);
%function [A B C G]=fttuck03(X,r1,r2,r3,W1,W2,W3,startmode,convlim,Ai,Bi,Ci);
%
%This TUCKER3 program calculates the W principal components in
%the problem Z = A(w1)*G(w1,w2,w3)*(B(w2)kronC(w3)) for wi=1,2,...,W(i);
%
%This particular procedure (fttuck03) normalizes A, B and C during
%the iterations.
%
% X        : Supermatrix (r1 x r2穜3) from cube (r1 x r2 x r3).
% r1,r2,r3 : Number of observations along each way.
% W1,W2,W3 : Number of factors.
% startmode: 1=Initiate A, B and C by TUCKER1/UPCA solutions.
%            2=Initiate A, B and C with Ai, Bi and Ci.
%            3=Random vectors (not a good choice)
%            In all cases G is computed from the derived X, A, B and C.
% convlim  : Max. difference between two successive sq. sums of misfits.
%            Default set to 1.0e-2 
% Ai,Bi,Ci : Initial guesses on A, B and C from the user.
% A,B,C	   : The found factors in columns. 
%            The dimensions will be A(r1 x W1) B(r2 x W2) C(r3 x W3).
% G        : Cube of interactions, the CORE
% Author   : Claus A. Andersson, May 1995 
% Copyright: Food Technology,
%            Royal Veterinary & Agricultural University
%            Copenhagen, Denmark
% E-mail   : [e-mail] 

%Initialize A, B and C
if startmode==1, 
  disp('Initializing A, B and C with TUCKER1/UPCA solutions');
  disp('Calling FTUPCA01.M')
  [A B C]=ftupca01(X,r1,r2,r3,W1,W2,W3);  
end;
if startmode==2,
  disp('Initializing A, B and C as defined by user.')
  A=Ai;
  B=Bi;
  C=Ci;
end; 
if startmode==3,
  disp('Initializing A, B and C by random vectors.')
  A=rand(r1,W1)+1;
  B=rand(r2,W2)+1;
  C=rand(r3,W3)+1;
end;

%Updating G - the core
tmp1=zeros(W1*W2,r3);
for t=1:r3,
  tmp=(A'*X(1:r1,(t-1)*r2+1:t*r2)*B)';
  tmp1(:,t)=tmp(:);
end;
tmp2=kron(A'*A,B'*B)\tmp1*C/(C'*C);
G=zeros(W1,W2*W3);
for t=1:W3,
  G(1:W1,(t-1)*W2+1:t*W2)=reshape(tmp2(:,t)',W2,W1)';
end;
clear tmp tmp1 tmp2 t

%Calculates the misfit of factors and core
diff=0;
for k=1:r3
  tmp1=G(:,1:W2)*C(k,1);
  for l=2:W3
    tmp1=tmp1+G(:,(l-1)*W2+1:l*W2)*C(k,l);
  end;
  diff=diff+sum(sum((X(:,(k-1)*r2+1:k*r2)-A*tmp1*B').^2));
end;
fprintf('The sq. sum of misfit is %12.8f at start. \n',diff);

%Iteration starts
if exist('convlim')==0,
  convlim=1e-2;
end;
if convlim<eps,
  convlim=1e-2;
end;
it=0;
converged=0;
olddiff=diff;
while converged==0,
  it=it+1;

  %Updating A
  tmp2=zeros(W1,r3*r2);
  for k=1:r3,
    tmp1=C(k,1)*G(:,1:W2)*B';
    for l=2:W3,
      tmp1=tmp1+C(k,l)*G(:,(l-1)*W2+1:l*W2)*B';
    end;
    tmp2(1:W1,(k-1)*r2+1:k*r2)=tmp1;
  end;
  A=X*tmp2'/(tmp2*tmp2');
  for t=1:W1,
    A(:,t)=A(:,t)/norm(A(:,t));
  end;

  %Updating B
  tmp1=zeros(r1*r3,W2);
  for k=1:r3,
    tmp2=A*C(k,1)*G(:,1:W2);
    for l=2:W3,
      tmp2=tmp2+A*C(k,l)*G(:,(l-1)*W2+1:l*W2);
    end;
    tmp1((k-1)*r1+1:k*r1,1:W2)=tmp2;
  end;
  Z=zeros(r1*r3,r2);
  for t=1:r3,
    Z((t-1)*r1+1:t*r1,1:r2)=X(1:r1,(t-1)*r2+1:t*r2);
  end;
  B=Z'*tmp1/(tmp1'*tmp1);
  for t=1:W2,
    B(:,t)=B(:,t)/norm(B(:,t));
  end;

  %Updating C
  Y=zeros(W3,W3);
  Z=zeros(r3,W3);
  for t=1:W3,
    for u=1:W3,
      Y(t,u)=trace(G(:,(t-1)*W2+1:t*W2)'*A'*A*G(:,(u-1)*W2+1:u*W2)*B'*B);
    end;
  end;
  for t=1:r3,
    for u=1:W3,
      Z(t,u)=trace(G(:,(u-1)*W2+1:u*W2)'*A'*X(:,(t-1)*r2+1:t*r2)*B);
    end;
  end;
  C=Z/Y;
  for t=1:W3,
    C(:,t)=C(:,t)/norm(C(:,t));
  end;

  %Updating G - the core
  tmp1=zeros(W1*W2,r3);
  for t=1:r3,
    tmp=(A'*X(1:r1,(t-1)*r2+1:t*r2)*B)';
    tmp1(:,t)=tmp(:);
  end;
  tmp2=kron(A'*A,B'*B)\tmp1*C/(C'*C);
  for t=1:W3,
    G(1:W1,(t-1)*W2+1:t*W2)=reshape(tmp2(:,t)',W2,W1)';
  end;

  %Calculates the misfit of factors and core
  diff=0;
  for k=1:r3
    tmp1=G(:,1:W2)*C(k,1);
    for l=2:W3
      tmp1=tmp1+G(:,(l-1)*W2+1:l*W2)*C(k,l);
    end;
    diff=diff+sum(sum((X(:,(k-1)*r2+1:k*r2)-A*tmp1*B').^2));
  end;
  fprintf('The sq. sum of misfit is %12.8f after %i iterations. \n',diff,it);

  %Checks if a minimum has been found
  if abs(olddiff-diff)<convlim, converged=1; end,
  olddiff=diff;

end;

disp('FTTUCK03: Execution ended.')