How to find the common eigenvectors of two matrices with distincts eigenvalues

Question

I am looking for finding or rather building common eigenvectors matrix X between 2 matrices A and B such as :

AX=aX with "a" the diagonal matrix corresponding to the eigenvalues

BX=bX with "b" the diagonal matrix corresponding to the eigenvalues

where A and B are square and diagonalizable matrices.

I took a look in a similar post but had not managed to conclude, i.e having valid results when I build the final wanted endomorphism F defined by : F = P D P^-1

I have also read the wikipedia topic and this interesting paper but couldn't have to extract methods pretty easy to implement.

Particularly, I am interested by the eig(A,B) Matlab function.

I tried to use it like this :

% Search for common build eigen vectors between FISH_sp and FISH_xc
[V,D] = eig(FISH_sp,FISH_xc);
% Diagonalize the matrix (A B^-1) to compute Lambda since we have AX=Lambda B X
[eigenv, eigen_final] = eig(inv(FISH_xc)*FISH_sp);
% Compute the final endomorphism : F = P D P^-1
FISH_final = V*eye(7).*eigen_final*inv(V)

But the matrix FISH_final don't give good results since I can do other computations from this matrix FISH_final (this is actually a Fisher matrix) and the results of these computations are not valid.

So surely, I must have done an error in my code snippet above. In a first time, I prefer to conclude in Matlab as if it was a prototype, and after if it works, look for doing this synthesis with MKL or with Python functions. Hence also tagging python.

How can I build these common eigenvectors and finding also the eigenvalues associated? I am a little lost between all the potential methods that exist to carry it out.

The screen capture below shows that the kernel of commutator has to be different from null vector :

EDIT 1: From maths exchange, one advices to use Singular values Decomposition (SVD) on the commutator [A,B], that is in Matlab doing by :

"If 𝑣 is a common eigenvector, then ‖(𝐴𝐵−𝐵𝐴)𝑣‖=0. The SVD approach gives you a unit-vector 𝑣 that minimizes ‖(𝐴𝐵−𝐵𝐴)𝑣‖ (with the constraint that ‖𝑣‖=1)"

So I extract the approximative eigen vectors V from :

[U,S,V] = svd(A*B-B*A)

Is there a way to increase the accuracy to minimize ‖(𝐴𝐵−𝐵𝐴)𝑣‖ as much as possible ?

IMPORTANT REMARK : Maybe some of you didn't fully understand my goal.

Concerning the common basis of eigen vectors, I am looking for a combination (vectorial or matricial) of V1 and V2, or directly using null operator on the 2 input Fisher marices, to build this new basis "P" in which, with others eigenvalues than known D1 and D2 (noted D1a and D2a), we could have :

F = P (D1a+D2a) P^-1

To compute the new Fisher matrix F, I need to know P, assuming that D1a and D2a are equal respectively to D1 and D2 diagonal matrices (coming from diagonalization of A and B matrices)

If I know common basis of eigen vectors P, I could deduce D1a and Da2 from D1 and D2, couldn't I ?

The 2 Fisher matrices are available on these links :

Matrix A

Matrix B

Selfsustaining

Have you looked at the documentation of eig? I don't see where eig(A,B) gives common eigenvectors/values of A and B. – Danyluk

Also I would ask Python and Matlab questions separately. They are different questions after all. – Danyluk

Argyll . Thanks for your answer. If you look at fr.mathworks.com/help/matlab/ref/eig.html#d122e338300 , it is indicated that [V,D] = eig(A,B) returns diagonal matrix D of generalized eigenvalues and full matrix V whose columns are the corresponding right eigenvectors, so that A*V = B*V*D. How to make the link to find, between the 2 matrices, a common eigen vectors basis from this relation A*V = B*V*D ? Regards – Selfsustaining

