0001 function [A,score,AA]=nt_mcca(C,N)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 if nargin<2; error('!'); end
0012 if size(C,1) ~= size(C,2); error('!'); end
0013 if size(C,1) ~= round(size(C,1)/N)*N; error('!'); end
0014
0015
0016 nblocks=size(C,1)/N;
0017 for iBlock=1:nblocks
0018 idx=(iBlock-1)*N + (1:N);
0019 CC=C(idx,idx);
0020 [V, S] = eig(CC) ;
0021 V=real(V); S=real(S);
0022 [E,idx2] = sort(diag(S)', 'descend');
0023 topcs=V(:,idx2);
0024 EXP=1-10^-12;
0025 E=E.^EXP;
0026 EE=(1./E); EE(find(E<=0))=0;
0027 A(idx,idx)=topcs*diag(sqrt(EE));
0028 end
0029 C=A'*C*A;
0030
0031
0032
0033 [V, S] = eig(C) ;
0034 V=real(V); S=real(S);
0035 [E, idx] = sort(diag(S)', 'descend') ;
0036 topcs = V(:,idx);
0037 A=A*topcs;
0038
0039
0040 C=topcs'*C*topcs;
0041 score=diag(C);
0042
0043
0044 if nargout>2;
0045 AA=[];
0046 for iBlock=1:nblocks
0047 AA{iBlock}=A(N*(iBlock-1)+(1:N),:);
0048 end
0049 end
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071 return
0072
0073
0074
0075
0076 figure(1); clf;
0077 x1=randn(10000,10); x2=randn(10000,10); x3=randn(10000,10);
0078 x=[x1,x2,x3];
0079 C=x'*x;
0080 [A,score,AA]=nt_mcca(C,10);
0081 subplot 131; nt_imagescc(A); title('mcca transform');
0082 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
0083 subplot 133; nt_imagescc(x'*(x*A)); title('crosscorr between raw & transf'); xlabel('transformed'); ylabel('raw');
0084 z=x*A;
0085 figure(11); clf;
0086 plot(mean(z.^2));
0087
0088
0089 figure(2); clf
0090 x1=randn(10000,10); x=[x1,x1,x1];
0091 C=x'*x;
0092 [A,score,AA]=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,score,AA]=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,score,AA]=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