Martin J. Gander
Section de Mathématiques
2-4 rue du Lièvre, CP 64
CH-1211 Genève
Martin.Gander(at)unige.ch   Fax: +41 22 379 11 76

Home Research Teaching Consulting Personal  
   
Seminars
Genève:
- Analyse numérique
McGill:
- CSE
- Applied Mathematics
Conference
Ascona 2013

Available projects
Postdoctoral:
- upon request
Doctoral degree:
- upon request
Master thesis:
- upon request

:: The FG-Algorithm
The FG-Algorithm is a very elegant iterative algorithm invented by Bernhard Flury and Walter Gautschi, which tries to simultaneously diagonalize a set of symmetric positive definite matrices, see

An Algorithm for Simultaneous Orthogonal Transformation of Several Positive Definite Symmetric Matrices to Nearly Diagonal Form, SIAM J. Sci. Stat. Comput. 7(1), pp. 169-184, 1986.

Working on the comments on selected papers by Walter Gautschi, I studied this paper, and decided to implement the algorithm in Matlab: here is the result
        
function B=F(A,n,tol);
% F F-step of the FG algorithm by Flury and Gautschi    
%   B=F(A); takes a structure A containing length(A) matrices of
%   the same size, and computes a matrix B which simultanously
%   diagonalizes these matrices as well as possible
p=size(A{1},1);
B=eye(p);                      % using identity fails in special cases:
% B=orth(rand(p));             % use rand if the method does not converge
e=1;                           % for the prereading scheme
while e>tol 
  Bold=B;
  % phi(A,B,n), for i=1:length(A), B'*A{i}*B, end; pause   % stepwise output
  for l=1:p-1
    for j=l+1:p
      for i=1:length(A)
        T{i}=[B(:,l)'*A{i}*B(:,l) B(:,l)'*A{i}*B(:,j);
              B(:,j)'*A{i}*B(:,l) B(:,j)'*A{i}*B(:,j)];
      end;      
      Q=Gold(T,n,tol);
      H=[B(:,l) B(:,j)]*Q;
      B(:,[l,j])=H;
    end;
  end; 
  e=norm(Bold-B);
end;

function Q=G(T,n,tol);
% G G-step of the FG algorithm by Flury and Gautschi
Q=eye(2); 
e=1;                                           % for the prereading scheme
while e>tol
  Qold=Q; U=zeros(2);
  for i=1:length(T)
    di1=Q(:,1)'*T{i}*Q(:,1);
    di2=Q(:,2)'*T{i}*Q(:,2);
    U=U+n(i)*(di1-di2)/(di1*di2)*T{i};
  end;  
  [Q,E]=eig(U);
  if abs(Q(1,1))==min(min(abs(Q)))             % construct unique rotation:
    Q=[Q(:,2) Q(:,1)];                         % chose cos to be maximal
  end;                            
  Q=[sign(Q(1,1))*Q(:,1) sign(Q(2,2))*Q(:,2)]; % and of positive sign 
  e=norm(Qold-Q);
end;

function r=phi(A,B,n)
r=1;
for i=1:length(A)
  C=B'*A{i}*B;
  r=r*(det(diag(diag(C)))/det(C))^n(i);
end;
       
      
In order to test the algorithm, one can use the commands
       
A{1}=[45 10  0  5   0   0
      10 45  5  0   0   0
       0  5 45 10   0   0
       5  0 10 45   0   0
       0  0  0  0 16.4 -4.8
       0  0  0  0 -4.8 13.6];
A{2}=[ 27.5 -12.5 -0.5 -4.5 -2.04  3.72
      -12.5  27.5 -4.5 -0.5  2.04 -3.72
       -0.5  -4.5 24.5 -9.5 -3.72 -2.04
       -4.5  -0.5 -9.5 24.5  3.72  2.04
       -2.04 2.04 -3.72 3.72 54.76 -4.68
        3.72 -3.72 -2.04 2.04 -4.68 51.24]
B=F(A,[1,1],0.0001)
       
        
which reproduce the example in the paper. The algorithm really finds the matrix that simultaneously diagonalizes the set of matrices, as one can see with the example below suggested by Ivan Graham during a conference in Urumqi in August 2012:
       
A{1}=-[-2 1 0 0
      1 -2 1 0
      0 1 -2 1
      0 0 1 -2];
A{2}=[4 1 0 0
      1 4 1 0
      0 1 4 1
      0 0 1 4];
B=F(A,[1,1],0.0001)
       
        
but now one gets an infinite loop with the initial guess B=eye(p) in the F algorithm: one needs to use the alternative random initial guess for the method to work.
   
Home Research Teaching Consulting Personal