cpfunc.m

这是一个三面阵列的数学工具软件

function [A,B,C,f,fp,iter,tripcos,ftiter,mintripcos,ttt]=cpfunc(X,n,m,p,r,ort1,ort2,ort3,start,conv,A,B,C,maxit);
% [A,B,C,f,fp,iter,tripcos,ftiter,mintripcos,ttt]=cpfunc(X,n,m,p,r,ort1,ort2,ort3,start,conv,A,B,C,maxit);
% program candecomp/parafac, minimizes sum(k) || X(k) - A D(k) B' ||
%
% Input: X: data
%        n,m,p: order
%        r: dimensionality of solution
%        ort1, ort2, ort3: constraints on A,B and C, resp., indicated by values
%         			1=No constraint
%						2=Orthogonality
%						3=Zero correlations
% 			start: 
%						0 = from SVD's
%						1 = random, orthonormalized
%						2 = user specified component matrices
%						3 = from Generalized eigenvalue decomposition, only based on first 2 slices
%			conv: convergence criterion, e.g., 1e-6
%        A,B,C: starting values for component matrices (only if start=2)
%        maxit: maximal number of iterations
% Output:
%        A,B,C = component matrices
%        f = loss function value
%			fp = fit percentage
%        iter= number of iterations
%        tripcos = minimal triple cosine between two components, across three component matrices
%        ftiter= function values and minimal triple cosine at every 10 iterations 
%	      ttt=computation time 
% uses permnew.m ssq.m tr.m

ftiter=zeros(maxit/10,2);
mintripcos=0;

