mfbox_rjd.m

toolbox for spm 5 for data, model free analysis

function [V,A]=mfbox_rjd(A,threshold)
% Copyright by Peter Gruber, 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
%
% joint diagonalization (possibly
% approximate) of REAL matrices.
% This function minimizes a joint diagonality criterion
% through n matrices of size m by m.
%
% Input :
% * the  n by nm matrix A is the concatenation of m matrices
%   with size n by n. We denote A = [ A1 A2 .... An ]
% * threshold is an optional small number (typically = 1.0e-8 see below).
%
% Output :
% * V is a n by n orthogonal matrix.
% * qDs is the concatenation of (quasi)diagonal n by n matrices:
%   qDs = [ D1 D2 ... Dn ] where A1 = V*D1*V' ,..., An =V*Dn*V'.
%
% The algorithm finds an orthogonal matrix V
% such that the matrices D1,...,Dn  are as diagonal as possible,
% providing a kind of `average eigen-structure' shared
% by the matrices A1 ,..., An.
% If the matrices A1,...,An do have an exact common eigen-structure
% ie a common othonormal set eigenvectors, then the algorithm finds it.
% The eigenvectors THEN are the column vectors of V
% and D1, ...,Dn are diagonal matrices.
% 
% The algorithm implements a properly extended Jacobi algorithm.
% The algorithm stops when all the Givens rotations in a sweep
% have sines smaller than 'threshold'.
% In many applications, the notion of approximate joint diagonalization
% is ad hoc and very small values of threshold do not make sense
% because the diagonality criterion itself is ad hoc.
% Hence, it is often not necessary to push the accuracy of
% the rotation matrix V to the machine precision.
% It is defaulted here to the square root of the machine precision.
%
% Author : Jean-Francois Cardoso. [e-mail]
% This software is for non commercial use only.
% It is freeware but not in the public domain.
% A version for the complex case is available
% upon request at [e-mail]
%-----------------------------------------------------
% Two References:
%
% The algorithm is explained in:
%
%@article{SC-siam,
%   HTML =	"ftp://sig.enst.fr/pub/jfc/Papers/siam_note.ps.gz",
%   author = "Jean-Fran\c{c}ois Cardoso and Antoine Souloumiac",
%   journal = "{SIAM} J. Mat. Anal. Appl.",
%   title = "Jacobi angles for simultaneous diagonalization",
%   pages = "161--164",
%   volume = "17",
%   number = "1",
%   month = jan,
%   year = {1995}}
%
%  The perturbation analysis is described in
%
%@techreport{PertDJ,
%   author = "{J.F. Cardoso}",
%   HTML =	"ftp://sig.enst.fr/pub/jfc/Papers/joint_diag_pert_an.ps",
%   institution = "T\'{e}l\'{e}com {P}aris",
%   title = "Perturbation of joint diagonalizers. Ref\# 94D027",
%   year = "1994" }

if (nargin<2), threshold = sqrt(eps); end;

[m,nm] = size(A);

V = eye(m);

B = [1,0,0;0,1,1;0,-i,i];
Bt = B';

encore = true;
while (encore)
    encore = false;
    for p=1:(m-1)
        for q=(p+1):m
            Ip = p:m:nm; Iq = q:m:nm;

            g = [A(p,Ip)-A(q,Iq);A(p,Iq);A(q,Ip)]; 
            [vcp,D] = eig(real(B*(g*g')*Bt));
            [la,K] = sort(diag(D));
            angles = vcp(:,K(3));
            if (angles(1)<0), angles = -angles; end
            c = sqrt(0.5+angles(1)/2);
            s = 0.5*complex(angles(2),-angles(3))/c; 

            if (abs(s)>threshold) %%% updates matrices A and V by a Givens rotation
                encore = true;
                pair = [p;q];
                G = [c,-conj(s);s,c];
                V(:,pair) = V(:,pair)*G;
                A(pair,:) = G'*A(pair,:);
                A(:,[Ip,Iq]) = [c*A(:,Ip)+s*A(:,Iq),-conj(s)*A(:,Ip)+c*A(:,Iq)];
            end
        end
    end
end

return