diff --git a/src/Quaternion.jl b/src/Quaternion.jl
index b5e6f303..9afc5881 100644
--- a/src/Quaternion.jl
+++ b/src/Quaternion.jl
@@ -293,3 +293,26 @@ function slerp(qa::Quaternion{T}, qb::Quaternion{T}, t::T) where {T}
         qa.v3 * ratio_a + qm.v3 * ratio_b,
     )
 end
+
+function _complex_representation(A::Matrix{Quaternion{T}}) where {T}
+    # convert a quatenion matrix to corresponding complex matrix
+    n = size(A, 1)
+    Ac = Matrix{Complex{T}}(undef, 2n, 2n)
+    Ac[1:n, 1:n] .= complex.(getfield.(A, :s), getfield.(A, :v1))
+    Ac[1:n, n+1:2n] .= complex.(getfield.(A, :v2), getfield.(A, :v3))
+    Ac[n+1:2n, 1:n] .= .- conj.(view(Ac, 1:n, n+1:2n))
+    Ac[n+1:2n, n+1:2n] .= conj.(view(Ac, 1:n, 1:n))
+    return Ac
+end
+
+function _quaternion_representation(A::Matrix{Complex{T}}) where {T}
+    n = round(Int, size(A, 1)/2)
+	
+    X = @view A[1:n, 1:n]
+    Y = @view A[1:n, n+1:2n]
+
+    return [Quaternion(X[i,j].re, X[i,j].im, Y[i,j].re, Y[i,j].im) for i in 1:n, j in 1:n]
+end
+
+# A quick way of doing the quaternionic matrix exponential
+exp(A::Matrix{Quaternion{T}}) where {T} = _quaternion_representation(exp(_complex_representation(A)))
diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl
index d0635f2e..6a9b630c 100644
--- a/test/test_Quaternion.jl
+++ b/test/test_Quaternion.jl
@@ -138,4 +138,13 @@ for _ in 1:100
         @test q ⊗ slerp(q1, q2, t) ≈ slerp(q ⊗ q1, q ⊗ q2, t)
         @test q ⊗ linpol(q1, q2, t) ≈ linpol(q ⊗ q1, q ⊗ q2, t)
     end
+    @testset "exp(::Matrix{Quaternion})" begin  # test matrix exponential
+        z = zeros(Quaternion{Float64}, 3, 3)
+        id = Matrix{Quaternion{Float64}}(I, 3, 3)
+        d = [Quaternion(randn(4)...) for i in 1:3]
+        expd = exp.(d)
+        D = exp(Array(Diagonal(d)))
+        @test exp(z) ≈ id
+        @test D ≈ Array(Diagonal(expd))
+    end
 end