parafac.m

多元统计分析是一种应用非常广泛的数据处理方法

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{i} =l;
      end
      if showfit~=-1
        disp(' Oops sorry - ended up with random instead')
      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 & FixMode(m2)==0
          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 strcmp(accel_pattern,'paired') % this chooses Bro standard two-steps version
        if it==1
            startiter=-acc+1;
            disp([' LS paired accel will be done every other time after iter ' ...
                int2str(startiter)])
        end

        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)));
            %G(it)=err;

            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
    elseif strcmp(accel_pattern,'shifted') % this chooses harshman's every-iter way
        if acc==do_acc-1;
            Load_o2=Factors;
            disp([' LS shifted accel will be done every time after iter ' int2str(it)])
        elseif acc==do_acc;
            Load_o1=Load_o2; % shift the one that was latest back to first
            acc=0;Load_o2=Factors; % latest factors become the newer set
            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)));
            %G(it)=err;

            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
    elseif strcmp(accel_pattern,'none')
        % do nothing
    else
        error('accel_pattern should be ''shifted'', ''paired'', or ''none''.')
    end


%   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=krb(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=krb(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
    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=krb(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=krb(L1,Z);
        end
        OldLoad=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
        L=pfls(Z,X,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));
      x(:)=X;
      X=x';
      x(:)=Weights;
      Weights=x';
    end
  end
  
  % POSTPROCES LOADINGS (ALL VARIANCE IN FIRST MODE)
  if ~any(FixMode)

    A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
    for i=2:ord
      B=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
      for ff=1:Fac
        A(:,ff)=A(:,ff)*norm(B(:,ff));
        B(:,ff)=B(:,ff)/norm(B(:,ff));
      end
      Factors(lidx(i,1):lidx(i,2))=B(:);
    end
    Factors(lidx(1,1):lidx(1,2))=A(:);