use symmetrize and symmetrize! where possible – and inplace of Y wher…

…e possible.

src/manifolds/Rotations.jl

@@ -397,7 +397,8 @@ to map it into the Lie algebra.
 """
 riemannian_Hessian(M::Rotations, p, G, H, X)
 function riemannian_Hessian!(::Rotations{N}, Y, p, G, H, X) where {N}
-    project!(SkewSymmetricMatrices(N), Y, p' * H - 1 / 2 .* X * (G' * p + p' * G))
+    symmetrize!(Y, G' * p)
+    project!(SkewSymmetricMatrices(N), Y, p' * H - X * Y)
     return Y
 end
 
@@ -418,7 +419,8 @@ The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by
 Weingarten(::Rotations, p, X, V)
 
 function Weingarten!(::Rotations, Y, p, X, V)
-    Y .= -1 / 2 * p * (V' * X - X' * V)
+    symmetrize!(Y, V' * X)
+    Y .= -p * Y
     return Y
 end
 

src/manifolds/StiefelCanonicalMetric.jl

@@ -214,12 +214,8 @@ function riemannian_Hessian!(
     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,
-    )
+    Gp = symmetrize(G' * p)
+    Z = symmetrize((I - p * p') * G * p')
+    project!(M, Y, p, H - X * Gp - Z * X)
     return Y
 end

src/manifolds/StiefelEuclideanMetric.jl

@@ -209,7 +209,8 @@ Compared to Eq. (5.6) also the metric conversion simplifies to the identity.
 riemannian_Hessian(M::Stiefel, p, G, H, X)
 
 function riemannian_Hessian!(M::Stiefel, Y, p, G, H, X)
-    project!(M, Y, p, H - 1 / 2 .* X * (G' * p + p' * G))
+    Z = symmetrize(G' * p)
+    project!(M, Y, p, H - X * Z)
     return Y
 end
 
@@ -232,6 +233,7 @@ The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by
 Weingarten(::Stiefel, p, X, V)
 
 function Weingarten!(::Stiefel, Y, p, X, V)
-    Y .= -X * p' * V - p * (X' * V + V' * X)
+    Z = symmetrize(X' * V)
+    Y .= -X * p' * V - p * Z
     return Y
 end

src/manifolds/SymmetricPositiveDefiniteAffineInvariant.jl

@@ -431,10 +431,10 @@ the symmetrization operator. Then the formula reads
 riemannian_Hessian(M::SymmetricPositiveDefinite, p, G, H, X)
 
 function riemannian_Hessian!(::SymmetricPositiveDefinite, Y, p, G, H, X)
-    # The following formulae even work for complex, H-symmetric positive definite matrices.
-    # the term within the seconds summands outer \sym:
-    Y .= X * 0.5 * (G' + G) * p
-    Y .= p * 1 / 2 * (H' + H) * p + 1 / 2 * (Y' + Y)
+    symmetrize!(Y, G)
+    symmetrize!(Y, X * Y * p)
+    Z = symmetrize(H)
+    Y .= p * Z * p + Y
     return Y
 end
 

src/manifolds/Unitary.jl

@@ -128,7 +128,8 @@ to map it into the Lie algebra.
 """
 riemannian_Hessian(M::UnitaryMatrices, p, G, H, X)
 function riemannian_Hessian!(M::UnitaryMatrices, Y, p, G, H, X)
-    project!(M, Y, p, p' * H - 1 / 2 .* X * (G' * p + p' * G))
+    symmetrize!(Y, G' * p)
+    project!(M, Y, p, p' * H - X * Y)
     return Y
 end
 
@@ -147,7 +148,8 @@ The formula is due to [AbsilMahonyTrumpf:2013](@cite) given by
 Weingarten(::UnitaryMatrices, p, X, V)
 
 function Weingarten!(::UnitaryMatrices, Y, p, X, V)
-    Y .= -1 / 2 * p * (V' * X - X' * V)
+    symmetrize!(Y, X' * V)
+    Y .= -p * Y
     return Y
 end
 

src/utils.jl

@@ -221,6 +221,25 @@ For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns
     return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...)
 end
 
+@doc raw"""
+    symmetrize!(Y,X)
+
+Given a quare matrix `X` compute `1/2 .* (X' + X)` in place of `Y`
+"""
+function symmetrize!(Y, X)
+    Y .= 1 / 2 .* (X' .+ X)
+    return Y
+end
+
+@doc raw"""
+    symmetrize(X)
+
+Given a quare matrix `X` compute `1/2 .* (X' + X)`.
+"""
+function symmetrize(X)
+    return 1 / 2 .* (X' .+ X)
+end
+
 @doc raw"""
     vec2skew!(X, v, k)