Diagonalization of quaternion matrices I would like to diagonalize a very large matrix that has the same property as quaternion matrices (in a sense that the matrix can be written as a linear combination of quaternions):
$$
H = \left(\begin{array}{cc} H_{11} & H_{12} \\ -H^\ast_{12} & H^\ast_{11} \end{array}\right),
\quad U^\dagger HU = \epsilon
$$
where $\epsilon$ is a diagonal matrix, $U$ is a unitary matrix and $U$ has the same symmetry
$$
U = \left(\begin{array}{cc} U_{11} & U_{12} \\ -U^\ast_{12} & U^\ast_{11} \end{array}\right)
$$
Note that there are $dim(\epsilon)/2$ sets of identical eivenvalues.
Of course it is trivial to diagonalize $H$ in a brute-force way to get right eivenvalues; however, for some reasons I need a symmetry-adapted set of eivenvectors as shown above. Is there a way to transform the matrix obtained by brute-force diagonalization to an symmetry-adapted matrix above? Although I am aware of some papers such as (J. Math. Phys. 46, 052106 (2005); http://dx.doi.org/10.1063/1.1896386), which suggest
$$
\left(\begin{array}{cc} V_{11} & V_{12} \\ V_{21} & V_{11} \end{array}\right)\to 
\left(\begin{array}{cc} U_{11} & U_{12} \\ -U^\ast_{12} & U^\ast_{11} \end{array}\right)
$$
where $V^\dagger H V = \epsilon$ and
$$
U_{11} = \frac{1}{2}(V_{11}+V_{22}^\ast),\quad U_{12}=\frac{1}{2}(V_{12}-V_{21}^\ast)
$$
it does not seem to preserve orthogonality of the eigenvectors.
 A: I am rather certain that you cannot modify a standard diagonalization to the structured diagonalization that you want.  I have seen such a method used, but it is expected to fail occasionally, especially for a large matrix.  What you expect are the eigenvalues to all have even multiplicity, and then you modify the basis of each eigenspace so that if [v;w] is an eigenvector, you also use [-conj(w); conj(v)].  (Pardon the Matlab notation.)  When numerical errors are involved, you will find it hard to figure which eigenvectors to group into even-multiplicity eigenspaces.
It seems your matrix is also Hermitian.  I suggest that you can use HAPACK, Fortran Software for (Skew-)Hamiltonian Eigenvalue Problems, available from http://www.tu-chemnitz.de/mathematik/hapack.  For a dense 10,000-by-10,000 matrix, this may work (give time and much RAM), as it takes advantage of one symmetry, skew-Hamiltonian in CS parlance.  However, the code probably ignores the fact that your matrix is Hermitian.  Note that Hermitian plus skew-Hamiltonian implies the quaternionic structure you have.
You can also use the Page / Van Loan algorithm to get to a block diagonal form and then deal with those separately.  This is discussed in "Topological insulators and C∗-algebras: Theory and numerical practice" by myself and Hastings, Annals of Physics Volume 326, Issue 7, July 2011, Pages 1699–1759.  You can grab some Matlab code associated to the paper "Computing a logarithm of a unitary matrix with general spectrum" Arxiv:1203.6151.  A newer version of the this paper exists and will soon be published, but the latest Matlab code is at http://hdl.handle.net/1928/23450.  
A: Sangwine and Le Bihan have developed an algorithm, implemented in MATLAB, that can be used for the diagonalization (or singular value decomposition) of matrices with quaternion elements.
A Jacobi algorithm is discussed in Computing the SVD of a quaternion matrix. The MATLAB toolbox is here.