Trying to remove the cobweb on my linear algebra, I don't think there is a link between A*V = B*V*D and common eigenvectors? I can write an answer highlighting Matlab's eigenvector methods. The answer will work on small matrices; otherwise I do not wish to devise an efficient algorithm on the spot. As far as I can tell, there is no standard numerical method to find common eigenvectors. If you are trying to understand Matlab, perhaps what I suggest would help. Otherwise, I think either you need to provide an algorithm or this question needs to be asked on a different stack exchange site. – Danyluk

So what is your situation? Do you understand how to get common eigenvectors conceptually and just looking to understand how to implement it on Matlab? Are your matrices supposed to be all diagonalizable? – Danyluk

@Danyluk I am looking for building common eigenvector on 2 small matrices (7x7), so your answer will be precious for me. Regards – Selfsustaining

Yes my matrices are diagonalizable ! – Selfsustaining

@Danyluk .and also, I tried a numerical solving for building a common eigenvectors under the form eigenv_sp*a+eigenv_xc*b with "a" and "b" 2 matrices to found (I mean all their elements). You can see this attempts on fr.mathworks.com/matlabcentral/answers/… – Selfsustaining

Let us continue this discussion in chat. – Selfsustaining

Do you need to build a basis of common eigenvectors (7 in your case) or just find as many common eigenvectors as you can? – Bedraggled

@Bedraggled . I need to build a basis of common eigenvectors (7 in my case) even if eigen values will change. Regards – Selfsustaining

Two things: 1) I don't think your indices are correct for either case. 2) m is not necessarily a single column. It can be empty, single column, or multiple columns. So V(:,i) = m should not work. You need to horizontally concatenate. For example, initialize with V=[]; then in each iteration, aggregate result with V=[V,m]. – Danyluk

I'm not sure but you may need to reword your question: "I have two square matrices (say I have two eigenvector matrices) and I want to find all pairs of columns that are linearly dependent". – Abb

@Abb . I don't understand well your comment : the main goal is to find or to build common eigen vectors for both square matrices, that is, forming a new basis in which I have : A = P D_a P^-1 and B = P D_b P^-1, where P is the passing matrix representing this common eigen vectors basis and D_a, D_b the eigen values: could you see clearer ? – Selfsustaining

In my understanding use [~,VA] = eig(A) and [~,VB] = eig(B) to find eigenvectors. If some of columns of VA are linearly dependent on some columns of VB they are common eigenvectors. – Abb

Answer 1

I don't think there is a built-in facility in Matlab for computing common eigenvalues of two matrices. I'll just outline brute force way and do it in Matlab in order to highlight some of its eigenvector related methods. We will assume the matrices A and B are square and diagonalizable.

Outline of steps:

1. Get eigenvectors/values for A and B respectively.

2. Group the resultant eigenvectors by their eigenspaces.

3. Check for intersection of the eigenspaces by checking linear dependency among the eigenvectors of A and B one pair eigenspaces at a time.

Matlab does provide methods for (efficiently) completing each step! Except of course step 3 involves checking linear dependency many many times, which in turn means we are likely doing unnecessary computation. Not to mention, finding common eigenvectors may not require finding all eigenvectors. So this is not meant to be a general numerical recipe.

How to get eigenvector/values

The syntax is

[V,D] = eig(A)

where D(i), V(:,i) are the corresponding eigenpairs.

Just be wary of numerical errors. In other words, if you check

tol=sum(abs(A*V(:,i)-D(i)*V(:,i)));

tol<n*eps should be true for some small n for a smallish matrix A but it's probably not true for 0 or 1.

Example:

>> A = gallery('lehmer',4);
>> [V,D] = eig(A);
>> sum(abs(A*V(:,1)-D(1)*V(:,1)))<eps
ans =
  logical
   0
>> sum(abs(A*V(:,1)-D(1)*V(:,1)))<10*eps
ans =
  logical
   1

How to group eigenvectors by their eigenspaces

In Matlab, eigenvalues are not automatically sorted in the output of [V,D] = eig(A). So you need to do that.

• Get diagonal entries of matrix: diag(D)

• Sort and keep track of the required permutation for sorting: [d,I]=sort(diag(D))

