mfbox_constraintica.m

toolbox for spm 5 for data, model free analysis

function [S,A,W]=mfbox_constraintica(X,varargin)
% perform ica with an additional constraint on the columns of the mixing matrix
%
% usage: [S,A,W]=mfbox_constraintica(X[,arg,val[...]])
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

error(nargchk(1,Inf,nargin));

maxiter = 150;
constraint = [];
delta = 0.05;
alpha = 0.5;
verbose = 0;
cortype = 'jade';
method = 'simult';
nummat = Inf;

orgs = size(X);
X = reshape(X,orgs(1),[]);

feig = 1;
leig = orgs(1);

l = nargin-1;
for i=1:2:l
    op = varargin{i};
    pa = varargin{i+1};
    switch lower(op)
        case 'constraint'
            constraint = pa(:);
        case 'method'
            method = pa;
        case 'maxiter'
            maxiter = pa;
        case 'alpha'
            alpha = pa;
        case 'delta'
            delta = pa;
        case 'firsteig'
            feig = pa;
        case 'lasteig'
            leig = pa;
        case 'nummat'
            nummat = pa;
        case 'icatype'
            cortype = pa;
        case 'verbose'
           p = find(strcmp(pa,{'off','on'}));
           verbose = p(1)-1;
       otherwise
           warning('Illegal parameter %s ',op);
    end
end

Xmean = mean(X,1);
for i=1:size(X,2), X(:,i) = X(:,i)-Xmean(i); end

[E,D,wNoise] = mfbox_pcamat(X,feig,leig);
[wMat,dewMat] = mfbox_whitenv(E,D,wNoise);
X = real(wMat*X);
[M,T] = size(X);

if (isempty(constraint))
    constraint = ones(size(dewMat,1),1);
