Sketch Submersion and Canonical Metric on Stiefel. Fix a few docs.

docs/src/manifolds/probabilitysimplex.md

@@ -1,7 +1,7 @@
 # The probability simplex
 
 ```@autodocs
-Modules = [Manifolds]
+Modules = [Manifolds, Base]
 Pages = ["manifolds/ProbabilitySimplex.jl"]
 Order = [:type, :function]
 Private=false

docs/src/manifolds/symmetricpositivedefinite.md

@@ -13,7 +13,7 @@ The manifold can be equipped with different metrics
 ## Common and metric independent functions
 
 ```@autodocs
-Modules = [Manifolds]
+Modules = [Manifolds, ManifoldsBase]
 Pages = ["manifolds/SymmetricPositiveDefinite.jl"]
 Order = [:function]
 Public=true

docs/src/misc/notation.md

@@ -31,7 +31,6 @@ Within the documented functions, the utf8 symbols are used whenever possible, as
 | ``\cdot^\mathrm{H}`` | Hermitian or conjugate transposed for both complex or quaternion matrices| |
 | ``\operatorname{Hess} f(p)`` | (Riemannian) Hessian of function ``f \colon T_p\mathcal{M} \to T_p\mathcal M`` (i.e. the 1-1-tensor form) at ``p \in \mathcal{M}`` | | |
 | ``\nabla^2 f(p)`` | (Euclidean) Hessian of function ``f`` in the embedding | `H` | |
-
 | ``e`` | identity element of a group | |
 | ``I_k`` | identity matrix of size ``k\times k`` | |
 | ``k`` | indices | ``i,j`` | |

src/Manifolds.jl

@@ -484,7 +484,7 @@ Base.in(p, M::AbstractManifold; kwargs...) = is_point(M, p, false; kwargs...)
     X ∈ TangentSpaceAtPoint(M,p)
 
 Check whether `X` is a tangent vector from (in) the tangent space $T_p\mathcal M$, i.e.
