parafac2.m

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

   if any(miss(:))
     MissingElements = 1;
     % Replace missing with mean over slab (not optimal but what the heck)
     % Iteratively they'll be replaced with model estimates
     x(find(miss)) = mean(x(find(~miss)));
     X{k} = x;
     MissNum = MissNum + prod(size(find(miss)));
     AllNum = AllNum + prod(size(x));
   end
end
if MissingElements
   if Options(5)==0
      PercMiss = 100*MissNum/AllNum;
      RoundedOf = .1*round(PercMiss*10);
      disp([' Missing data handled by EM   : ',num2str(RoundedOf),'%'])
   end
end
clear x

% Initialize by ten small runs
if nargin<5
   if initi==0
      if Options(5)==0
         disp([' Use best of ',num2str(NumRep)]) 
         disp(' initially fitted models')
      end
      Opt = Options;
      Opt = Options(1)/20;
      Opt(2) = NumItInRep; % Max NumItInRep iterations
      Opt(3) = 1;  % Init with SVD
      Opt(4) = 0;
      Opt(5) = 1;
      [A,H,C,P,bestfit]=parafac2(X,F,Constraints,Opt);
      AllFit = bestfit;
      for i = 2:NumRep
         Opt(3) = 2;   % Init with random
         [a,h,c,p,fit]=parafac2(X,F,Constraints,Opt);
         AllFit = [AllFit fit];
         if fit<bestfit
            A=a;H=h;C=c;P=p;
            bestfit = fit;
         end
      end
      AddiOutput.AllFit = AllFit;
      if Options(5)==0
         for ii=1:length(AllFit)
            disp([' Initial Model Fit            : ',num2str(AllFit(ii))])
         end
      end
      % Initialize by SVD
   elseif initi==1
      if Options(5)==0
         disp(' SVD based initialization')
      end
      XtX=X{1}*X{1}';
      for k = 2:K
         XtX = XtX + X{k}*X{k}';
      end
      [A,s,v]=svd(XtX,0);  
      A=A(:,1:F);
      C=ones(K,F)+randn(K,F)/10;
      H = eye(F);
   elseif initi==2
      if Options(5)==0
         disp(' Random initialization')
      end
      A = rand(I,F);
      C = rand(K,F);
      H = eye(F);
   else
      error(' Options(2) wrongly specified')
   end
end

if initi~=1
   XtX=X{1}*X{1}'; % Calculate for evaluating fit (but if initi = 1 it has been calculated)
   for k = 2:K
      XtX = XtX + X{k}*X{k}';
   end
end  
fit    = sum(diag(XtX));
oldfit = fit*2;
fit0   = fit;
it     = 0;
Delta = 1;

if Options(5)==0
   disp(' ')
   disp(' Fitting model ...')
   disp(' Loss-value      Iteration     %VariationExpl')
end

