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.


I still can't answer question 1, but I can offer a fix for question 2 by offering a completely different proof of the theorem (this is essentially the same as here, but this was the only reference I could find.

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).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.