• Identify repeating elements in d: [~,ia,~]=unique(d,'stable')

ia(i) tells you the beginning index of the ith eigenspace. So you can expect d(ia(i):ia(i+1)-1) to be identical eigenvalues and thus the eigenvectors belonging to the ith eigenspace are the columns W(:,ia(i):ia(i+1)-1) where W=V(:,I). Of course, for the last one, the index is ia(end):end

The last step happens to be answered here in true generality. Here, unique is sufficient at least for small A.

(Feel free to ask a separate question on how to do this whole step of "shuffling columns of one matrix based on another diagonal matrix" efficiently. There are probably other efficient methods using built-in Matlab functions.)

For example,

>> A=[1,2,0;1,2,2;3,6,1];
>> [V,D] = eig(A),
V =
         0         0    0.7071
    1.0000   -0.7071         0
         0    0.7071   -0.7071
D =
     3     0     0
     0     5     0
     0     0     3
>> [d,I]=sort(diag(D));
>> W=V(:,I),
W =
         0    0.7071         0
    1.0000         0   -0.7071
         0   -0.7071    0.7071
>> [~,ia,~]=unique(d,'stable'),
ia =
     1
     3

which makes sense because the 1st eigenspace is the one with eigenvalue 3 comprising of span of column 1 and 2 of W, and similarly for the 2nd space.

How to get linear intersect of (the span of) two sets

To complete the task of finding common eigenvectors, you do the above for both A and B. Next, for each pair of eigenspaces, you check for linear dependency. If there is linear dependency, the linear intersect is an answer.

There are a number of ways for checking linear dependency. One is to use other people's tools. Example: https://www.mathworks.com/matlabcentral/fileexchange/32060-intersection-of-linear-subspaces

One is to get the RREF of the matrix formed by concatenating the column vectors column-wise.

Let's say you did the computation in step 2 and arrived at V1, D1, d1, W1, ia1 for A and V2, D2, d2, W2, ia2 for B. You need to do

for i=1:numel(ia1)
    for j=1:numel(ia2)
         check_linear_dependency(col1,col2);
    end
end

where col1 is W1(:,ia1(i):ia1(i+1)-1) as mentioned in step 2 but with the caveat for the last space and similarly for col2 and by check_linear_dependency we mean the followings. First we get RREF:

[R,p] = rref([col1,col2]);

You are looking for, first, rank([col1,col2])<size([col1,col2],2). If you have computed rref anyway, you already have the rank. You can check the Matlab documentation for details. You will need to profile your code for selecting the more efficient method. I shall refrain from guess-estimating what Matlab does in rank(). Although whether doing rank() implies doing the work in rref can make a good separate question.

In cases where rank([col1,col2])<size([col1,col2],2) is true, some rows don't have leading 1s and I believe p will help you trace back to which columns are dependent on which other columns. And you can build the intersect from here. As usual, be alert of numerical errors getting in the way of == statements. We are getting to the point of a different question -- ie. how to get linear intersect from rref() in Matlab, so I am going to leave it here.

There is yet another way using fundamental theorem of linear algebra (*sigh at that unfortunate naming):

