function [A,B,C,D]=QALS(X,F,A,B,C,D); % Plain Vanilla Quadrilinear ALS % X = IxJxKxL rank-F 4-way array % Nikos Sidiropoulos % UVA, Sept. 1999 % uses krp.m, which implements the Khatri-Rao product % Magic numbers, control ALS loop below: SMALLNUMBER = 10^8*eps; MAXNUMITER = 100; [I, J, K, L]=size(X); % initial estimates: if (nargin < 6) % disp('Using random initial estimates ...'); A = randn(I,F); B = randn(J,F); C = randn(K,F); D = randn(L,F); end % Quadrilinear ALS: % first construct the four different unfolded data matrices: % Used in the A-update step: UA = []; for j=1:J, for k=1:K, UA = [UA; squeeze(X(:,j,k,:)).']; end end % Used in the B-update step: UB = []; for k=1:K, for l=1:L, UB = [UB; squeeze(X(:,:,k,l))]; end end % Used in the C-update step: UC = []; for l=1:L, for i=1:I, UC = [UC; squeeze(X(i,:,:,l))]; end end % Used in the D-update step: UD = []; for i=1:I, for j=1:J, UD = [UD; squeeze(X(i,j,:,:))]; end end % compute current fit: fit = 0; for k=1:K, for l=1:L, model(:,:,k,l) = A*diag(C(k,:))*diag(D(l,:))*B.'; fit = fit + norm(squeeze(X(:,:,k,l))-squeeze(model(:,:,k,l)),'fro')^2; end end %fprintf('fit = %12.10f\n',fit); fitold = 2*fit; fitinit = fit; it = 0; allfits = []; while abs((fit-fitold)/fitold) > SMALLNUMBER & it < MAXNUMITER & fit > 10^5*eps it=it+1; fitold=fit; % re-estimate A: A = (pinv(krp(B,krp(C,D)))*UA).'; % re-estimate B: B = (pinv(krp(C,krp(D,A)))*UB).'; % re-estimate C: C = (pinv(krp(D,krp(A,B)))*UC).'; % re-estimate D: D = (pinv(krp(A,krp(B,C)))*UD).'; % compute new fit: fit = 0; for k=1:K, for l=1:L, model(:,:,k,l) = A*diag(C(k,:))*diag(D(l,:))*B.'; fit = fit + norm(squeeze(X(:,:,k,l))-squeeze(model(:,:,k,l)),'fro')^2; end end allfits = [allfits; fit]; plot(allfits); drawnow; fprintf('fit = %12.10f\n',fit); end % while loop % end of algorithm function AkrpB = krp(A,B); [I F] = size(A); [J F1] = size(B); if (F1 ~= F) disp('krp.m: column dimensions do not match!!! - exiting matlab'); exit; end AkrpB = []; for f=1:F, AkrpB = [AkrpB kron(A(:,f),B(:,f))]; end