I've come up with the following piece of Python code (using the library Sympy):



```
def double_diagonalize(m1, m2):
    V, _ = (m1.T * m2).diagonalize()
    U, _ = (m1 * m2.T).diagonalize()
    return U, V

```

What I've found is that given many (but not all) random pairs of $
n \times n$ matrices $M\_+$ and $M\_-$, it produces $U$ and $V$ such that both $U^{-1} M\_+ V^{-T}$ and $U^T M\_- V$ evaluate to diagonal matrices.


My question is, is this result known? And by what name? It looks similar to [Generalised Schur Decomposition](https://en.wikipedia.org/wiki/Schur_decomposition#Generalized_Schur_decomposition) except that the matrices $U$ and $V$ are not necessarily unitary, and the (pair of) "normal forms" $U^{-1} M\_+ V^{-T}$ and $U^T M\_- V$ are usually fully diagonal instead of merely upper triangular.


It might be related to this as well: <https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix#Generalized_eigenvalue_problem>



Double-diagonalisation of nxn matrices?
It's the solution of the [product eigenvalue problem](https://doi.org/10.1137/S0036144504443110) $M\_1^T M\_2$, once you transpose the second relation.


