# Double-diagonalisation of nxn matrices?

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

It's the solution of the product eigenvalue problem $$M_1^T M_2$$, once you transpose the second relation.