parafac.m

强大的多维工具箱.应用在Matlab中,可分析多纬数据结构.直接安装.

  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{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 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
    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
        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(:);
  end
  % APPLY SIGN CONVENTION IF NO FIXED MODES
  %  FixMode=1
  if ~any(FixMode)&~(any(const==2)|any(const==3))
    Sign = ones(1,Fac);
    for i=ord:-1:2
      A=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
      Sign2=ones(1,Fac);
      for ff=1:Fac
        [out,sig]=max(abs(A(:,ff)));
        Sign(ff) = Sign(ff)*sign(A(sig,ff));
        Sign2(ff) = sign(A(sig,ff));
      end
      A=A*diag(Sign2);
      Factors(lidx(i,1):lidx(i,2))=A(:);
    end
    A=reshape(Factors(lidx(1,1):lidx(1,2)),DimX(1),Fac);
    A=A*diag(Sign);
    Factors(lidx(1,1):lidx(1,2))=A(:);
  end

  % Check if nonneg_obeyed
  for i=1:ord
    if const(i)==2|const(i)==3
      A=reshape(Factors(lidx(i,1):lidx(i,2)),DimX(i),Fac);
      if any(A(:))<0
        nonneg_obeyed=0;
      end
    end
  end

  % EVALUATE SOFAR
  % Convert to new format
  clear ff,id1 = 0;
  for i = 1:length(DimX)