% Iterative part
while abs(fit-oldfit)>oldfit*ConvCrit & it<MaxIt & fit>1000*eps
    oldfit = fit;
    it   = it + 1;
    
    % Update P
    for k = 1:K
      Qk       = X{k}'*(A*diag(C(k,:))*H');
      P{k}     = Qk*psqrt(Qk'*Qk);
    %  [u,s,v]  = svd(Qk.');P{k}  = v(:,1:F)*u(:,1:F)';
      Y(:,:,k) = X{k}*P{k};
    end
    
    % Update A,H,C using PARAFAC-ALS
    [A,H,C,ff]=parafac(reshape(Y,I,F*K),[I F K],F,1e-4,[Constraints(1) ConstB Constraints(2)],A,H,C,5);
    [fit,X] = pf2fit(X,A,H,C,P,K,MissingElements,MissingOnes);
      
    % Print interim result
    if rem(it,ShowFit)==0|it == 1
       if Options(5)==0
          fprintf(' %12.10f       %g        %3.4f \n',fit,it,100*(1-fit/fit0));
          subplot(2,2,1)
          plot(A),title('First mode')
          subplot(2,2,2)
          plot(C),title('Third mode')
          subplot(2,2,3)
          plot(P{1}*H),title('Second mode (only first k-slab shown)')
          drawnow
       end
    end

end

if rem(it,ShowFit)~=0 %Show final fit if not just shown
   if Options(5)==0
      fprintf(' %12.10f       %g        %3.4f \n',fit,it,100*(1-fit/fit0));
   end
end



function [fit,X]=pf2fit(X,A,H,C,P,K,MissingElements,MissingOnes);

   % Calculate fit and impute missing elements from model

   fit = 0;
   for k = 1:K
     M   = A*diag(C(k,:))*(P{k}*H)';
     % if missing values replace missing elements with model estimates
     if nargout == 2 
       if any(MissingOnes{k})
         x=X{k};
         x(find(MissingOnes{k})) = M(find(MissingOnes{k}));
         X{k} = x;
       end
     end
     fit = fit + sum(sum(abs (X{k} - M ).^2));
   end


function X = psqrt(A,tol)

   % Produces A^(-.5) even if rank-problems

   [U,S,V] = svd(A,0);
   if min(size(S)) == 1
     S = S(1);
   else
     S = diag(S);
   end
   if (nargin == 1)
     tol = max(size(A)) * S(1) * eps;
   end
   r = sum(S > tol);
   if (r == 0)
     X = zeros(size(A'));
   else
     S = diag(ones(r,1)./sqrt(S(1:r)));
     X = V(:,1:r)*S*U(:,1:r)';
  end
  
  
function [A,B,C,fit,it] = parafac(X,DimX,Fac,crit,Constraints,A,B,C,maxit,DoLineSearch);

% Complex PARAFAC-ALS
% Fits the PARAFAC model Xk = A*Dk*B.' + E
% where Dk is a diagonal matrix holding the k'th
% row of C.
%
% Uses on-the-fly projection-compression to speed up 
% the computations. This requires that the first mode 
% is the largest to be effective
% 
% INPUT
% X          : Data
% DimX       : Dimension of X
% Fac        : Number of factors
% OPTIONAL INPUT
% crit       : Convergence criterion (default 1e-6)
% Constraints: [a b c], if e.g. a=0 => A unconstrained, a=1 => A nonnegative
% A,B,C      : Initial parameter values
%
% I/O
% [A,B,C,fit,it]=parafac(X,DimX,Fac,crit,A,B,C);
%
% Copyright 1998
% Rasmus Bro
% KVL, Denmark, [e-mail]

% Initialization
if nargin<9
  maxit   = 2500;      % Maximal number of iterations
end
showfit = pi;         % Show fit every 'showfit'th iteration (set to pi to avoid)

if nargin<4
  crit=1e-6;
end

if crit==0
  crit=1e-6;
end

I = DimX(1);
J = DimX(2);
K = DimX(3);

InitWithRandom=0;
if nargin<8
   InitWithRandom=1;
end
if nargin>7 & size(A,1)~=I
  InitWithRandom=1;
end

if nargin<5
   ConstA = 0;ConstB = 0;ConstC = 0;
else
   ConstA = Constraints(1);ConstB = Constraints(2);ConstC = Constraints(3);
end

if InitWithRandom

  if I<Fac
    A = rand(I,Fac);
  else
    A = orth(rand(I,Fac));
  end
  if J<Fac
    B = rand(J,Fac);
  else
    B = orth(rand(J,Fac));
  end
  if K<Fac
    C = rand(K,Fac);
  else
    C = orth(rand(K,Fac));
  end
end

SumSqX = sum(sum(abs(X).^2));
fit    = SumSqX;
fit0   = fit;
fitold = 2*fit;
it     = 0;
Delta  = 5;

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

   % Do line-search
   if rem(it+2,2)==-1
      [A,B,C,Delta]=linesrch(X,DimX,A,B,C,Ao,Bo,Co,Delta);
   end
   
   Ao=A;Bo=B;Co=C;
   % Update A
   Xbc=0;
   for k=1:K
     Xbc = Xbc + X(:,(k-1)*J+1:k*J)*conj(B*diag(C(k,:)));
   end
   if ConstA == 0 % Unconstrained
      A = Xbc*pinv((B'*B).*(C'*C)).';
   elseif ConstA == 1 % Nonnegativity, requires reals
      Aold = A;
      for i = 1:I
         ztz = (B'*B).*(C'*C);
         A(i,:) = fastnnls(ztz,Xbc(i,:)')';
      end
      if any(sum(A)<100*eps*I)
         A = .99*Aold+.01*A; % To prevent a matrix with zero columns
      end
   elseif ConstA == 2 % Orthogonality
      A = Xbc*(Xbc'*Xbc)^(-.5);
   elseif ConstA == 3 % Unimodality
      A = unimodalcrossproducts((B'*B).*(C'*C),Xbc',A);
   end

   % Project X down on orth(A) - saves time if first mode is large
   [Qa,Ra]=qr(A,0);
   x=Qa'*X;

   % Update B
   if ConstB == 10 % Procrustes
      B = eye(Fac);
   else
      Xac=0;
      for k=1:K
         Xac = Xac + x(:,(k-1)*J+1:k*J).'*conj(Ra*diag(C(k,:)));
      end
      if ConstB == 0 % Unconstrained
         B = Xac*pinv((Ra'*Ra).*(C'*C)).';
      elseif ConstB == 1 % Nonnegativity, requires reals
         Bold = B;
         for j = 1:J
            ztz = (Ra'*Ra).*(C'*C);
            B(j,:) = fastnnls(ztz,Xac(j,:)')';
         end
         if any(sum(B)<100*eps*J)
            B = .99*Bold+.01*B; % To prevent a matrix with zero columns
         end
      end
   end
  
    % Update C
    if ConstC == 0 % Unconstrained
       ab=pinv((Ra'*Ra).*(B'*B));
       for k=1:K 
          C(k,:) = (ab*diag(Ra'* x(:,(k-1)*J+1:k*J)*conj(B))).';
       end
    elseif ConstC == 1  % Nonnegativity, requires reals
       Cold = C;
       ztz = (Ra'*Ra).*(B'*B);
       for k = 1:K
          xab = diag(Ra'* x(:,(k-1)*J+1:k*J)*B);
          C(k,:) = fastnnls(ztz,xab)';
       end
       if any(sum(C)<100*eps*K)
          C = .99*Cold+.01*C; % To prevent a matrix with zero columns
       end
    elseif ConstC == 2 % Orthogonality
       Z=(Ra'*Ra).*(B'*B);
       Y=[];
       for k=1:K
          d=diag(Ra'*x(:,(k-1)*J+1:k*J)*B)'; 
          Y=[Y;d];
       end;
       [P,D,Q]=svd(Y,0);
       C=P*Q';
    elseif ConstC == 3 % Unimodality
       xab = [];
       for k = 1:K
          xab = [xab diag(Ra'* x(:,(k-1)*J+1:k*J)*B)];
       end
       C = unimodalcrossproducts((Ra'*Ra).*(B'*B),xab,C);
    elseif ConstC == 10 % GPA => Isotropic scaling factor
       ab=(Ra'*Ra).*(B'*B);
       ab = pinv(ab(:));
       C(1,:) = 1;
       for k=2:K 
          yy = [];
          yyy = diag(Ra'* x(:,(k-1)*J+1:k*J)*conj(B)).';
          for f=1:Fac
             yy = [yy;yyy(:)];
          end
          C(k,:) = ab*yy;
       end
    end
      
    % Calculating fit. Using orthogonalization instead
   %fit=0;for k=1:K,residual=X(:,(k-1)*J+1:k*J)-A*diag(C(k,:))*B.';fit=fit+sum(sum((abs(residual).^2)));end
   [Qb,Rb]=qr(B,0);
   [Z,Rc]=qr(C,0);
   fit=SumSqX-sum(sum(abs(Ra*ppp(Rb,Rc).').^2));
   
   if rem(it,showfit)==0
      fprintf(' %12.10f       %g        %3.4f \n',fit,it,100*(1-fit/fit0));
   end
end

% ORDER ACCORDING TO VARIANCE
Tuck     = diag((A'*A).*(B'*B).*(C'*C));
[out,ID] = sort(Tuck);
A        = A(:,ID);
if ConstB ~= 10 % Else B is eye
   B        = B(:,ID);
end
C        = C(:,ID);
% NORMALIZE A AND C (variance in B)
if ConstB ~= 10 % Then B is eye
   for f=1:Fac,normC(f) = norm(C(:,f));end
   for f=1:Fac,normA(f) = norm(A(:,f));end
   B        = B*diag(normC)*diag(normA);  
   A        = A*diag(normA.^(-1));
   C        = C*diag(normC.^(-1));
   
   % APPLY SIGN CONVENTION
   SignA = sign(sum(sign(A))+eps);
   SignC = sign(sum(sign(C))+eps);
   A = A*diag(SignA);
   C = C*diag(SignC);
   B = B*diag(SignA)*diag(SignC);
end

function [NewA,NewB,NewC,DeltaMin] = linesrch(X,DimX,A,B,C,Ao,Bo,Co,Delta);

dbg=0;

if nargin<5