 Martin J. Gander Section de Mathématiques Université de Genève Rue du Conseil-Général 9, CP 64 1211 Genève 4, Suisse 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