From 6ac1e839dd4183d1f6480f14467158d5a02471e3 Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 00:15:03 +0200
Subject: [PATCH 1/6] propose the fix for #522.

---
 src/manifolds/GeneralUnitaryMatrices.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index bce4582e0e..9521671c1a 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -612,7 +612,7 @@ Return the dimension of the manifold orthogonal matrices and of the manifold of
 manifold_dimension(::GeneralUnitaryMatrices{n,ℝ}) where {n} = div(n * (n - 1), 2)
 @doc raw"""
     manifold_dimension(M::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices})
-    
+
 Return the dimension of the manifold of special unitary matrices.
 ```math
 \dim_{\mathrm{SU}(n)} = n^2-1.
@@ -666,16 +666,17 @@ end
 
 Orthogonally project the tangent vector ``X ∈ 𝔽^{n × n}``, ``\mathbb F ∈ \{\mathbb R, \mathbb C\}``
 to the tangent space of `M` at `p`,
-and change the representer to use the corresponding Lie algebra, i.e. we compute
+where `X` is already assumed to be in a Lie algebra represemtation, that is,
+the projection is perfomed onto the Lie algebra. We compute
 
 ```math
-    \operatorname{proj}_p(X) = \frac{pX-(pX)^{\mathrm{T}}}{2},
+    \operatorname{proj}_p(X) = \frac{X-X^{\mathrm{T}}}{2},
 ```
 """
 project(::GeneralUnitaryMatrices, p, X)
 
 function project!(::GeneralUnitaryMatrices{n,𝔽}, Y, p, X) where {n,𝔽}
-    project!(SkewHermitianMatrices(n, 𝔽), Y, p \ X)
+    project!(SkewHermitianMatrices(n, 𝔽), Y, X)
     return Y
 end
 

From 6c873351307ae11f693f02367cb77204c7b0631f Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 07:40:59 +0200
Subject: [PATCH 2/6] proposes a variant that can do _both_.

