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

Answer: 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:  Get eigenvectors/values for A and B respectively.  Group the resultant eigenvectors by their eigenspaces.  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.