Unify use of an atol when comparing to zero (#692)

* Unify atols in isapproxes. * fix dispatch route for direction vector transport in Circle. * Make the atol more generic/tolerant * update the literature entry for the Manifolds.jl paper --------- Co-authored-by: Mateusz Baran <[email protected]>
JuliaManifolds · Dec 25, 2023 · deede80 · deede80 · kellertuer · Dec 25, 2023
1 parent b2d83b3
commit deede80
Show file tree

Hide file tree

Showing 43 changed files with 558 additions and 310 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,6 +5,15 @@ 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.9] – 2023-12-24
+
+### Fixed
+
+* introduced a nonzero `atol` for all point and vector checks that compre to zero.
+  This makes those checks a bit more relaxed by default and resolves [#630](https://github.com/JuliaManifolds/Manifolds.jl/issues/630).
+* `default_estimation_method(M, f)` is deprecated, use `default_approximation_method(M, f)` for your specific method `f` on the manifold `M`.
+* `AbstractEstimationMethod` is deprecated, use `AbstractApproximationMethod` instead.
+
 ## [0.9.8] - 2023-11-17
 
 ### Fixed

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.9.8"
+version = "0.9.9"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -42,7 +42,7 @@ ManifoldsRecipesBaseExt = ["Colors", "RecipesBase"]
 ManifoldsTestExt = "Test"
 
 [compat]
-BoundaryValueDiffEq = "4, 5"
+BoundaryValueDiffEq = "4, 5.6.1"
 Colors = "0.12"
 Distributions = "0.22.6, 0.23, 0.24, 0.25"
 Einsum = "0.4"
@@ -51,7 +51,7 @@ HybridArrays = "0.4"
 Kronecker = "0.4, 0.5"
 LinearAlgebra = "1.6"
 ManifoldDiff = "0.3.7"
-ManifoldsBase = "0.15.0"
+ManifoldsBase = "0.15.6"
 Markdown = "1.6"
 MatrixEquations = "2.2"
 OrdinaryDiffEq = "6.31"

diff --git a/README.md b/README.md
@@ -53,17 +53,19 @@ If you have any questions regarding the Manifolds.jl ecosystem feel free to reac
 
 ## Citation
 
-If you use `Manifolds.jl` in your work, please cite the following
+If you use `Manifolds.jl` in your work, please cite the following open access article
 
 ```biblatex
 @article{AxenBaranBergmannRzecki:2023,
-    AUTHOR     = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki},
-    EPRINT     = {2021.08777},
-    EPRINTTYPE = {arXiv},
-    JOURNAL    = {AMS Transactions on Mathematical Software},
-    NOTE       = {accepted for publication},
-    TITLE      = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds},
-    YEAR       = {2023}
+    AUTHOR    = {Axen, Seth D. and Baran, Mateusz and Bergmann, Ronny and Rzecki, Krzysztof},
+    ARTICLENO = {33},
+    DOI       = {10.1145/3618296},
+    JOURNAL   = {ACM Transactions on Mathematical Software},
+    MONTH     = {dec},
+    NUMBER    = {4},
+    TITLE     = {Manifolds.Jl: An Extensible Julia Framework for Data Analysis on Manifolds},
+    VOLUME    = {49},
+    YEAR      = {2023}
 }
 ```
 

diff --git a/docs/src/references.bib b/docs/src/references.bib
@@ -89,13 +89,15 @@ @book{AyJostLeSchwachhoefer:2017
     PUBLISHER = {Springer Cham}
 }
 @article{AxenBaranBergmannRzecki:2023,
-    AUTHOR     = {Seth D. Axen and Mateusz Baran and Ronny Bergmann and Krzysztof Rzecki},
-    EPRINT     = {2021.08777},
-    EPRINTTYPE = {arXiv},
-    JOURNAL    = {AMS Transactions on Mathematical Software},
-    NOTE       = {accepted for publication},
-    TITLE      = {Manifolds.jl: An Extensible {J}ulia Framework for Data Analysis on Manifolds},
-    YEAR       = {2023}
+    AUTHOR    = {Axen, Seth D. and Baran, Mateusz and Bergmann, Ronny and Rzecki, Krzysztof},
+    ARTICLENO = {33},
+    DOI       = {10.1145/3618296},
+    JOURNAL   = {ACM Transactions on Mathematical Software},
+    MONTH     = {dec},
+    NUMBER    = {4},
+    TITLE     = {Manifolds.Jl: An Extensible Julia Framework for Data Analysis on Manifolds},
+    VOLUME    = {49},
+    YEAR      = {2023}
 }
 #
 # B

diff --git a/src/Manifolds.jl b/src/Manifolds.jl
@@ -48,6 +48,7 @@ import ManifoldsBase:
     check_vector,
     copy,
     copyto!,
+    default_approximation_method,
     default_inverse_retraction_method,
     default_retraction_method,
     default_vector_transport_method,
@@ -149,19 +150,14 @@ import ManifoldsBase:
     submanifold_component,
     submanifold_components,
     vector_space_dimension,
-    vector_transport_along, # just specified in Euclidean - the next 5 as well
-    vector_transport_along_diff,
-    vector_transport_along_project,
+    vector_transport_along,           # just specified in Euclidean - the next 5 as well
     vector_transport_along!,
-    vector_transport_along_diff!,
-    vector_transport_along_project!,
+    vector_transport_along_diff!,     # For consistency these are imported, but for now not
+    vector_transport_along_project!,  # overwritten with new definitons.
     vector_transport_direction,
-    vector_transport_direction_diff,
     vector_transport_direction!,
     vector_transport_direction_diff!,
     vector_transport_to,
-    vector_transport_to_diff,
-    vector_transport_to_project,
     vector_transport_to!,
     vector_transport_to_diff!,
     vector_transport_to_project!, # some overwrite layer 2
@@ -186,6 +182,9 @@ import ManifoldDiff:
     riemannian_Hessian,
     riemannian_Hessian!
 
+import Statistics: mean, mean!, median, median!, cov, var
+import StatsBase: mean_and_var
+
 using Base.Iterators: repeated
 using Distributions
 using Einsum: @einsum
@@ -198,6 +197,7 @@ using ManifoldsBase:
     ℝ,
     ℂ,
     ℍ,
+    AbstractApproximationMethod,
     AbstractBasis,
     AbstractDecoratorManifold,
     AbstractInverseRetractionMethod,
@@ -225,20 +225,26 @@ using ManifoldsBase:
     CotangentSpaceType,
     CoTFVector,
     CoTVector,
+    CyclicProximalPointEstimation,
     DefaultBasis,
     DefaultOrthogonalBasis,
     DefaultOrthonormalBasis,
     DefaultOrDiagonalizingBasis,
     DiagonalizingBasisData,
     DiagonalizingOrthonormalBasis,
     DifferentiatedRetractionVectorTransport,
+    EfficientEstimator,
     EmbeddedManifold,
     EmptyTrait,
     EuclideanMetric,
     ExponentialRetraction,
+    ExtrinsicEstimation,
     Fiber,
     FiberType,
     FVector,
+    GeodesicInterpolation,
+    GeodesicInterpolationWithinRadius,
+    GradientDescentEstimation,
     InverseProductRetraction,
     IsIsometricEmbeddedManifold,
     IsEmbeddedManifold,
@@ -295,6 +301,7 @@ using ManifoldsBase:
     VectorSpaceFiber,
     VectorSpaceType,
     VeeOrthogonalBasis,
+    WeiszfeldEstimation,
     @invoke_maker,
     _euclidean_basis_vector,
     combine_allocation_promotion_functions,
@@ -607,8 +614,9 @@ include("deprecated.jl")
 export test_manifold
 export test_group, test_action
 
-#
+# Abstract main types
 export CoTVector, AbstractManifold, AbstractManifoldPoint, TVector
+# Manifolds
 export AbstractSphere, AbstractProjectiveSpace
 export Euclidean,
     ArrayProjectiveSpace,
@@ -748,9 +756,10 @@ export AbstractInverseRetractionMethod,
     ShootingInverseRetraction,
     SoftmaxInverseRetraction
 # Estimation methods for median and mean
-export AbstractEstimationMethod,
+export AbstractApproximationMethod,
     GradientDescentEstimation,
     CyclicProximalPointEstimation,
+    EfficientEstimator,
     GeodesicInterpolation,
     GeodesicInterpolationWithinRadius,
     ExtrinsicEstimation
@@ -788,9 +797,10 @@ export ×,
     convert,
     complex_dot,
     decorated_manifold,
-    default_vector_transport_method,
+    default_approximation_method,
     default_inverse_retraction_method,
     default_retraction_method,
+    default_vector_transport_method,
     det_local_metric,
     differential_canonical_project,
     differential_canonical_project!,

diff --git a/src/deprecated.jl b/src/deprecated.jl
@@ -1 +1,6 @@
 
+@deprecate default_estimation_method(M::AbstractManifold, f) default_approximation_method(
+    M,
+    f,
+)
+@deprecate ExtrinsicEstimation() ExtrinsicEstimation(EfficientEstimator())
diff --git a/src/groups/addition_operation.jl b/src/groups/addition_operation.jl
@@ -83,8 +83,14 @@ function inv_diff!(::AdditionGroupTrait, G::AbstractDecoratorManifold, Y, p, X)
     return Y
 end
 
-function is_identity(::AdditionGroupTrait, G::AbstractDecoratorManifold, q; kwargs...)
-    return isapprox(G, q, zero(q); kwargs...)
+function is_identity(
+    ::AdditionGroupTrait,
+    G::AbstractDecoratorManifold,
+    q::T;
+    atol::Real=sqrt(prod(representation_size(G))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
+    return isapprox(G, q, zero(q); atol=atol, kwargs...)
 end
 # resolve ambiguities
 function is_identity(

diff --git a/src/groups/group_action.jl b/src/groups/group_action.jl
@@ -116,7 +116,7 @@ where element `a` is acting on `p`, with respect to the group element.
 
 Let ``\mathcal G`` be the group acting on manifold ``\mathcal M`` by the action `A`.
 The action is of element ``g ∈ \mathcal G`` on a point ``p ∈ \mathcal M``.
-The differential transforms vector `X` from the tangent space at `a ∈ \mathcal G`, 
+The differential transforms vector `X` from the tangent space at `a ∈ \mathcal G`,
 ``X ∈ T_a \mathcal G`` into a tangent space of the manifold ``\mathcal M``.
 When action on element `p` is written as ``\mathrm{d}τ^p``, with the specified left or right
 convention, the differential transforms vectors
@@ -193,7 +193,7 @@ end
         A::AbstractGroupAction,
         pts,
         p,
-        mean_method::AbstractEstimationMethod = GradientDescentEstimation(),
+        mean_method::AbstractApproximationMethod = GradientDescentEstimation(),
     )
 
 Calculate an action element ``a`` of action `A` that is the mean element of the orbit of `p`
@@ -210,7 +210,7 @@ function center_of_orbit(
     A::AbstractGroupAction,
     pts::AbstractVector,
     q,
-    mean_method::AbstractEstimationMethod=GradientDescentEstimation(),
+    mean_method::AbstractApproximationMethod=GradientDescentEstimation(),
 )
     alignments = map(p -> optimal_alignment(A, q, p), pts)
     return mean(base_group(A), alignments, mean_method)

diff --git a/src/groups/group_operation_action.jl b/src/groups/group_operation_action.jl
@@ -199,7 +199,7 @@ function center_of_orbit(
     A::GroupOperationAction,
     pts::AbstractVector,
     p,
-    mean_method::AbstractEstimationMethod,
+    mean_method::AbstractApproximationMethod,
 )
     μ = mean(A.group, pts, mean_method)
     return inverse_apply(switch_direction_and_side(A), p, μ)

diff --git a/src/groups/special_linear.jl b/src/groups/special_linear.jl
@@ -49,9 +49,15 @@ function check_point(G::SpecialLinear, p; kwargs...)
     return nothing
 end
 
-function check_vector(G::SpecialLinear, p, X; kwargs...)
+function check_vector(
+    G::SpecialLinear,
+    p,
+    X::T;
+    atol::Real=sqrt(prod(representation_size(G))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     trX = tr(inverse_translate_diff(G, p, p, X, LeftForwardAction()))
-    if !isapprox(trX, 0; kwargs...)
+    if !isapprox(trX, 0; atol=atol, kwargs...)
         return DomainError(
             trX,
             "The matrix $(X) does not lie in the tangent space of $(G) at $(p), since " *

diff --git a/src/manifolds/CenteredMatrices.jl b/src/manifolds/CenteredMatrices.jl
@@ -35,9 +35,14 @@ zero.
 
 The tolerance for the column sums of `p` can be set using `kwargs...`.
 """
-function check_point(M::CenteredMatrices, p; kwargs...)
+function check_point(
+    M::CenteredMatrices,
+    p::T;
+    atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     m, n = get_parameter(M.size)
-    if !isapprox(sum(p, dims=1), zeros(1, n); kwargs...)
+    if !isapprox(sum(p, dims=1), zeros(1, n); atol=atol, kwargs...)
         return DomainError(
             p,
             string(
@@ -56,9 +61,15 @@ Check whether `X` is a tangent vector to manifold point `p` on the
 sum to zero and its values are from the correct [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system).
 The tolerance for the column sums of `p` and `X` can be set using `kwargs...`.
 """
-function check_vector(M::CenteredMatrices, p, X; kwargs...)
+function check_vector(
+    M::CenteredMatrices,
+    p,
+    X::T;
+    atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     m, n = get_parameter(M.size)
-    if !isapprox(sum(X, dims=1), zeros(1, n); kwargs...)
+    if !isapprox(sum(X, dims=1), zeros(1, n); atol=atol, kwargs...)
         return DomainError(
             X,
             "The vector $(X) is not a tangent vector to $(p) on $(M), since its columns do not sum to zero.",

diff --git a/src/manifolds/CholeskySpace.jl b/src/manifolds/CholeskySpace.jl
@@ -28,10 +28,15 @@ it's size fits the manifold, it is a lower triangular matrix and has positive
 entries on the diagonal.
 The tolerance for the tests can be set using the `kwargs...`.
 """
-function check_point(M::CholeskySpace, p; kwargs...)
+function check_point(
+    M::CholeskySpace,
+    p::T;
+    atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     cks = check_size(M, p)
     cks === nothing || return cks
-    if !isapprox(norm(strictlyUpperTriangular(p)), 0.0; kwargs...)
+    if !isapprox(norm(strictlyUpperTriangular(p)), 0.0; atol=atol, kwargs...)
         return DomainError(
             norm(UpperTriangular(p) - Diagonal(p)),
             "The point $(p) does not lie on $(M), since it strictly upper triangular nonzero entries",
@@ -54,8 +59,14 @@ after [`check_point`](@ref)`(M,p)`, `X` has to have the same dimension as `p`
 and a symmetric matrix.
 The tolerance for the tests can be set using the `kwargs...`.
 """
-function check_vector(M::CholeskySpace, p, X; kwargs...)
-    if !isapprox(norm(strictlyUpperTriangular(X)), 0.0; kwargs...)
+function check_vector(
+    M::CholeskySpace,
+    p,
+    X;
+    atol::Real=sqrt(prod(representation_size(M)) * eps(float(eltype(p)))),
+    kwargs...,
+)
+    if !isapprox(norm(strictlyUpperTriangular(X)), 0.0; atol=atol, kwargs...)
         return DomainError(
             norm(UpperTriangular(X) - Diagonal(X)),
             "The matrix $(X) is not a tangent vector at $(p) (represented as an element of the Lie algebra) since it is not lower triangular.",

diff --git a/src/manifolds/Circle.jl b/src/manifolds/Circle.jl
@@ -78,8 +78,14 @@ check_vector(::Circle{ℝ}, ::Any...; ::Any...)
 function check_vector(M::Circle{ℝ}, p, X; kwargs...)
     return nothing
 end
-function check_vector(M::Circle{ℂ}, p, X; kwargs...)
-    if !isapprox(abs(complex_dot(p, X)), 0.0; kwargs...)
+function check_vector(
+    M::Circle{ℂ},
+    p,
+    X::T;
+    atol::Real=sqrt(eps(real(float(number_eltype(T))))),
+    kwargs...,
+) where {T}
+    if !isapprox(abs(complex_dot(p, X)), 0; atol=atol, kwargs...)
         return DomainError(
             abs(complex_dot(p, X)),
             "The value $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
@@ -160,7 +166,7 @@ function Base.exp(M::Circle{ℂ}, p::Number, X::Number, t::Number)
 end
 
 exp!(::Circle{ℝ}, q, p, X) = (q .= sym_rem(p + X))
-exp!(::Circle{ℝ}, q, p, X, t::Number) = (q .= sym_rem(p + t * X))
+exp!(::Circle{ℝ}, q, p, X, t::Number) = (q .= sym_rem(p[] + t * X[]))
 function exp!(M::Circle{ℂ}, q, p, X)
     θ = norm(M, p, X)
     q .= cos(θ) * p + usinc(θ) * X
@@ -538,6 +544,11 @@ function parallel_transport_to!(M::Circle{ℂ}, Y, p, X, q)
     return Y
 end
 
+# dispatch before allocation
+function _vector_transport_direction(M::Circle, p, X, d, ::ParallelTransport)
+    return parallel_transport_to(M, p, X, exp(M, p, d))
+end
+
 """
     volume_density(::Circle, p, X)