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.



Diagonalization of quaternion matrices
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>. 


