I'd like to go between $n \times m$ matrices with quaternionic entries to (using, it would seem, the Pauli matrix conversion scheme) to $2 n \times 2 m$ matrices with complex entries, and vice versa. I'd appreciate some efficient coding to accomplish this pair of tasks.
Say we have an $n \times m$ matrix the $n m$ entries of which are quaternions. We want to generate a $2 n \times 2 m$ matrix through the rule
Quaternion[a_, b_, c_, d_] -> (a IdentityMatrix - b PauliMatrix - c PauliMatrix - d PauliMatrix)
So, it seems that we have to appropriately reshape the resultant array (using, presumably,
ArrayReshape) to the desired $2 n \times 2 m$ form.
Then, what off-hand seems more challenging, we want to be able to perform the inverse operation (going from the $2 n \times 2 m$ form to the [quaternionic-entry] $n \times m$ form]). Such an inverse operation would assume that the matrix to be transformed back possesses the indicated special (Pauli matrix) structure.
As some background to this problem, what I am really interested in accomplishing is finding the appropriate way to generate random quaternionic $4 \times 4$ density matrices with respect to Hilbert-Schmidt measure. The manner of accomplishing this for random complex $4 \times 4$ density matrices and random real $4 \times 4$ density matrices is presented in
(see eq. (1) and bot. p. 9 there). It seems that I may have to convert the $4 \times 4$ quaternionic-entry matrices to $8 \times 8$ form and use the developed rule for complex-entry matrices.
Some background references to the general quaternionic/Pauli matrix topic are (the amazingly detailed)
which has a very explicit (eqs. (2)-(5)) straightforward treatment pertaining to the issue raised.
As well as wikipedia,