tucker.m

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

function [Factors,G,ExplX,Xm]=tucker(X,DimX,Fac,Options,ConstrF,ConstrG,Factors,G);
%function [Factors,G,ExplX,Xm]=tucker(X,DimX,Fac[,Options[,ConstrF,[ConstrG[,Factors[,G]]]]]);
%
% $ Version 1.12 $ Date 14. Nov. 1999 $ Not compiled $
%
% Change: True LS unimodality now supported.
%
% This algorithm requires access to:
% 'fnipals' 'gsm' 'inituck' 'calcore'
% 'nmodel' 'nonneg' 'setopts' 'misssum'
% 'missmean' 't3core'
%
% See also:
% 'parafac' 'maxvar3' 'maxdia3' 'maxswd3'
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. Furthermore, the
% code can not be made part of anything but the 'N-way Toolbox'.
% In case of doubt, contact the holder of the copyrights.
%
% Claus A. Andersson
% Chemometrics Group, Food Technology
% Department of Food and Dairy Science
% Royal Veterinary and Agricultutal University
% Rolighedsvej 30, DK-1958 Frederiksberg, Denmark
% E-mail: [e-mail]
%
% ---------------------------------------------------------           
%             The general N-way Tucker model
% ---------------------------------------------------------
%    
% [Factors,G,ExplX,Xm]=tucker(X,DimX,Fac[,Options[,ConstrF,[ConstrG[,Factors[,G]]]]]);
% [Factors,G,ExplX,Xm]=tucker(X,DimX,Fac);
%
% X        : The unfolded data.
% DimX     : Row-vector describing the dimensions of X.
% Fac      : Row-vector describing the number of factors
%            in each of the N modes. A '-1' (minus one)
%            will tell the algorithm not to estimate factors
%            for this mode, yielding a Tucker2 model.
%            Ex. [3 2 4]
% Options  : See parafac.
% ConstrF  : Constraints that must apply to 'Factors'.
%            Define a row-vector of size N that describes how
%            each mode should be treated.
%            '0' orthogonality (default)
%            '1' non-negativity
%            '2' unconstrained
%            '3' keep unchanged from users input 'Factors'.
%            '4' unimodality and non-negativity.
%            E.g.: [0 2 1] yields ortho in first mode, uncon in the second
%            and non-neg in the third mode.
%            Note: The algorithm uses random values if there are no
%            non-negative components in the iteration intermediates. Thus,
%            if non-negativity is applied, the iterations may be
%            non-monotone in minor sequences.
% ConstrG  : Constraints that must apply to 'G'.
%            '[]' or '0' will not constrain the elements of 'G'.
%            To define what core elements should be allowed, give a core that
%            is 1 (one) on all active positions and zero elsewhere - this boolean
%            core array must have the same dimensions as defined by 'Fac'.
%            If 'ConstrG'==3 the core is kept constant as defined by input. If the 
%            core is fixed, the intermediate outputs will no longer correspond to the
%            sum of squared core entries (would be constants)
%            they will reflect some iteration intermediates.
% Factors  : A row-vector containing the solutions.
% G        : Core array that matches the dimensions defined by 'Fac'.
% ExplX    :
% Xf       : Models predictions of the missing elements
%
% This algorithm applies to the general N-way case, so
% the unfolded X can have any number of dimensions. The
% principles of 'projections' and 'systematic unfolding 
% methodology (SUM)' are used in this algorithm to provide
% a fast approach - also for larger data arrays. This
% algorithm can handle missing values if denoted
% by NaN's. It can also be used to make TUCKER2 models by
% properly setting the elements of 'Fac' to -1.
%
% Note: When estimating a Tucker model on data using non-orthogonal factors,
%       the sum of square of the core may differ between models of the
%       same dataset. This is in order since the factors may
%       thus be correlated. However, the expl. var. should always be the same.
%      

format long
format compact
dbg=0;

