# Finding the symplectic matrix in Williamson's theorem

tl;dr: How do I construct the symplectic matrix in Williamson's theorem?

I am interested in a constructive proof/version of Williamson's theorem in symplectic linear algebra. Maybe I'm just missing a simple step, so here is what I know:

Let us fix the symplectic form $J=\begin{pmatrix} 0_n & 1_n \\ -1_n & 0_n\end{pmatrix}$. Recall:

Theorem: Let $M\in\mathbb{R}^{2n\times 2n}$ be a positive-definite matrix, then there exists a symplectic matrix $S\in Sp(2n,\mathbb{R})$ and a diagonal matrix $D\in\mathbb{R}^{n\times n}$ such that $$M=S^T\tilde{D}S$$ where $\tilde{D}=\operatorname{diag}(D,D)$ is diagonal.

My basic interest is in how to construct the $S$ (i.e. write a matlab program that does this). This theorem is cited and used in various cases, however one cannot find many proofs despite the original. Here is the one proof I found (the most important step for me is unjustified).

Let us denote by $\langle \cdot,\cdot \rangle_M:=\langle \cdot, M\cdot \rangle$ the symmetric bilinear form defined by $M$ and by $\sigma(\cdot,\cdot)=\langle \cdot, J\cdot \rangle$ the symplectic form. The theorem asserts that there is an $M$-orthogonal and symplectic basis. The failure of the basis to be $M$-orthonormal is then given by $D$.

In order to find such a basis, it seems a good idea to look at $JM$. As this is a real matrix, which is antisymmetric with respect to $\langle\cdot,\cdot\rangle_M$, its eigenvalues come in pairs $\pm i\lambda_j$ ($j=1,\ldots n$) and the corresponding eigenvectors in pairs $e_j\pm if_j$ with real $e_j,f_j$ (the proof in the link considers the eigenvalues/eigenvectors of $M^{-1}J$ instead).

The claim in the proof is now that $\{e_j,f_j\}_j$ forms an $M$-orthonormal basis.

If this is the case, one can easily see that $\delta_{jk}=\langle e_j,Me_k\rangle=-\lambda_k\sigma(e_j,f_k)$ and $0=\langle e_j,Mf_k\rangle=\lambda_k\sigma(e_j,e_k)$ and similarly for $f_k$, hence the basis is indeed $M$-orthogonal and symplectic after normalization by $\sqrt{\lambda_j}$. The matrix $S$ should then consist of the normalized $f_j$ and $e_j$ as columns.

Q1: Why is $\{e_j,f_j\}_j$ an orthonormal basis w.r.t to $M$? Matlab suggests that most of the times, it is, but I seem to miss something fundamental, because I don't see it. Why most of the times? Eigenvalue multiplicities seem to play a role here:

Q2: If $\lambda_j\neq \lambda_k$ always, then Q1 seems to be true. If there are some eigenvalue multiplicities, this seems to be wrong, probably because of the non-uniqueness of eigenvectors. How can I fix this? I suspect that some Gram-Schmidt wrt. $M$ is necessary, but since I can't answer Q1, I can't see whether this does the trick.

Without loss of generality, we can write $S=M^{-1/2}K\overline{D}^{1/2}$ with $K\in O(2n)$, since then $S^TMS=\overline{D}$ by construction. We already know how to find $\overline{D}$ (spectrum of $JM$), so if we can find $K$ such that $S$ is symplectic, we are done.
Write: $$S^TJS=J \quad \Leftrightarrow \quad K^TM^{-1/2}JM^{-1/2}K=\overline{D}^{-1/2}J\overline{D}^{-1/2}$$
The right hand side is of the form $$\overline{D}^{-1/2}J\overline{D}^{-1/2}=\begin{pmatrix}{} 0_n & D^{-1} \\ -D^{-1} & 0_n \end{pmatrix}$$ and $M^{-1/2}JM^{-1/2}$ is antisymmetric, too. Up to a change of basis, we see that the right hand side is the real diagonalization of the antisymmetric left hand side (as in the spectral theory of skew-symmetric matrices).
This means we can find $K$ by finding the real diagonalization of $M^{-1/2}JM^{-1/2}$ in the other basis, which can be done e.g. by the (real) Schur decomposition (which coincidentally computes $D$ as well).