Home > NoiseTools > nt_mcca.m

nt_mcca

PURPOSE ^

[A,score,AA]=nt_mcca(C,N) - multiple cca

SYNOPSIS ^

function [A,score,AA]=nt_mcca(C,N,Nkeep)

DESCRIPTION ^

[A,score,AA]=nt_mcca(C,N) - multiple cca

  A: transform matrix
  score: commonality score (ranges from 1 to N^2)
  AA: array of subject-specific MCCA transform matrices
 
  C: covariance matrix of aggregated data sets
  N: number of channels of each data set
  Nkeep: number of components to keep (for orthogonal transforms)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [A,score,AA]=nt_mcca(C,N,Nkeep)
0002 %[A,score,AA]=nt_mcca(C,N) - multiple cca
0003 %
0004 %  A: transform matrix
0005 %  score: commonality score (ranges from 1 to N^2)
0006 %  AA: array of subject-specific MCCA transform matrices
0007 %
0008 %  C: covariance matrix of aggregated data sets
0009 %  N: number of channels of each data set
0010 %  Nkeep: number of components to keep (for orthogonal transforms)
0011 
0012 if nargin<3; Nkeep=[]; end
0013 if nargin<2; error('!'); end
0014 if size(C,1) ~= size(C,2); error('!'); end
0015 if size(C,1) ~= round(size(C,1)/N)*N; error('!'); end
0016 
0017 % sphere by blocks
0018 nblocks=size(C,1)/N;
0019 for iBlock=1:nblocks
0020     idx=(iBlock-1)*N + (1:N);
0021     CC=C(idx,idx);
0022     [V, S] = eig(CC) ;  
0023     V=real(V); S=real(S);
0024     [E,idx2] = sort(diag(S)', 'descend');
0025     topcs=V(:,idx2);
0026     EXP=1-10^-12; 
0027     E=E.^EXP; % break symmetry when x and y perfectly correlated (otherwise cols of x*A and y*B are not orthogonal)
0028     EE=(1./E); EE(find(E<=0))=0;
0029     A(idx,idx)=topcs*diag(sqrt(EE));
0030 end
0031 C=A'*C*A;
0032 
0033 
0034 % final PCA
0035 [V, S] = eig(C) ;
0036 V=real(V); S=real(S);
0037 [E, idx] = sort(diag(S)', 'descend') ;
0038 topcs = V(:,idx);
0039 A=A*topcs;
0040 %A=A(:,1:N);
0041 
0042 C=topcs'*C*topcs;
0043 score=diag(C);
0044 
0045 
0046 if nargout>2;
0047     AA=[];
0048     for iBlock=1:nblocks
0049         AA{iBlock}=A(N*(iBlock-1)+(1:N),:);
0050     end
0051 end
0052 
0053 if nargout>3;
0054     if isempty(Nkeep); error('must specify Nkeep'); end
0055     AAA=[];
0056     for iBlock=1:nblocks
0057         % covariance of subject's data
0058         idx=(iBlock-1)*N + (1:N);
0059         C11=C(idx,idx);
0060         % covariance of selected MCCA components
0061         tmp=A(:,1:Nkeep);
0062         C22=tmp'*C*tmp;
0063         % cross covariance between subject's data and transformed data
0064         C12=C(idx,:)*tmp; clear tmp
0065         C21=C12';
0066         % CCA:
0067         [tmp]=nt_cca([],[],[],[C11,C12;C21,C22],N);
0068         AAA{iBlock}=tmp;
0069     end
0070 end
0071 
0072 
0073 return
0074 
0075 % test code
0076 
0077 % 3 uncorrelated data sets
0078 figure(1); clf;
0079 x1=randn(10000,10); x2=randn(10000,10); x3=randn(10000,10); 
0080 x=[x1,x2,x3];
0081 C=x'*x;
0082 A=nt_mcca(C,10);
0083 subplot 131; nt_imagescc(A); title('mcca transform');
0084 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
0085 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
0086 
0087 
0088 % 3 identical data sets
0089 figure(2); clf
0090 x1=randn(10000,10); x=[x1,x1,x1]; 
0091 C=x'*x; 
0092 A=nt_mcca(C,10);
0093 subplot 131; nt_imagescc(A); title('mcca transform');
0094 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
0095 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
0096 
0097 % 3 data sets with shared parts
0098 figure(3); clf
0099 x1=randn(10000,5); x2=randn(10000,5); x3=randn(10000,5); x4=randn(10000,5); 
0100 x=[x2,x1,x3,x1,x4,x1];
0101 C=x'*x; 
0102 A=nt_mcca(C,10);
0103 subplot 131; nt_imagescc(A); title('mcca transform');
0104 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
0105 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
0106 
0107 % 3 data sets with parts shared 2 by 2
0108 figure(4); clf
0109 x1=randn(10000,5); x2=randn(10000,5); x3=randn(10000,5); x4=randn(10000,5); 
0110 x=[x2,x1,x3,x1,x3,x4];
0111 C=x'*x; 
0112 A=nt_mcca(C,10);
0113 subplot 131; nt_imagescc(A); title('mcca transform');
0114 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
0115 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
0116 
0117 
0118 
0119 
0120 
0121

Generated on Thu 30-Nov-2017 17:26:18 by m2html © 2005