Dear Edward,
> We met briefly at the Brain Mapping meeting in San Antonio where we
> discussed the procedure for a nonlinear PCA analysis. Could you please
> post the code for that analysis and variables I should be aware of in
> working through it. Thanks in advance for your assistance.
My pleasure - here it is. We have not as yet written this up or evaluated
its limits but you are welcome to see how it works.
With best wishes - Karl
===========================================================================
function [U,V,V0,V1,V2] = spm_nlpca(Y,n,q)
% Nonlinear PCA - based on a second order model
% FORLAT [U,V,V0,V1,V2] = spm_nlpca(Y,n,q);
% Y - {m k} data matrix
% n - number of sources to model
% q - number of components to explain
% U - {m n} 1st order variates
% V - {m p} 2nd order variates
% V0 - {1 k} 0th order modes
% V1 - {n k} 1st order modes
% V2 - {p k} 2nd order modes
%__________________________________________________________________________
%
% Orthogonal decomposition into zeroth, first and second order components
%
% Y = ones(m,1)*V0 + U*V1 + V*V2
%
% where V = U(:,i).*U(:,j) and phi(U) = 2./(1 + exp(-U)) - 1;
%
% subject to the constraint that U are orthogonal
%__________________________________________________________________________
% dimension reduction
%--------------------------------------------------------------------------
[X S W] = svd(Y,0);
W = W(:,[1:q]);
W = W/std(Y*W(:,1));
X = Y*W;
W = pinv(W);
% minimization
%--------------------------------------------------------------------------
OPT = foptions;
OPT(14) = 512*q*n;
F = fmins('pca_update',eye(q,n),OPT,[],X);
% output arguments
%--------------------------------------------------------------------------
[H,U,V,V0,V1,V2] = pca_update(F,X);
V0 = V0*W;
V1 = V1*W;
V2 = V2*W;
% total SSQ (adjusted for zeroth order effects)
%--------------------------------------------------------------------------
Q = sum(diag(S).^2) - sum(V0'.^2);
% grphical output
%--------------------------------------------------------------------------
subplot(2,2,1)
bar(sum(U.^2).*sum(V1'.^2)/Q)
xlabel('source')
ylabel('proportion of variance')
title('1st order effects')
axis square
subplot(2,2,2)
bar(sum((V - U*(pinv(U)*V)).^2).*sum(V2'.^2)/Q)
xlabel('interactions')
title('2nd order effects')
axis square
subplot(4,2,6)
imagesc(spm_en(V1'))
title('1st order modes')
ylabel('observation')
axis square
subplot(4,2,8)
imagesc(spm_en(V2'))
title('2nd order modes')
xlabel('source/interactions')
ylabel('observation')
axis square
subplot(4,2,5)
imagesc(spm_en(U))
xlabel('sources')
ylabel('observation')
title('expression of sources')
axis square
subplot(4,2,7)
plot(U(:,1),U(:,2),'.')
xlabel('1st source')
ylabel('2nd source')
axis square
grid on
===========================================================================
function [H,U,V,V0,V1,V2] = pca_update(F,X)
%__________________________________________________________________________
% 1st order nodes
%--------------------------------------------------------------------------
[m k] = size(X);
[k n] = size(F);
p = n*(n - 1)/2;
% Lateral weights
%--------------------------------------------------------------------------
[e v] = eig(F'*X'*X*F);
C = sqrt(v)*e';
L = inv(C)*diag(diag(C));
U = X*F*L;
% 2nd order nodes
%--------------------------------------------------------------------------
V = zeros(m,p);
q = 1;
for i = 1:n
for j = (i + 1):n
V(:,q) = phi(U(:,i).*U(:,j));
q = q + 1;
end
end
% LSQ estimates
%--------------------------------------------------------------------------
Q = [ones(m,1) U V];
P = pinv(Q)*X;
V0 = P(1,:);
V1 = P([1:n] + 1,:);
V2 = P([1:p] + 1 + n,:);
% objective
%--------------------------------------------------------------------------
R = X - Q*P;
H = sum(sum(R.^2));
function [y] = phi(x)
y = 2./(1 + exp(-x)) - 1;
===========================================================================
function [X] = spm_en(X)
% Euclidean normalization
% FORMAT [X] = spm_en(X);
% X - matrix
%___________________________________________________________________________
%
% spm_en performs a Euclidean normalization setting column-wise sum of
% squares to unity
for i = 1:size(X,2)
X(:,i) = X(:,i)/sqrt(sum(X(:,i).^2));
end
===========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|