null( [null(col1.').' ; null(col2.').'] )

The formula I got from here. I think ftla is why it should work. If that's not why or if you want to be sure that the formula works (which you probably should), please ask a separate question. Just beware that purely math questions should go on a different stackexchange site.

---

Now I guess we are done!

---

EDIT 1:

Let's be extra clear with how ia works with an example. Let's say we named everything with a trailing 1 for A and 2 for B. We need

for i=1:numel(ia1)
    for j=1:numel(ia2)
        if i==numel(ia1)
            col1 = W1(:,ia1(end):end);
        else
            col1 = W1(:,ia1(i):ia1(i+1)-1);
        end
        if j==numel(ia2)
            col2 = W2(:,ia2(j):ia2(j+1)-1);
        else
            col2 = W2(:,ia2(end):end);
        end
        check_linear_dependency(col1,col2);
    end
end

---

EDIT 2:

I should mention the observation that common eigenvectors should be those in the nullspace of the commutator. Thus, perhaps null(A*B-B*A) yields the same result.

But still be alert of numerical errors. With the brute force method, we started with eigenpairs with low tol (see definition in earlier sections) and so we already verified the "eigen" part in the eigenvectors. With null(A*B-B*A), the same should be done as well.

Of course, with multiple methods at hand, it's good idea to compare results across methods.

Danyluk

Thanks a lot for your detailled answer, I am going to implement your method ! – Selfsustaining

@youpilat13: Good luck! – Danyluk

@youpilat13: btw, aren't common eigenvectors exactly those in the nullspace of the commutator? In other words, in Matlab, null(A*B-B*A). May be worth comparing the brute force result with something like that. – Danyluk

Hi, coud you take a look in my EDIT 1 please at my script to see how to build the common eigenvectors ? (the eigen values will also change I guess easily). Regards – Selfsustaining

@youpilat13: They look right. check_linear_dependency() is meant to be a custom function for you to build. If the use of RREF doesn't immediately come to use, try the last null( [null(col1.').' ; null(col2.').'] ) method first. – Danyluk

@youpilat13: And don't forget the last index is different. So you need to do something with if i==numel(ia1) and if j==numel(ia2) – Danyluk

