function [A,B,C,FIT,IT]=comfac(X,Fac,Options,DoComp,CompIt,Init);
%COMFAC Algorithm for fitting the complex-valued PARAFAC model
%
% See e.g. Rasmus Bro, N. D. Sidiropoulos, and G. B. Giannakis. A Fast
% Least Squares Algorithm for Separating Trilinear Mixtures.
% Proc. ICA99 � Int. Workshop on Independent Component Analysis
% and Blind Signal Separation. Jan. 11�15, Aussois, France:289-294, 1999.
%
% N. D. Sidiropoulos, G. B. Giannakis, and Rasmus Bro. Blind PARAFAC
% Receivers for DS-CDMA Systems. IEEE Transactions on Signal Processing March, 2000.
%
% The algorithm works by first compressing the data using a Tucker3 models. Subsequently the
% PARAFAC model is fitted to the compressed array, either initialized with DTLD (~ESPRIT) or
% with PARAFAC-ALS. The solution is de-compressed to the original space, and a few safe-guard
% PARAFAC-ALS steps are performed. Further optimization of this algorithm is possible, e.g.
% if the data are very large, if the profiles are very correlated or if the noise is huge. This
% has not been pursued here, in order to have a generally applicable algorithm
%
% INPUT
% X : I x J x K data three-way array
% Fac : Number of factors in PARAFAC model
%
% OPTIONAL INPUTS
% Options : A 1x2 vector
% (1) : Mode to compress to dimension two in
% DTLD if smallest mode-dimension is more
% than two. Default is the smallest dimension
% (2) : Number of extra components in Tucker3
% compression model. Default if no given
% is one
%
% ADVANCED OPTIONS
% DoComp : 0 => No compression
% 1 => Compression
% CompIt : CompIt => Max number of iterations in Compression
% Init : 0 => GRAM/DTLD
% 1 => Ten x PARAFAC with 10 iterations started from RandOrth
%
%
% I/O
% [A,B,C,FIT,IT]=comfac(X,Fac,Options,DoComp,CompIt,Init);
%
% Or short : [A,B,C]=comfac(X,Fac);
%
%
% 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. In case of doubt,
% contact the holder of the copyrights.
%
% Copyright 1998
% Rasmus Bro
% KVL, Denmark, [email protected]
% &
% Nikos Sidiropoulos
% Univ. Minnesota, [email protected]
%
% 31/3/00 RB fixed bug when dimensions not suitable for DTLD
crit = 1e-5; % Criterion used for fitting with Alternating LS algorithm
if length(size(X))~=3
error(' The data must be held in a three-way array')
end
DimX = size(X);
X = reshape(X,DimX(1),prod(DimX(2:3)));
dbg=1;
DeafultExtraCompInTucker=2;
if exist('DoComp')~=1|isempty(DoComp)
DoComp = 1;
end
if exist('CompIt')~=1|isempty(CompIt)
CompIt = 3;
end
if exist('Init')~=1|isempty(Init)
Init = 1;
end
if length(find(DimX>Fac))<2 % GRAM/ESPRIT cannot be used
ReplaceTLDwithRand = 1;
else
ReplaceTLDwithRand = 0;
end
if dbg
% home
disp(' ')
disp([' Fitting ',num2str(Fac),'-comp. PARAFAC model using COMFAC ...'])
disp([' Array size: ',num2str(DimX(1)),' x ',num2str(DimX(2)),' x ',num2str(DimX(3))])
end
if nargin<3|isempty(Options)
Options=[0 DeafultExtraCompInTucker];
end
if Options(1)==0
[out,CompToTwo] = min(DimX);
else
CompToTwo = Options(1);
end
if length(Options)<2
Options(2)=DeafultExtraCompInTucker; % Number of additional components in compression as compared to model
end
W = Options(2)+Fac;
W = [W W W];
for i=1:3
if W(i)>DimX(i)
W(i)=DimX(i);
end
end
% Display things
if DoComp & dbg
disp([' Array compressed to size ',num2str(W(1)),' x ',num2str(W(2)),' x ',num2str(W(3))])
end
% Compression
if DoComp
[Ut,Vt,Zt,Gt,fit]=tucker3(X,DimX,W,CompIt);
else
Gt = X;
W = DimX;
Ut=speye(DimX(1));
Vt=speye(DimX(2));
Zt=speye(DimX(3));
end
% Initialize PARAFAC using DTLD
if Init == 0
if ~ReplaceTLDwithRand
disp(' Initializing using direct trilinear decomposition')
[Adtld,Bdtld,Cdtld]=cdtld(Gt,W,Fac,CompToTwo);
else
disp(' Initializing using random values because tld cannot be used due to the sizes')
[Adtld,Bdtld,Cdtld,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
end
elseif Init == 1
disp(' Initializing using best of 11 initial small runs')
if ~ReplaceTLDwithRand
[Adtld,Bdtld,Cdtld]=cdtld(Gt,W,Fac,CompToTwo);
else
[Adtld,Bdtld,Cdtld,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
end
fitout=sum(sum(abs( Gt-Adtld*ppp(Bdtld,Cdtld).').^2));
for rep=1:10
[adtld,bdtld,cdtldvar,fit]=cparafac(Gt,W,Fac,crit,0,0,0,10);
if fit<fitout
Adtld=adtld;Bdtld=bdtld;Cdtld=cdtldvar;fitout=fit;
end
end
end
% Fit PARAFAC model in compressed space
if DoComp
disp(' Fitting PARAFAC in compressed space')
[Acomp,Bcomp,Ccomp,fit2,it2]=cparafac(Gt,W,Fac,crit,Adtld,Bdtld,Cdtld);
% Decompress to original domain
disp(' Transforming interim solution to original space')
Ainit = Ut*Acomp;
Binit = Vt*Bcomp;
Cinit = Zt*Ccomp;
else
Ainit = Adtld;
Binit = Bdtld;
Cinit = Cdtld;
it2=0;
fit2=NaN;
end
% Fit PARAFAC model in raw space
disp(' Fitting PARAFAC in original space')
[A,B,C,FIT,IT]=cparafac(X,DimX,Fac,crit,Ainit,Binit,Cinit);
disp(' Algorithm converged')
function [A,B,C,fit,it]=cparafac(X,DimX,Fac,crit,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)
% 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, [email protected]
% Initialization
if nargin<8
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<7
InitWithRandom=1;
end
if nargin>6 & size(A,1)~=I
InitWithRandom=1;
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 = sum(sum(abs(X-A*ppp(B,C).')));
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
A = Xbc*inv((B'*B).*(C'*C)).';
% Project X down on orth(A) - saves time if first mode is large
[Qa,Ra]=qr(A,0);
x=Qa'*X;
% Update B
Xac=0;
for k=1:K
Xac = Xac + x(:,(k-1)*J+1:k*J).'*conj(Ra*diag(C(k,:)));
end
B = Xac*inv((Ra'*Ra).*(C'*C)).';
% Update C
ab=inv((Ra'*Ra).*(B'*B));
for k=1:K
C(k,:) = (ab*diag(Ra'* x(:,(k-1)*J+1:k*J)*conj(B))).';
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
评论1