The formal way to simultaneously diagonalise two Hermitian operators that commute? Given two Hermitian operators $A$ and $B$, such that $[A,B] = 0$, if one [or both] of these operators are degenerate, how does one define a formal way of going about simultaneously diagonalising both of them? Preferably, if possible, with some example of the method.
Previous chat discussions:

*

*I asked this question on the h-bar a while back, and while I got a good response, I think the nature of brevity of the chat meant that when I got to thinking about it more, I realised I still wasn't entirely sure of the whole thing, plus I figured it would be nice to put this question on the site as I personally couldn't find it when I needed it.


*Discussed in the hbar chat room here.


*Previous conversation in chat about the subject.
 A: For each pair (eigenspace for $A$, eigenspace for $B$), pick a basis for their intersection. (In particular, if an intersection is 0-dimensional, it will not have/contribute with any basis vectors.) Finally concatenate all the bases.
A: The simplest method is to first diagonalise $A$. Then consider in turn each eigenvalue $\lambda$ and a basis of the associated eigenspace $\mathcal{E}_\lambda$: $\newcommand{\ket}[1]{|{#1}\rangle}(\ket{1},\ket{2},\cdots,\ket{n})$. You then construct the matrix $\newcommand{\bra}[1]{\langle{#1}|}M_B=(\bra{i}B\ket{j})_{1\le i,j \le n}$. There is one such matrix for each eigenvalue $\lambda$, just to be crystal clear, but I will refrain from tagging a $\lambda$ to $M_B$ to keep the notation readable. Finally you diagonalise $M_B$. This will gives you eigenvalues $\mu_1, \cdots, \mu_n$ (not necessarily distinct) and eigenvectors $u_i$, which are column vectors, 
$$u_i=\begin{pmatrix} u_{i1}\\ \vdots \\u_{in} \end{pmatrix},$$
such that $M_Bu_i=\mu_iu_i$. Then you construct $\ket{i'}=\sum_j u_{ij}\ket{j}$, and now $(\ket{1'}, \cdots, \ket{n'})$ are eigenvectors for both $A$, all for the eigenvalue $\lambda$, and for $B$, for the respective eigenvalues $\mu_1,\cdots,\mu_n$. 
Of course, if $\lambda$ was not degenerate in the first place, i.e. $n=1$, then there is nothing to do! If it is degenerate, most often, $n$ will be small, at least way smaller than the dimension of the eigenproblem for $A$, so the diagonalisation of $M_B$ will be comparatively easy. Then, again, don't forget you have to do this for each eigenvalue $\lambda$ of $A$. Of course, you could start by diagonalising $B$ instead: just do what looks simpler.
There are way more efficient numerical methods note, which would make a big difference for large matrices, but for your bread and butter quantum system, the method I highlighted should be tractable.
A: I suppose that the Hilbert space $H$ is finite dimensional
and I indicate by $H_\lambda$ the eigenspace of $A$ with eigenvalue $\lambda$ and some dimension $d_\lambda\geq 1$. 
You know that $$H_\lambda \perp H_{\lambda'}\quad \mbox{if $\lambda\neq \lambda'$}\tag{1}$$ and $$\oplus_\lambda H_\lambda= H\tag{2}\:.$$
The fundamental idea is that every $H_\lambda$ is invariant under the action of $B$, i.e.,  $$B(H_\lambda) \subset H_\lambda\:.$$ This is  because $A(Bx)= BAx= B\lambda x= \lambda (Bx)$ if $x\in H_\lambda$. 
As a consequence you can 
(a) restrict $B$ to $H_\lambda$  and, noticing that $B|_{H_\lambda} : H_\lambda \to H_\lambda$ is still Hermitian as you easily prove,
(b)  find an orthonormal basis of $H_\lambda$ made of  eigenvectors $\{x^{(\lambda)}_n\}_{n=1,\ldots d_\lambda}$ of $B$ with corresponding eigenvalues $\mu^{(\lambda)}_n$.
Varying both


*

*$\lambda$ in the set of eigenvalues of $A$ and 

*$n=1, \ldots , d_\lambda$, 


due to (1) and (2) the set of all unit mutually orthogonal vectors $x^{(\lambda)}_n$ form an orthonormal basis of the entire  $H$. 
This basis is made of simultaneous eigenvectors of $A$ and $B$ because $$Ax^{(\lambda)}_n =\lambda x^{(\lambda)}_n$$ and $$Bx^{(\lambda)}_n= \mu^{(\lambda)}_n x^{(\lambda)}_n\:.$$
Obviously it may happen that $\mu^{(\lambda)}_n =\mu^{(\lambda')}_m$ for some $\lambda \neq \lambda'$.
A: The next simplest method (particularly useful if you only need numerical result) is to consider the operator 
$$
{\cal O}=\alpha A+\beta B
$$
with $\alpha$ and $\beta$ chosen so that ${\cal O}$ has no repeated eigenvalues.  Then the eigenvectors of ${\cal O}$ are simultaneous eigenvectors of $A$ and $B$.
