We consider the quaternion linear system \(AX = B\) for the unknown matrix X, where A, B are given \(n\times n\) , \(n\times s\) matrices with quaternion entries, motivated by applications that arise from fields such as quantum mechanics and signal processing. Our primary concern is the large-scale setting when n is large so that direct solutions are not feasible. We describe a block Krylov subspace method for the iterative solution of these quaternion linear systems. One difference compared to usual block Krylov subspace methods over complex Euclidean spaces is that the multiplication of quaternion scalars is not commutative. We describe a block quaternion Arnoldi process, taking noncommutativity features of quaternions into account, to generate an orthonormal basis for the quaternion Krylov space \(\text {blockspan} \{ R_0, A R_0, \dots , A^k R_0 \}\) , where \(R_0 = B - A X_0\) and \(X_0\) is an initial guess for the solution. Then the best solution of \(AX = B\) in the least-squares sense is sought in the generated Krylov space. We explain how these least-squares problems over quaternion Krylov spaces can be solved efficiently by means of Householder reflectors. Most notably, we analyze rigorously the convergence of the proposed block quaternion GMRES approach when A is diagonalizable, and in the more general setting when A is not necessarily diagonalizable by making use of the Jordan form of A. Finally, we report numerical results that confirm the validity of the deduced theoretical convergence results, in particular illustrate that the proposed block quaternion Krylov subspace method converges quickly when A has clustered eigenvalues.
