Solving Generalized eigenvectors and eigen values to get a common basis

Question

I am looking for solving a generalized eigenvectors and eigen value problem in Matlab. For this, I have tested 2 methods.

1. if Generalized problem is formulated as :

Then, we could multiply by B^(-1) on each side, such as :

So, from a theorical point of view, it is the simple and classical eigen values problem.

Finally, in Matlab, I simply did, with A=FISH_sp and B=FISH_xc :

[Phi, Lambda] = eig(inv(FISH_xc)*FISH_sp);

But results are not correct when I make after a simple Fisher synthesis (constraints are too bad and also making appear nan values. I don't know why I don't get the same results than the second ones below.

2. the second method comes from the following paper.

To summarize, the algorithm used is described on page 7. I have followed all the steps of this algorithm and it seems to give better results when I make the Fisher synthesis.

Here the interested part (sorry, I think Latex is not available on stakoverflow) :

Here my little Matlab script for this method :

% Diagonalize A = FISH_sp and B = Fish_xc
[V1,D1] = eig(FISH_sp);
[V2,D2] = eig(FISH_xc);

% Applying each  step of algorithm 1 on page 7
phiB_bar = V2*(D2.^(0.5)+1e-10*eye(7))^(-1);
barA = inv(phiB_bar)*FISH_sp*phiB_bar;
[phiA, vA] = eig(barA);
Phi = phiB_bar*phiA;

So at the end, I find phi eigenvectors matrix (phi) and lambda diagonal matrix (D1).

3. Now, I would like to do the link between this generalized problem and the eventual common eigenvectors between A and B matrices (respectively Fish_sp and Fish_xc). Is there a way to perform this?

Indeed, what I have done up to now is to to find a parallel relation between A*Phi and B*Phi, linked by Lambda diagonal matrix. Maybe, we could arrange this relation such that :

A*Phi'=Phi'*Lambda_A'

and

B*Phi'=Phi'*Lambda_B'

4. From a numerical point of view, why don't I get the same results between the method in 1) and the method in 3) ? I mean about the Phi eigen vectors matrix and the Lambda diagonal matrix.

However, this is the same formulation.

EDIT :

I have wrong results if I want to say that phi diagonalizes both A=FISH_sp and B=FISH_xc matrices.

Indeed, by doing :

% Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = inv(FISH_GCsp);
COV_XC_first = inv(FISH_XC);
COV_GCsp = COV_GCsp_first(1:N,1:N);
COV_XC = COV_XC_first(1:N,1:N);
% Invert to get Fisher matrix
FISH_sp = inv(COV_GCsp);
FISH_xc = inv(COV_XC);
% Diagonalize
[V1,D1] = eig(FISH_sp);
[V2,D2] = eig(FISH_xc);

% Build phi matrix
% V2 corresponds to eigen vectors of FISH_xc
phiB_bar = V2*diag(diag(D2.^(-0.5)));
% DEBUG : check identity matrix => OK, Identity matrix found !
id = (phiB_bar')*FISH_xc*phiB_bar
% phi matrix
barA = (phiB_bar')*FISH_sp*phiB_bar
[phiA, vA] = eig(barA);
phi = phiB_bar*phiA;

% Check eigen values : OK, columns of eigenvalues found !
FISH_sp*V1./V1
% Check eigen values : OK, columns of eigenvalues found !
FISH_xc*V2./V2

% Check if phi diagolize FISH_sp : NOT OK, not identical eigenvalues 
FISH_sp*phi./phi
% Check if phi diagolize FISH_sp : NOT OK, not identical eigenvalues 
FISH_xc*phi./phi

So, I don't find that matrix of eigenvectors Phi diagonalizes A and B since the eigenvalues expected are not columns of identical values.

By the way, I find the eigenvalues D1 and D2 coming from :

[V1,D1] = eig(FISH_sp);
[V2,D2] = eig(FISH_xc);

% Check eigen values : OK, columns of eigenvalues D1 found !
FISH_sp*V1./V1
% Check eigen values : OK, columns of eigenvalues D2 found !
FISH_xc*V2./V2

How could I fix this wrong result (I am talking about the ratios :

FISH_sp*phi./phi
FISH_xc*phi./phi

which don't give same values for a given column of FISH_sp and FISH_xc) ) ?? In the paper, they say that phi diagonalizes A=FISH_sp and B=FISH_xc but I can't reproduce it.

If someone could see where is my error ...

Answer 1

You defined barA = inv(phiB_bar)*FISH_sp*phiB_bar. From Eq. (39) in the manuscript it looks like it should be barA = transpose(phiB_bar)*FISH_sp*phiB_bar instead.

Also, your mehtod 1 fails when B is singular (an inverse does not exist). MATLAB's eig(A,B) however should handle also singular Bs if memory serves me well.