From d30c34a3a5e756e083e0164a67a9c00094910aa7 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 18 Aug 2023 14:29:12 +0200 Subject: [PATCH] Introduce the Riemannian Hessian (#30) --- Project.toml | 4 ++-- docs/Project.toml | 3 ++- docs/make.jl | 51 ++++++++++++++++++++++++++++++++++----- docs/src/backends.md | 38 +++++++++++++++++++++++++++++ docs/src/internals.md | 11 +++++++++ docs/src/library.md | 53 +---------------------------------------- docs/src/references.bib | 23 ++++++++++++++++++ docs/src/references.md | 4 ++++ src/ManifoldDiff.jl | 5 +++- src/differentials.jl | 7 +----- src/riemannian_diff.jl | 32 +++++++++++++++++++++++++ test/test_gradients.jl | 26 ++++++++++++++++++-- 12 files changed, 187 insertions(+), 70 deletions(-) create mode 100644 docs/src/backends.md create mode 100644 docs/src/internals.md create mode 100644 docs/src/references.bib create mode 100644 docs/src/references.md diff --git a/Project.toml b/Project.toml index 3ea51c4..b0331bf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ManifoldDiff" uuid = "af67fdf4-a580-4b9f-bbec-742ef357defd" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann "] -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" diff --git a/docs/Project.toml b/docs/Project.toml index 6170d36..bb62070 100644 --- a/docs/Project.toml +++ b/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" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 05a2103..adee529 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/backends.md b/docs/src/backends.md new file mode 100644 index 0000000..d60253c --- /dev/null +++ b/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] +``` \ No newline at end of file diff --git a/docs/src/internals.md b/docs/src/internals.md new file mode 100644 index 0000000..74702cf --- /dev/null +++ b/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 +``` \ No newline at end of file diff --git a/docs/src/library.md b/docs/src/library.md index b544b37..529efcf 100644 --- a/docs/src/library.md +++ b/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 -``` diff --git a/docs/src/references.bib b/docs/src/references.bib new file mode 100644 index 0000000..48400b9 --- /dev/null +++ b/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}, +} \ No newline at end of file diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 0000000..1bfafad --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,4 @@ +# Literature + +```@bibliography +``` \ No newline at end of file diff --git a/src/ManifoldDiff.jl b/src/ManifoldDiff.jl index 6307d36..44f8386 100644 --- a/src/ManifoldDiff.jl +++ b/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 diff --git a/src/differentials.jl b/src/differentials.jl index d80464a..61ac604 100644 --- a/src/differentials.jl +++ b/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, diff --git a/src/riemannian_diff.jl b/src/riemannian_diff.jl index 82e89ef..6ea1391 100644 --- a/src/riemannian_diff.jl +++ b/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 diff --git a/test/test_gradients.jl b/test/test_gradients.jl index 55480a5..261a8ca 100644 --- a/test/test_gradients.jl +++ b/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