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

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.

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

Answer 1

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.

Answer 2

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.