diff --git a/Project.toml b/Project.toml
index 44da9a80cc..42d39cabf2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.8.19"
+version = "0.8.20"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ManifoldsBase = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MatrixEquations = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
+Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
@@ -31,9 +32,10 @@ Einsum = "0.4"
 Graphs = "1.4"
 HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
-MatrixEquations = "2.2"
 ManifoldsBase = "0.13.13"
+MatrixEquations = "2.2"
 Plots = "1"
+Quaternions = "0.5"
 RecipesBase = "1.1"
 RecursiveArrayTools = "2"
 Requires = "0.5, 1"
diff --git a/docs/src/manifolds/generalunitary.md b/docs/src/manifolds/generalunitary.md
index 6dcdfb5bd4..d640a47551 100644
--- a/docs/src/manifolds/generalunitary.md
+++ b/docs/src/manifolds/generalunitary.md
@@ -7,7 +7,7 @@ Both [`OrthogonalMatrices`](@ref) and [`UnitaryMatrices`](@ref) are quite simila
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/Orthogonal.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ## Unitary Matrices
@@ -15,13 +15,14 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/Unitary.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
+
 ## [Common functions](@id generalunitarymatrices)
 
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["manifolds/GeneralUnitaryMatrices.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
diff --git a/docs/src/manifolds/group.md b/docs/src/manifolds/group.md
index a35e28880f..7cbe4f08bf 100644
--- a/docs/src/manifolds/group.md
+++ b/docs/src/manifolds/group.md
@@ -40,7 +40,7 @@ As a concrete wrapper for manifolds (e.g. when the manifold per se is a group ma
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/GroupManifold.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Generic Operations
@@ -52,7 +52,7 @@ For groups based on an addition operation or a group operation, several default
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/addition_operation.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 #### Multiplication Operation
@@ -60,7 +60,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/multiplication_operation.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Circle group
@@ -68,7 +68,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/circle_group.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### General linear group
@@ -76,7 +76,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/general_linear.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Heisenberg group
@@ -84,7 +84,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/heisenberg.jl"]
-Order = [:type,:function]
+Order = [:constant, :type, :function]
 ```
 
 ### (Special) Orthogonal and (Special) Unitary group
@@ -97,7 +97,7 @@ many common functions, these are also implemented on a common level.
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/general_unitary_groups.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 #### Orthogonal group
@@ -105,7 +105,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/orthogonal.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 #### Special orthogonal group
@@ -113,7 +113,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/special_orthogonal.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 #### Special unitary group
@@ -121,7 +121,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/special_unitary.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 #### Unitary group
@@ -129,7 +129,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/unitary.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Power group
@@ -137,7 +137,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/power_group.jl"]
-Order = [:type, :function, :constant]
+Order = [:constant, :type, :function]
 ```
 
 ### Product group
@@ -145,7 +145,7 @@ Order = [:type, :function, :constant]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/product_group.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Semidirect product group
@@ -153,7 +153,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/semidirect_product_group.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Special Euclidean group
@@ -161,7 +161,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/special_euclidean.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Special linear group
@@ -169,7 +169,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/special_linear.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ### Translation group
@@ -177,7 +177,7 @@ Order = [:type, :function]
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["groups/translation_group.jl"]
-Order = [:type, :function]
+Order = [:constant, :type, :function]
 ```
 
 ## Group actions
diff --git a/docs/src/misc/notation.md b/docs/src/misc/notation.md
index c343c0355d..c0020f5bb1 100644
--- a/docs/src/misc/notation.md
+++ b/docs/src/misc/notation.md
@@ -27,7 +27,7 @@ Within the documented functions, the utf8 symbols are used whenever possible, as
 | ``\gamma`` | a geodesic | ``\gamma_{p;q}``, ``\gamma_{p,X}`` | connecting two points ``p,q`` or starting in ``p`` with velocity ``X``. |
 | ``∇ f(p)`` | gradient of function ``f \colon \mathcal{M} \to \mathbb{R}`` at ``p \in \mathcal{M}`` | | |
 | ``\circ`` | a group operation | |
-| ``\cdot^\mathrm{H}`` | Hermitian or conjugate transposed| |
+| ``\cdot^\mathrm{H}`` | Hermitian or conjugate transposed for both complex or quaternion matrices| |
 | ``e`` | identity element of a group | |
 | ``I_k`` | identity matrix of size ``k\times k`` | |
 | ``k`` | indices | ``i,j`` | |
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 4f572c5be9..2b2ab1e923 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -259,6 +259,7 @@ using ManifoldsBase:
     trait
 using Markdown: @doc_str
 using MatrixEquations: lyapc
+using Quaternions
 using Random
 using RecipesBase
 using RecipesBase: @recipe, @series
diff --git a/src/groups/unitary.jl b/src/groups/unitary.jl
index 50806c47be..761e7667ea 100644
--- a/src/groups/unitary.jl
+++ b/src/groups/unitary.jl
@@ -1,32 +1,35 @@
 @doc raw"""
-     Unitary{n,𝔽} = GeneralUnitaryMultiplicationGroup{n,ℂ,AbsoluteDeterminantOneMatrices}
+     Unitary{n,𝔽} = GeneralUnitaryMultiplicationGroup{n,𝔽,AbsoluteDeterminantOneMatrices}
 
-The group of unitary matrices ``\mathrm{U}(n)``.
+The group of unitary matrices ``\mathrm{U}(n, 𝔽)``, either complex (when 𝔽=ℂ) or quaternionic
+(when 𝔽=ℍ)
 
-The group consists of all points ``p ∈ \mathbb C^{n × n}`` where ``p^\mathrm{H}p = pp^\mathrm{H} = I``.
+The group consists of all points ``p ∈ 𝔽^{n × n}`` where ``p^\mathrm{H}p = pp^\mathrm{H} = I``.
 
 The tangent spaces are if the form
 
 ```math
-T_p\mathrm{U}(x) = \bigl\{ X \in \mathbb C^{n×n} \big| X = pY \text{ where } Y = -Y^{\mathrm{H}} \bigr\}
+T_p\mathrm{U}(n) = \bigl\{ X \in 𝔽^{n×n} \big| X = pY \text{ where } Y = -Y^{\mathrm{H}} \bigr\}
 ```
 
 and we represent tangent vectors by just storing the [`SkewHermitianMatrices`](@ref) ``Y``,
-or in other words we represent the tangent spaces employing the Lie algebra ``\mathfrak{u}(n)``.
+or in other words we represent the tangent spaces employing the Lie algebra ``\mathfrak{u}(n, 𝔽)``.
+
+Quaternionic unitary group is isomorphic to the compact symplectic group of the same dimension.
 
 # Constructor
 
-    Unitary(n)
+    Unitary(n, 𝔽::AbstractNumbers=ℂ)
 
-Construct ``\mathrm{U}(n)``.
+Construct ``\mathrm{U}(n, 𝔽)``.
 See also [`Orthogonal(n)`](@ref) for the real-valued case.
 """
-const Unitary{n} = GeneralUnitaryMultiplicationGroup{n,ℂ,AbsoluteDeterminantOneMatrices}
+const Unitary{n,𝔽} = GeneralUnitaryMultiplicationGroup{n,𝔽,AbsoluteDeterminantOneMatrices}
 
-Unitary(n) = Unitary{n}(UnitaryMatrices(n))
+Unitary(n, 𝔽::AbstractNumbers=ℂ) = Unitary{n,𝔽}(UnitaryMatrices(n, 𝔽))
 
 @doc raw"""
-    exp_lie(G::Unitary{2}, X)
+    exp_lie(G::Unitary{2,ℂ}, X)
 
 Compute the group exponential map on the [`Unitary(2)`](@ref) group, which is
 
@@ -35,15 +38,19 @@ Compute the group exponential map on the [`Unitary(2)`](@ref) group, which is
 ```
 
 where ``θ = \frac{1}{2} \sqrt{4\det(X) - \operatorname{tr}(X)^2}``.
- """
-exp_lie(::Unitary{2}, X)
+"""
+exp_lie(::Unitary{2,ℂ}, X)
+
+function exp_lie(::Unitary{1,ℍ}, X::Number)
+    return exp(X)
+end
 
 function exp_lie!(::Unitary{1}, q, X)
-    q[1] = exp(X[1])
+    q[] = exp(X[])
     return q
 end
 
-function exp_lie!(::Unitary{2}, q, X)
+function exp_lie!(::Unitary{2,ℂ}, q, X)
     size(X) === (2, 2) && size(q) === (2, 2) || throw(DomainError())
     @inbounds a, d = imag(X[1, 1]), imag(X[2, 2])
     @inbounds b = (X[2, 1] - X[1, 2]') / 2
@@ -69,7 +76,11 @@ function exp_lie!(G::Unitary, q, X)
 end
 
 function log_lie!(::Unitary{1}, X, p)
-    X[1] = log(p[1])
+    X[] = log(p[])
+    return X
+end
+function log_lie!(::Unitary{1}, X::AbstractMatrix, p::AbstractMatrix)
+    X[] = log(p[])
     return X
 end
 function log_lie!(G::Unitary, X, p)
@@ -78,6 +89,13 @@ function log_lie!(G::Unitary, X, p)
     return X
 end
 
+identity_element(::Unitary{1,ℍ}) = Quaternion(1.0)
+
+function log_lie(::Unitary{1}, q::Number)
+    return log(q)
+end
+
 Base.inv(::Unitary, p) = adjoint(p)
 
-show(io::IO, ::Unitary{n}) where {n} = print(io, "Unitary($(n))")
+show(io::IO, ::Unitary{n,ℂ}) where {n} = print(io, "Unitary($(n))")
+show(io::IO, ::Unitary{n,ℍ}) where {n} = print(io, "Unitary($(n), ℍ)")
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index 5ab50843ba..7b6a851fe0 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -170,6 +170,8 @@ function default_estimation_method(
     return GeodesicInterpolationWithinRadius(π / 2 / √2)
 end
 
+embed(::GeneralUnitaryMatrices, p) = p
+
 @doc raw"""
     embed(M::GeneralUnitaryMatrices{n,𝔽}, p, X)
 
diff --git a/src/manifolds/Unitary.jl b/src/manifolds/Unitary.jl
index 61b8759895..d8fd294984 100644
--- a/src/manifolds/Unitary.jl
+++ b/src/manifolds/Unitary.jl
@@ -1,8 +1,9 @@
 
 @doc raw"""
-    const UnitaryMatrices{n} = AbstarctUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices}
+    const UnitaryMatrices{n,𝔽} = AbstarctUnitaryMatrices{n,𝔽,AbsoluteDeterminantOneMatrices}
 
-The manifold ``U(n)`` of ``n×n`` complex matrices such that
+The manifold ``U(n,𝔽)`` of ``n×n`` complex matrices (when 𝔽=ℂ) or quaternionic matrices
+(when 𝔽=ℍ) such that
 
 ``p^{\mathrm{H}}p = \mathrm{I}_n,``
 
@@ -17,15 +18,37 @@ The tangent spaces are given by
     \bigr\}
 ```
 
-But note that tangent vectors are represented in the Lie algebra, i.e. just using ``Y`` in the representation above.
+But note that tangent vectors are represented in the Lie algebra, i.e. just using ``Y`` in
+the representation above.
 
 # Constructor
-     UnitaryMatrices(n)
+    
+    UnitaryMatrices(n, 𝔽::AbstractNumbers=ℂ)
 
 see also [`OrthogonalMatrices`](@ref) for the real valued case.
 """
-const UnitaryMatrices{n} = GeneralUnitaryMatrices{n,ℂ,AbsoluteDeterminantOneMatrices}
+const UnitaryMatrices{n,𝔽} = GeneralUnitaryMatrices{n,𝔽,AbsoluteDeterminantOneMatrices}
 
-UnitaryMatrices(n::Int) = UnitaryMatrices{n}()
+UnitaryMatrices(n::Int, 𝔽::AbstractNumbers=ℂ) = UnitaryMatrices{n,𝔽}()
 
-show(io::IO, ::UnitaryMatrices{n}) where {n} = print(io, "UnitaryMatrices($(n))")
+check_size(::UnitaryMatrices{1,ℍ}, p::Number) = nothing
+check_size(::UnitaryMatrices{1,ℍ}, p, X::Number) = nothing
+
+embed(::UnitaryMatrices{1,ℍ}, p::Number) = SMatrix{1,1}(p)
+
+embed(::UnitaryMatrices{1,ℍ}, p, X::Number) = SMatrix{1,1}(X)
+
+function exp(::UnitaryMatrices{1,ℍ}, p, X::Number)
+    return p * exp(X)
+end
+
+function log(::UnitaryMatrices{1,ℍ}, p::Number, q::Number)
+    return log(conj(p) * q)
+end
+
+project(::UnitaryMatrices{1,ℍ}, p) = normalize(p)
+
+project(::UnitaryMatrices{1,ℍ}, p, X) = (X - conj(X)) / 2
+
+show(io::IO, ::UnitaryMatrices{n,ℂ}) where {n} = print(io, "UnitaryMatrices($(n))")
+show(io::IO, ::UnitaryMatrices{n,ℍ}) where {n} = print(io, "UnitaryMatrices($(n), ℍ)")
diff --git a/test/groups/general_unitary_groups.jl b/test/groups/general_unitary_groups.jl
index a19cef73d4..8002f14317 100644
--- a/test/groups/general_unitary_groups.jl
+++ b/test/groups/general_unitary_groups.jl
@@ -42,6 +42,44 @@ include("group_utils.jl")
         end
     end
 
+    @testset "Quaternionic Unitary Group" begin
+        QU1 = Unitary(1, ℍ)
+        @test repr(QU1) == "Unitary(1, ℍ)"
+
+        @test identity_element(QU1) === Quaternion(1.0)
+
+        p = QuaternionF64(
+            0.4815296357756736,
+            0.6041613272484806,
+            -0.2322369798903669,
+            0.5909181717450419,
+            true,
+        )
+        X = Quaternion(0.0, 0, 0, 1)
+        q = exp(QU1, p, X)
+        X2 = log(QU1, p, q)
+        @test isapprox(QU1, p, X, X2)
+        q2 = exp_lie(QU1, X)
+        X3 = log_lie(QU1, q2)
+        @test isapprox(QU1, p, X, X3)
+
+        q3 = Ref(Quaternion(0.0))
+        exp_lie!(QU1, q3, X)
+        @test isapprox(QU1, q3[], q2)
+        X4 = Ref(Quaternion(0.0))
+        log_lie!(QU1, X4, q2)
+        @test isapprox(QU1, identity_element(QU1), X3, X4[])
+
+        X5 = fill(Quaternion(0.0), 1, 1)
+        log_lie!(QU1, X5, fill(q2, 1, 1))
+        @test isapprox(QU1, identity_element(QU1), X3, X5[])
+
+        @test inv(QU1, p) == conj(p)
+
+        @test project(QU1, p, Quaternion(1.0, 2.0, 3.0, 4.0)) ===
+              Quaternion(0.0, 2.0, 3.0, 4.0)
+    end
+
     @testset "Special Unitary Group" begin
         SU2 = SpecialUnitary(2)
         @test repr(SU2) == "SpecialUnitary(2)"
diff --git a/test/manifolds/symplecticstiefel.jl b/test/manifolds/symplecticstiefel.jl
index ba48f8fb5e..36f01dbaec 100644
--- a/test/manifolds/symplecticstiefel.jl
+++ b/test/manifolds/symplecticstiefel.jl
@@ -22,16 +22,14 @@ end
     @testset "Real" begin
         SpSt_6_4 = SymplecticStiefel(2 * 3, 2 * 2)
 
-        p_6_4 = Array{Float64}(
-            [
-                0 0 -5 -1
-                0 0 9 -2
-                0 0 -2 1
-                -2 -9 -3 6
-                -3 -13 -21 9
-                -8 -36 18 -6
-            ],
-        )
+        p_6_4 = Array{Float64}([
+            0 0 -5 -1
+            0 0 9 -2
+            0 0 -2 1
+            -2 -9 -3 6
+            -3 -13 -21 9
+            -8 -36 18 -6
+        ])
         q_6_4 = Array{Float64}([
             0 0 -3 1
             0 0 8 -3
diff --git a/test/manifolds/unitary_matrices.jl b/test/manifolds/unitary_matrices.jl
index c65ec2bbdd..25c920e180 100644
--- a/test/manifolds/unitary_matrices.jl
+++ b/test/manifolds/unitary_matrices.jl
@@ -1,5 +1,7 @@
 include("../utils.jl")
 
+using Quaternions
+
 @testset "Orthogonal Matrices" begin
     M = OrthogonalMatrices(3)
     @test repr(M) == "OrthogonalMatrices(3)"
@@ -39,3 +41,42 @@ end
     X2 = log(M, p, r)
     @test isapprox(M, p, X, X2)
 end
+
+@testset "Quaternionic Unitary Matrices" begin
+    M = UnitaryMatrices(1, ℍ)
+    @test repr(M) == "UnitaryMatrices(1, ℍ)"
+
+    # wrong length of size
+    @test_throws DomainError is_point(M, zeros(2, 2), true)
+
+    # Determinant not one
+    pF2 = [1im 1.0; 0.0 -1im]
+    @test_throws DomainError is_point(M, pF2, true)
+    p = QuaternionF64(
+        0.4815296357756736,
+        0.6041613272484806,
+        -0.2322369798903669,
+        0.5909181717450419,
+        true,
+    )
+
+    @test is_point(M, fill(p, 1, 1))
+    @test is_point(M, p)
+
+    @test_throws DomainError is_vector(M, p, zeros(2, 2), true)
+    # not skew
+    @test_throws DomainError is_vector(M, p, Quaternion(1, 0, 0, 0), true)
+    X = Quaternion(0.0, 0, 0, 1)
+    @test is_vector(M, p, X)
+
+    pu = QuaternionF64(
+        -0.2178344173900564,
+        -0.4072721993877449,
+        -2.2484219560115712,
+        -0.4718862793239344,
+        false,
+    )
+    q = project(M, pu)
+    @test is_point(M, q)
+    @test isapprox(q, normalize(pu))
+end