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, nikos@ece.umn.edu % and % Rasmus Bro % KVL, Denmark, rb@kvl.dk % 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