-the [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`.
+the [`TangentSpaceAtPoint`](@ref) at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  `M`.
 This method uses [`is_vector`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/functions.html#ManifoldsBase.is_vector) deactivating the error throw option.
 """
 function Base.in(X, TpM::TangentSpaceAtPoint; kwargs...)

src/manifolds/FixedRankMatrices.jl

@@ -566,10 +566,15 @@ function retract_polar!(
 end
 
 @doc raw"""
-    riemannian_Hessian(M::FixedRankMatrices, p, g, H, X)
+    Y = riemannian_Hessian(M::FixedRankMatrices, p, G, H, X)
+    riemannian_Hessian!(M::FixedRankMatrices, Y, p, G, H, X)
+
+Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the
+Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`,
+where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,.
 
 The Riemannian Hessian can be computed as stated in Remark 4.1 [Nguyen:2023](@cite)
-or Section 2.3 [Vandereycken:2013](@cite), that B. Vandereycken adopted for Manopt (Matlab).
+or Section 2.3 [Vandereycken:2013](@cite), that B. Vandereycken adopted for [Manopt (Matlab)](https://www.manopt.org/reference/manopt/manifolds/fixedrank/fixedrankembeddedfactory.html).
 """
 riemannian_Hessian(M::FixedRankMatrices, p, G, H, X)
 function riemannian_Hessian(::FixedRankMatrices, Y, p, G, H, X)

src/manifolds/GrassmannStiefel.jl

@@ -298,7 +298,7 @@ end
     riemannian_Hessian(M::Grassmann, p, G, H, X)
 
 The Riemannian Hessian can be computed by adopting Eq. (6.6) [Nguyen:2023](@cite),
-where we use for the [`EuclideanMetric`](@ref) ``α_0=α_1=1`` in their formula.
+where we use for the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric) ``α_0=α_1=1`` in their formula.
 Let ``\nabla f(p)`` denote the Euclidean gradient `G`,
 ``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`. Then the formula reads
 

src/manifolds/HyperbolicHyperboloid.jl

@@ -434,18 +434,22 @@ function parallel_transport_to!(::Hyperbolic, Y, p, X, q)
 end
 
 @doc raw"""
-    riemannian_Hessian(M::Hyperbolic, p, G, H, X)
+    Y = riemannian_Hessian(M::Hyperbolic, p, G, H, X)
+    riemannian_Hessian!(M::Hyperbolic, Y, p, G, H, X)
 
-The Riemannian Hessian can be computed by adopting Remark 4.1 in [Nguyen:2023](@cite).
-Let ``\nabla f(p)`` denote the Euclidean gradient `G`,
-``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`, and
-``\mathbf{g} = \mathbf{g}^{-1} = \operatorname{diag}(1,...,1,-1)``.
-The formula reads
+Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the
+Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`,
+where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,.
+
+Let ``\mathbf{g} = \mathbf{g}^{-1} = \operatorname{diag}(1,...,1,-1)``.
+Then using Remark 4.1 [Nguyen:2023](@cite) the formula reads
 
 ```math
-    \operatorname{Hess}f(p)[X]
-    =
-    \operatorname{proj}_{T_p\mathcal M}\bigl( \mathbf{g}^-1\nabla^2f(p)[X] + X⟨p,\mathbf{g}^{-1}∇f(p)⟩_p\bigr).
+\operatorname{Hess}f(p)[X]
+=
+\operatorname{proj}_{T_p\mathcal M}\bigl(
+    \mathbf{g}^{-1}\nabla^2f(p)[X] + X⟨p,\mathbf{g}^{-1}∇f(p)⟩_p
+\bigr).
 ```
 """
 riemannian_Hessian(M::Hyperbolic, p, G, H, X)

src/manifolds/StiefelCanonicalMetric.jl

@@ -186,3 +186,52 @@ function inverse_retract!(
     @views mul!(X, qfact.U, LV[1:(2k), 1:k])
     return X
 end
+
+@doc raw"""
+    Y = riemannian_Hessian(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, G, H, X)
+    riemannian_Hessian!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, Y, p, G, H, X)
+
+Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the
+Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`,
+where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,.
+
+Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), for the [`CanonicalMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric)
+``α_0=1, α_1=\frac{1}{2}`` in their formula. The formula reads
+
+```math
+    \operatorname{Hess}f(p)[X]
+    =
+    \operatorname{proj}_{T_p\mathcal M}\Bigl(
+        ∇^2f(p)[X] - \frac{1}{2} X \bigl( (∇f(p))^{\mathrm{H}}p + p^{\mathrm{H}}∇f(p)\bigr)
+        - \frac{1}{2} \bigl( P ∇f(p) p^{\mathrm{H}} + p ∇f(p))^{\mathrm{H}} P)X
+    \Bigr),
+```
+where ``P = I-pp^{\mathrm{H}}``.
+
+Compared to Eq. (5.6) we have ``α_0 = 1`` and ``α_1 = \frac{1}{2}``.
+"""
+riemannian_Hessian(
+    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
+    p,
+    G,
+    H,
+    X,
+) where {n,k}
+
+function riemannian_Hessian!(
+    ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
+    Y,
+    p,
+    G,
+    H,
+    X,
+) where {n,k}
+    project!(
+        M,
+        Y,
+        p,
+        H - 1 / 2 .* X * (G' * p + p' * G) -
+        1 / 2 * ((I - p * p') * G * p' + p * G' * (I - p * p')) * X,
+    )
+    return Y
+end

src/manifolds/StiefelEuclideanMetric.jl

@@ -186,12 +186,15 @@ function retract_project!(M::Stiefel, q, p, X, t::Number)
 end
 
 @doc raw"""
-    riemannian_Hessian(M::Stiefel, p, G, H, X)
+    Y = riemannian_Hessian(M::Stiefel, p, G, H, X)
+    riemannian_Hessian!(M::Stiefel, Y, p, G, H, X)
 
