Merge cbd8de9 into ab49d04

NEWS.md

@@ -5,6 +5,16 @@ 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-10-xx
+
+### Added
+
+- Functions `inv_diff` and `inv_diff!` that correspond to differentials of group inversion.
+
+### 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

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},

ext/ManifoldsTestExt/tests_group.jl

@@ -16,6 +16,7 @@ using Manifolds:
         test_invariance = false,
         test_lie_bracket=false,
         test_adjoint_action=false,
+        test_inv_diff=false,
         diff_convs = [(), (LeftForwardAction(),), (RightBackwardAction(),)],
     )
 
@@ -41,6 +42,7 @@ function test_group(
     test_invariance=false,
     test_lie_bracket=false,
     test_adjoint_action=false,
+    test_inv_diff=false,
     diff_convs=[(), (LeftForwardAction(),), (RightBackwardAction(),)],
     test_log_from_identity=false,
     test_exp_from_identity=false,
@@ -300,6 +302,15 @@ 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.@test isapprox(
+            inv_diff(G, inv(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))
@@ -710,15 +721,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 +735,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

src/Manifolds.jl

@@ -989,6 +989,8 @@ export adjoint_action,
     identity_element!,
     inv,
     inv!,
+    inv_diff,
+    inv_diff!,
     inverse_apply,
     inverse_apply!,
     inverse_apply_diff,

src/groups/addition_operation.jl

@@ -55,6 +55,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

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,

src/groups/group.jl

@@ -529,6 +529,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,

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,
+            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,
+            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)

src/groups/multiplication_operation.jl

@@ -23,6 +23,8 @@ 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
 
@@ -124,6 +126,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)

src/groups/power_group.jl

@@ -128,6 +128,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
 

src/groups/product_group.jl

@@ -102,6 +102,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(

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},