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.
|