---
 src/groups/general_unitary_groups.jl    | 21 +++++++++++++++++++++
 src/manifolds/GeneralUnitaryMatrices.jl |  7 +++----
 2 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index ba22702444..619850d259 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -258,6 +258,27 @@ function log_lie!(
     return project!(G, X, Identity(G), X)
 end
 
+@doc raw"""
+     project(M::Orthogonal{n}, p, X)
+     project(M::SpecialUnitary{n}, p, X)
+     project(M::SpecialOrthogonal{n}, p, X)
+     project(M::SpecialUnitary{n}, p, X)
+
+Orthogonally project the tangent vector ``X ∈ 𝔽^{n × n}``, ``\mathbb F ∈ \{\mathbb R, \mathbb C\}``
+to the tangent space of `M` at `p`,
+where `X` is already assumed to be in a Lie algebra represemtation, that isand change the representer to use the corresponding Lie algebra, i.e. we compute
+
+```math
+    \operatorname{proj}_p(X) = \frac{X-X^{\mathrm{T}}}{2},
+```
+"""
+project(::GeneralUnitaryMultiplicationGroup, p, X)
+
+function project!(::GeneralUnitaryMultiplicationGroup{n,𝔽}, Y, p, X) where {n,𝔽}
+    project!(SkewHermitianMatrices(n, 𝔽), Y, X)
+    return Y
+end
+
 function Random.rand!(G::GeneralUnitaryMultiplicationGroup, pX; kwargs...)
     rand!(G.manifold, pX; kwargs...)
     return pX
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index 9521671c1a..6342fe6fa4 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -666,17 +666,16 @@ end
 
 Orthogonally project the tangent vector ``X ∈ 𝔽^{n × n}``, ``\mathbb F ∈ \{\mathbb R, \mathbb C\}``
 to the tangent space of `M` at `p`,
-where `X` is already assumed to be in a Lie algebra represemtation, that is,
-the projection is perfomed onto the Lie algebra. We compute
+and change the representer to use the corresponding Lie algebra, i.e. we compute
 
 ```math
-    \operatorname{proj}_p(X) = \frac{X-X^{\mathrm{T}}}{2},
+    \operatorname{proj}_p(X) = \frac{pX-(pX)^{\mathrm{T}}}{2},
 ```
 """
 project(::GeneralUnitaryMatrices, p, X)
 
 function project!(::GeneralUnitaryMatrices{n,𝔽}, Y, p, X) where {n,𝔽}
-    project!(SkewHermitianMatrices(n, 𝔽), Y, X)
+    project!(SkewHermitianMatrices(n, 𝔽), Y, p \ X)
     return Y
 end
 

From bfb6377ac609add5f4001e651aa3385c393511eb Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 08:01:30 +0200
Subject: [PATCH 3/6] fix project for general unitary.

---
 src/manifolds/GeneralUnitaryMatrices.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index 6342fe6fa4..7c447f9424 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -128,13 +128,13 @@ end
     check_vector(M::GeneralUnitaryMatrices{n,𝔽}, p, X; kwargs... )
 
 Check whether `X` is a tangent vector to `p` on the [`UnitaryMatrices`](@ref)
-space `M`, i.e. after [`check_point`](@ref)`(M,p)`, `X` has to be skew symmetric (Hermitian)
-and orthogonal to `p`.
+space `M`, i.e. after [`check_point`](@ref)`(M,p)`,
+``p^{-1}X`` has to be skew symmetric (Hermitian).
 
 The tolerance for the last test can be set using the `kwargs...`.
 """
 function check_vector(M::GeneralUnitaryMatrices{n,𝔽}, p, X; kwargs...) where {n,𝔽}
-    return check_point(SkewHermitianMatrices(n, 𝔽), X; kwargs...)
+    return check_point(SkewHermitianMatrices(n, 𝔽), p \ X; kwargs...)
 end
 
 @doc raw"""
@@ -669,13 +669,14 @@ to the tangent space of `M` at `p`,
 and change the representer to use the corresponding Lie algebra, i.e. we compute
 
 ```math
-    \operatorname{proj}_p(X) = \frac{pX-(pX)^{\mathrm{T}}}{2},
+    \operatorname{proj}_p(X) = p\frac{p^{-1}X-(p^{-1}X)^{\mathrm{T}}}{2},
 ```
 """
 project(::GeneralUnitaryMatrices, p, X)
 
 function project!(::GeneralUnitaryMatrices{n,𝔽}, Y, p, X) where {n,𝔽}
     project!(SkewHermitianMatrices(n, 𝔽), Y, p \ X)
+    mul!(Y, p, Y)
     return Y
 end
 

From 330dc8b8afc8d3af6cedae1e9b441e5002748c41 Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 10:24:38 +0200
Subject: [PATCH 4/6] Trying a few next steps, it is a little tricky to unify.

---
 src/groups/general_unitary_groups.jl    |  97 +++++++++++++++++++++
 src/manifolds/GeneralUnitaryMatrices.jl | 109 ++++--------------------
 src/manifolds/Orthogonal.jl             |  15 +++-
 src/manifolds/Rotations.jl              |   9 ++
 src/manifolds/Unitary.jl                |  16 ++--
 test/manifolds/rotations.jl             |   8 +-
 6 files changed, 149 insertions(+), 105 deletions(-)

diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index 619850d259..87ffdf0f8f 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -44,6 +44,8 @@ end
 
 decorated_manifold(G::GeneralUnitaryMultiplicationGroup) = G.manifold
 
+embed!(G::GeneralUnitaryMultiplicationGroup, Y, p, X) = copyto!(G, Y, p, X)
+
 @doc raw"""
      exp_lie(G::Orthogonal{2}, X)
      exp_lie(G::SpecialOrthogonal{2}, X)
@@ -178,6 +180,101 @@ function log(G::GeneralUnitaryMultiplicationGroup, ::Identity{MultiplicationOper
     return log_lie(G, q)
 end
 
+@doc raw"""
+    log(M::Orthogonal, p, X)
+    log(M::SpecialOrthogonal, p, X)
+    log(M::SpecialUnitary, p, X)
+    log(M::Unitary, p, X)
+
+Compute the logarithmic map, that is, since the resulting ``X`` is represented in the Lie algebra,
+
+```
+log_p q = \log(p^{\mathrm{H}q)
+```
+which is projected onto the skew symmetric matrices for numerical stability.
+"""
+log(::GeneralUnitaryMultiplicationGroup, p, q)
+
+@doc raw"""
+    log(M::GeneralUnitaryMultiplicationGroup, p, q)
+
+Compute the logarithmic map on the [`Rotations`](@ref) manifold
+`M` which is given by
+
+```math
+\log_p q = \operatorname{log}(p^{\mathrm{T}}q)
+```
+
+where $\operatorname{Log}$ denotes the matrix logarithm. For numerical stability,
+the result is projected onto the set of skew symmetric matrices.
+
+For antipodal rotations the function returns deterministically one of the tangent vectors
+that point at `q`.
+"""
+log(::GeneralUnitaryMultiplicationGroup{n,ℝ}, ::Any...) where {n}
+function ManifoldsBase.log(M::GeneralUnitaryMultiplicationGroup{2,ℝ}, p, q)
+    U = transpose(p) * q
+    @assert size(U) == (2, 2)
+    @inbounds θ = atan(U[2], U[1])
+    return get_vector(M, p, θ, DefaultOrthogonalBasis())
+end
+function log!(::GeneralUnitaryMultiplicationGroup{n,ℝ}, X, p, q) where {n}
+    U = transpose(p) * q
+    X .= real(log_safe(U))
+    return project!(SkewSymmetricMatrices(n), X, p, X)
+end
+function log!(M::GeneralUnitaryMultiplicationGroup{2,ℝ}, X, p, q)
+    U = transpose(p) * q
+    @assert size(U) == (2, 2)
+    @inbounds θ = atan(U[2], U[1])
+    return get_vector!(M, X, p, θ, DefaultOrthogonalBasis())
+end
+function log!(M::GeneralUnitaryMultiplicationGroup{3,ℝ}, X, p, q)
+    U = transpose(p) * q
+    cosθ = (tr(U) - 1) / 2
+    if cosθ ≈ -1
+        eig = eigen_safe(U)
+        ival = findfirst(λ -> isapprox(λ, 1), eig.values)
+        inds = SVector{3}(1:3)
+        ax = eig.vectors[inds, ival]
+        return get_vector!(M, X, p, π * ax, DefaultOrthogonalBasis())
+    end
+    X .= U ./ usinc_from_cos(cosθ)
+    return project!(SkewSymmetricMatrices(3), X, p, X)
+end
+function log!(::GeneralUnitaryMultiplicationGroup{4,ℝ}, X, p, q)
+    U = transpose(p) * q
+    cosα, cosβ = Manifolds.cos_angles_4d_rotation_matrix(U)
+    α = acos(clamp(cosα, -1, 1))
+    β = acos(clamp(cosβ, -1, 1))
+    if α ≈ 0 && β ≈ π
+        A² = Symmetric((U - I) ./ 2)
+        P = eigvecs(A²)
+        E = similar(U)
+        fill!(E, 0)
+        @inbounds begin
+            E[2, 1] = -β
+            E[1, 2] = β
+        end
+        copyto!(X, P * E * transpose(P))
+    else
+        det(U) < 0 && throw(
+            DomainError(
+                "The logarithm is not defined for $p and $q with a negative determinant of p'q) ($(det(U)) < 0).",
+            ),
+        )
+        copyto!(X, real(Manifolds.log_safe(U)))
+    end
+    return project!(SkewSymmetricMatrices(4), X, p, X)
+end
+
+function log!(::GeneralUnitaryMultiplicationGroup{n,𝔽}, X, p, q) where {n,𝔽}
+    log_safe!(X, adjoint(p) * q)
+    project!(SkewHermitianMatrices(n, 𝔽), X, q)
+    X .= p * X
+    return X
+end
+
 function log!(
     G::GeneralUnitaryMultiplicationGroup,
     X,
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index 7c447f9424..f665c327a4 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -175,19 +175,14 @@ embed(::GeneralUnitaryMatrices, p) = p
 @doc raw"""
     embed(M::GeneralUnitaryMatrices{n,𝔽}, p, X)
 
-Embed the tangent vector `X` at point `p` in `M` from
-its Lie algebra representation (set of skew matrices) into the
-Riemannian submanifold representation
-
-The formula reads
-```math
-X_{\text{embedded}} = p * X
+Embed the tangent vector `X` at point `p` in `M` is assumed to be of the form
+``X = pY``, where  ``Y``is skew symmetric, so the embedding is the identity.
 ```
 """
 embed(::GeneralUnitaryMatrices, p, X)
 
-function embed!(::GeneralUnitaryMatrices, Y, p, X)
-    return mul!(Y, p, X)
+function embed!(G::GeneralUnitaryMatrices, Y, p, X)
+    return copyto!(G, Y, p, X)
 end
 
 @doc raw"""
@@ -510,95 +505,24 @@ inner(::GeneralUnitaryMatrices, p, X, Y) = dot(X, Y)
     log(M::OrthogonalMatrices, p, X)
     log(M::UnitaryMatrices, p, X)
 
-Compute the logarithmic map, that is, since the resulting ``X`` is represented in the Lie algebra,
+Compute the logarithmic map, that is,
 
 ```
-log_p q = \log(p^{\mathrm{H}q)
+log_p q = p\log(p^{\mathrm{H}q)
 ```
-which is projected onto the skew symmetric matrices for numerical stability.
 """
 log(::GeneralUnitaryMatrices, p, q)
 
-@doc raw"""
-    log(M::Rotations, p, q)
-
-Compute the logarithmic map on the [`Rotations`](@ref) manifold
-`M` which is given by
-
-```math
-\log_p q = \operatorname{log}(p^{\mathrm{T}}q)
-```
-
-where $\operatorname{Log}$ denotes the matrix logarithm. For numerical stability,
-the result is projected onto the set of skew symmetric matrices.
-
-For antipodal rotations the function returns deterministically one of the tangent vectors
-that point at `q`.
-"""
-log(::GeneralUnitaryMatrices{n,ℝ}, ::Any...) where {n}
-function ManifoldsBase.log(M::GeneralUnitaryMatrices{2,ℝ}, p, q)
-    U = transpose(p) * q
-    @assert size(U) == (2, 2)
-    @inbounds θ = atan(U[2], U[1])
-    return get_vector(M, p, θ, DefaultOrthogonalBasis())
-end
-function log!(::GeneralUnitaryMatrices{n,ℝ}, X, p, q) where {n}
-    U = transpose(p) * q
-    X .= real(log_safe(U))
-    return project!(SkewSymmetricMatrices(n), X, p, X)
-end
-function log!(M::GeneralUnitaryMatrices{2,ℝ}, X, p, q)
-    U = transpose(p) * q
-    @assert size(U) == (2, 2)
-    @inbounds θ = atan(U[2], U[1])
-    return get_vector!(M, X, p, θ, DefaultOrthogonalBasis())
-end
-function log!(M::GeneralUnitaryMatrices{3,ℝ}, X, p, q)
-    U = transpose(p) * q
-    cosθ = (tr(U) - 1) / 2
-    if cosθ ≈ -1
-        eig = eigen_safe(U)
-        ival = findfirst(λ -> isapprox(λ, 1), eig.values)
-        inds = SVector{3}(1:3)
-        ax = eig.vectors[inds, ival]
-        return get_vector!(M, X, p, π * ax, DefaultOrthogonalBasis())
-    end
-    X .= U ./ usinc_from_cos(cosθ)
-    return project!(SkewSymmetricMatrices(3), X, p, X)
-end
-function log!(::GeneralUnitaryMatrices{4,ℝ}, X, p, q)
-    U = transpose(p) * q
-    cosα, cosβ = Manifolds.cos_angles_4d_rotation_matrix(U)
-    α = acos(clamp(cosα, -1, 1))
-    β = acos(clamp(cosβ, -1, 1))
-    if α ≈ 0 && β ≈ π
-        A² = Symmetric((U - I) ./ 2)
-        P = eigvecs(A²)
-        E = similar(U)
-        fill!(E, 0)
-        @inbounds begin
-            E[2, 1] = -β
-            E[1, 2] = β
-        end
-        copyto!(X, P * E * transpose(P))
-    else
-        det(U) < 0 && throw(
-            DomainError(
-                "The logarithm is not defined for $p and $q with a negative determinant of p'q) ($(det(U)) < 0).",
-            ),
-        )
-        copyto!(X, real(Manifolds.log_safe(U)))
-    end
-    return project!(SkewSymmetricMatrices(4), X, p, X)
+function ManifoldsBase.log(::GeneralUnitaryMatrices{2,ℝ}, p, q)
+    return p * log(GeneralUnitaryMultiplicationGroup{n,2}(), p, q)
 end
-
-function log!(::GeneralUnitaryMatrices{n,𝔽}, X, p, q) where {n,𝔽}
-    log_safe!(X, adjoint(p) * q)
-    project!(SkewHermitianMatrices(n, 𝔽), q, q)
-    return q
+function log!(::GeneralUnitaryMatrices{n,ℝ}, X, p, q) where {n,𝔽}
+    log!(GeneralUnitaryMultiplicationGroup{m,𝔽}(), X, p, q)
+    X .= p * X
+    return X
 end
 
-norm(::GeneralUnitaryMatrices, p, X) = norm(X)
+norm(::GeneralUnitaryMatrices, p, X) = norm(p \ X)
 
 @doc raw"""
     manifold_dimension(M::Rotations)
@@ -610,6 +534,7 @@ Return the dimension of the manifold orthogonal matrices and of the manifold of
 ```
 """
 manifold_dimension(::GeneralUnitaryMatrices{n,ℝ}) where {n} = div(n * (n - 1), 2)
+
 @doc raw"""
     manifold_dimension(M::GeneralUnitaryMatrices{n,ℂ,DeterminantOneMatrices})
 
@@ -676,7 +601,7 @@ project(::GeneralUnitaryMatrices, p, X)
 
 function project!(::GeneralUnitaryMatrices{n,𝔽}, Y, p, X) where {n,𝔽}
     project!(SkewHermitianMatrices(n, 𝔽), Y, p \ X)
-    mul!(Y, p, Y)
+    Y .= p * Y
     return Y
 end
 
@@ -712,14 +637,14 @@ This is also the default retraction on these manifolds.
 retract(::GeneralUnitaryMatrices{n,𝔽}, ::Any, ::Any, ::QRRetraction) where {n,𝔽}
 
 function retract_qr!(::GeneralUnitaryMatrices{n,𝔽}, q::AbstractArray{T}, p, X) where {n,𝔽,T}
-    A = p + p * X
+    A = p + X
     qr_decomp = qr(A)
     d = diag(qr_decomp.R)
     D = Diagonal(sign.(d .+ convert(T, 0.5)))
     return copyto!(q, qr_decomp.Q * D)
 end
 function retract_polar!(M::GeneralUnitaryMatrices{n,𝔽}, q, p, X) where {n,𝔽}
-    A = p + p * X
+    A = p + X
     return project!(M, q, A; check_det=false)
 end
 
diff --git a/src/manifolds/Orthogonal.jl b/src/manifolds/Orthogonal.jl
index 0c7b31ce0b..e2f41d9d6c 100644
--- a/src/manifolds/Orthogonal.jl
+++ b/src/manifolds/Orthogonal.jl
@@ -3,7 +3,20 @@
 
 The manifold of (real) orthogonal matrices ``\mathrm{O}(n)``.
 
-    OrthogonalMatrices(n)
+The tangent space at ``p`` are all matrices of the form ``X = pY``, where ``Y```
+is skew symmetric.
+
+If you prefer to use the representation of tangent vectors in the Lie algebra, see
+the [`Orthogonal`](@ref)`(n)` group.
+
+# Constructor
+
+   OrthogonalMatrices(n)
+
+Generate the ``n × x`` orthognoal matrices manifold.
+
+## See also
+[`Rotations`](@ref), [`UnitaryMatrices`](@ref), [`Orthogonal`](@ref).
 """
 const OrthogonalMatrices{n} = GeneralUnitaryMatrices{n,ℝ,AbsoluteDeterminantOneMatrices}
 
diff --git a/src/manifolds/Rotations.jl b/src/manifolds/Rotations.jl
index 332ee80969..0346186176 100644
--- a/src/manifolds/Rotations.jl
+++ b/src/manifolds/Rotations.jl
@@ -4,11 +4,20 @@
 The manifold of rotation matrices of sice ``n × n``, i.e.
 real-valued orthogonal matrices with determinant ``+1``.
 
+The tangent space at ``p`` are all matrices of the form ``X = pY``, where ``Y```
+is skew symmetric.
+
+If you prefer to use the representation of tangent vectors in the Lie algebra, see
+the [`SpecialOrthogonal`](@ref)`(n)` group.
+
 # Constructor
 
     Rotations(n)
 
 Generate the manifold of ``n × n`` rotation matrices.
+
+## See also
+[`OrthogonalMatrices`](@ref), [`SpecialOrthogonal`](@ref).
 """
 const Rotations{n} = GeneralUnitaryMatrices{n,ℝ,DeterminantOneMatrices}
 
diff --git a/src/manifolds/Unitary.jl b/src/manifolds/Unitary.jl
index 54dc783d1b..2d2306eaf2 100644
--- a/src/manifolds/Unitary.jl
+++ b/src/manifolds/Unitary.jl
@@ -18,11 +18,11 @@ 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.
+If you prefer to use the representation of tangent vectors in the Lie algebra, see
+the [`Unitary`](@ref)`(n)` group.
 
 # Constructor
-    
+
     UnitaryMatrices(n, 𝔽::AbstractNumbers=ℂ)
 
 see also [`OrthogonalMatrices`](@ref) for the real valued case.
@@ -39,7 +39,7 @@ 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)
+    return p * exp(p \ X)
 end
 
 function get_coordinates_orthonormal(::UnitaryMatrices{1,ℍ}, p, X, ::QuaternionNumbers)
