parafac.m

多维数据分析,有nPLS,PARAFAC,TURKER等

end

% Display convergence criterion
if showfit~=-1
   disp([' The convergence criterion is ',num2str(crit)]) 
end

% Define loading as one ((r1*r2*r3*...*r7)*Fac x 1) vector [A(:);B(:);C(:);...].
% The i'th loading goes from lidx(i,1) to lidx(i,2)
lidx=[1 DimX(1)*Fac];
for i=2:ord
   lidx=[lidx;[lidx(i-1,2)+1 sum(DimX(1:i))*Fac]];
end

% Check if weighted regression required
if size(Weights,1)==size(X,1)&prod(size(Weights))/size(X,1)==size(X,2)
    Weights = reshape(Weights,size(Weights,1),prod(size(Weights))/size(X,1));
   if showfit~=-1
      disp(' Given weights will be used for weighted regression')
   end
   DoWeight=1;
else
   if showfit~=-1
      disp(' No weights given')
   end
   DoWeight=0;
end

% Make idx matrices if missing values
if any(isnan(X(:)))
   MissMeth=1;
else
   MissMeth=0;
end
if MissMeth
   id=sparse(find(isnan(X)));
   idmiss2=sparse(find(~isnan(X)));
   if showfit~=-1
      disp([' ', num2str(100*(length(id)/prod(DimX))),'% missing values']);
      disp(' Expectation maximization will be used for handling missing values')
   end
   SSX=sum(sum(X(idmiss2).^2)); % To be used for evaluating the %var explained
   % If weighting to zero should be used
   % Replace missing with mean values or model estimates initially
   %Chk format ok.
   dimisok = 1;
   if length(OldLoad)==length(DimX)
     for i=1:length(DimX)
       if ~all(size(OldLoad{i})==[DimX(i) Fac])
         dimisok = 0;
       end
     end
   else
     dimisok = 0;
   end
   if dimisok
      model=nmodel(OldLoad);
      model = reshape(model,DimX);
      X(id)=model(id);
   else
      meanX=mean(X(find(~isnan(X))));
      meanX=mean(meanX);
      X(id)=meanX*ones(size(id));
   end
else
   if showfit~=-1
      disp(' No missing values')
   end
   SSX=sum(sum(X.^2)); % To be used for evaluating the %var explained
end

% Check if weighting is tried used together with unimodality or orthogonality
if any(const==3)|any(const==1)
   if DoWeight==1
      disp(' ')
      disp(' Weighting is not possible together with unimodality and orthogonality.')
      disp(' It can be done using majorization, but has not been implemented here')
      disp(' Please contact the authors for further information')
      error
   end
end

% Acceleration
acc=-5;     
do_acc=1;   % Do acceleration every do_acc'th time
acc_pow=2;  % Extrapolate to the iteration^(1/acc_pow) ahead
acc_fail=0; % Indicate how many times acceleration have failed 
max_fail=4; % Increase acc_pow with one after max_fail failure
if showfit~=-1
   disp(' Line-search acceleration scheme initialized')
end

% Find initial guesses for the loadings if no initial values are given

% Use old loadings
if length(OldLoad)==ord % Use old values
   if showfit~=-1
      disp(' Using old values for initialization')
   end
   Factors=OldLoad;
   % Use DTLD
elseif Init==0
   if min(DimX)>1&ord==3&MissMeth==0
       if sum(DimX<Fac)<2
         if showfit~=-1
           disp(' Using direct trilinear decomposition for initialization')
         end
         try 
           [A,B,C]=dtld(reshape(X,DimX),Fac);
         catch
           A = rand(DimX(1),Fac);B = rand(DimX(2),Fac);C = rand(DimX(3),Fac);
         end
       else
         if showfit~=-1
           disp(' Using random values for initialization')
         end
         for i=1:length(DimX)
           Factors{i}=rand(DimX(i),Fac);
         end
         A = Factors{1};B=Factors{2};C = Factors{3};
       end
        A=real(A);B=real(B);C=real(C);
      % Check for signs and reflect if appropriate
      for f=1:Fac
         if sign(sum(A(:,f)))<0
            if sign(sum(B(:,f)))<0
               B(:,f)=-B(:,f);
               A(:,f)=-A(:,f);
            elseif sign(sum(C(:,f)))<0
               C(:,f)=-C(:,f);
               A(:,f)=-A(:,f);
            end
         end
         if sign(sum(B(:,f)))<0
            if sign(sum(C(:,f)))<0
               C(:,f)=-C(:,f);
               B(:,f)=-B(:,f);
            end
         end
      end
      Factors{1}=A;Factors{2}=B;Factors{3}=C;
      
    else
      if showfit~=-1
        disp(' Using singular values for initialization')
      end
      try 
        Factors=ini(reshape(X,DimX),Fac,2);
      catch
        Factors=[];
        for i=1:length(DimX);
          l = rand(DimX(i),Fac);
          Factors = [Factors;l(:)];
        end          
      end
    end

    % Use SVD 
 elseif Init==1
   if all(DimX>=Fac)
     if showfit~=-1
       disp(' Using singular values for initialization')
     end
     try 
       Factors=ini(reshape(X,DimX),Fac,2);
     catch
       Factors=[];
       for i=1:length(DimX);
         l = rand(DimX(i),Fac);
         Factors = [Factors;l(:)];
       end          
     end

   else
     if showfit~=-1
       disp(' Using random values for initialization')
     end
     for i=1:length(DimX)
       Factors{i}=rand(DimX(i),Fac);
     end
   end
   
   % Use random (orthogonal)
