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)
      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:

- nt_cca [A,B,R]=nt_cca(x,y,shifts,C,m,thresh,demeanflag) - canonical correlation
- nt_imagescc nt_imagescc - plot image with symmetric scaling
- nt_mcca [A,score,AA]=nt_mcca(C,N) - multiple cca

This function is called by:

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

## 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)
    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,score,AA]=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('crosscorr between raw & transf'); xlabel('transformed'); ylabel('raw');
    0086 z=x*A;
    0087 figure(11); clf;
    0088 plot(mean(z.^2));
    0089
    0090 % 3 identical data sets
    0091 figure(2); clf
    0092 x1=randn(10000,10); x=[x1,x1,x1];
    0093 C=x'*x;
    0094 [A,score,AA]=nt_mcca(C,10);
    0095 subplot 131; nt_imagescc(A); title('mcca transform');
    0096 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
    0097 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
    0098
    0099 % 3 data sets with shared parts
    0100 figure(3); clf
    0101 x1=randn(10000,5); x2=randn(10000,5); x3=randn(10000,5); x4=randn(10000,5);
    0102 x=[x2,x1,x3,x1,x4,x1];
    0103 C=x'*x;
    0104 [A,score,AA]=nt_mcca(C,10);
    0105 subplot 131; nt_imagescc(A); title('mcca transform');
    0106 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
    0107 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
    0108
    0109 % 3 data sets with parts shared 2 by 2
    0110 figure(4); clf
    0111 x1=randn(10000,5); x2=randn(10000,5); x3=randn(10000,5); x4=randn(10000,5);
    0112 x=[x2,x1,x3,x1,x3,x4];
    0113 C=x'*x;
    0114 [A,score,AA]=nt_mcca(C,10);
    0115 subplot 131; nt_imagescc(A); title('mcca transform');
    0116 subplot 132; nt_imagescc(A'*C*A); title('covariance of transformed data');
    0117 subplot 133; nt_imagescc(x'*(x*A)); title('cross correlation'); xlabel('transformed'); ylabel('raw');
    0118