@@ -61,12 +61,12 @@ Base.isapprox(::UnitaryMatrices{1,ℍ}, x, y; kwargs...) = isapprox(x[], y[]; kw
 Base.isapprox(::UnitaryMatrices{1,ℍ}, p, X, Y; kwargs...) = isapprox(X[], Y[]; kwargs...)
 
 function log(::UnitaryMatrices{1,ℍ}, p::Number, q::Number)
-    return log(conj(p) * q)
+    return p * log(conj(p) * q)
 end
 
 @doc raw"""
     manifold_dimension(M::UnitaryMatrices{n,ℂ}) where {n}
-    
+
 Return the dimension of the manifold unitary matrices.
 ```math
 \dim_{\mathrm{U}(n)} = n^2.
@@ -75,7 +75,7 @@ Return the dimension of the manifold unitary matrices.
 manifold_dimension(::UnitaryMatrices{n,ℂ}) where {n} = n^2
 @doc raw"""
     manifold_dimension(M::UnitaryMatrices{n,ℍ})
-    
+
 Return the dimension of the manifold unitary matrices.
 ```math
 \dim_{\mathrm{U}(n, ℍ)} = n(2n+1).
@@ -87,7 +87,7 @@ Manifolds.number_of_coordinates(::UnitaryMatrices{1,ℍ}, ::AbstractBasis{ℍ})
 
 project(::UnitaryMatrices{1,ℍ}, p) = normalize(p)
 
-project(::UnitaryMatrices{1,ℍ}, p, X) = (X - conj(X)) / 2
+project(::UnitaryMatrices{1,ℍ}, p, X) = p * (X - conj(X)) / 2
 
 function Random.rand(M::UnitaryMatrices{1,ℍ}; vector_at=nothing)
     if vector_at === nothing
diff --git a/test/manifolds/rotations.jl b/test/manifolds/rotations.jl
index fde041624f..167492144f 100644
--- a/test/manifolds/rotations.jl
+++ b/test/manifolds/rotations.jl
@@ -200,12 +200,12 @@ include("../utils.jl")
     @testset "Convert from Lie algebra representation of tangents to Riemannian submanifold representation" begin
         M = Manifolds.Rotations(3)
         p = project(M, collect(reshape(1.0:9.0, (3, 3))))
-        x = [[0, -1, 3] [1, 0, 2] [-3, -2, 0]]
-        @test is_vector(M, p, x, true)
-        @test embed(M, p, x) == p * x
+        X = [[0, -1, 3] [1, 0, 2] [-3, -2, 0]]
+        @test is_vector(M, p, X, true)
+        @test embed(M, p, X) == X
         Y = zeros((3, 3))
         embed!(M, Y, p, x)
-        @test Y == p * x
+        @test Y == X
         @test Y ≈ p * (p'Y - Y'p) / 2
     end
     @testset "Edge cases of Rotations" begin

From c496be2d6e83e5ce8ab857bfa76954b47c9c560e Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 15:05:44 +0200
Subject: [PATCH 5/6] Further work towards the different representation.

---
 src/groups/general_unitary_groups.jl    | 101 +++++++++++++++++++++
 src/manifolds/GeneralUnitaryMatrices.jl | 112 ++----------------------
 2 files changed, 109 insertions(+), 104 deletions(-)

diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index 87ffdf0f8f..7b088e274c 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -46,6 +46,107 @@ decorated_manifold(G::GeneralUnitaryMultiplicationGroup) = G.manifold
 
 embed!(G::GeneralUnitaryMultiplicationGroup, Y, p, X) = copyto!(G, Y, p, X)
 
+
+@doc raw"""
+    exp(M::Rotations, p, X)
+    exp(M::OrthogonalMatrices, p, X)
+    exp(M::UnitaryMatrices, p, X)
+
+Compute the exponential map, that is, since ``X`` is represented in the Lie algebra,
+
+```
+exp_p(X) = p\mathrm{e}^X
+```
+
+For different sizes, like ``n=2,3,4`` there is specialised implementations
+
+The algorithm used is a more numerically stable form of those proposed in
+[^Gallier2002] and [^Andrica2013].
+
+[^Gallier2002]:
+    > Gallier J.; Xu D.; Computing exponentials of skew-symmetric matrices
+    > and logarithms of orthogonal matrices.
+    > International Journal of Robotics and Automation (2002), 17(4), pp. 1-11.
+    > [pdf](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205).
+
+[^Andrica2013]:
+    > Andrica D.; Rohan R.-A.; Computing the Rodrigues coefficients of the
+    > exponential map of the Lie groups of matrices.
+    > Balkan Journal of Geometry and Its Applications (2013), 18(2), pp. 1-2.
+    > [pdf](https://www.emis.de/journals/BJGA/v18n2/B18-2-an.pdf).
+
+"""
+exp(::GeneralUnitaryMultiplicationGroup, p, X)
+
+function exp!(M::GeneralUnitaryMultiplicationGroup, q, p, X)
+    return copyto!(M, q, p * exp(X))
+end
+
+function exp(M::GeneralUnitaryMultiplicationGroup{2,ℝ}, p::SMatrix, X::SMatrix)
+    θ = get_coordinates(M, p, X, DefaultOrthogonalBasis())[1]
+    sinθ, cosθ = sincos(θ)
+    return p * SA[cosθ -sinθ; sinθ cosθ]
+end
+function exp!(M::GeneralUnitaryMultiplicationGroup{2,ℝ}, q, p, X)
+    @assert size(q) == (2, 2)
+    θ = get_coordinates(M, p, X, DefaultOrthogonalBasis())[1]
+    sinθ, cosθ = sincos(θ)
+    return copyto!(q, p * SA[cosθ -sinθ; sinθ cosθ])
+end
+function exp!(M::GeneralUnitaryMultiplicationGroup{3,ℝ}, q, p, X)
+    θ = norm(M, p, X) / sqrt(2)
+    if θ ≈ 0
+        a = 1 - θ^2 / 6
+        b = θ / 2
+    else
+        a = sin(θ) / θ
+        b = (1 - cos(θ)) / θ^2
+    end
+    pinvq = I + a .* X .+ b .* (X^2)
+    return copyto!(q, p * pinvq)
+end
+function exp!(::GeneralUnitaryMultiplicationGroup{4,ℝ}, q, p, X)
+    T = eltype(X)
+    α, β = angles_4d_skew_sym_matrix(X)
+    sinα, cosα = sincos(α)
+    sinβ, cosβ = sincos(β)
+    α² = α^2
+    β² = β^2
+    Δ = β² - α²
+    if !isapprox(Δ, 0; atol=1e-6)  # Case α > β ≥ 0
+        sincα = sinα / α
+        sincβ = β == 0 ? one(T) : sinβ / β
+        a₀ = (β² * cosα - α² * cosβ) / Δ
+        a₁ = (β² * sincα - α² * sincβ) / Δ
+        a₂ = (cosα - cosβ) / Δ
+        a₃ = (sincα - sincβ) / Δ
+    elseif α == 0 # Case α = β = 0
+        a₀ = one(T)
+        a₁ = one(T)
+        a₂ = T(1 / 2)
+        a₃ = T(1 / 6)
+    else  # Case α ⪆ β ≥ 0, α ≠ 0
+        sincα = sinα / α
+        r = β / α
+        c = 1 / (1 + r)
+        d = α * (α - β) / 2
+        if α < 1e-2
+            e = @evalpoly(α², T(1 / 3), T(-1 / 30), T(1 / 840), T(-1 / 45360))
+        else
+            e = (sincα - cosα) / α²
+        end
+        a₀ = (α * sinα + (1 + r - d) * cosα) * c
+        a₁ = ((3 - d) * sincα - (2 - r) * cosα) * c
+        a₂ = (sincα - (1 - r) / 2 * cosα) * c
+        a₃ = (e + (1 - r) * (e - sincα / 2)) * c
+    end
+
+    X² = X * X
+    X³ = X² * X
+    pinvq = a₀ * I + a₁ .* X .+ a₂ .* X² .+ a₃ .* X³
+    return copyto!(q, p * pinvq)
+end
+
 @doc raw"""
      exp_lie(G::Orthogonal{2}, X)
      exp_lie(G::SpecialOrthogonal{2}, X)
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index f665c327a4..bcc9bd48da 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -134,7 +134,7 @@ space `M`, i.e. after [`check_point`](@ref)`(M,p)`,
 The tolerance for the last test can be set using the `kwargs...`.
 """
 function check_vector(M::GeneralUnitaryMatrices{n,𝔽}, p, X; kwargs...) where {n,𝔽}
-    return check_point(SkewHermitianMatrices(n, 𝔽), p \ X; kwargs...)
+    return check_point(SkewHermitianMatrices(n, 𝔽), transpose(p) * X; kwargs...)
 end
 
 @doc raw"""
@@ -185,104 +185,8 @@ function embed!(G::GeneralUnitaryMatrices, Y, p, X)
     return copyto!(G, Y, p, X)
 end
 
-@doc raw"""
-    exp(M::Rotations, p, X)
-    exp(M::OrthogonalMatrices, p, X)
-    exp(M::UnitaryMatrices, p, X)
-
-Compute the exponential map, that is, since ``X`` is represented in the Lie algebra,
-
-```
-exp_p(X) = p\mathrm{e}^X
-```
-
-For different sizes, like ``n=2,3,4`` there is specialised implementations
-
-The algorithm used is a more numerically stable form of those proposed in
-[^Gallier2002] and [^Andrica2013].
-
-[^Gallier2002]:
-    > Gallier J.; Xu D.; Computing exponentials of skew-symmetric matrices
-    > and logarithms of orthogonal matrices.
-    > International Journal of Robotics and Automation (2002), 17(4), pp. 1-11.
-    > [pdf](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.35.3205).
-
-[^Andrica2013]:
-    > Andrica D.; Rohan R.-A.; Computing the Rodrigues coefficients of the
-    > exponential map of the Lie groups of matrices.
-    > Balkan Journal of Geometry and Its Applications (2013), 18(2), pp. 1-2.
-    > [pdf](https://www.emis.de/journals/BJGA/v18n2/B18-2-an.pdf).
-
-"""
-exp(::GeneralUnitaryMatrices, p, X)
-
-function exp!(M::GeneralUnitaryMatrices, q, p, X)
-    return copyto!(M, q, p * exp(X))
-end
-
-function exp(M::GeneralUnitaryMatrices{2,ℝ}, p::SMatrix, X::SMatrix)
-    θ = get_coordinates(M, p, X, DefaultOrthogonalBasis())[1]
-    sinθ, cosθ = sincos(θ)
-    return p * SA[cosθ -sinθ; sinθ cosθ]
-end
-function exp!(M::GeneralUnitaryMatrices{2,ℝ}, q, p, X)
-    @assert size(q) == (2, 2)
-    θ = get_coordinates(M, p, X, DefaultOrthogonalBasis())[1]
-    sinθ, cosθ = sincos(θ)
-    return copyto!(q, p * SA[cosθ -sinθ; sinθ cosθ])
-end
-function exp!(M::GeneralUnitaryMatrices{3,ℝ}, q, p, X)
-    θ = norm(M, p, X) / sqrt(2)
-    if θ ≈ 0
-        a = 1 - θ^2 / 6
-        b = θ / 2
-    else
-        a = sin(θ) / θ
-        b = (1 - cos(θ)) / θ^2
-    end
-    pinvq = I + a .* X .+ b .* (X^2)
-    return copyto!(q, p * pinvq)
-end
-function exp!(::GeneralUnitaryMatrices{4,ℝ}, q, p, X)
-    T = eltype(X)
-    α, β = angles_4d_skew_sym_matrix(X)
-    sinα, cosα = sincos(α)
-    sinβ, cosβ = sincos(β)
-    α² = α^2
-    β² = β^2
-    Δ = β² - α²
-    if !isapprox(Δ, 0; atol=1e-6)  # Case α > β ≥ 0
-        sincα = sinα / α
-        sincβ = β == 0 ? one(T) : sinβ / β
-        a₀ = (β² * cosα - α² * cosβ) / Δ
-        a₁ = (β² * sincα - α² * sincβ) / Δ
-        a₂ = (cosα - cosβ) / Δ
-        a₃ = (sincα - sincβ) / Δ
-    elseif α == 0 # Case α = β = 0
-        a₀ = one(T)
-        a₁ = one(T)
-        a₂ = T(1 / 2)
-        a₃ = T(1 / 6)
-    else  # Case α ⪆ β ≥ 0, α ≠ 0
-        sincα = sinα / α
-        r = β / α
-        c = 1 / (1 + r)
-        d = α * (α - β) / 2
-        if α < 1e-2
-            e = @evalpoly(α², T(1 / 3), T(-1 / 30), T(1 / 840), T(-1 / 45360))
-        else
-            e = (sincα - cosα) / α²
-        end
-        a₀ = (α * sinα + (1 + r - d) * cosα) * c
-        a₁ = ((3 - d) * sincα - (2 - r) * cosα) * c
-        a₂ = (sincα - (1 - r) / 2 * cosα) * c
-        a₃ = (e + (1 - r) * (e - sincα / 2)) * c
-    end
-
-    X² = X * X
-    X³ = X² * X
-    pinvq = a₀ * I + a₁ .* X .+ a₂ .* X² .+ a₃ .* X³
-    return copyto!(q, p * pinvq)
+function exp!(M::GeneralUnitaryMatrices{n,𝔽,S}, q, p, X) where {n,𝔽,S}
+    return exp!(GeneralUnitaryMultiplicationGroup{n,𝔽,S}(M), q, p, adjoint(p)*X)
 end
 
 @doc raw"""
@@ -513,16 +417,16 @@ log_p q = p\log(p^{\mathrm{H}q)
 """
 log(::GeneralUnitaryMatrices, p, q)
 
-function ManifoldsBase.log(::GeneralUnitaryMatrices{2,ℝ}, p, q)
-    return p * log(GeneralUnitaryMultiplicationGroup{n,2}(), p, q)
+function ManifoldsBase.log(M::GeneralUnitaryMatrices{2,ℝ,S}, p, q) where {S}
+    return p * log(GeneralUnitaryMultiplicationGroup{2,ℝ,S}(M), p, q)
 end
-function log!(::GeneralUnitaryMatrices{n,ℝ}, X, p, q) where {n,𝔽}
-    log!(GeneralUnitaryMultiplicationGroup{m,𝔽}(), X, p, q)
+function log!(M::GeneralUnitaryMatrices{n,𝔽,S}, X, p, q) where {n,𝔽,S}
+    log!(GeneralUnitaryMultiplicationGroup{n,𝔽,S}(M), X, p, q)
     X .= p * X
     return X
 end
 
-norm(::GeneralUnitaryMatrices, p, X) = norm(p \ X)
+norm(::GeneralUnitaryMatrices, p, X) = norm(adjoint(p) * X)
 
 @doc raw"""
     manifold_dimension(M::Rotations)

From a64c1a68a13ff492b1ae5c2c4111785ce6be3087 Mon Sep 17 00:00:00 2001
From: Ronny Bergmann <git@ronnybergmann.net>
Date: Sat, 13 Aug 2022 15:08:45 +0200
Subject: [PATCH 6/6] formatter.

---
 src/groups/general_unitary_groups.jl    | 1 -
 src/manifolds/GeneralUnitaryMatrices.jl | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index 7b088e274c..becce68b0e 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -46,7 +46,6 @@ decorated_manifold(G::GeneralUnitaryMultiplicationGroup) = G.manifold
 
 embed!(G::GeneralUnitaryMultiplicationGroup, Y, p, X) = copyto!(G, Y, p, X)
 
-
 @doc raw"""
     exp(M::Rotations, p, X)
     exp(M::OrthogonalMatrices, p, X)
diff --git a/src/manifolds/GeneralUnitaryMatrices.jl b/src/manifolds/GeneralUnitaryMatrices.jl
index bcc9bd48da..ac121a9665 100644
--- a/src/manifolds/GeneralUnitaryMatrices.jl
+++ b/src/manifolds/GeneralUnitaryMatrices.jl
@@ -186,7 +186,7 @@ function embed!(G::GeneralUnitaryMatrices, Y, p, X)
 end
 
 function exp!(M::GeneralUnitaryMatrices{n,𝔽,S}, q, p, X) where {n,𝔽,S}
-    return exp!(GeneralUnitaryMultiplicationGroup{n,𝔽,S}(M), q, p, adjoint(p)*X)
+    return exp!(GeneralUnitaryMultiplicationGroup{n,𝔽,S}(M), q, p, adjoint(p) * X)
 end
 
 @doc raw"""