« Generalized simultaneous diagonalization » of a pair of symmetric, non-commuting, positive semi-definite matrices I hope my question is trivial for some of you but for the time being I’m lost somewhere between the generalized eigenproblem, simultaneous diagonalization of quadratic forms, simultaneous SVD, generalized SVD, etc., none of them matching my problem.
Let $A$  and $B$  be two symmetric, positive semi-definite (but not positive-definite) matrices in  ${\mathbb{R}^{n \times n}}$. $\left[ {A,B} \right] \ne 0$ . Both of them are diagonalizable, one of them is diagonal. Find a pair of non-singular matrices $P$  and $Q$  such as
$\left\{ \begin{gathered}
  PAQ = {D_1} \hfill \\
  PBQ = {D_2} \hfill \\ 
\end{gathered}  \right.$
${D_1}$  and ${D_2}$  diagonal. No other property of  $P$  and $Q$   is required.
 A: From the theoretical point of view: there is a Jordan-like canonical form for pairs of Hermitian matrices $(A,B)$ under the equivalence relation $(A,B) \sim (PAP^T, PBP^T)$ (with a nonsingular $P$), see https://www.sciencedirect.com/science/article/pii/0024379576900215 . This is known today in my area as "(a variant of) the even Kronecker canonical form"; see e.g. my paper https://doi.org/10.1137/120861679 which contains a self-contained statement for Hermitian-anti-Hermitian (which is the same up to multiplication by $i$ of one of them).
It is easy to check that only a few of these blocks can be used to produce positive semidefinite $A$ and $B$: if I am not wrong, if $A,B$ are spd then you can only get in the canonical form size-1 blocks corresponding to purely real eigenvalues (E2), size-1 blocks corresponding to infinite eigenvalues (E3), and size-1 singular blocks (E4, corresponding to common kernel of the two matrices). All these blocks are size-1: in particular, this means that this decomposition will produce diagonal $PAP^T,PBP^T$, and hence it is the decomposition that you want. I.e., the form that you ask for exists for all pairs $(A,B)$, and you can obtain it even if you restrict to transformations with $Q=P^T$ (which is nice, since they preserve symmetry and definiteness).
Numerically, you can first identify the common kernel of $A$ and $B$ and then change basis to put it in evidence, to reduce the problem to
$$
\begin{bmatrix}A_{11} & 0\\ 0 & 0\end{bmatrix}, \begin{bmatrix}B_{11} & 0\\ 0 & 0\end{bmatrix}
$$ 
and then apply a generalized eigensolver (e.g., Matlab's eig(A11, B11)) to $A_{11}$ and $B_{11}$, which is now a nonsingular problem (i.e., $\det(A+Bz)$ is not identically zero) in view of the previous discussion.