if nargin==0,
   help('tucker.m');
   error(['Error calling ''tucker.m''. Since no input arguments were given, the ''help'' command was initiated.'])
   return;
end;
if nargin<3,
   help('tucker.m');
   error(['Error calling ''tucker.m''. At least three (3) input arguments must be given. Read the text above.'])
   return;
end;
if size(Fac,2)==1,
   help('tucker.m');
   error(['Error calling ''tucker.m''. ''Fac'' must be a row-vector.'])
end;    

% Initialize system variables
N=size(Fac,2);
Fac_orig=Fac;
finda=find(Fac==-1);
if ~isempty(finda),
   Fac(finda)=zeros(size(finda));
end;
FIdx0=cumsum([1 DimX(1:N-1).*Fac(1:N-1)]);
FIdx1=cumsum([DimX.*Fac]);
pmore=30;
pout=0;
Xm=[];
MissingExist=any(isnan(X(:)));
if MissingExist,
   IdxIsNans=find(isnan(X));
end;
SSX=misssum(misssum(X.^2));

if exist('Options'),
   Options_=Options;
else
   Options_=[0];
end;
load noptiot3.mat;
i=find(Options_);
Options(i)=Options_(i);
if isnan(Options(5)),
   prlvl = 0;
else 
   prlvl = 1;
end;
Options12=Options(1);
Options11=Options12*10;
Options21=Options(2);
Options31=Options(3);
Options41=Options(4);
Options51=Options(5);
Options61=Options(6);
Options71=Options(7);
Options81=Options(8);
Options91=Options(9);
Options101=Options(10);

if ~exist('ConstrF'),
   ConstrF=[];
end;
if isempty(ConstrF),
   ConstrF=zeros(size(DimX));
end;
if ConstrF==0 ,
   ConstrF=zeros(size(DimX));
end;

if ~exist('ConstrG')
   ConstrG=[];
end;
if isempty(ConstrG),
   ConstrG=0;
end;

if ~exist('Factors'),
   Factors=[];
end;

if ~exist('G'),
   G=[];
end;

%Give a status/overview
if prlvl>0,
   fprintf('\n\n');
   fprintf('=================   RESUME  &  PARAMETERS   ===================\n');
   fprintf('Array                 : %i-way array with dimensions (%s)\n',N,int2str(DimX));
   if any(Fac==0),
      fprintf('Model                 : (%s) TUCKER2 model\n',int2str(Fac));
   else
      fprintf('Model                 : (%s) TUCKER3 model\n',int2str(Fac));
   end;   
end

%Mth initialization
txt1=str2mat('derived by SVD (orthogonality constrained).');
txt1=str2mat(txt1,'derived by NIPALS (orthogonality constrained).');
txt1=str2mat(txt1,'derived by Gram-Schmidt (orthogonality constrained).');
txt1=str2mat(txt1,'This mode is not compressed/calculated, i.e., TUCKER2 model.');
txt1=str2mat(txt1,'derived by non-negativity least squares.');
txt1=str2mat(txt1,'derived by unconstrained simple least squares.');
txt1=str2mat(txt1,'unchanged, left as defined in input ''Factors''.');
txt1=str2mat(txt1,'derived by unimodality constrained regression.');
MethodO=1;
for k=1:N,
   UpdateCore(k)=1;
   if ConstrF(k)==0,
      if Fac(k)>0,
         if 0<DimX(k) & DimX(k)<=180,
            Mth(k)=1;
         end;
         if 180<DimX(k) & DimX(k)<=Inf,
            Mth(k)=3;
         end;
         if Fac(k)<=6 & 180<DimX(k),
            Mth(k)=2;
         end;
      end;
      UpdateWithPinv(k)=1; %Update with the L-LS-P-w/Kron approach
      CalcOrdinar(k)=1;
   end;
   if ConstrF(k)==1,
      Mth(k)=5; %nonneg
      MethodO=2; %use the flexible scheme
      UpdateCore(k)=1; %Update the core in this mode
      CalcOrdinar(k)=1;
   end;
   if ConstrF(k)==2,
      Mth(k)=6; %uncon
      MethodO=2; %use the flexible scheme
      UpdateCore(k)=1; %Update the core in this mode
      CalcOrdinar(k)=1;
   end;
   if ConstrF(k)==3,
      Mth(k)=7; %unchanged
      MethodO=2;
      UpdateCore(k)=1; %Update the core in this mode
      CalcOrdinar(k)=1;
   end;
   if ConstrF(k)==4,
      Mth(k)=8; %unimod
      MethodO=2; %use the flexible scheme
      UpdateCore(k)=1; %Update the core in this mode
      CalcOrdinar(k)=1;
   end;   
   if Fac_orig(k)==-1
      Mth(k)=4;
      UpdateCore(k)=0; %Do not update core for this mode
      CalcOrdinar(k)=1;
   end;
   if Options91>=1,
      if prlvl>0,
         if Mth(k)~=4,
            fprintf('Mode %i                : %i factors %s\n',k,Fac(k),txt1(Mth(k),:));
         else
            fprintf('Mode %i                : %s\n',k,txt1(Mth(k),:));
         end;
      end;
   end;
