# Matlab - Compute approximative common eigenvectors basis between two matrices as a function of tolerance

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.

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, I mean for a tolerance as small as possible ?

Are there alternative methods or routines to perform this minimization of commutator combined with the vector $$v$$, that is ‖(𝐴𝐵−𝐵𝐴)𝑣‖ ?

I saw there is another function called rref which can accept a tolerance parameter but :

1. What's the difference with singular values decomposition svd
2. Which criterion could I apply for a pertinent choice of tolerance value

The 2 matrices to find approximative common eigen vectors matrix are available here :

A matrix

B matrix

Anyone could try to apply a function Matlab appropriate to find a basis of common eigen vectors or write a small Matlab script for this ? Even approximative basis would be enough, everything depends of the tolerance that I am ready to accept but currently I don't know how to introduce this tolerance parameter with SVD algorithm.

UPDATE 1: among different methods that I have tried to use, anyone could explain the method of Pool Variance ? If it can be easy to implement. I remember that it is consisted by taking the half of each diagonalised Fisher matrices and sum them to come back after into final Covariance space, i.e by just applying : Cov = P Fisher_diagonal_sum P^-1 ?

With this method, I get interesting results, but from a theorical point of view, impossible for me to justify the principle of this Pool variance matrix (that is, by taking the half of diagonal Fisher matrices values and come back to Covariance) : why the half ?

Any suggestion/track/clue help is welcome

• @FedericoPoloni. Sorry, I just fix them. Regards Jan 16 at 11:21
• You are writing a lot, but it's a bit confuse. I suggest you first define what you mean by "approximate computation of common eigenvectors"? Approximate in which measure, i.e. what is tooptimize? The link you give to Wikipedia considers the generalized eigenvalue problem, which is a different thing. So maybe beforehand you should first state what you really want -- is it to solve the generalized Eigenvalue problem? For this, you don't need common eigenvectors, and also no approximations of them. Jan 21 at 14:25
• @davidhigh . ok, by talking about an approximate computation of common eigenvectors, I mean that given the fact that my 2 matrices don't commute, there is no analytical solution. So that's why I would like an approximative solution (depending of the tolerance I want). If I manage to get this approximative basis, I can do the sum of 2 eigen values of each matrix in the same basis (Fisher matrices synthesis and I know already these eigen values). That's the main goal : find an approximative with a tolerance fixed (that can be adapted of course). Regards Jan 22 at 7:17
• Ok, and approximate in which sense? Take the Euclidian basis vectors, they are an approximate eigenbasis for any two non-commuting matrices. Not the right approximation? Why not? Jan 22 at 18:07
• @davidhigh By approximative, I mean that the eigenvectors of first matrix (passing matrix) looks like to the eigenvectors of second matrix . You will tell me : up to which point ? that's the problem when I say "looks like" or "approximative" or even "similar". Only a tolerance factor will allow me to choose or not to consider the 2 matrices as "similar" or "approximative similar". But I can't see for the moment how to introduce this "tolerance" factor in my code and in the building of this pseudo common eigenvectors matrix, I say "pseudo" since there are no analytical solutions. Jan 22 at 18:40

This problem is known as 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.

The simplest method I know of seeks a unitary matrix $$U$$ that minimizes the sum of squares of the off-diagonal elements of $$U^HAU$$ and $$U^HBU$$ (the $$H$$ superscript denotes the Hermitian transpose). The matrix $$U$$ is composed as a product of Givens matrices just like in the Jacobi diagonalization algorithm. Implementations in several languages, including matlab, can be found here: http://www2.iap.fr/users/cardoso/jointdiag.html

If you relax the constraint that the diagonalizating matrix is unitary, you may be able to reduce the cost function further, though the optimization becomes harder. See for example: https://www.eng.tau.ac.il/~arie/Files/tsp02.pdf.

• Hello @AmitHochman. Thank you for your answer. I don't know if I should consider a unitary matrix U or not. Which criterion should indicate me the right choice to make ? As a starting point, I thought that, for combining the 2 matrices (I mean cross-correlations), I should add each eigen values of the 2 Fisher matrix : this summing only makes sense if we are in the same eigenvectors basis when doing this summation. Up to this, it seems correct. But if inverse again the summed diagonal matrix to get back in individual parameters basis, constraints given by fiinal Fisher matrix are not good. Jan 28 at 19:40
• I wonder if you think that my approach is correct. Jan 28 at 19:40
• I would try the orthogonal methods first. Only if you can’t find an orthonormal basis that meets your requirements try the more involved non-orthogonal methods. Moreover, since your matrices are symmetric, an orthonormal basis makes sense. Also, you wrote in math.stackexchange.com/questions/3922971/… that you want an orthogonal matrix. Jan 28 at 22:12

As its documentation suggests in multiple places, rref is "mainly of academic interest" (read: "used only to explain Gaussian elimination to undergrads"), and is not a serious competitor of SVD-based algorithms in terms of stability. I recommend against it.

• @youpilat13 Nullspace detection with the SVD comes naturally with a tolerance parameter. But the choice of this tolerance is application-dependent: in general, the residual that you can expect in svd(A*B-B*A) will depend on the magnitude of the error that you have in your input matrices. Jan 16 at 11:40
• @youpilat13 Ah, the XY problem at play again. You should have asked "how do I find a common basis of eigenvectors, for a pair of matrices that has one". That's a different problem. Jan 16 at 16:52
• Anyhow, I don't know of any careful perturbation theory, but in practice a quick and dirty way to find a common basis of eigenvectors is eig(x * A + y * B), where x,y are random scalars. But that should be a separate question I guess. :) Jan 18 at 9:14
• @youpilat13 No, I mean the V output of [V, D] = eig(x * A + y * B) is your basis of common eigenvectors, for almost all choices of x,y. No need to do anything more than that one-liner. Jan 18 at 12:46
• @youpilat13 Please ask a new question at this point... :) Jan 18 at 13:04