-The Riemannian Hessian can be computed by adopting Eq. (5.6) [Nguyen:2023](@cite),
-where we use for the [`EuclideanMetric`](@ref) ``α_0=α_1=1`` in their formula.
-Let ``\nabla f(p)`` denote the Euclidean gradient `G`,
-``\nabla^2 f(p)[X]`` the Euclidean Hessian `H`. Then the formula reads
+Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the
+Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`,
+where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,.
+
+Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), where we use for the [`EuclideanMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric)
+``α_0=α_1=1`` in their formula. Then the formula reads
 
 ```math
     \operatorname{Hess}f(p)[X]

src/manifolds/StiefelSubmersionMetric.jl

@@ -257,6 +257,56 @@ function log!(
     return inverse_retract!(M, X, p, q, inverse_retraction)
 end
 
+@doc raw"""
+    Y = riemannian_Hessian(M::MetricManifold{ℝ,Stiefel{n,k,ℝ}, StiefelSubmersionMetric},, p, G, H, X)
+    riemannian_Hessian!(MetricManifold{ℝ,Stiefel{n,k,ℝ}, StiefelSubmersionMetric},, Y, p, G, H, X)
+
+Compute the Riemannian Hessian ``\operatorname{Hess} f(p)[X]`` given the
+Euclidean gradient ``∇ f(\tilde p)`` in `G` and the Euclidean Hessian ``∇^2 f(\tilde p)[\tilde X]`` in `H`,
+where ``\tilde p, \tilde X`` are the representations of ``p,X`` in the embedding,.
+
+Here, we adopt Eq. (5.6) [Nguyen:2023](@cite), for the [`CanonicalMetric`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds.html#ManifoldsBase.EuclideanMetric)
+``α_0=1, α_1=\frac{1}{2}`` in their formula. The formula reads
+
+```math
+    \operatorname{Hess}f(p)[X]
+    =
+    \operatorname{proj}_{T_p\mathcal M}\Bigl(
+        ∇^2f(p)[X] - \frac{1}{2} X \bigl( (∇f(p))^{\mathrm{H}}p + p^{\mathrm{H}}∇f(p)\bigr)
+        - \frac{2α+1}{2(α+1)} \bigl( P ∇f(p) p^{\mathrm{H}} + p ∇f(p))^{\mathrm{H}} P)X
+    \Bigr),
+```
+where ``P = I-pp^{\mathrm{H}}``.
+
+Compared to Eq. (5.6) we have that their ``α_0 = 1``and ``\alpha_1 =  \frac{2α+1}{2(α+1)} + 1``.
+"""
+riemannian_Hessian(
+    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric},
+    p,
+    G,
+    H,
+    X,
+) where {n,k}
+
+function riemannian_Hessian!(
+    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},<:StiefelSubmersionMetric},
+    Y,
+    p,
+    G,
+    H,
+    X,
+) where {n,k}
+    α = metric(M).α
+    project!(
+        M,
+        Y,
+        p,
+        H - 1 / 2 .* X * (G' * p + p' * G) -
+        (2 * α + 1) / (2 * (α + 1)) * ((I - p * p') * G * p' + p * G' * (I - p * p')) * X,
+    )
+    return Y
+end
+
 # StiefelFactorization code
 # Note: intended only for internal use
 

src/manifolds/VectorBundle.jl

@@ -85,7 +85,7 @@ TangentSpaceAtPoint(M::AbstractManifold, p) = VectorSpaceAtPoint(TangentBundleFi
 """
     TangentSpace(M::AbstractManifold, p)
 
-Return a [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) representing tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`.
+Return a [`TangentSpaceAtPoint`](@ref) representing tangent space at `p` on the [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold) `M`.
 """
 TangentSpace(M::AbstractManifold, p) = VectorSpaceAtPoint(TangentBundleFibers(M), p)
 
@@ -603,7 +603,7 @@ end
 @doc raw"""
     injectivity_radius(M::TangentSpaceAtPoint)
 
-Return the injectivity radius on the [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) `M`, which is $∞$.
+Return the injectivity radius on the [`TangentSpaceAtPoint`](@ref) `M`, which is $∞$.
 """
 injectivity_radius(::TangentSpaceAtPoint) = Inf
 
@@ -726,7 +726,7 @@ end
 """
     is_flat(::TangentSpaceAtPoint)
 
-Return true. [`TangentSpaceAtPoint`](@ref Manifolds.TangentSpaceAtPoint) is a flat manifold.
+Return true. [`TangentSpaceAtPoint`](@ref) is a flat manifold.
 """
 is_flat(::TangentSpaceAtPoint) = true
 """