nt_quad2square

[tosquare,D]=nt_quad2square(toquad,order) - quadratic to squared linear component
 
   tosquares: matrix to linear components with closest square
   D: scores

   toquad: vector that defines quadratic component from cross products
   order: order of cross products, 'colwise', 'diagwise' (default)

0001 function [tosquares,D]=nt_quad2square(toquad,order)
0002 %[tosquare,D]=nt_quad2square(toquad,order) - quadratic to squared linear component
0003 %
0004 %   tosquares: matrix to linear components with closest square
0005 %   D: scores
0006 %
0007 %   toquad: vector that defines quadratic component from cross products
0008 %   order: order of cross products, 'colwise', 'diagwise' (default)
0009 
0010 if nargin<2; order='diagwise'; end
0011 if nargin<1; error('!'); end
0012 
0013 if size(toquad,2)~=1; error('toquad should be column vector'); end
0014 
0015 nquads=size(toquad,1);
0016 nchans=floor(sqrt(nquads*2));
0017 if nquads ~= nchans*(nchans+1)/2; 
0018     [nchans nquads]
0019     error('unexpected size for toquad');
0020 end
0021 
0022 ii=1;
0023 A=zeros(nchans);
0024 switch order
0025     case 'colwise'
0026         for iRow=1:nchans;
0027             for iCol=1:iRow
0028                 if iRow==iCol;
0029                     A(iRow,iCol)=toquad(ii);
0030                 else
0031                     A(iRow,iCol)=toquad(ii)/2;
0032                     A(iCol,iRow)=toquad(ii)/2;
0033                 end
0034                 ii=ii+1;
0035             end
0036         end
0037     case 'diagwise'
0038         for iDiag=0:nchans-1
0039             for kk=1:(nchans-iDiag)
0040                 iRow=kk+iDiag;
0041                 iCol=kk;
0042                 if iRow==iCol;
0043                     A(iRow,iCol)=toquad(ii);
0044                 else
0045                     A(iRow,iCol)=toquad(ii)/2;
0046                     A(iCol,iRow)=toquad(ii)/2;
0047                 end
0048                 ii=ii+1;
0049             end
0050         end
0051     otherwise
0052         error('unexpected order');
0053 end
0054 
0055 % eigenvectors & values
0056 [V,D]=eig(A);
0057 D=diag(D);
0058 [~,idx]=sort(abs(D),'descend'); 
0059 D=D.^2/sum(D.^2);
0060 V=V(:,idx);
0061 D=D(idx);
0062 tosquares=V;