end;

UserFactors=1;
if isempty(Factors),
   UserFactors=0;
end;

usefacinput=0;
if MissingExist,
   if ~UserFactors
      [i j]=find(isnan(X));
      mnx=missmean(X)/3;
      mny=missmean(X')/3;
      n=size(i,1);
      for k=1:n,
         i_=i(k);
         j_=j(k);
         X(i_,j_) = mny(i_) + mnx(j_);
      end;
      mnz=(missmean(mnx)+missmean(mny))/2;
      p=find(isnan(X));
      X(p)=mnz;
   else
      usefacinput=1;
      Xm=nmodel(Factors,G,DimX,Fac_orig);
      X(IdxIsNans)=Xm(IdxIsNans);
   end;
   SSMisOld=sum(sum( X(IdxIsNans).^2 ));
   SSMis=SSMisOld;
end;

% Initialize the Factors by some method
UserFactors=1;
if isempty(Factors),
   Factors=inituck(X,DimX,Fac_orig,2,[]);
   UserFactors=0;
end;

% Initialize the core
Core_uncon=0;
Core_nonneg=0;
Core_cmplex=0;
Core_const=0;
G_cons=ConstrG;
if all(ConstrG(:)==1),
   ConstrG=0;
end;
if ConstrG==0,
   G=calcore(X,DimX,Fac_orig,Factors,[],1,MissingExist);   
   Core_uncon=1;
elseif prod(size(ConstrG)==[Fac(1) prod(Fac(2:N))]),
   tmpM2=1;
   for k=1:N;
      if Mth(k)==4,
         tmpM1=eye(DimX(k));
      else
         tmpM1=reshape(Factors(FIdx0(k):FIdx1(k)),DimX(k),Fac(k));
      end;
      tmpM2=ckron(tmpM2,tmpM1);
   end
   w=ConstrG(:);
   fwz=find(w==1);
   fwnn=find(w==2);
   fwfix=find(w==3);
   G=zeros(prod(Fac),1);
   G(fwz)=tmpM2(:,fwz)\X(:);
   enda=size(Fac,2); %!
   G=reshape(G,Fac(1),prod(Fac(2:enda)));
   Core_cmplex=1;
   MethodO=2;
   UpdateCore=2*ones(1,N);
elseif ConstrG==2,
   G=calcore(X,DimX,Fac_orig,Factors,[],1,MissingExist);   
   G(find(G<0))=0;
   Core_nonneg=1;
elseif ConstrG==3,
   Core_const=1;
end;   

if prlvl>0
   fprintf('Type of algorithm     : ');
   if MethodO==1,
      fprintf('Orthogonal projections.\n');
   elseif MethodO==2,
      fprintf('Flexible scheme.\n');
   end;
   fprintf('Core                  : ');
   if Core_cmplex==1 & Core_nonneg==0,
      fprintf('Unconstrained composite/sparse core.\n');
   elseif Core_cmplex==1 & Core_nonneg==1,
      fprintf('Non-negativity constrained and composite/sparse core.\n');
   elseif Core_const==1,
      fprintf('Fixed core.\n');
   else
      fprintf('Full unconstrained core.\n');
   end;
end;

if prlvl>0,
   if MissingExist,
      if Options91>=1,
         fprintf('Missing data          : Yes, 2 active loops (expectation maximization).\n');
         fprintf('                        %i values (%.2f%%) out of %i are NaNs/missing.\n',prod(size(IdxIsNans)),100*prod(size(IdxIsNans))/prod(size(X)),prod(size(X)));
         if usefacinput==0,
            fprintf('                        Missing values initialized from column and row means.\n');
         else
            fprintf('                        Missing values initialized from model based on the given input.\n');
         end;
         fprintf('Convergence crit. 1   : %.5g (relative) sum of sq. core elements (corrected for missing values).\n',Options12);
         fprintf('Convergence crit. 2   : %.5g (relative) sum of sq. of pred. missing values.\n',Options11);
         fprintf('Iteration limit       : %i is the maximum number of overall iterations.\n',Options61);
      end;
   else
      if Options91>=1,
         fprintf('Missing data          : No, 1 active loop.\n');
         fprintf('Convergence crit. 1   : %.5g (relative) sum of sq. core elements.\n',Options12);
         fprintf('Iteration limit       : %i is the maximum number of overall iterations.\n',Options61);
      end;
   end;
   fprintf('\n');
   if MissingExist,
      str1=' Iter. 1  |  Corrected sum of   |  Sum of sq. miss.  |   Expl.  ';
      str2='     #    |  sq. core elements  |       values       |  var. [%]';
      fprintf('%s\n',str1);
      fprintf('%s\n\n',str2);
   else
      str1=' Iter. 1  |       Sum of        |   Expl.  ';
      str2='     #    |  sq. core elements  |  var. [%]';
      fprintf('%s\n',str1);
      fprintf('%s\n\n',str2);
   end;
end;
Conv_true=0;
if MethodO==1, %Can use the faster projection technique
   SSGOld=0;
   Converged2=0;
   it1=0;
   itlim1=0;
   t0=clock;
   while ~Converged2, 
      Converged1=0;
      while ~Converged1,
         it1=it1+1;   
         % Iterate over the modes
         for c=1:N,
            %Compress the data by projections
            if Mth(c)~=4,
               CurDimX=DimX;
               RedData=X;
               for k=1:N;
                  if k~=c,
                     if Mth(k)~=4,
                        kthFactor=reshape(Factors(FIdx0(k):FIdx1(k)),DimX(k),Fac(k));
                        RedData=kthFactor'*RedData;
                        CurDimX(k)=Fac(k);
                     else
                        RedData=RedData;
                     end,
                  end,
                  if k~=N,
                     newi=CurDimX(k+1);
                     newj=prod(CurDimX)/newi;
                  else
                     newi=CurDimX(1);
                     newj=prod(CurDimX)/newi;
                  end;
                  RedData=reshape(RedData',newi,newj);
               end;
               %Reshape to the proper unfolding
               for k=1:(c-1);
                  if k~=c,
                     newi=CurDimX(k+1);
                     newj=prod(CurDimX)/CurDimX(k+1);
                  else,
                     newi=CurDimX(1);
                     newj=prod(CurDimX)/CurDimX(1);
                  end;
                  RedData=reshape(RedData',newi,newj);
               end;
               %Find a basis in the projected space
               %...using the robust SVD
               if Mth(c)==1,
                  if MissingExist,
                     [U S V]=svd(RedData',0);
                     cthFactor=V(:,1:Fac(c));
                  else
                     [U S V]=svd(RedData',0);
                     cthFactor=V(:,1:Fac(c));