Introduce the Riemannian Hessian (#30)

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldDiff"
 uuid = "af67fdf4-a580-4b9f-bbec-742ef357defd"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -30,7 +30,7 @@ ChainRulesCore = "1"
 DoubleFloats = ">= 0.9.2"
 ForwardDiff = "0.10"
 Manifolds = "0.8.43"
-ManifoldsBase = "0.13.29, 0.14"
+ManifoldsBase = "0.14.10"
 RecursiveArrayTools = "2"
 Requires = "1"
 StaticArrays = "1"

docs/Project.toml

@@ -2,6 +2,7 @@
 ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -14,5 +15,5 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 [compat]
 Documenter = "0.27"
 ManifoldDiff = "0.3"
+Manifolds = "0.8"
 ManifoldsBase = "0.14"
-Manifolds = "0.8"

docs/make.jl

@@ -1,13 +1,52 @@
-using ForwardDiff, ReverseDiff, FiniteDifferences, Zygote
-using Manifolds, ManifoldsBase, ManifoldDiff, Documenter
+# (a) if docs is not the current active environment, switch to it
+# (from https://github.com/JuliaIO/HDF5.jl/pull/1020/) 
+if Base.active_project() != joinpath(@__DIR__, "Project.toml")
+    using Pkg
+    Pkg.activate(@__DIR__)
+    Pkg.develop(PackageSpec(; path = (@__DIR__) * "/../"))
+    Pkg.resolve()
+    Pkg.instantiate()
+end
 
+using ManifoldDiff, ManifoldsBase, ManifoldDiff
+using Documenter, DocumenterCitations
+using FiniteDiff, ForwardDiff, ReverseDiff, FiniteDifferences, Zygote
+
+bib = CitationBibliography(joinpath(@__DIR__, "src", "references.bib"); style = :alpha)
 makedocs(
-    # for development, we disable prettyurls
-    format = Documenter.HTML(prettyurls = false, assets = ["assets/favicon.ico"]),
-    modules = [ManifoldDiff],
+    bib;
+    format = Documenter.HTML(
+        prettyurls = get(ENV, "CI", nothing) == "true",
+        assets = ["assets/favicon.ico"],
+    ),
+    modules = [
+        ManifoldDiff,
+        isdefined(Base, :get_extension) ?
+        Base.get_extension(ManifoldDiff, :ManifoldDiffFiniteDiffExt) :
+        ManifoldDiff.ManifoldDiffFiniteDiffExt,
+        isdefined(Base, :get_extension) ?
+        Base.get_extension(ManifoldDiff, :ManifoldDiffFiniteDifferencesExt) :
+        ManifoldDiff.ManifoldDiffFiniteDifferencesExt,
+        isdefined(Base, :get_extension) ?
+        Base.get_extension(ManifoldDiff, :ManifoldDiffForwardDiffExt) :
+        ManifoldDiff.ManifoldDiffForwardDiffExt,
+        isdefined(Base, :get_extension) ?
+        Base.get_extension(ManifoldDiff, :ManifoldDiffReverseDiffExt) :
+        ManifoldDiff.ManifoldDiffReverseDiffExt,
+        isdefined(Base, :get_extension) ?
+        Base.get_extension(ManifoldDiff, :ManifoldDiffZygoteExt) :
+        ManifoldDiff.ManifoldDiffZygoteExt,
+    ],
     authors = "Seth Axen, Mateusz Baran, Ronny Bergmann, and contributors.",
     sitename = "ManifoldDiff.jl",
-    pages = ["Home" => "index.md", "Library" => "library.md"],
+    strict = Documenter.except(:autodocs_block),
+    pages = [
+        "Home" => "index.md",
+        "Backends" => "backends.md",
+        "Library of functions" => "library.md",
+        "Internals" => "internals.md",
+        "Literature" => "references.md",
+    ],
 )
 deploydocs(
     repo = "github.com/JuliaManifolds/ManifoldDiff.jl.git",

docs/src/backends.md

@@ -0,0 +1,38 @@
+# Differentiation backends
+
+```@docs
+set_default_differential_backend!
+default_differential_backend
+```
+
+## EmbeddedDiff
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["embedded_diff.jl"]
+Order = [:type, :function, :constant]
+```
+
+## ForwardDiff.jl
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["forward_diff.jl"]
+Order = [:type, :function, :constant]
+```
+
+## FiniteDiff.jl
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["finite_diff.jl"]
+Order = [:type, :function, :constant]
+```
+
+## FiniteDifferenes.jl
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["finite_differences.jl"]
+Order = [:type, :function, :constant]
+```

docs/src/internals.md

@@ -0,0 +1,11 @@
+## Internal functions
+
+```@docs
+ManifoldDiff.AbstractDiffBackend
+ManifoldDiff.CurrentDiffBackend
+ManifoldDiff._current_default_differential_backend
+ManifoldDiff._hessian
+ManifoldDiff._jacobian
+ManifoldDiff._gradient
+ManifoldDiff._derivative
+```

docs/src/library.md

@@ -15,14 +15,7 @@ Private = true
 
 ```@autodocs
 Modules = [ManifoldDiff]
-Pages = ["adjoint_differentials.jl"]
-Order = [:type, :function, :constant]
-Private = true
-```
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["differentials.jl"]
+Pages = ["adjoint_differentials.jl","differentials.jl"]
 Order = [:type, :function, :constant]
 Private = true
 ```
@@ -58,47 +51,3 @@ Modules = [ManifoldDiff]
 Pages = ["riemannian_diff.jl"]
 Order = [:type, :function, :constant]
 ```
-
-## Manifold-specific specializations
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["manifolds.jl"]
-Order = [:type, :function, :constant]
-```
-
-## Differentiation backends
-
-### EmbeddedDiff
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["embedded_diff.jl"]
-Order = [:type, :function, :constant]
-```
-
-### ForwardDiff.jl
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["forward_diff.jl"]
-Order = [:type, :function, :constant]
-```
-
-### FiniteDifferenes.jl
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["finite_differences.jl"]
-Order = [:type, :function, :constant]
-```
-
-## Internal functions
-
-```@autodocs
-Modules = [ManifoldDiff]
-Pages = ["ManifoldDiff.jl"]
-Order = [:type, :function, :constant]
-Private = true
-Public=false
-```

docs/src/references.bib

@@ -0,0 +1,23 @@
+@book{Boumal:2023,
+    TITLE     = {An Introduction to Optimization on Smooth Manifolds},
+    AUTHOR    = {Boumal, Nicolas},
+    YEAR      = {2023},
+    MONTH     = mar,
+    EDITION   = {First},
+    PUBLISHER = {Cambridge University Press},
+    DOI       = {10.1017/9781009166164},
+    URLDATE   = {2023-07-13},
+    ABSTRACT  = {Optimization on Riemannian manifolds-the result of smooth geometry and optimization merging into one elegant modern framework-spans many areas of science and engineering, including machine learning, computer vision, signal processing, dynamical systems and scientific computing. This text introduces the differential geometry and Riemannian geometry concepts that will help students and researchers in applied mathematics, computer science and engineering gain a firm mathematical grounding to use these tools confidently in their research. Its charts-last approach will prove more intuitive from an optimizer's viewpoint, and all definitions and theorems are motivated to build time-tested optimization algorithms. Starting from first principles, the text goes on to cover current research on topics including worst-case complexity and geodesic convexity. Readers will appreciate the tricks of the trade for conducting research and for numerical implementations sprinkled throughout the book.},
+    ISBN      = {978-1-00-916616-4}
+}
+
+@article{Zimmermann:2020,
+    AUTHOR        = {Zimmermann, Ralf},
+    TITLE         = {Hermite Interpolation and Data Processing Errors on Riemannian Matrix Manifolds},
+    JOURNAL       = {SIAM Journal on Scientific Computing},
+    VOLUME        = {42},
+    NUMBER        = {5},
+    PAGES         = {A2593-A2619},
+    YEAR          = {2020},
+    DOI           = {10.1137/19M1282878},
+}

docs/src/references.md

@@ -0,0 +1,4 @@
+# Literature
+
+```@bibliography
+```

src/ManifoldDiff.jl

@@ -14,6 +14,7 @@ using ManifoldsBase:
     TangentSpaceType,
     PowerManifoldNested,
     PowerManifoldNestedReplacing,
+    Weingarten,
     allocate_result,
     get_iterator,
     _write,
@@ -203,7 +204,7 @@ function __init__()
     # There is likely no way to set defaults without Requires.jl
     @require FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" begin
         if default_differential_backend() === NoneDiffBackend()
-            set_default_differential_backend!(FiniteDifferencesBackend())
+            set_default_differential_backend!(FiniteDifferencesBackend()) #This expects a method in the inner ()?
         end
     end
 
@@ -240,4 +241,6 @@ include("forward_diff.jl")
 include("reverse_diff.jl")
 include("zygote.jl")
 
+export riemannian_gradient, riemannian_gradient!, riemannian_Hessian, riemannian_Hessian!
+export set_default_differential_backend!, default_differential_backend
 end # module

src/differentials.jl

@@ -132,17 +132,12 @@ function differential_exp_argument_lie_approx! end
     )
 
 Approximate the differential of the inverse retraction `invretr` using a finite difference
-formula (see Eq. (16) in [^Zimmermann2019]):
+formula (see Eq. (16) in [Zimmermann:2020](@cite)
 ```math
 \frac{\operatorname{retr}^{-1}_q(\operatorname{retr}_p(hX)) - \operatorname{retr}^{-1}_q(\operatorname{retr}_p(-hX))}{2h}
 ```
 where ``h`` is the finite difference step `h`, ``\operatorname{retr}^{-1}`` is the inverse
 retraction `invretr` and ``\operatorname{retr}`` is the retraction `retr`.
-
-[^Zimmermann2019]:
-    > R. Zimmermann, “Hermite interpolation and data processing errors on Riemannian matrix
-    > manifolds,” arXiv:1908.05875 [cs, math], Sep. 2019,
-    > Available: http://arxiv.org/abs/1908.05875
 """
 function differential_inverse_retract_argument_fd_approx(
     M::AbstractManifold,

src/riemannian_diff.jl

@@ -335,3 +335,35 @@ function riemannian_gradient!(
     change_representer!(M, X, embedding_metric, p, X)
     return X
 end
+
+@doc raw"""
+    riemannian_Hessian(M, p, eG, eH, X)
+    riemannian_Hessian!(M, Y, p, eG, eH, X)
+
+Convert the Euclidean Hessian `eH=```\operatorname{Hess} \tilde f(p) [X]``
+of a function ``f \colon \mathcal M \to \mathbb R``, which is the restriction of ``\tilde f`` to ``\mathcal M``,
+given additionally the (Euclidean) gradient ``\operatorname{grad} \tilde f(p)``.
+
+The Riemannian Hessian is then computed by
+
+```math
+\operatorname{Hess} f(p)[X]
+= \operatorname{proj}_{T_p\mathcal M}\bigl(\operatorname{Hess} \tilde f(p)[X])
++ \mathcal W_p\Bigl( X, \operatorname{proj}_{N_p\mathcal M}\bigl( \operatorname{grad} \tilde f (p) \bigr) \Bigr),
+```
+
+where ``N_p\mathcal M`` denotes the normal space, i.e. the orthogonal complement of the tangent space in the embedding,
+and ``\mathcal W_p`` denotes the [`Weingarten`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions/#ManifoldsBase.Weingarten-Tuple{AbstractManifold,%20Any,%20Any,%20Any}) map. See [Boumal:2023](@cite) for more details
+
+The function is inspired by `ehess2rhess` in the [Matlab package Manopt](https://manopt.org).
+"""
+function riemannian_Hessian(M::AbstractManifold, p, eG, eH, X)
+    Y = zero_vector(M, p)
+    riemannian_Hessian!(M, Y, p, eG, eH, X)
+    return Y
+end
+function riemannian_Hessian!(M::AbstractManifold, Y, p, eG, eH, X)
+    project!(M, Y, p, eH) #first term - project the Euclidean Hessian
+    Y .+= Weingarten(M, p, X, eG - project(M, p, eG))
+    return Y
+end

test/test_gradients.jl

@@ -1,5 +1,6 @@
-using Manifolds
-using ManifoldsBase
+using Manifolds, ManifoldDiff, ManifoldsBase, Test
+import ManifoldsBase: Weingarten!
+Weingarten!(::Sphere, Y::Vector, p::Vector, X::Vector, V::Vector) = (Y .= -X * (p' * V))
 
 @testset "Gradients" begin
     M = Sphere(2)
@@ -27,4 +28,25 @@ using ManifoldsBase
         @test U == V
         @test U == W
     end
+    @testset "Riemannian Gradient / Hessian Conversion" begin
+        M = Sphere(2)
+        E = ℝ^3
+        A = [1.0 0.1 0.2; 0.1 2.0 0.3; 0.2 0.3 3.0]
+        ef(E, p) = 1 / 2 * p'A * p
+        grad_ef(E, p) = A * p
+        grad_f(M, p) = A * p - (p'A * p) * p
+        Hess_ef(E, p, X) = A * X
+        Hess_f(M, p, X) = A * X - (p' * A * X) .* p - (p' * A * p) .* X
+
+        p = [1.0, 0.0, 0.0]
+        X = [0.0, 0.1, 0.2]
+
+        Y1 = riemannian_gradient(M, p, grad_ef(E, p))
+        Y2 = grad_f(M, p)
+        @test isapprox(M, p, Y1, Y2)
+
+        Z1 = riemannian_Hessian(M, p, grad_ef(E, p), Hess_ef(M, p, X), X)
+        Z2 = Hess_f(M, p, X)
+        @test isapprox(M, p, Z1, Z2)
+    end
 end