[nn,mp]=size(X);
if nn~=n | mp~=m*p, 
  beep
  fprintf(' Specified size of array (%g x %g x %g) does not correspond to actual size\n',n,m,p');
  keyboard
end;

t0=clock;
ssx=ssq(X);

if start==0 
  % rational starts via eigendecompositions, or random if r is too big
  if n>=r,[K,L]=ed(X*X');A=K(:,1:r);else A=orth(rand(r,r)-.5);A=A(1:n,:);end;
  Z=permnew(X,n,m,p);		% yields m x p x n array
  if m>=r,[K,L]=ed(Z*Z');B=K(:,1:r);else B=orth(rand(r,r)-.5);B=B(1:m,:);end;
  Z=permnew(Z,m,p,n);		% yields p x n x m
  if p>=r,[K,L]=ed(Z*Z');C=K(:,1:r);else C=orth(rand(r,r)-.5);C=C(1:p,:);end;
end;    

if start==1
  if n>=r,A=orth(rand(n,r)-.5);else A=orth(rand(r,r)-.5);A=A(1:n,:);end;
  if m>=r,B=orth(rand(m,r)-.5);else B=orth(rand(r,r)-.5);B=B(1:m,:);end;
  if p>=r,C=orth(rand(p,r)-.5);else C=orth(rand(r,r)-.5);C=C(1:p,:);end;
end;

if start==3
  disp(' only based on two slabs !! ')

  [k,l]=ed((A'*X(:,m+1:2*m)*B)\(A'*X(:,1:m)*B)) ; 
  T1=eye(r);T2=eye(r);
  [rhok,argk]=complpol(k);
  [rhol,argl]=complpol(diag(l));
  j=1;
  while j<=r
    if abs(imag(l(j,j)))<.00000001  % real eigenvalue
      if abs(imag(k(1,j)))>=.00000001 % complex eigenvector
        T1(j,j)=exp(i*argk(1,j));
      end;
    end;
    if abs(imag(l(j,j)))>=.00000001  % j-th en j+1-th are complex eigenvalues
      c=argk(1,j)+argk(1,j+1);
      T1(j,j)=exp(i*c/2);
      T1(j+1,j+1)=exp(i*c/2);
      T2(j:j+1,j:j+1)=[1 1;i -i];
      j=j+1;
    end;
    j=j+1;
  end;
  k=real(k/T1/T2)
  l=T2*T1*l/T1/T2
  l=real(diag(diag(l)));
  C(2,:)=ones(1,r);
  C(1,:)=diag(l)';
  A=X(:,m+1:2*m)*B*k;
  B=B/k'
end;

H=zeros(r,r^2);		% superidentity 3-way array
for ii=1:r	
  H(ii,(ii-1)*r+ii)=1;
end;

H1=permnew(H,r,r,r);
H1=permnew(B*H1,m,r,r);
H1=permnew(C*H1,p,r,m);  % gives r x m*p array  H*C'xB'
f=ssq(X-A*H1);
fprintf(' Function value at Start is %12.8f \n',f);
fold=f+2*conv*f;
iter=0;
BB=B'*B;CC=C'*C;

while (fold-f>conv*f | iter<2) & f>conv^2 & iter<maxit

  fold=f;  
  
  % Update A
  % first compute XF = X H C'xB' en FF' (implicitly)

  Z1=permnew(X,n,m,p);
  Z1=permnew(B'*Z1,r,p,n);
  Z1=permnew(C'*Z1,r,n,r);   % gives n x r^2 array
  XF=Z1*H';
  if ort1==1
    FF=BB.*CC;
    A=XF/FF;
  end;
  if ort1==2
    [P,D,Q]=svd(XF,0);
    A=P*Q';
  end;
  if ort1==3
    FF=BB.*CC;
    [P,D,Q]=svd(XF-ones(n,1)*mean(XF),0);
    A=P*Q'+ones(n,1)*mean(XF)/FF;
  end;
  AA=A'*A;

  % Update B
  % first compute XF = X H A'xC' en FF' (implicitly)

  Z=permnew(X,n,m,p);
  Z1=permnew(Z,m,p,n);
  Z1=permnew(C'*Z1,r,n,m);
  Z1=permnew(A'*Z1,r,m,r);   % gives m x r^2 array
  XF=Z1*H';
  if ort2==1
    FF=AA.*CC;
    B=XF/FF;
  end;
  if ort2==2
    [P,D,Q]=svd(XF,0);
    B=P*Q';
  end;
  if ort2==3
    FF=AA.*CC;
    [P,D,Q]=svd(XF-ones(m,1)*mean(XF),0);
    B=P*Q'+ones(m,1)*mean(XF)/FF;
  end;
  BB=B'*B;

  % Update C
  % first compute XF = X H B'xA' en FF' (implicitly)

  Z=permnew(Z,m,p,n);
  Z1=permnew(Z,p,n,m);
  Z1=permnew(A'*Z1,r,m,p);
  Z1=permnew(B'*Z1,r,p,r);   % gives p x r^2 array
  XF=Z1*H';
  if ort3==1
    FF=AA.*BB;
    C=XF/FF;
  end;
  if ort3==2
    [P,D,Q]=svd(XF,0);
    C=P*Q';
  end;
  if ort3==3
    FF=AA.*BB;
    [P,D,Q]=svd(XF-ones(p,1)*mean(XF),0);
    C=P*Q'+ones(p,1)*mean(XF)/FF;
  end;
  CC=C'*C; 

  % Evaluate

  if ort3==1,f=ssx-tr(CC*FF);end;
  if ort3~=1,
    H1=permnew(H,r,r,r);
    H1=permnew(B*H1,m,r,r);
    H1=permnew(C*H1,p,r,m);  % gives r x m*p array  H*C'xB'
    f=ssq(X-A*H1);
  end;

  iter=iter+1;
  if rem(iter,10)==0
    tripcos=min(min(phi(A,A).*phi(B,B).*phi(C,C)));
    if iter==10,mintripcos=tripcos;end;
    if tripcos<mintripcos,mintripcos=tripcos;end;
    if rem(iter,1000)==0,fprintf(' Minimal Triple cosine = %12.8f \n',tripcos);end;
    ftiter(iter/10,:)=[f tripcos];
  end;
  if rem(iter,50)==0,fprintf('f= %12.8f after %g iters ; diff.=%12.8f \n',f,iter,(fold-f));end;
end;
ttt=etime(clock,t0);
fp=100-100*f/ssq(X);
tripcos=min(min(phi(A,A).*phi(B,B).*phi(C,C)));
if iter<10,mintripcos=tripcos;end;
fprintf(' Function value is %12.8f after %g iterations \n',f,iter);
fprintf(' Fit percentage is %6.2f %%\n',fp);
fprintf(' Procedure used %g seconds \n',ttt);