Simultaneous diagonalization of two positive semi-definite matrices Let matrices $A, B$ be two $n \times n $ positive semi-definite matrices ; they can be represented in the following form 
 $$A=\sum_{i=1}^{n} \psi_{i}p_{i}p_{i}^{T}=P\Psi P^{T}, \quad B=\sum_{i=1}^{n} \phi_{i}q_{i}q_{i}^{T}=Q \Phi Q^{T} $$
where $\{p_{i}, \cdots, p_{n} \}$ and $\{q_{i}, \cdots, q_{n} \}$ are standard orthogonal sets in $\mathcal{R}^{n}$. There always exist a matrix $U$  diagonalizing simultaneously $A,B$, i.e.,
$$U^{T}AU=\Lambda_{1}, \quad U^{T}BU=\Lambda_{2}$$,
where $\Lambda_{1}, \Lambda_{2}$ are diagonal matrices.  My question is how to find the matrices $\Lambda_{1}, \Lambda_{2}$ form the eigenvalues and eigenvectors of $A, B$, that is to say, find the relationship between $(U, \Lambda_{1}, \Lambda_{2})$ and $(P,\Psi, Q, \Phi)$.
 A: There is an easy proof when at least one of the two matrices is positive definite, linked to in the other answer, but it does not work when both matrices are singular. (You can't factorise a singular matrix as a product of non-singular matrices!)
Here is a proof that does not assume positive definiteness. Let $V = \{ v \in \mathbb{R}^n : v^\dagger A v = 0 \}$. By the Cauchy–Schwarz inequality, $V = \ker A$; in particular, $V$ is a vector subspace. Let $U$ be any vector subspace satisfying the following conditions:

*

*$U \cap V = \{ 0 \}$.

*$U + V = \mathbb{R}^n$.

*If $u \in U$ and $v \in V$, then $u^\dagger B v = 0$.

Such a $U$ exists: let $W = \{ u \in \mathbb{R}^n : \forall v \in V . u^\dagger B v = 0 \}$, choose a basis of $V \cap W$, extend it to a basis of $W$, and then discard the basis of $V \cap W$ to obtain a basis of $U$; by construction, $U + V = V + W$, and to see that $V + W = \mathbb{R}^n$, choose a basis of $V$ that is orthogonal with respect to $B$ and use the usual projection formula. (The basis vectors of $V$ that are null with respect to $B$ may be ignored.)
Now, observe that $A$ defines a positive definite quadratic form on $U$, so we may find a basis of $U$ that is orthonormal with respect to $A$ and orthogonal with respect to $B$. On the other hand, $A$ is null on $V$, and there is a basis of $V$ that is orthogonal with respect to $B$. Putting these together yields a basis of $\mathbb{R}^n$ diagonalising $A$ and $B$ simultaneously.
Here is an outline algorithm to compute the required basis numerically:

*

*Obtain a basis $M_V$ of $V$ by using e.g. SVD or eigendecomposition to find the null space of $A$.

*Obtain a basis $M_W$ of $W$ by using e.g. SVD to find the null space of $M_V^\dagger B$.

*Obtain a basis $M_U$ of $U$ by using e.g. QR decomposition of $\begin{pmatrix} M_V & M_W \end{pmatrix}$; alternatively, if $M_V$ is an orthonormal basis of $V$, then one can use e.g. SVD to find an orthonormal basis $M_K$ of $(I - M_V M_V^\dagger) M_W$ then use SVD or QR decomposition to find the column space of $M_W (I - M_K M_K^\dagger)$.

*Obtain an invertible matrix $R$ such that $R^\dagger M_U^\dagger A M_U R = I$ by using e.g. Cholesky decomposition of $M_U^\dagger A M_U$.

*Obtain a unitary matrix $Q$ such that $Q^\dagger R^\dagger M_U^\dagger B M_U R Q$ is a diagonal matrix by using e.g. eigendecomposition of $R^\dagger M_U^\dagger B M_U R$.

*Obtain an invertible matrix $P$ such that $P^\dagger M_V^\dagger B M_V P$ is a diagonal matrix by using e.g. eigendecomposition of $M_V^\dagger B M_V$.

The required basis is then $M = \begin{pmatrix} M_U R Q &  M_V P \end{pmatrix}$. Indeed, by construction:
$$M^\dagger A M = \begin{pmatrix} I & 0 \\ 0 & 0 \end{pmatrix}$$
$$M^\dagger B M = \begin{pmatrix} Q^\dagger R^\dagger M_U^\dagger B M_U R Q & 0 \\ 0 & P^\dagger M_V^\dagger B M_V P \end{pmatrix}$$
I would be interested to hear if there is a simpler / more efficient algorithm...
A: See this: A property of positive definite matrices
I believe a possible singularity does not change much the proof for positive definite matrices.
