4

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 :

generalized problem

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

other formulation

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.

  1. 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) :

description of generalized eigenvectors and eigenvalues

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

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

  • I only quickly read through half of the question, but I can already see an error in your definition `phiB_bar = V2*(D2.^(0.5)+1e-10*eye(7)^(-1));` in your second method: from the formula from the paper that should be `phiB_bar = V2*inv(D2.^(0.5)+1e-10*eye(7));`, shouldn't it? – Matteo V Jan 28 '21 at 08:44
  • @MatteoV . Thanks for this suggestion, you are right : I have replace by `phiB_bar = V2*(D2.^(0.5)+1e-10*eye(7))^(-1);` But the `phi` doesn't still appear as a eigenvector of both `A=FISH_sp` and `B=FISH_xp`. It is frustrating... –  Jan 28 '21 at 11:21
  • Just did a quick test with two 5x5 random matrices, and the method 1 seems to be correct (`norm(A*phi - B*phi*lambda) = 1.933967216188714e-14` in my case), while for method 2 it returns `norm(A*phi - B*phi*lambda) = 2.228744851244737`. Maybe there are some underlying assumptions in the paper (disclaimer: I haven't read it)? – Matteo V Jan 28 '21 at 11:56
  • `[V, D] = eig(A, B)` should solve the Generalized Eigenvalue Problem `A*V = B*V*D`. Is this not working for your problem? For questions about the paper math.stackexchange.com may be the better place to ask. – pH 74 Jan 28 '21 at 12:08
  • @pH74 yes, I can use directly the computation with `[V, D] = eig(A, B)` but I would like to do the link with a common eigenevctors basis between A and B. One suggested me to look a the problem of **joint diagonalization**, and it has two variants: orthogonal, in which the basis vectors are orthonormal, and non-orthogonal, which is harder to solve, but which may be more appropriate to your application. You can see this advise on : https://scicomp.stackexchange.com/questions/36670/matlab-compute-approximative-common-eigenvectors-basis-between-two-matrices-as . –  Jan 28 '21 at 13:02
  • But I don't know how to make the link between this joint diagonalization and the generalized eigenvalues problem. –  Jan 28 '21 at 13:02

1 Answers1

0

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.

Luk
  • 36
  • 1
  • this error has too little impact on the constraints but glad that you found it, it's already that fixed. About the matrix B, I handle only with B not singular since I use often the inverse of B (which is equal to the covariance of a probe represented by B), so no worries about that. Best regards –  Feb 02 '21 at 15:10