elseif Init==2
   if showfit~=-1
      disp(' Using orthogonal random for initialization')
   end
   Factors=ini(reshape(X,DimX),Fac,1);
   
elseif Init==3
   error(' Initialization option set to three has been changed to 10')
   
   % Use several small ones of the above
elseif Init==10
   if showfit~=-1
      disp(' Using several small runs for initialization')
   end
   Opt=Options;
   Opt(5) = NaN;
   Opt(6) = NumbIteraInitia;
   Opt(2) = 0;
   ERR=[];
   [Factors,it,err] = parafac(reshape(X,DimX),Fac,Opt,const,[],[],Weights);
   ERR = [ERR;err];
   Opt(2) = 1;
   [F,it,Err] = parafac(reshape(X,DimX),Fac,Opt,const,[],[],Weights);
   ERR=[ERR;Err];
   if Err<err
      Factors=F;
      err=Err;
   end
   Opt(2)=2;
   for rep=1:3
      [F,it,Err]=parafac(reshape(X,DimX),Fac,Opt,const,[],[],Weights);
      ERR=[ERR;Err];
      if Err<err
         Factors=F;
         err=Err;
      end
   end
   if showfit~=-1
      disp(' ')
      disp(' Obtained fit-values')
      disp([' Method   Fit'])
      disp([' DTLD     ',num2str(ERR(1))])
      disp([' SVD      ',num2str(ERR(2))])
      disp([' RandOrth ',num2str(ERR(3))])
      disp([' RandOrth ',num2str(ERR(4))])
      disp([' RandOrth ',num2str(ERR(5))])
   end
else
   error(' Problem in PARAFAC initialization - Not set correct')
end

% Check for signs and reflect if appropriate

for f=1:Fac
    for m=1:ord-1
        if sign(sum(Factors{m}(:,f)<0)) & FixMode(m)==0
            contin=1;
            for m2 = m+1:ord
                if contin
                    if sign(sum(Factors{m2}(:,f)<0))
                        Factors{m}(:,f)=-Factors{m}(:,f);
                        Factors{m2}(:,f)=-Factors{m2}(:,f);
                        contin=0;
                    end
                end
            end
        end
    end
end
% Convert to old format
if iscell(Factors)
    ff = [];
    for f=1:length(Factors)
        ff=[ff;Factors{f}(:)];
    end
    Factors = ff;
end
% ALTERNATING LEAST SQUARES
err=SSX;
f=2*crit;
it=0;
connew=2;conold=1; % for missing values
ConstraintsNotRight = 0; % Just to ensure that iterations are not stopped if constraints are not yet fully imposed

if showfit~=-1
   disp(' ')
   disp(' Sum-of-Squares   Iterations  Explained')
   disp(' of residuals                 variation')
end

while (((f>crit) | (norm(connew-conold)/norm(conold)>MissConvCrit) | ConstraintsNotRight) & it<maxit)|~ nonneg_obeyed
   conold=connew; % for missing values
   it=it+1;
   acc=acc+1; 
   if acc==do_acc;
      Load_o1=Factors;
   end
   if acc==do_acc+1;
      acc=0;Load_o2=Factors;
      Factors=Load_o1+(Load_o2-Load_o1)*(it^(1/acc_pow));
      % Convert to new format
      clear ff,id1 = 0;
      for i = 1:length(DimX) 
         id2 = sum(DimX(1:i).*Fac);ff{i} = reshape(Factors(id1+1:id2),DimX(i),Fac);id1 = id2;
      end
      model=nmodel(ff);
      model = reshape(model,DimX(1),prod(DimX(2:end)));
      
      if MissMeth
         connew=model(id);
         errX=X-model;
         if DoWeight==0
            nerr=sum(sum(errX(idmiss2).^2));
         else
            nerr=sum(sum((Weights(idmiss2).*errX(idmiss2)).^2));
         end
      else
         if DoWeight==0
            nerr=sum(sum((X-model).^2));
         else
            nerr=sum(sum((X.*Weights-model.*Weights).^2));
         end
      end
      if nerr>err
         acc_fail=acc_fail+1;
         Factors=Load_o2;
         if acc_fail==max_fail,
            acc_pow=acc_pow+1+1;
            acc_fail=0;
            if showfit~=-1
               disp(' Reducing acceleration');
            end
         end
      else
         if MissMeth
            X(id)=model(id);
         end
      end
   end
   
   if DoWeight==0
      for ii=ord:-1:1
         if ii==ord;
            i=1;
         else
            i=ii+1;
         end
         idd=[i+1:ord 1:i-1];
         l_idx2=lidx(idd,:);
         dimx=DimX(idd);
         if ~FixMode(i)
            L1=reshape(Factors(l_idx2(1,1):l_idx2(1,2)),dimx(1),Fac);
            if ord>2
               L2=reshape(Factors(l_idx2(2,1):l_idx2(2,2)),dimx(2),Fac);
               Z=kr(L2,L1);
            else
               Z = L1;
            end
            for j=3:ord-1
               L1=reshape(Factors(l_idx2(j,1):l_idx2(j,2)),dimx(j),Fac);
               Z=kr(L1,Z);
            end
            ZtZ=Z'*Z;
            ZtX=Z'*X';
            OldLoad=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
            L=pfls(ZtZ,ZtX,DimX(i),const(i),OldLoad,DoWeight,Weights);
            Factors(lidx(i,1):lidx(i,2))=L(:);
         end
         x=zeros(prod(DimX([1:ii-1 ii+1:ord])),DimX(ii));  % Rotate X so the current last mode is the first
         x(:)=X;
         X=x';
      end

   else