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 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