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