else
    constraint = constraint(:)-repmat(mean(constraint(:)),length(constraint),1);
    constraint = constraint/sqrt(constraint'*constraint);
end
lc = length(constraint);
la = size(dewMat,1);

switch cortype
    case 'jade'
        [e,C,CM] = mfbox_cumindependence(X,false,nummat);
    case 'sobi'
        [e,C,CM] = mfbox_sobimats(reshape(X',[orgs(2:end),N]),maskcolor,nummat);
end

delta = delta*sum(C(:));

dM = ones(size(dewMat,1),0);
for i=1:ceil(lc/4):(la-lc+1), dM(i:(i+lc-1),end+1) = constraint; end

Ro = eye(M);
A = normalize(eye(M)+randn(M));
W = A^(-1);
i = 1;
cpi = zeros(1,maxiter);
conl = zeros(1,maxiter); 
corl = zeros(1,maxiter); 

intv = floor(maxiter/20);

while (i<maxiter)
    switch method
        case 'simple'
            if (alpha<1)
                dW = findmincovgrad(CM,W,-1,(plotting&&(mod(i,intv)==1)));
                dW = (1-alpha)*dW;
            else
                dW = 0;
            end
            if (alpha>0)
                dA = findmaxconstraintgrad(dM,dewMat*A,-1,(plotting&&(mod(i,intv)==1)));
                dA = alpha*A'*dewMat'*dA*A'; % pullback + invertieren
            else
                dA = 0;
            end
            R = normalize(W+(dW+dA)*10^-3);
            W = R;

            if (alpha>0), A = W^(-1); end
            cpi(i) = norm(R-Ro,'fro');
            Ro = R;
        case 'seesaw'
            if (alpha<1)
                [dW,dWl] = findmincovgrad(CM,W,0,(plotting&&(mod(i,intv)==1)));
                nW = (1-alpha)*normalize(W+dW);
                corl(i) = [1,0]*dWl(:);
            else
                nW = 0;
            end

            W = normalize(nW);
            A = W^(-1);

            if (alpha>0)
                [dA,dAl] = findmaxconstraintgrad(dM,dewMat*A,0,(plotting&&(mod(i,intv)==1)));
                nA = alpha*normalize(A+wMat*dA)^(-1); % kein pullback !
                conl(i) = [1,0]*dAl(:);
            else
                nA = 0;
            end

            A = normalize(nA);
            W = A^(-1);
            R = W;

            if (cond(W)>1/eps)
                if (verbose), fprintf('Oops almost singular ...\n'); end
                W = W+rand(M);
                break;
            end

            cpi(i) = norm(R-Ro,'fro');
            Ro = R;
        case 'simult'
            if (alpha<1)
                [dW,dWl] = findmincovgrad(CM,W,1,(plotting&&(mod(i,intv)==1)));
                nW = (1-alpha)*normalize(W+dW);
                corl(i) = [1,0]*dWl(:);
            else
                nW = 0;
            end

            if (alpha>0)
                [dA,dAl] = findmaxconstraintgrad(dM,dewMat*A,1,(plotting&&(mod(i,intv)==1)));
                nA = alpha*normalize(A+wMat*dA)^(-1); % kein pullback, weil neue matrix und nicht gradient
                conl(i) = [1,0]*dAl(:);
            else
                nA = 0;
            end

            R = normalize(nA+nW);
            W = R;

            if (cond(W)>1/eps)
                if (verbose), fprintf('Oops almost singular ...\n'); end
                break;
                W = W+rand(M);
            end

            if (alpha>0), A = W^(-1); end
            cpi(i) = norm(R-Ro,'fro');
            Ro = R;
        case 'combinedsimult'
            if (alpha<1), dW = findmincovgrad(CM,W,-1,(plotting&&(mod(i,intv)==1))); dW = (1-alpha)*dW;
            else dW = 0;
            end

            if (alpha>0), dA = findmaxconstraintgrad(dM,dewMat*A,-1,(plotting&&(mod(i,intv)==1))); dA = alpha*A'*dewMat'*dA*A'; % pullback + invertieren
            else dA = 0;
            end

            dW = dA+dW; % combine

            if (alpha<1)
                [dW,dWl] = findmincovgrad(CM,W,dW,(plotting&&(mod(i,intv)==1)));
                nW = (1-alpha)*normalize(W+dW);
                corl(i) = [1,0]*dWl(:);
            else
                nW = 0;
            end

            if (alpha>0)
                [dA,dAl] = findmaxconstraintgrad(dM,dewMat*A,dewMat*W'*dW*W',(plotting&&(mod(i,intv)==1)));
                nA = alpha*normalize(A+wMat*dA)^(-1); % kein pullback, weil neue matrix und nicht gradient
                conl(i) = [1,0]*dAl(:);
            else
                nA = 0;
            end

            R = nA+nW;
            W = normalize(R);

            if (cond(W)>1/eps)
                if (verbose), fprintf('Oops almost singular ...\n'); end
                W = W+rand(M);
                break;
            end

            if (alpha>0), A = W^(-1); end
            cpi(i) = norm(R-Ro,'fro');
            Ro = R;
    end

    if (cpi(i)<delta), break; end
    if (verbose&&(mod(i,intv)==0))
        fprintf('Iteration %d (change %f correlation cost %f constraint cost %f)\n',i,cpi(i),corl(i),conl(i));
    end
    i = i+1;
end

if (alpha==0), A = W^(-1); end

cpi = cpi(1:(i-1));
corl = corl(1:(i-1));
conl = conl(1:(i-1));

if (i==maxiter), warning('No convergence after %d iterations',i); end

S = W*X+repmat(W*wMat*Xmean,1,T);
W = W*wMat; 
A = dewMat*A;

return


function [M]=normalize(M)

sM = sum(M.^2,2);
M = M./repmat(sqrt(sM),1,size(M,2));


function [G,l,J]=findmaxconstraintgrad(Md,W,f)

if (nargin<3), f = 0; end % normalize gradient (0 not, 1 full)

NUMRES = 10; % use that many residues at once
sW = size(W);
V = zeros(sW);

smd = size(Md,2);
if (smd==0)
    G = V; l = 0;
    return
end

docor = false;
if (smd==1)
    if (sum((Md(:)-1).^2)~=0), docor = true; end
end

for i=1:smd
    imd = find(~isnan(Md(:,i)));
    d = Md(imd,i); w = W(imd,:);
    dd = d*d';

    if (docor), o = eye(length(d))-1/length(d); w = o*w;
    else o = 1;
    end

    x = sum(w.^2,1)'; dix = diag(1./x);
    xd = (w'*d).^2; dxd = diag(xd);

    V(imd,:) = V(imd,:)+o*(dd*w-w*dix*dxd)*dix;
end
V = V*W'*W; % nat. gradient

if (length(f)>1)
    V = f;
elseif (f>0)
    X = W/norm(W(:),'fro');
    V = V-f*X*(V(:)'*X(:));
elseif (f<0)
    G = V;
    return
end

% V = V/norm(V,'fro');

MP = zeros(smd,size(W,2),6);
%DMP = zeros(smd,size(W,2),8);

for i=1:smd
    imd = find(~isnan(Md(:,i)));
    d = Md(imd,i); w = W(imd,:); v = V(imd,:);

    if (docor), o = eye(length(d))-1/length(d); w = o*w; v = o*v;
    end

    ww = sum(w.^2,1)'; wddw = (w'*d).^2;
    wv = sum(w.*v,1)'; wddv = (v'*d).*(w'*d);
    vv = sum(v.^2,1)'; vddv = (v'*d).^2;

    MP(i,:,:) = [vv,2*wv,ww,vddv,2*wddv,wddw];
end

MP = reshape(MP,[],6); %cost is sum_i pol(MP(i,4:6))/pol(MP(i,1:3))

if (exist('OCTAVE_HOME')==5)
    RPE = []; K = 0;
    for i=1:size(MP,1)
        [r,p,k,e] = residue(MP(i,4:6),MP(i,1:3));
        RPE(end+1:end+2,1:3) = [r(:),p(:),e(:)];
        K = polyadd(K,k);
    end
else % matlab
    RPE = []; K = 0;
    for i=1:size(MP,1)
        [r,p,k] = residue(MP(i,4:6),MP(i,1:3));
        e = ones(size(p));
        for j=2:length(p), if (p(j-1)==p(j)), e(j) = e(j-1)+1; end, end
        RPE(end+1:end+2,1:3) = [r(:),p(:),e(:)];
        K = polyadd(K,k);
    end
end
% diff it
K = polyderiv(K);
RPE(:,1) = -RPE(:,1).*RPE(:,3);
RPE(:,3) = RPE(:,3)+1;

[t,p] = sort(RPE(:,2)); % sort position of the residue
RPE = RPE(p,:);

l = 0; % TODO: faster
nrpe = size(RPE,1);
nv = ceil(max(1,nrpe-NUMRES)/2);
NN = cell(1,nv); ZZ = cell(1,nv);
for v=1:2:max(1,nrpe-NUMRES)
    p = RPE(v:min(v+NUMRES,nrpe),:);
    Z = 0; N = 1;
    for i=1:size(p,1)
        if (p(i,3)>1), j = 1; for k=1:p(i,3), j = conv(j,[1,-p(i,2)]); end
        else j = [1,-p(i,2)]; end
        Z = polyadd(conv(Z,j),p(i,1)*N);
        N = conv(N,j);
    end
    if (any(K~=0)), P = conv(N,K); P = polyadd(P,Z);
    else P = Z; end
    NN{(v+1)/2} = N'; 
    ZZ{(v+1)/2} = P';
    l = [l;roots(P)];
end

l = real(l);
%l = real(l((imag(l)).^2<eps));

t = cumprod(repmat(l(:)',2,1));
t = [t(2:-1:1,:);ones(1,size(t,2))];

J = 0;
for i=1:100:size(MP,1)
    j = i:min(i+99,size(MP,1));
    J = J+sum(((MP(j,4:6)*t)./(MP(j,1:3)*t)),1);
end


[ll,ld] = max(J);

G = l(ld)*V;
l = J([1,ld]);


function c=polyadd(a,b)
la = length(a);
lb = length(b);
if (la==lb)
    c = a+b;
elseif (la<lb)
    c = b;
    c((lb-la+1):lb) = c((lb-la+1):lb)+a;
else
    c = a;
    c((la-lb+1):la) = c((la-lb+1):la)+b;
end


function [G,l,J]=findmincovgrad(MC,W,f)

if (nargin<3), f = 0; end % normalize gradient (0 not, 1 full)

I = zeros(1,5);
V = zeros(size(W));
CW = cell(1,size(MC,3));
P = cell(1,size(MC,3));
Q = cell(1,size(MC,3));
QP = cell(1,size(MC,3));

for i=1:size(MC,3)
    C = squeeze(MC(:,:,i));
    CW{i} = C*W';
    P{i} = W*CW{i}; Q{i} = P{i}+diag(-diag(P{i}));
    QP{i} = Q{i}*P{i};
    V = V+(QP{i}*W);
%    V = V+(diag(diag(P{i}).^(-2))*QP{i}*W);
end

if (length(f)>1)
    V = f;
elseif (f>0)
    X = W/norm(W(:),'fro');
    V = V-f*X*(V(:)'*X(:));
elseif (f<0)
    G = V;
    return
end

% V = V/norm(V,'fro');

for i=1:size(MC,3)
    C = squeeze(MC(:,:,i));
    R = V*CW{i}; S = R+diag(-diag(R));
    T = V*C*V'; U = T+diag(-diag(T)); 

    I = I+[trace(T*U),...
        2*trace(R*U+T*S), ...
        trace(P{i}*U+(R+R')*(S+S')+T*Q{i}), ...
        2*trace(P{i}*S+R*Q{i}), ...
        trace(QP{i})];
end

J = (I(1:4).*(4:-1:1))/(4*I(1)); % diff it (TODO simply calculate diff'd version)

l = [0;(-(sqrt(3)*j)/2-1/2)*(sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3)+ ...
    (((sqrt(3)*j)/2-1/2)*(J(2)^2-3*J(3)))/(9*(sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3))-J(2)/3; ...
    ((sqrt(3)*j)/2-1/2)*(sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3)+ ...
    ((-(sqrt(3)*j)/2-1/2)*(J(2)^2-3*J(3)))/(9*(sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3))-J(2)/3; ...
    (sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3)+ ...
    (J(2)^2-3*J(3))/(9*(sqrt(27*J(4)^2+(4*J(2)^3-18*J(2)*J(3))*J(4)+4*J(3)^3-J(2)^2*J(3)^2)/(6*sqrt(3))-(27*J(4)-9*J(2)*J(3)+2*J(2)^3)/54)^(1/3))-J(2)/3];

l = real(l(imag(l).^2<eps));
J = polyval(I,l);


[ll,ld] = min(J);

G = l(ld)*V;
l = J([1,ld]);