[spacedome]

I found this approach to work well for linear systems, here is some rough code I used for the representation map (note the jmag/kmag functions are part of my implementation, not sure what the equivalent is for Quaternions.jl).

  function cmatrix(Q::AbstractVecOrMat{Quaternion{T}}) where {T}
      [complex.( real.(Q),  imag.(Q)) complex.( jmag.(Q), kmag.(Q));
       complex.(-jmag.(Q),  kmag.(Q)) complex.( real.(Q), -imag.(Q))]
  end

  function qmatrix(C::AbstractMatrix{Complex{T}}) where {T}
      n, m = size(C)
      quat.(C[1:n÷2, 1:m÷2], C[1:n÷2, m÷2+1:m])
  end