0001 function [A,score,AA]=nt_mcca(C,N,Nkeep)
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
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;
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
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
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
0058 idx=(iBlock-1)*N + (1:N);
0059 C11=C(idx,idx);
0060
0061 tmp=A(:,1:Nkeep);
0062 C22=tmp'*C*tmp;
0063
0064 C12=C(idx,:)*tmp; clear tmp
0065 C21=C12';
0066
0067 [tmp]=nt_cca([],[],[],[C11,C12;C21,C22],N);
0068 AAA{iBlock}=tmp;
0069 end
0070 end
0071
0072
0073 return
0074
0075
0076
0077
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
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
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
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