It is pretty difficult for me to understand : have I to put null( [null(col1.').' ; null(col2.').'] ) inside the double loop ? this would be clearer if you modified my script that I put. Thanks – Selfsustaining

the instruction null(A*B-B*A) returns a column vector -0.0085 -0.0048 -0.2098 0.9776 -0.0089 -0.0026 0.0109 . Does it mean that it is possible to find common eigen vectors ? – Selfsustaining

How to circumvent the issue about the last index is different. So you need to do something with if i==numel(ia1) and if j==numel(ia2) ? I tried by setting num(ia1)-1 and num(ia2)-1 in limits of double loop, but this wat, I have only 6x6 matrices and not 7x7 as expected. – Selfsustaining

@youpilat13: Yes, null( [null(col1.').' ; null(col2.').'], according to the link, can be literally how you implement check_linear_dependency(). You do need to account for the last index issue. You want to think through what ia means. (And I know it's more about getting used to Matlab.) ia(i) is the beginning index for the ith eigenspace. So it follows that ia(i+1) is the beginning index for the i+1th eigenspace. Thus, ia(i+1)-1 is the ending index for the ith eigenspace. However, that is only true until you hit the last eigenspace. – Danyluk

@youpilat13: You do not want to discard last eigenspace. You need to instead split into: if it is the last eigenspace, col1=something; else col1=something else. – Danyluk

@youpilat13: re what you got from null(A*B-B*A) in your example. It appears to me what I said to you about commutator is true. But you need to further examine the mathematics part yourself or on math stackexchange. Next, do visit the documentation on null(). If you get a single column vector, then nullspace of the commutator is 1-dimensional as far as Matlab can tell. With that result, you can think about how to verify that vector is indeed an eigenvector of A and B down to a small tolerance. – Danyluk

Thanks for your help. You can see in my EDIT 2 my attempt to apply what you suggested but I get an error. If you could see what's wrong, I would be grateful. – Selfsustaining

I have a new version of the code snippet in EDIT3 concerning the intersection of eigen spaces. The results gives elements with zero values, I am opened to any suggestions/fix/clue/track is welcome. By the way, I wonder if a double loop is necessary and not only a single ? Regards – Selfsustaining

@youpilat13: You get the error because your conditional expression is incorrect. Revisit what needs to be done: the last eigenspace for A needs to be read differently; the last eigenspace for B needs to be read differently. So it is one if block for each of col1 and col2. And it is in the definition of col1 and col2. m and V don't need to be in the if block. With your effort as the jumping point, I can put what I expect that will work in my answer. Hope that works for you. – Danyluk

@youpilat13: As for double loop vs single loop, if you mean a single loop that will be equivalent as the double I showed, yes. Otherwise, no. The brute force idea is about checking intersect of each pair of eigenspaces. For each i,j combination, the eigenspaces in questions are distinct. You wouldn't want to miss any combination. – Danyluk

isn't there a little error in your code snippet ? i.e col2 = W2(:,ia2(i):ia2(j+1)-1)instead of col2 = W2(:,ia2(i):ia2(i+1)-1) ? – Selfsustaining

@youpilat13: Yes! – Danyluk

Thanks ! as you saw, I have launched a bounty to draw attention on this computation of common basis of eigen vectors. Just 2 questions : 1) given the fact I get always V = 7x0 empty double matrix returned, I would like how to include inside the double loop the null(A*B-B*A) which returns the colum vector (here row for convenience) = [ 0.0085, -0.0048, -0.2098, 0.9776, -0.0089, -0.0026, 0.0109 ] but I would need to introduce i and j indices. – Selfsustaining

2) How to include the tolerance parameter tol inside this double loop ? I would like to include this : % Diagonalize [V1,D1] = eig(FISH_sp); [V2,D2] = eig(FISH_xc); V1V2 = V1'*V2; [j k] = find(abs(abs(V1V2)-1)<tol); % show that the jth A eigenvector and kth B eigenvector are proportional V1(:,j(1))./V2(:,k(1)) ; V1(:,j(2))./V2(:,k(2)) with for example a tol = 0.1 . Regards – Selfsustaining

Could you take a look please at my EDIT 1 ? I would be grateful. I did tests with the research of common eigen vector(s) but results are strange with null(AB-BA) and I don't understand why. Regards – Selfsustaining

@youpilat13: I would make your questions focused and ask them one at a time. When the question shifts to "how to verify a vector is eigenvector numerically", the question is that. It has shifted away from the original question. And you need to be aware of the adverse consequence for the people who try to answer your question, who in turn need to process the materials and consider the details -- they are not your real life advisors obviously; but also those who may get into similar situation and can benefit from your question and its answers. De-focused questions won't be helpful in the future. – Danyluk

@youpilat13: And taking away potential benefit in question/answers by not focusing them erodes away some motivation for others to review the question to begin with. That's something I think you need to keep in mind when engaging in SO questions/answers. – Danyluk

Could you post please a simple Python script that summarizes all the different parts of you method mentionned in your answer. Sorry, I have difficulties to do a concise code from all of what you indicate. I would be grateful if you could do this. Regards – Selfsustaining

Answer 2

I suspect this is rather a delicate matter.

First off, mathematically, A and B are simultaneously diagonalisable iff they commute, that is iff

A*B - B*A == 0 (often A*B-B*A is written [A,B])
(for if A*X = X*a and B*X = X*b with a, b diagonal then
A = X*a*inv(X), B = X*b*inv(X) 
[A,B] = X*[a,b]*inv(X) = 0 since a and b, being diagonal, commute)

So I'd say the first thing to check is that your A and B do commute, and here is the first awkward issue: since [A,B] as computed is unlikely to be all zeroes due to rounding error, you'll need to decide if [A,B] being non-zero is just due to rounding error or if, actually, A and B don't commute.

Now suppose x is an eigenvector of A, with eigenvalue e. Then

A*B*x = B*A*x = B*e*x = e*B*x

And so we have, mathematically, two possibilities: either Bx is 0, or Bx is also an eigenvector of A with eigenvalue e.

A nice case is when all the elements of a are different, that is when each eigenspace of A is one dimensional. In that case: if AX = Xa for diagonal a, then BX = Xb for diagonal b (which you'll need to compute). If you diagonalize A, and all the eigenvalues are sufficiently different, then you can assume each eigenspace is of dimension 1, but what does 'sufficiently' mean? Another delicate question, alas. If two computed eigenvalues are very close, are the eigenvalues different or is the difference rounding error? Anyway, to compute the eigenvalues of b for each eigenvector x of A compute Bx. If ||Bx|| is small enough compared to ||x|| then the eigenvalue of B is 0, otherwise it's

x'*B*x/(x'*x) 

In the general case, some of the eigenspaces may have dimension greater than 1. The one dimensional eigen spaces can be dealt with as above, but more computation is required for the higher dimensional ones.

Suppose that m eigenvectors x[1].. x[m] of A correspond to the eigenvalue e. Since A and B commute, it is easy to see that B maps the space spanned by the xs to itself. Let C be the mxm matrix

C[i,j] = x[j]'*B*x[i]

Then C is symmetric and so we can diagonalize it, ie find orthogonal V and diagonal c with

C = V'*c*V

If we define

y[l] = Sum{ k | V[l,k]*x[k] } l=1..m

then a little algebra shows that y[l] is an eigenvector of B, with eigenvalue c[l]. Moreover, since each x[i] is an eigenvector of A with the same eigenvalue e, each y[l] is also an eigenvector of A with eigenvector e.

So all in all, I think a method would be:

1. Compute [A,B] and if its really not 0, give up
2. Diagonalise A, and sort the eigenvalues to be increasing (and sort the eigenvectors!)
3. Identify the eigenspaces of A. For the 1 dimensional spaces the corresponding eigenvector of A is an eigenvector of B, and all you need to compute is the eigenvalue of B. For higher dimensional ones, proceed as in the previous paragraph.

A relatively expensive (in computational effort) but reasonably reliable way to test whether the commutator is zero would be to compute the svd of the commutator and take is largest singular value, c say, and also to take the largest singular value (or largest absolute value of the eigenvalues) a of A and b of B. Unless c is a lot smaller (eg 1e-10 times) the lesser of a and b, you should conclude the commutator is not zero.

Bedraggled

Two matrices don't have to be simultaneously diagonalizable in order to have some common eigenvectors. It's only true that the common ones are in the nullspace of the commutator. I pointed that out to the OP below my question. As usual, more involved math discussions would happen on a different site. – Danyluk

@Argyll, true but in the comments tp the question the OP answered my comment by saying he did require and eigenbasis. – Bedraggled

Oh ok. One of the links he provided points to only examining common eigenvectors. – Danyluk

@Bedraggled . If I do : >> null(FISH_sp*FISH_xc-FISH_xc*FISH_sp) , I get the result : ans = -0.0085 -0.0048 -0.2098 0.9776 -0.0089 -0.0026 0.0109 what could I conclude about this ? How can I build 7 columns vectors like this to form a passing matrix P representing all the common basis of eigenvectors ? Regards – Selfsustaining

@youpilat13 Its a bit difficult to say from this. Is that saying that the null space of the commutator is one dimensional? If so then that is evidence that the commutator is not zero and that therefore you cannot find a basis of common eigenvectors. Another check would be to look at the singular values (from the svd) of the commutator and also of the original matrices. Unless the commutators singular values are very much smaller than the singular values of the originals, it is likely that your problem has no solution – Bedraggled

@Bedraggled . I think that I may introduce a tolerance to consider the commutator [A,B] near to 0, what do you think about this ? – Selfsustaining

@youpilat13 I've added a final paragraph to my answer – Bedraggled

@Bedraggled you can see the results in my EDIT 2. what should I conclude ? – Selfsustaining

@Bedraggled . Could you write please a simple Python code that makes the synthesis of all the different parts of you method mentionned in your answer. Sorry, I have difficulties to do a brief code from all of what you indicate. I would be grateful if you could do this. Regards – Selfsustaining

@youpilat13 Alas no, for I do not know python – Bedraggled

@Bedraggled . Or maybe you know Matlab ? – Selfsustaining

Answer 3

Finding approximate common eigenvector of two matrices is known as finding common principal components (e.g. B. Flury). They can be found by either optimizing loglikelihood or least squares cost function. In particular, for the above two matrices, the use of LS is a better choice as Matrix B is problematic.

Cholent