diff --git a/NEWS.md b/NEWS.md
index b0b19478a5..173678d256 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,12 +5,17 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.9.4] - 2023-11-04
+## [0.9.4] - 2023-11-06
 
 ### Added
 
+- Functions `inv_diff`, `inv_diff!`, `adjoint_inv_diff` and `adjoint_inv_diff!` that correspond to differentials and pullbacks of group inversion.
 - Julia 1.10-rc CI workflow.
 
+### Fixed
+
+- Fixed issue with incorrect implementation of `apply_diff_group` in `GroupOperationAction` with left backward and right forward action [#669](https://github.com/JuliaManifolds/Manifolds.jl/issues/669).
+
 ## [0.9.3] - 2023-10-28
 
 ### Added
diff --git a/Project.toml b/Project.toml
index 60c1d39781..97d08c3128 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.9.3"
+version = "0.9.4"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/docs/src/references.bib b/docs/src/references.bib
index f313545904..ecd9632754 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -338,6 +338,19 @@ @article{GaoSonAbsilStykel:2021
     TITLE   = {Riemannian Optimization on the Symplectic Stiefel Manifold},
     JOURNAL = {SIAM Journal on Optimization}
 }
+@inproceedings{Giles:2008,
+	address = {Berlin, Heidelberg},
+	series = {Lecture {Notes} in {Computational} {Science} and {Engineering}},
+	title = {Collected {Matrix} {Derivative} {Results} for {Forward} and {Reverse} {Mode} {Algorithmic} {Differentiation}},
+	isbn = {978-3-540-68942-3},
+	doi = {10.1007/978-3-540-68942-3_4},
+	booktitle = {Advances in {Automatic} {Differentiation}},
+	publisher = {Springer},
+	author = {Giles, Mike B.},
+	editor = {Bischof, Christian H. and Bücker, H. Martin and Hovland, Paul and Naumann, Uwe and Utke, Jean},
+	year = {2008},
+	pages = {35--44},
+}
 @incollection{GuiguiMaignantTroouvePennec:2021,
     DOI       = {10.1007/978-3-030-80209-7_12},
     YEAR      = {2021},
diff --git a/ext/ManifoldsTestExt/tests_general.jl b/ext/ManifoldsTestExt/tests_general.jl
index 54b37d7cd4..872cb91733 100644
--- a/ext/ManifoldsTestExt/tests_general.jl
+++ b/ext/ManifoldsTestExt/tests_general.jl
@@ -153,7 +153,7 @@ function test_manifold(
             ) for i in 1:n
         ]
     end
-    Test.@testset "dimension" begin
+    Test.@testset "dimension" begin # COV_EXCL_LINE
         Test.@test isa(manifold_dimension(M), expected_dimension_type)
         Test.@test manifold_dimension(M) ≥ 0
     end
@@ -180,7 +180,7 @@ function test_manifold(
         end
     end
 
-    Test.@testset "is_point" begin
+    Test.@testset "is_point" begin # COV_EXCL_LINE
         for pt in pts
             atol = is_point_atol_multiplier * find_eps(pt)
             Test.@test is_point(M, pt; atol=atol)
@@ -303,7 +303,7 @@ function test_manifold(
         mutating=is_mutating,
     )
 
-    Test.@testset "(inverse &) retraction tests" begin
+    Test.@testset "(inverse &) retraction tests" begin # COV_EXCL_LINE
         for (p, X) in zip(pts, tv)
             epsx = find_eps(p)
             point_atol = is_point_atol_multiplier * find_eps(p)
@@ -325,7 +325,7 @@ function test_manifold(
                 end
                 Test.@test is_point(M, new_pt; atol=point_atol)
                 (test_inplace && is_mutating) &&
-                    Test.@testset "inplace test for retract!" begin
+                    Test.@testset "inplace test for retract!" begin # COV_EXCL_LINE
                         p2 = copy(M, p)
                         X2 = copy(M, p, X)
                         q = retract(M, p2, X2, retr_method)
@@ -367,7 +367,7 @@ function test_manifold(
         end
     end
 
-    Test.@testset "atlases" begin
+    Test.@testset "atlases" begin # COV_EXCL_LINE
         if !isempty(test_atlases)
             Test.@test get_default_atlas(M) isa AbstractAtlas{ℝ}
         end
@@ -411,7 +411,7 @@ function test_manifold(
         end
     end
 
-    Test.@testset "basic linear algebra in tangent space" begin
+    Test.@testset "basic linear algebra in tangent space" begin # COV_EXCL_LINE
         for (p, X) in zip(pts, tv)
             Test.@test isapprox(M, p, 0 * X, zero_vector(M, p); atol=find_eps(pts[1]))
             Test.@test isapprox(M, p, 2 * X, X + X)
@@ -421,7 +421,7 @@ function test_manifold(
     end
 
     test_tangent_vector_broadcasting &&
-        Test.@testset "broadcasted linear algebra in tangent space" begin
+        Test.@testset "broadcasted linear algebra in tangent space" begin # COV_EXCL_LINE
             for (p, X) in zip(pts, tv)
                 Test.@test isapprox(M, p, 3 * X, 2 .* X .+ X)
                 Test.@test isapprox(M, p, -X, X .- 2 .* X)
@@ -469,7 +469,7 @@ function test_manifold(
     !(
         default_retraction_method === nothing ||
         default_inverse_retraction_method === nothing
-    ) && Test.@testset "vector transport" begin
+    ) && Test.@testset "vector transport" begin # COV_EXCL_LINE
         tvatol = is_tangent_atol_multiplier * find_eps(pts[1])
         X1 = inverse_retract(M, pts[1], pts[2], default_inverse_retraction_method)
         X2 = inverse_retract(M, pts[1], pts[3], default_inverse_retraction_method)
@@ -499,7 +499,7 @@ function test_manifold(
             vector_transport_retractions,
             vector_transport_inverse_retractions,
         )
-            Test.@testset "vector transport method $(vtm)" begin
+            Test.@testset "vector transport method $(vtm)" begin # COV_EXCL_LINE
                 tvatol = is_tangent_atol_multiplier * find_eps(pts[1])
                 X1 = inverse_retract(M, pts[1], pts[2], irtr_m)
                 X2 = inverse_retract(M, pts[1], pts[3], irtr_m)
@@ -531,7 +531,7 @@ function test_manifold(
                         vector_transport_to!(M, v1t1_m, pts[1], X1, pts32, vtm)
                         Test.@test isapprox(M, pts32, v1t1, v1t1_m; atol=tvatol)
                         test_inplace &&
-                            Test.@testset "inplace test for vector_transport_to!" begin
+                            Test.@testset "inplace test for vector_transport_to!" begin # COV_EXCL_LINE
                                 X1a = copy(M, pts[1], X1)
                                 Xt = vector_transport_to(M, pts[1], X1, pts32, vtm)
                                 vector_transport_to!(M, X1a, pts[1], X1a, pts32, vtm)
@@ -559,7 +559,7 @@ function test_manifold(
     end
 
     for btype in basis_types_vecs
-        Test.@testset "Basis support for $(btype)" begin
+        Test.@testset "Basis support for $(btype)" begin # COV_EXCL_LINE
             p = pts[1]
             b = get_basis(M, p, btype)
             Test.@test isa(b, CachedBasis)
@@ -713,7 +713,7 @@ function test_manifold(
         end
     end
 
-    Test.@testset "number_eltype" begin
+    Test.@testset "number_eltype" begin # COV_EXCL_LINE
         for (p, X) in zip(pts, tv)
             Test.@test number_eltype(X) == number_eltype(p)
             p = retract(M, p, X, default_retraction_method)
@@ -822,7 +822,7 @@ function test_manifold(
         end
     end
 
-    Test.@testset "tangent vector distributions" begin
+    Test.@testset "tangent vector distributions" begin # COV_EXCL_LINE
         for tvd in tvector_distributions
             supp = Manifolds.support(tvd)
             Test.@test supp isa
@@ -864,9 +864,10 @@ function test_parallel_transport(
 )
     length(P) < 2 &&
         error("The Parallel Transport test set requires at least 2 points in P")
-    Test.@testset "Test Parallel Transport" begin
+    Test.@testset "Test Parallel Transport" begin # COV_EXCL_LINE
         along && @warn "parallel transport along test not yet implemented"
-        Test.@testset "To (a point)" begin # even with to =false this displays no tests
+        Test.@testset "To (a point)" begin # COV_EXCL_LINE
+            # even with to =false this displays no tests
             if to
                 for i in 1:(length(P) - 1)
                     p = P[i]
@@ -890,7 +891,7 @@ function test_parallel_transport(
                 end
             end
         end
-        Test.@testset "(Tangent Vector) Direction" begin
+        Test.@testset "(Tangent Vector) Direction" begin # COV_EXCL_LINE
             if direction
                 for i in 1:(length(P) - 1)
                     p = P[i]
diff --git a/ext/ManifoldsTestExt/tests_group.jl b/ext/ManifoldsTestExt/tests_group.jl
index 10ead50598..1b6023ff44 100644
--- a/ext/ManifoldsTestExt/tests_group.jl
+++ b/ext/ManifoldsTestExt/tests_group.jl
@@ -3,19 +3,90 @@ using Base: IdentityUnitRange
 using Manifolds:
     LeftForwardAction, LeftBackwardAction, RightForwardAction, RightBackwardAction
 
+translate_diff_id(G, q, X, conv) = translate_diff(G, q, identity_element(G), X, conv)
+
+"""
+    _move(G::AbstractManifold, q, X, ::LeftSide)
+
+The left translation ``τ_L(g, X) = gX``.
+"""
+_move(G::AbstractManifold, q, X, ::LeftSide) =
+    translate_diff_id(G, q, X, (LeftAction(), LeftSide()))
+"""
+    _move(G::AbstractManifold, q, X, ::RightSide)
+
+The right translation ``τ_R(g, X) = Xg``.
+"""
+_move(G::AbstractManifold, q, X, ::RightSide) =
+    translate_diff_id(G, q, X, (RightAction(), RightSide()))
+
+"""
+    test_inv_diff_fn(G::AbstractManifold, q, X, conv)
+
+Test the differential of the inverse on a Lie group.
+Denote this inverse by ``I(g) := g^{-1}``.
+If the left and right transports are ``τ_L(g,X) := gX``
+and ``τ_R(g,X) := Xg`` respectively, then
+```math
+⟨DI(g), τ_L(g,X)⟩ = -τ_R(g^{-1}, X)
+```
+and
+``` math
+⟨DI(g), τ_R(g,X)⟩ = -τ_L(g^{-1}, X)
+```
+"""
+function test_inv_diff_fn(G::AbstractManifold, q, X, conv)
+    q_ = inv(G, q)
+    Test.@test isapprox(
+        TangentSpace(G, q_),
+        inv_diff(G, q, _move(G, q, X, conv)),
+        -_move(G, q_, X, switch_side(conv)),
+    )
+end
+
+_get_side(::LeftAction) = RightSide()
+_get_side(::RightAction) = LeftSide()
+
+_transporter(G::AbstractManifold, q, X, dir) = _move(G, q, X, _get_side(dir))
+
+@doc raw"""
+    test_apply_diff_group_fn(A::AbstractGroupAction{TAD}, q, X, p) where {TAD}
+
+This should hold for *any* group action ``A`` on any manifold.
+If you define ``π_p(g) := A(g, p)`` for ``g ∈ \mathcal G`` and ``p ∈ \mathcal M``,
+and define, for ``X in Alg(G)``,
+    ``τ_R(g, X) := Xg`` (the right translation),
+and ``τ_L(g, X) := gX`` (the left translation), then we have the identity:
+```math
+⟨Dπ_{q}(g), τ(g, X)⟩ = ⟨Dπ_{A(g, q)}(1), X⟩
+```
+where, for a *left* action, ``τ`` is the *right* translation,
+and for a *right* action, ``τ`` is the *left* translation.
+"""
+function test_apply_diff_group_fn(A::AbstractGroupAction{TAD}, q, X, p) where {TAD}
+    G = base_group(A)
+    p_ = apply(A, q, p)
+    X1 = apply_diff_group(A, q, _transporter(G, q, X, TAD()), p)
+    X2 = apply_diff_group(A, identity_element(G), X, p_)
+    Test.@test isapprox(TangentSpace(G, p_), X1, X2)
+end
+
 """
     test_group(
         G,
         g_pts::AbstractVector,
-        X_pts::AbstractVector = [],
-        Xe_pts::AbstractVector = [];
-        atol = 1e-10,
-        test_mutating = true,
-        test_exp_lie_log = true,
-        test_diff = false,
-        test_invariance = false,
-        test_lie_bracket=false,
-        test_adjoint_action=false,
+        X_pts::AbstractVector=[],
+        Xe_pts::AbstractVector=[];
+        atol::Real=1e-10,
+        test_mutating::Bool=true,
+        test_exp_lie_log::Bool=true,
+        test_diff::Bool=false,
+        test_invariance::Bool=false,
+        test_lie_bracket::Bool=false,
+        test_adjoint_action::Bool=false,
+        test_inv_diff::Bool=false,
+        test_adjoint_inv_diff::Bool=false,
+        test_apply_diff_group::Bool=false,
         diff_convs = [(), (LeftForwardAction(),), (RightBackwardAction(),)],
     )
 
@@ -33,18 +104,21 @@ function test_group(
     g_pts::AbstractVector,
     X_pts::AbstractVector=[],
     Xe_pts::AbstractVector=[];
-    atol=1e-10,
-    test_mutating=true,
-    test_exp_lie_log=true,
-    test_one_arg_identity_element=true,
-    test_diff=false,
-    test_invariance=false,
-    test_lie_bracket=false,
-    test_adjoint_action=false,
+    atol::Real=1e-10,
+    test_mutating::Bool=true,
+    test_exp_lie_log::Bool=true,
+    test_one_arg_identity_element::Bool=true,
+    test_diff::Bool=false,
+    test_invariance::Bool=false,
+    test_lie_bracket::Bool=false,
+    test_adjoint_action::Bool=false,
+    test_inv_diff::Bool=false,
+    test_adjoint_inv_diff::Bool=false,
     diff_convs=[(), (LeftForwardAction(),), (RightBackwardAction(),)],
-    test_log_from_identity=false,
-    test_exp_from_identity=false,
-    test_vee_hat_from_identity=false,
+    test_log_from_identity::Bool=false,
+    test_exp_from_identity::Bool=false,
+    test_vee_hat_from_identity::Bool=false,
+    test_apply_diff_group::Bool=false,
 )
     e = Identity(G)
 
@@ -300,6 +374,21 @@ function test_group(
         end
     end
 
+    test_inv_diff && Test.@testset "Differential of inverse" begin # COV_EXCL_LINE
+        Test.@test isapprox(inv_diff(G, e, Xe_pts[1]), -Xe_pts[1]; atol=atol)
+        Test.@testset "test_inv_diff" for side in [LeftSide(), RightSide()]
+            test_inv_diff_fn(G, g_pts[1], X_pts[1], side)
+        end # COV_EXCL_LINE
+    end
+    test_adjoint_inv_diff && Test.@testset "Differential of inverse" begin # COV_EXCL_LINE
+        Test.@test isapprox(adjoint_inv_diff(G, e, Xe_pts[1]), -Xe_pts[1]; atol=atol)
+        Test.@test isapprox(
+            adjoint_inv_diff(G, g_pts[1], inv_diff(G, g_pts[1], X_pts[1])),
+            X_pts[1];
+            atol=atol,
+        )
+    end
+
     test_exp_lie_log && Test.@testset "group exp/log properties" begin
         Test.@testset "e = exp(0)" begin
             X = log_lie(G, Identity(G))
@@ -317,7 +406,7 @@ function test_group(
             end
         end
 
-        Test.@testset "X = log(exp(X))" begin
+        Test.@testset "X = log(exp(X))" begin # COV_EXCL_LINE
             for X in Xe_pts
                 g = exp_lie(G, X)
                 Test.@test is_point(G, g; atol=atol)
@@ -517,6 +606,19 @@ function test_group(
         end
     end
 
+    test_apply_diff_group && Test.@testset "apply_diff_group" begin
+        Test.@testset "$dir, $side" for side in [LeftSide(), RightSide()],
+            dir in [LeftAction(), RightAction()]
+
+            test_apply_diff_group_fn(
+                GroupOperationAction(G, (dir, side)),
+                g_pts[1],
+                X_pts[1],
+                g_pts[2],
+            )
+        end # COV_EXCL_LINE
+    end
+
     return nothing
 end
 
@@ -710,15 +812,6 @@ function test_action(
                 Test.@test is_vector(M, am, aX, true; atol=atol)
                 Test.@test is_vector(M, ainvm, ainvv, true; atol=atol)
             end
-
-            a12 = compose(A, a_pts[1], a_pts[2])
-            a2m = apply(A, a_pts[2], m)
-            a12X = apply_diff(A, a12, m, X)
-            a2X = apply_diff(A, a_pts[2], m, X)
-            Test.@test isapprox(M, a2m, apply_diff(A, a_pts[1], a2m, a2X), a12X; atol=atol)
-
-            Test.@test isapprox(M, m, apply_diff(A, e, m, X), X; atol=atol)
-            Test.@test isapprox(M, m, inverse_apply_diff(A, e, m, X), X; atol=atol)
         end
 
         test_mutating_action && Test.@testset "mutating" begin
@@ -733,22 +826,6 @@ function test_action(
                     Test.@test is_vector(M, am, aX, true; atol=atol)
                     Test.@test is_vector(M, ainvm, ainvv, true; atol=atol)
                 end
-
-                a12 = compose(A, a_pts[1], a_pts[2])
-                a2m = apply(A, a_pts[2], m)
-                a12m = apply(A, a12, m)
-                a12X, a2X, a1_a2X = allocate(X), allocate(X), allocate(X)
-                Test.@test apply_diff!(A, a12X, a12, m, X) === a12X
-                Test.@test apply_diff!(A, a2X, a_pts[2], m, X) === a2X
-                Test.@test apply_diff!(A, a1_a2X, a_pts[1], a2m, a2X) === a1_a2X
-                Test.@test isapprox(M, a12m, a1_a2X, a12X; atol=atol)
-
-                eX = allocate(X)
-                Test.@test apply_diff!(A, eX, e, m, X) === eX
-                Test.@test isapprox(M, m, eX, X; atol=atol)
-                eX = allocate(X)
-                Test.@test inverse_apply_diff!(A, eX, e, m, X) === eX
-                Test.@test isapprox(M, m, eX, X; atol=atol)
             end
         end
     end
diff --git a/src/Manifolds.jl b/src/Manifolds.jl
index 2c9a5315cb..d5f551a0fb 100644
--- a/src/Manifolds.jl
+++ b/src/Manifolds.jl
@@ -953,6 +953,8 @@ export adjoint_action,
     adjoint_action!,
     adjoint_apply_diff_group,
     adjoint_apply_diff_group!,
+    adjoint_inv_diff,
+    adjoint_inv_diff!,
     affine_matrix,
     apply,
     apply!,
@@ -989,6 +991,8 @@ export adjoint_action,
     identity_element!,
     inv,
     inv!,
+    inv_diff,
+    inv_diff!,
     inverse_apply,
     inverse_apply!,
     inverse_apply_diff,
diff --git a/src/groups/addition_operation.jl b/src/groups/addition_operation.jl
index 2c0f110d47..d4a8920e7c 100644
--- a/src/groups/addition_operation.jl
+++ b/src/groups/addition_operation.jl
@@ -28,6 +28,20 @@ function adjoint_action!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y,
     return copyto!(G, Y, p, X)
 end
 
+"""
+    adjoint_inv_diff(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X)
+
+Compute the value of pullback of additive matrix inversion ``p ↦ -p`` at ``X``, i.e. ``-X``.
+"""
+function adjoint_inv_diff(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X)
+    return -X
+end
+function adjoint_inv_diff!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X)
+    Y .= X
+    Y .*= -1
+    return Y
+end
+
 identity_element(::AdditionGroupTrait, G::AbstractDecoratorManifold, p::Number) = zero(p)
 
 function identity_element!(::AdditionGroupTrait, G::AbstractDecoratorManifold, p)
@@ -55,6 +69,20 @@ function inv!(
     return q
 end
 
+"""
+    inv_diff(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X)
+
+Compute the value of differential of additive matrix inversion ``p ↦ -p`` at ``X``, i.e. ``-X``.
+"""
+function inv_diff(::AdditionGroupTrait, G::AbstractDecoratorManifold, p, X)
+    return -X
+end
+function inv_diff!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X)
+    Y .= X
+    Y .*= -1
+    return Y
+end
+
 function is_identity(::AdditionGroupTrait, G::AbstractDecoratorManifold, q; kwargs...)
     return isapprox(G, q, zero(q); kwargs...)
 end
diff --git a/src/groups/general_unitary_groups.jl b/src/groups/general_unitary_groups.jl
index 804e6eecea..077497c0fd 100644
--- a/src/groups/general_unitary_groups.jl
+++ b/src/groups/general_unitary_groups.jl
@@ -306,7 +306,8 @@ function translate_diff!(
     X,
     ::RightForwardAction,
 )
-    return copyto!(G, Y, X)
+    copyto!(G, Y, X)
+    return Y
 end
 function translate_diff!(
     G::GeneralUnitaryMultiplicationGroup,
@@ -316,7 +317,8 @@ function translate_diff!(
     X,
     ::LeftBackwardAction,
 )
-    return copyto!(G, Y, p * X * inv(G, p))
+    copyto!(G, Y, p * X * inv(G, p))
+    return Y
 end
 function translate_diff!(
     G::GeneralUnitaryMultiplicationGroup,
diff --git a/src/groups/group.jl b/src/groups/group.jl
index b9c6608b20..ed99cf2ab7 100644
--- a/src/groups/group.jl
+++ b/src/groups/group.jl
@@ -463,6 +463,28 @@ function adjoint_action!(
     return Y
 end
 
+@doc raw"""
+    adjoint_inv_diff(G::AbstractDecoratorManifold, p, X)
+
+Compute the value of pullback of inverse ``p^{-1} ∈ \mathcal{G}`` of an element
+``p ∈ \mathcal{G}`` at tangent vector `X` at ``p^{-1}``. The result is a tangent vector
+at ``p``.
+"""
+adjoint_inv_diff(G::AbstractDecoratorManifold, p)
+
+@trait_function adjoint_inv_diff(G::AbstractDecoratorManifold, p, X)
+function adjoint_inv_diff(
+    ::TraitList{<:IsGroupManifold},
+    G::AbstractDecoratorManifold,
+    p,
+    X,
+)
+    Y = allocate_result(G, inv_diff, X, p)
+    return adjoint_inv_diff!(G, Y, p, X)
+end
+
+@trait_function adjoint_inv_diff!(G::AbstractDecoratorManifold, Y, p, X)
+
 function ManifoldDiff.differential_exp_argument_lie_approx!(
     M::AbstractManifold,
     Z,
@@ -529,6 +551,22 @@ function inv!(
     return e
 end
 
+@doc raw"""
+    inv_diff(G::AbstractDecoratorManifold, p, X)
+
+Compute the value of differential of inverse ``p^{-1} ∈ \mathcal{G}`` of an element
+``p ∈ \mathcal{G}`` at tangent vector `X` at `p`. The result is a tangent vector at ``p^{-1}``.
+"""
+inv_diff(G::AbstractDecoratorManifold, p)
+
+@trait_function inv_diff(G::AbstractDecoratorManifold, p, X)
+function inv_diff(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p, X)
+    Y = allocate_result(G, inv_diff, X, p)
+    return inv_diff!(G, Y, p, X)
+end
+
+@trait_function inv_diff!(G::AbstractDecoratorManifold, Y, p, X)
+
 function Base.copyto!(
     ::TraitList{IsGroupManifold{O}},
     ::AbstractDecoratorManifold,
diff --git a/src/groups/group_operation_action.jl b/src/groups/group_operation_action.jl
index 34e6de888c..2a7d06557d 100644
--- a/src/groups/group_operation_action.jl
+++ b/src/groups/group_operation_action.jl
@@ -65,12 +65,35 @@ end
 
 function adjoint_apply_diff_group(A::GroupOperationAction, a, X, p)
     G = base_group(A)
-    return inverse_translate_diff(G, a, p, X, reverse_direction_and_side(A))
+    if direction_and_side(A) === LeftForwardAction() ||
+       direction_and_side(A) === RightBackwardAction()
+        return inverse_translate_diff(G, a, p, X, reverse_direction_and_side(A))
+    else
+        return inverse_translate_diff(
+            G,
+            p,
+            a,
+            adjoint_inv_diff(G, apply(A, a, p), X),
+            (direction(A), switch_side(action_side(A))),
+        )
+    end
 end
 
 function adjoint_apply_diff_group!(A::GroupOperationAction, Y, a, X, p)
     G = base_group(A)
-    return inverse_translate_diff!(G, Y, a, p, X, reverse_direction_and_side(A))
+    if direction_and_side(A) === LeftForwardAction() ||
+       direction_and_side(A) === RightBackwardAction()
+        return inverse_translate_diff!(G, Y, a, p, X, reverse_direction_and_side(A))
+    else
+        return inverse_translate_diff!(
+            G,
+            Y,
+            p,
+            a,
+            adjoint_inv_diff(G, apply(A, a, p), X),
+            (direction(A), switch_side(action_side(A))),
+        )
+    end
 end
 
 apply(A::GroupOperationAction, a, p) = translate(A.group, a, p, direction_and_side(A))
@@ -95,13 +118,65 @@ function apply_diff!(A::GroupOperationAction, Y, a, p, X)
     return translate_diff!(A.group, Y, a, p, X, direction_and_side(A))
 end
 
+@doc raw"""
+    apply_diff_group(A::GroupOperationAction, a, X, p)
+
+Compute differential of [`GroupOperationAction`](@ref) `A` with respect to group element
+at tangent vector `X`:
+
+````math
+(\mathrm{d}τ^p) : T_{a} \mathcal G → T_{τ_a p} \mathcal G
+````
+
+There are four cases:
+* left action from the left side: ``L_a: p ↦ a \circ p``, where
+````math
+(\mathrm{d}L_a) : T_{a} \mathcal G → T_{a \circ p} \mathcal G.
+````
+* right action from the left side: ``L'_a: p ↦ a^{-1} \circ p``, where
+````math
+(\mathrm{d}L'_a) : T_{a} \mathcal G → T_{a^{-1} \circ p} \mathcal G.
+````
+* right action from the right side: ``R_a: p ↦ p \circ a``, where
+````math
+(\mathrm{d}R_a) : T_{a} \mathcal G → T_{p \circ a} \mathcal G.
+````
+* left action from the right side: ``R'_a: p ↦ p \circ a^{-1}``, where
+````math
+(\mathrm{d}R'_a) : T_{a} \mathcal G → T_{p \circ a^{-1}} \mathcal G.
+````
+
+"""
 function apply_diff_group(A::GroupOperationAction, a, X, p)
     G = base_group(A)
-    return translate_diff(G, p, a, X, reverse_direction_and_side(A))
+    if direction_and_side(A) === LeftForwardAction() ||
+       direction_and_side(A) === RightBackwardAction()
+        return translate_diff(G, p, a, X, reverse_direction_and_side(A))
+    else
+        return translate_diff(
+            G,
+            p,
+            a,
+            inv_diff(G, a, X),
+            (direction(A), switch_side(action_side(A))),
+        )
+    end
 end
 function apply_diff_group!(A::GroupOperationAction, Y, a, X, p)
     G = base_group(A)
-    return translate_diff!(G, Y, p, a, X, reverse_direction_and_side(A))
+    if direction_and_side(A) === LeftForwardAction() ||
+       direction_and_side(A) === RightBackwardAction()
+        return translate_diff!(G, Y, p, a, X, reverse_direction_and_side(A))
+    else
+        return translate_diff!(
+            G,
+            Y,
+            p,
+            a,
+            inv_diff(G, a, X),
+            (direction(A), switch_side(action_side(A))),
+        )
+    end
 end
 
 function inverse_apply_diff(A::GroupOperationAction, a, p, X)
diff --git a/src/groups/multiplication_operation.jl b/src/groups/multiplication_operation.jl
index de21633021..b828650287 100644
--- a/src/groups/multiplication_operation.jl
+++ b/src/groups/multiplication_operation.jl
@@ -23,9 +23,47 @@ Base.:\(p, ::Identity{MultiplicationOperation}) = inv(p)
 Base.:\(::Identity{MultiplicationOperation}, p) = p
 Base.:\(e::Identity{MultiplicationOperation}, ::Identity{MultiplicationOperation}) = e
 
+Base.inv(e::Identity{MultiplicationOperation}) = e
+
 LinearAlgebra.det(::Identity{MultiplicationOperation}) = true
 LinearAlgebra.adjoint(e::Identity{MultiplicationOperation}) = e
 
+@doc raw"""
+    adjoint_inv_diff(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, X)
+
+Compute the value of differential of matrix inversion ``p ↦ p^{-1}`` at ``X``.
+When tangent vectors are represented in Lie algebra in a left-invariant way, the formula
+reads ``-p^\mathrm{T}X(p^{-1})^\mathrm{T}``. For matrix groups with ambient space tangent
+vectors, the formula would read ``-(p^{-1})^\mathrm{T}X(p^{-1})^\mathrm{T}``. See the
+section about matrix inverse in [Giles:2008](@cite).
+"""
+function adjoint_inv_diff(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, X)
+    return -(p' * X * inv(G, p)')
+end
+function adjoint_inv_diff(
+    ::MultiplicationGroupTrait,
+    G::AbstractDecoratorManifold,
+    p::AbstractArray{<:Number,0},
+    X::AbstractArray{<:Number,0},
+)
+    p_inv = inv(p[])
+    return -(p[] * X * p_inv)
+end
+
+function adjoint_inv_diff!(
+    ::MultiplicationGroupTrait,
+    G::AbstractDecoratorManifold,
+    Y,
+    p,
+    X,
+)
+    p_inv = inv(p)
+    Z = X * p_inv'
+    mul!(Y, p', Z)
+    Y .*= -1
+    return Y
+end
+
 function identity_element!(
     ::MultiplicationGroupTrait,
     G::AbstractDecoratorManifold,
@@ -124,6 +162,35 @@ function inv!(
     return q
 end
 
+"""
+    inv_diff(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, X)
+
+Compute the value of differential of matrix inversion ``p ↦ p^{-1}`` at ``X``.
+When tangent vectors are represented in Lie algebra in a left-invariant way, the formula
+reads ``-pXp^{-1}``. For matrix groups with ambient space tangent vectors, the formula would
+read ``-p^{-1}Xp^{-1}``. See the section about matrix inverse in [Giles:2008](@cite).
+"""
+function inv_diff(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, X)
+    return -(p * X * inv(G, p))
+end
+function inv_diff(
+    ::MultiplicationGroupTrait,
+    G::AbstractDecoratorManifold,
+    p::AbstractArray{<:Number,0},
+    X::AbstractArray{<:Number,0},
+)
+    p_inv = inv(p[])
+    return -(p[] * X * p_inv)
+end
+
+function inv_diff!(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, Y, p, X)
+    p_inv = inv(p)
+    Z = X * p_inv
+    mul!(Y, p, Z)
+    Y .*= -1
+    return Y
+end
+
 compose(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, p, q) = p * q
 
 function compose!(::MultiplicationGroupTrait, G::AbstractDecoratorManifold, x, p, q)
diff --git a/src/groups/power_group.jl b/src/groups/power_group.jl
index d454e7e63e..d1bf597f81 100644
--- a/src/groups/power_group.jl
+++ b/src/groups/power_group.jl
@@ -61,6 +61,29 @@ end
     end
 end
 
+function adjoint_inv_diff!(G::PowerGroup, Y, p, X)
+    GM = G.manifold
+    rep_size = representation_size(GM.manifold)
+    for i in get_iterator(GM)
+        adjoint_inv_diff!(
+            GM.manifold,
+            _write(GM, rep_size, Y, i),
+            _read(GM, rep_size, p, i),
+            _read(GM, rep_size, X, i),
+        )
+    end
+    return Y
+end
+function adjoint_inv_diff!(G::PowerGroupNestedReplacing, Y, p, X)
+    GM = G.manifold
+    N = GM.manifold
+    rep_size = representation_size(N)
+    for i in get_iterator(GM)
+        Y[i...] = adjoint_inv_diff(N, _read(GM, rep_size, p, i), _read(GM, rep_size, X, i))
+    end
+    return Y
+end
+
 function identity_element!(G::PowerGroup, p)
     GM = G.manifold
     N = GM.manifold
@@ -128,6 +151,29 @@ function inv!(
     return q
 end
 
+function inv_diff!(G::PowerGroup, Y, p, X)
+    GM = G.manifold
+    rep_size = representation_size(GM.manifold)
+    for i in get_iterator(GM)
+        inv_diff!(
+            GM.manifold,
+            _write(GM, rep_size, Y, i),
+            _read(GM, rep_size, p, i),
+            _read(GM, rep_size, X, i),
+        )
+    end
+    return Y
+end
+function inv_diff!(G::PowerGroupNestedReplacing, Y, p, X)
+    GM = G.manifold
+    N = GM.manifold
+    rep_size = representation_size(N)
+    for i in get_iterator(GM)
+        Y[i...] = inv_diff(N, _read(GM, rep_size, p, i), _read(GM, rep_size, X, i))
+    end
+    return Y
+end
+
 # lower level methods are added instead of top level ones to not have to deal
 # with `Identity` disambiguation
 
diff --git a/src/groups/product_group.jl b/src/groups/product_group.jl
index 571f108ed6..7565c5f925 100644
--- a/src/groups/product_group.jl
+++ b/src/groups/product_group.jl
@@ -40,6 +40,18 @@ end
     end
 end
 
+function adjoint_inv_diff!(G::ProductGroup, Y, p, X)
+    M = G.manifold
+    map(
+        adjoint_inv_diff!,
+        M.manifolds,
+        submanifold_components(G, Y),
+        submanifold_components(G, p),
+        submanifold_components(G, X),
+    )
+    return Y
+end
+
 function identity_element(G::ProductGroup)
     M = G.manifold
     return ArrayPartition(map(identity_element, M.manifolds))
@@ -102,6 +114,18 @@ function inv!(G::ProductGroup, q, p)
 end
 inv!(::ProductGroup, q::Identity{ProductOperation}, ::Identity{ProductOperation}) = q
 
+function inv_diff!(G::ProductGroup, Y, p, X)
+    M = G.manifold
+    map(
+        inv_diff!,
+        M.manifolds,
+        submanifold_components(G, Y),
+        submanifold_components(G, p),
+        submanifold_components(G, X),
+    )
+    return Y
+end
+
 _compose(G::ProductGroup, p, q) = _compose(G.manifold, p, q)
 function _compose(M::ProductManifold, p::ArrayPartition, q::ArrayPartition)
     return ArrayPartition(
diff --git a/src/groups/rotation_translation_action.jl b/src/groups/rotation_translation_action.jl
index ea706e3639..1883198be6 100644
--- a/src/groups/rotation_translation_action.jl
+++ b/src/groups/rotation_translation_action.jl
@@ -186,7 +186,8 @@ function apply_diff!(
     p,
     X,
 )
-    return mul!(Y, a.x[2], X)
+    mul!(Y, a.x[2], X)
+    return Y
 end
 function apply_diff!(
     ::RotationTranslationActionOnVector{LeftAction},
diff --git a/test/groups/circle_group.jl b/test/groups/circle_group.jl
index b11672ab6f..cd6ce39923 100644
--- a/test/groups/circle_group.jl
+++ b/test/groups/circle_group.jl
@@ -46,6 +46,8 @@ using Manifolds:
             test_invariance=true,
             test_lie_bracket=true,
             test_adjoint_action=true,
+            test_inv_diff=true,
+            test_adjoint_inv_diff=true,
         )
     end
 
@@ -65,6 +67,8 @@ using Manifolds:
             test_invariance=true,
             test_lie_bracket=true,
             test_adjoint_action=true,
+            test_inv_diff=true,
+            test_adjoint_inv_diff=true,
         )
     end
 
diff --git a/test/groups/general_linear.jl b/test/groups/general_linear.jl
index b3e76768fb..24fe8026d6 100644
--- a/test/groups/general_linear.jl
+++ b/test/groups/general_linear.jl
@@ -111,6 +111,7 @@ using NLsolve
                 test_invariance=true,
                 test_lie_bracket=true,
                 test_adjoint_action=true,
+                test_apply_diff_group=true,
             )
             test_manifold(
                 G,
diff --git a/test/groups/general_unitary_groups.jl b/test/groups/general_unitary_groups.jl
index d393c07819..7168d59651 100644
--- a/test/groups/general_unitary_groups.jl
+++ b/test/groups/general_unitary_groups.jl
@@ -19,7 +19,7 @@ include("group_utils.jl")
             X2 = log_lie(On, p)
             e = Identity(MultiplicationOperation())
             # They are not yet inverting, p2 is on the “other half”
-            # but taking the log again should also be X ahain
+            # but taking the log again should also be X again
             @test isapprox(On, e, X2, X)
             @test log_lie(On, e) == zeros(n, n)
         end
diff --git a/test/groups/group_operation_action.jl b/test/groups/group_operation_action.jl
index 3992a90c5c..f99542588b 100644
--- a/test/groups/group_operation_action.jl
+++ b/test/groups/group_operation_action.jl
@@ -111,18 +111,63 @@ using Manifolds:
         test_switch_direction=true,
     )
 
+    m = m_pts[1]
+    X = X_pts[1]
+    e = identity_element(G)
+
+    @testset "apply_diff" begin
+        @test isapprox(M, m, apply_diff(A_left_fwd, e, m, X), X)
+        @test isapprox(M, m, inverse_apply_diff(A_left_fwd, e, m, X), X)
+        @test isapprox(M, m, apply_diff(A_right_back, e, m, X), X)
+        @test isapprox(M, m, inverse_apply_diff(A_right_back, e, m, X), X)
+
+        @test isapprox(M, m, apply_diff(A_left_back, e, m, X), X)
+        @test isapprox(M, m, inverse_apply_diff(A_left_back, e, m, X), X)
+        @test isapprox(M, m, apply_diff(A_right_fwd, e, m, X), X)
+        @test isapprox(M, m, inverse_apply_diff(A_right_fwd, e, m, X), X)
+
+        eX = allocate(X)
+        @test apply_diff!(A_left_fwd, eX, e, m, X) === eX
+        @test isapprox(M, m, eX, X)
+        eX = allocate(X)
+        @test inverse_apply_diff!(A_left_fwd, eX, e, m, X) === eX
+        @test isapprox(M, m, eX, X)
+    end
+
     @testset "apply_diff_group" begin
-        @test apply_diff_group(A_left_fwd, a_pts[1], X_pts[1], m_pts[1]) ≈
-              translate_diff(G, m_pts[1], a_pts[1], X_pts[1], RightBackwardAction())
-        Y = similar(X_pts[1])
-        apply_diff_group!(A_left_fwd, Y, a_pts[1], X_pts[1], m_pts[1])
-        @test Y ≈ translate_diff(G, m_pts[1], a_pts[1], X_pts[1], RightBackwardAction())
-
-        @test adjoint_apply_diff_group(A_left_fwd, a_pts[1], X_pts[1], m_pts[1]) ≈
-              inverse_translate_diff(G, a_pts[1], m_pts[1], X_pts[1], RightBackwardAction())
-        Y = similar(X_pts[1])
-        adjoint_apply_diff_group!(A_left_fwd, Y, a_pts[1], X_pts[1], m_pts[1])
-        @test Y ≈
-              inverse_translate_diff(G, a_pts[1], m_pts[1], X_pts[1], RightBackwardAction())
+        @test apply_diff_group(A_left_fwd, a_pts[1], X, m) ≈
+              translate_diff(G, m, a_pts[1], X, RightBackwardAction())
+        Y = similar(X)
+        apply_diff_group!(A_left_fwd, Y, a_pts[1], X, m)
+        @test Y ≈ translate_diff(G, m, a_pts[1], X, RightBackwardAction())
+
+        @test adjoint_apply_diff_group(A_left_fwd, a_pts[1], X, m) ≈
+              inverse_translate_diff(G, a_pts[1], m, X, RightBackwardAction())
+        Y = similar(X)
+        adjoint_apply_diff_group!(A_left_fwd, Y, a_pts[1], X, m)
+        @test Y ≈ inverse_translate_diff(G, a_pts[1], m, X, RightBackwardAction())
+
+        @test adjoint_apply_diff_group(A_right_fwd, a_pts[1], X, m) ≈
+              inverse_translate_diff(
+            G,
+            m,
+            a_pts[1],
+            adjoint_inv_diff(G, a_pts[1], X),
+            RightBackwardAction(),
+        )
+
+        Y = similar(X)
+        adjoint_apply_diff_group!(A_right_fwd, Y, a_pts[1], X, m)
+        @test Y ≈ inverse_translate_diff(
+            G,
+            m,
+            a_pts[1],
+            -a_pts[1]' * X / a_pts[1]',
+            RightBackwardAction(),
+        )
+
+        @test apply_diff_group(A_right_fwd, a_pts[1], X, m) ≈
+              -m \ (a_pts[1] * X / a_pts[1]) * m
+        @test apply_diff_group(A_left_back, a_pts[1], X, m) ≈ -a_pts[1] * X / a_pts[1]
     end
 end
diff --git a/test/groups/groups_general.jl b/test/groups/groups_general.jl
index 39a445722c..6e3cbb4d8d 100644
--- a/test/groups/groups_general.jl
+++ b/test/groups/groups_general.jl
@@ -18,7 +18,10 @@ using Manifolds:
         @test repr(eg) === "Identity(NotImplementedOperation)"
         @test adjoint(eg) == eg
         @test number_eltype(eg) == Bool
+        @test !is_group_manifold(NotImplementedManifold())
         @test !is_group_manifold(NotImplementedManifold(), NotImplementedOperation())
+        @test !has_biinvariant_metric(NotImplementedManifold())
+        @test !has_invariant_metric(NotImplementedManifold(), LeftForwardAction())
         @test is_identity(G, eg) # identity transparent
         @test_throws MethodError identity_element(G) # but for a NotImplOp there is no concrete id.
         @test isapprox(G, eg, eg)
@@ -215,6 +218,7 @@ using Manifolds:
         @test one(ge) === ge
         @test transpose(ge) === ge
         @test det(ge) == 1
+        @test inv(ge) === ge
         @test ge * p ≈ p
         @test p * ge ≈ p
         @test ge * ge === ge
@@ -287,6 +291,35 @@ using Manifolds:
         @test e_mul * e_add === e_add
         @test mul!(e_mul, e_mul, e_mul) === e_mul
     end
+
+    @testset "Issue #669" begin
+        G = SpecialOrthogonal(3)
+
+        id = identity_element(G)
+        X = hat(G, Identity(G), [1.0, 2.0, 3.0])
+        function apply_at_id(X, d, s)
+            A = GroupOperationAction(G, (d, s))
+            return apply_diff_group(A, id, X, id)
+        end
+        function apply_at_id!(Y, X, d, s)
+            A = GroupOperationAction(G, (d, s))
+            return apply_diff_group!(A, Y, id, X, id)
+        end
+
+        _get_sign(::LeftAction, ::LeftSide) = 1 # former case
+        _get_sign(::LeftAction, ::RightSide) = -1 # new case
+        _get_sign(::RightAction, ::LeftSide) = -1 # new case
+        _get_sign(::RightAction, ::RightSide) = 1 # former case
+        Y = similar(X)
+
+        for d in [LeftAction(), RightAction()]
+            for s in [LeftSide(), RightSide()]
+                @test apply_at_id(X, d, s) ≈ _get_sign(d, s) * X
+                apply_at_id!(Y, X, d, s)
+                @test Y ≈ _get_sign(d, s) * X
+            end
+        end
+    end
 end
 
 struct NotImplementedAction <: AbstractGroupAction{LeftAction} end
diff --git a/test/groups/power_group.jl b/test/groups/power_group.jl
index 30525c8466..7277cbe3a2 100644
--- a/test/groups/power_group.jl
+++ b/test/groups/power_group.jl
@@ -36,6 +36,8 @@ include("group_utils.jl")
             test_exp_from_identity=true,
             test_vee_hat_from_identity=true,
             test_adjoint_action=true,
+            test_inv_diff=true,
+            test_adjoint_inv_diff=true,
         )
 
         X = log_lie(G, pts[1])
diff --git a/test/groups/product_group.jl b/test/groups/product_group.jl
index 9dc00ba148..05e251e1a6 100644
--- a/test/groups/product_group.jl
+++ b/test/groups/product_group.jl
@@ -51,6 +51,8 @@ using RecursiveArrayTools
         test_exp_from_identity=true,
         test_log_from_identity=true,
         test_vee_hat_from_identity=true,
+        test_inv_diff=true,
+        test_adjoint_inv_diff=true,
     )
     @test isapprox(
         G,
diff --git a/test/groups/rotation_translation_action.jl b/test/groups/rotation_translation_action.jl
index 7c2c8fb6a3..2bb61eb137 100644
--- a/test/groups/rotation_translation_action.jl
+++ b/test/groups/rotation_translation_action.jl
@@ -69,6 +69,30 @@ include("group_utils.jl")
             apply_diff_group!(A_left, Y, Identity(G), a_X_pts[1], m_pts[1])
             @test Y ≈ a_X_pts[1].x[2] * m_pts[1]
         end
+
+        @testset "apply_diff" begin
+            @test apply_diff(A_left, Identity(G), m_pts[1], X_pts[1]) ≈ X_pts[1]
+            Y = similar(X_pts[1])
+            apply_diff!(A_left, Y, Identity(G), m_pts[1], X_pts[1])
+            @test Y ≈ X_pts[1]
+
+            @test apply_diff(A_left, a_pts[1], m_pts[1], X_pts[1]) ≈
+                  a_pts[1].x[2] * X_pts[1]
+            Y = similar(X_pts[1])
+            apply_diff!(A_left, Y, a_pts[1], m_pts[1], X_pts[1])
+            @test Y ≈ a_pts[1].x[2] * X_pts[1]
+
+            @test apply_diff(A_right, Identity(G), m_pts[1], X_pts[1]) ≈ X_pts[1]
+            Y = similar(X_pts[1])
+            apply_diff!(A_right, Y, Identity(G), m_pts[1], X_pts[1])
+            @test Y ≈ X_pts[1]
+
+            @test apply_diff(A_right, a_pts[1], m_pts[1], X_pts[1]) ≈
+                  a_pts[1].x[2] \ X_pts[1]
+            Y = similar(X_pts[1])
+            apply_diff!(A_right, Y, a_pts[1], m_pts[1], X_pts[1])
+            @test Y ≈ a_pts[1].x[2] \ X_pts[1]
+        end
     end
 end
 
diff --git a/test/groups/special_orthogonal.jl b/test/groups/special_orthogonal.jl
index f7cc93c295..5ac62be630 100644
--- a/test/groups/special_orthogonal.jl
+++ b/test/groups/special_orthogonal.jl
@@ -59,6 +59,9 @@ using Manifolds: LeftForwardAction, RightBackwardAction
                 test_exp_from_identity=true,
                 test_log_from_identity=true,
                 test_vee_hat_from_identity=true,
+                test_inv_diff=true,
+                test_adjoint_inv_diff=true,
+                test_apply_diff_group=true,
             )
 
             @testset "log_lie edge cases" begin
diff --git a/test/groups/translation_group.jl b/test/groups/translation_group.jl
index e88704ff72..457f4bcaa0 100644
--- a/test/groups/translation_group.jl
+++ b/test/groups/translation_group.jl
@@ -42,6 +42,8 @@ using Manifolds: LeftForwardAction, RightBackwardAction
                 test_exp_from_identity=true,
                 test_log_from_identity=true,
                 test_vee_hat_from_identity=true,
+                test_inv_diff=true,
+                test_adjoint_inv_diff=true,
             )
         end
 
diff --git a/test/manifolds/embedded_torus.jl b/test/manifolds/embedded_torus.jl
index bee4aed68c..a5ccf63e36 100644
--- a/test/manifolds/embedded_torus.jl
+++ b/test/manifolds/embedded_torus.jl
@@ -92,8 +92,14 @@ using BoundaryValueDiffEq
     sol_log = Manifolds.solve_chart_log_bvp(M, p0x, a2, A, (0, 0))
     @test sol_log(0.0)[1:2] ≈ p0x
     @test sol_log(1.0)[1:2] ≈ a2
-    @test norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]) ≈
-          Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0))
+    # a test randomly failed here on Julia 1.6 once for no clear reason?
+    # so I bumped tolerance considerably
+    bvp_atol = VERSION < v"1.7" ? 2e-3 : 1e-15
+    @test isapprox(
+        norm(M, A, (0, 0), p0x, sol_log(0.0)[3:4]),
+        Manifolds.estimate_distance_from_bvp(M, p0x, a2, A, (0, 0));
+        atol=bvp_atol,
+    )
 
     @test Manifolds.IntegratorTerminatorNearChartBoundary().check_chart_switch_kwargs ===
           NamedTuple()