tals.m

多维数据处理：MATLAB源程序用于处理多维数据

function [estA,estH,estST]=TALS(X,M,estA,estH,estST,MAXNUMITER);

% Plain Vanilla Trilinear ALS
% X = KxNxP rank-M 3-way array

% Copyright, 1998 -
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. In case of doubt,
% contact the holder of the copyrights.
%
% Copyright 1998
% Nikos Sidiropoulos
% Univ. Minnesota, [e-mail]
% and 
% Rasmus Bro
% KVL, Denmark, [e-mail]

% Magic numbers, control TALS loop below:

SMALLNUMBER = eps;

[K,N,P]=size(X);

% initial estimates:

if (nargin < 5)
 disp('Using random initial estimates ...');
 estA = randn(K,M);
 estH = randn(P,M);
 estS = randn(M,N);
else
 estS = estST.';
end

if (nargin < 6)
 MAXNUMITER = 100000;
end

% Trilinear ALS:

% first construct the three different unfolded data matrices:

Xtilde = [];
ap1 = zeros(K,N);
for p=1:1:P,
 for k=1:1:K,
  for n=1:1:N,
   ap1(k,n) = X(k,n,p);
  end
 end
 Xtilde = [Xtilde; ap1];
end

Ytilde = [];
ap2 = zeros(N,P);
for k=1:1:K,
 for n=1:1:N,
  for p=1:1:P,
   ap2(n,p) = X(k,n,p);
  end
 end
 Ytilde = [Ytilde; ap2];
end

Ztilde = [];
ap3 = zeros(P,K);
for n=1:1:N,
 for p=1:1:P,
  for k=1:1:K,
   ap3(p,k) = X(k,n,p);
  end
 end
 Ztilde = [Ztilde; ap3];
end

% compute current fit:
fit = 0;
for p=1:P,
  fit = fit + norm(X(:,:,p) - estA*diag(estH(p,:))*estS,'fro')^2;
end
fprintf('fit = %12.10f\n',fit);
fitold = 2*fit;
fitinit = fit;
it     = 0;


while abs((fit-fitold)/fitold) > SMALLNUMBER & it < MAXNUMITER & fit > 10*eps
 it=it+1;
 fitold=fit;

 % re-estimate A:

 mix = [];
 for n=1:1:N,
  mix = [mix; estH*diag(estS(:,n))];
 end
 if (rank(mix) == M)
  estA = (pinv(mix)*Ztilde).';
 end

 % re-estimate H:

 mix = [];
 for k=1:1:K,
  mix = [mix; estS.'*diag(estA(k,:))];
 end
 if (rank(mix) == M)
  estH = (pinv(mix)*Ytilde).';
 end

 % re-estimate S:

 mix = [];
 for p=1:1:P,
  mix = [mix; estA*diag(estH(p,:))];
 end
 if (rank(mix) == M)
  estS = pinv(mix)*Xtilde;
 end

 % compute new fit:

 fit = 0;
 for p=1:P,
  fit = fit + norm(X(:,:,p) - estA*diag(estH(p,:))*estS,'fro')^2;
 end
 fprintf('fit = %12.10f\n',fit);
 
end % while loop 

estST = estS.';

% end of algorithm