Unify use of an atol when comparing to zero

NEWS.md

@@ -5,6 +5,13 @@ 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-14
+
+### 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).
+
 ## [0.9.8] - 2023-11-17
 
 ### Fixed

Project.toml

@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.9.8"
+version = "0.9.9"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/groups/addition_operation.jl

@@ -83,8 +83,14 @@
     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=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(

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

src/groups/special_linear.jl

@@ -49,9 +49,15 @@
     return nothing
 end
 
-function check_vector(G::SpecialLinear, p, X; kwargs...)
+function check_vector(
+    G::SpecialLinear,
+    p,
+    X::T;
+    atol=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 " *

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=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=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.",

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=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=sqrt(prod(representation_size(M)) * eps(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.",

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=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.",
@@ -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)
 

src/manifolds/Elliptope.jl

@@ -85,10 +85,16 @@ zero diagonal.
 The tolerance for the base point check and zero diagonal can be set using the `kwargs...`.
 Note that symmetric of $X$ holds by construction an is not explicitly checked.
 """
-function check_vector(M::Elliptope, q, Y; kwargs...)
+function check_vector(
+    M::Elliptope,
+    q,
+    Y::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     X = q * Y' + Y * q'
     n = diag(X)
-    if !all(isapprox.(n, 0.0; kwargs...))
+    if !all(isapprox.(n, 0.0; atol=atol, kwargs...))
         return DomainError(
             n,
             "The vector $(X) is not a tangent to a point on $(M) (represented py $(q) and $(Y), since its diagonal is nonzero.",

src/manifolds/EmbeddedTorus.jl

@@ -52,9 +52,9 @@
 The method checks if the vector `X` is orthogonal to the vector normal to the torus,
 see [`normal_vector`](@ref). Absolute tolerance can be set using `atol`.
 """
-function check_vector(M::EmbeddedTorus, p, X; atol=eps(eltype(p)), kwargs...)
+function check_vector(M::EmbeddedTorus, p, X; atol=eps(float(eltype(p))), kwargs...)
     dot_nX = dot(normal_vector(M, p), X)
-    if !isapprox(dot_nX, 0; atol, kwargs...)
+    if !isapprox(dot_nX, 0; atol=atol, kwargs...)
         return DomainError(dot_nX, "The vector $(X) is not tangent to $(p) from $(M).")
     end
     return nothing

src/manifolds/EssentialManifold.jl

@@ -125,8 +125,14 @@ exp(::EssentialManifold, ::Any...)
 
 get_iterator(::EssentialManifold) = Base.OneTo(2)
 
-function _isapprox(M::EssentialManifold, p, q; kwargs...)
-    return isapprox(distance(M, p, q), 0.0; kwargs...)
+function _isapprox(
+    M::EssentialManifold,
+    p,
+    q::T;
+    atol=eps(real(float(number_eltype(number_eltype(T))))),
+    kwargs...,
+) where {T}
+    return isapprox(distance(M, p, q), 0.0; atol=atol, kwargs...)
 end
 
 """

src/manifolds/FixedRankMatrices.jl

@@ -324,15 +324,21 @@ Check whether the tangent [`UMVTVector`](@ref) `X` is from the tangent space of
 [`FixedRankMatrices`](@ref) `M`, i.e. that `v.U` and `v.Vt` are (columnwise) orthogonal to `x.U` and `x.Vt`,
 respectively, and its dimensions are consistent with `p` and `X.M`, i.e. correspond to `m`-by-`n` matrices of rank `k`.
 """
-function check_vector(M::FixedRankMatrices, p::SVDMPoint, X::UMVTVector; kwargs...)
+function check_vector(
+    M::FixedRankMatrices,
+    p::SVDMPoint,
+    X::UMVTVector;
+    atol=sqrt(prod(representation_size(M)) * eps(float(eltype(p.U)))),
+    kwargs...,
+)
     m, n, k = get_parameter(M.size)
-    if !isapprox(X.U' * p.U, zeros(k, k); kwargs...)
+    if !isapprox(X.U' * p.U, zeros(k, k); atol=atol, kwargs...)
         return DomainError(
             norm(X.U' * p.U - zeros(k, k)),
             "The tangent vector $(X) is not a tangent vector to $(p) on $(M) since v.U'x.U is not zero. ",
         )
     end
-    if !isapprox(X.Vt * p.Vt', zeros(k, k); kwargs...)
+    if !isapprox(X.Vt * p.Vt', zeros(k, k); atol=atol, kwargs...)
         return DomainError(
             norm(X.Vt * p.Vt - zeros(k, k)),
             "The tangent vector $(X) is not a tangent vector to $(p) on $(M) since v.V'x.V is not zero.",

src/manifolds/HyperbolicHyperboloid.jl

@@ -35,8 +35,14 @@ function check_point(M::Hyperbolic, p; kwargs...)
     return nothing
 end
 
-function check_vector(M::Hyperbolic, p, X; kwargs...)
-    if !isapprox(minkowski_metric(p, X), 0.0; kwargs...)
+function check_vector(
+    M::Hyperbolic,
+    p,
+    X::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
+    if !isapprox(minkowski_metric(p, X), 0; atol=atol, kwargs...)
         return DomainError(
             abs(minkowski_metric(p, X)),
             "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal (with respect to the Minkowski inner product) in the embedding.",

src/manifolds/KendallsPreShapeSpace.jl

@@ -37,9 +37,14 @@ representation_size(M::KendallsPreShapeSpace) = get_parameter(M.size)
 Check whether `p` is a valid point on [`KendallsPreShapeSpace`](@ref), i.e. whether
 each row has zero mean. Other conditions are checked via embedding in [`ArraySphere`](@ref).
 """
-function check_point(M::KendallsPreShapeSpace, p; atol=sqrt(eps(eltype(p))), kwargs...)
+function check_point(
+    M::KendallsPreShapeSpace,
+    p;
+    atol=sqrt(eps(float(eltype(p)))),
+    kwargs...,
+)
     for p_row in eachrow(p)
-        if !isapprox(mean(p_row), 0; atol, kwargs...)
+        if !isapprox(mean(p_row), 0; atol=atol, kwargs...)
             return DomainError(
                 mean(p_row),
                 "The point $(p) does not lie on the $(M) since one of the rows does not have zero mean.",
@@ -55,9 +60,15 @@ end
 Check whether `X` is a valid tangent vector on [`KendallsPreShapeSpace`](@ref), i.e. whether
 each row has zero mean. Other conditions are checked via embedding in [`ArraySphere`](@ref).
 """
-function check_vector(M::KendallsPreShapeSpace, p, X; atol=sqrt(eps(eltype(X))), kwargs...)
+function check_vector(
+    M::KendallsPreShapeSpace,
+    p,
+    X;
+    atol=sqrt(eps(float(eltype(X)))),
+    kwargs...,
+)
     for X_row in eachrow(X)
-        if !isapprox(mean(X_row), 0; atol, kwargs...)
+        if !isapprox(mean(X_row), 0; atol=atol, kwargs...)
             return DomainError(
                 mean(X_row),
                 "The vector $(X) is not a tangent vector to $(p) on $(M), since one of the rows does not have zero mean.",

src/manifolds/MultinomialDoublyStochastic.jl

@@ -60,10 +60,15 @@ end
 Checks whether `p` is a valid point on the [`MultinomialDoubleStochastic`](@ref)`(n)` `M`,
 i.e. is a  matrix with positive entries whose rows and columns sum to one.
 """
-function check_point(M::MultinomialDoubleStochastic, p; kwargs...)
+function check_point(
+    M::MultinomialDoubleStochastic,
+    p::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     n = get_parameter(M.size)[1]
     r = sum(p, dims=2)
-    if !isapprox(norm(r - ones(n, 1)), 0.0; kwargs...)
+    if !isapprox(norm(r - ones(n, 1)), 0.0; atol=atol, kwargs...)
         return DomainError(
             r,
             "The point $(p) does not lie on $M, since its rows do not sum up to one.",
@@ -78,9 +83,15 @@ Checks whether `X` is a valid tangent vector to `p` on the [`MultinomialDoubleSt
 This means, that `p` is valid, that `X` is of correct dimension and sums to zero along any
 column or row.
 """
-function check_vector(M::MultinomialDoubleStochastic, p, X; kwargs...)
+function check_vector(
+    M::MultinomialDoubleStochastic,
+    p,
+    X::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     r = sum(X, dims=2) # check for stochastic rows
-    if !isapprox(norm(r), 0.0; kwargs...)
+    if !isapprox(norm(r), 0.0; atol=atol, kwargs...)
         return DomainError(
             r,
             "The matrix $(X) is not a tangent vector to $(p) on $(M), since its rows do not sum up to zero.",

src/manifolds/ProbabilitySimplex.jl

@@ -129,8 +129,14 @@ after [`check_point`](@ref check_point(::ProbabilitySimplex, ::Any))`(M,p)`,
 `X` has to be of same dimension as `p` and its elements have to sum to one.
 The tolerance for the last test can be set using the `kwargs...`.
 """
-function check_vector(M::ProbabilitySimplex, p, X; kwargs...)
-    if !isapprox(sum(X), 0.0; kwargs...)
+function check_vector(
+    M::ProbabilitySimplex,
+    p,
+    X::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
+    if !isapprox(sum(X), 0.0; atol=atol, kwargs...)
         return DomainError(
             sum(X),
             "The vector $(X) is not a tangent vector to $(p) on $(M), since its elements do not sum up to 0.",

src/manifolds/ProjectiveSpace.jl

@@ -126,8 +126,14 @@ Check whether `X` is a tangent vector in the tangent space of `p` on the
 tangent space of the embedding and that the Frobenius inner product
 $⟨p, X⟩_{\mathrm{F}} = 0$.
 """
-function check_vector(M::AbstractProjectiveSpace, p, X; kwargs...)
-    if !isapprox(dot(p, X), 0; kwargs...)
+function check_vector(
+    M::AbstractProjectiveSpace,
+    p,
+    X::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
+    if !isapprox(dot(p, X), 0; atol=atol, kwargs...)
         return DomainError(
             dot(p, X),
             "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not" *

src/manifolds/SPDFixedDeterminant.jl

@@ -82,8 +82,14 @@ and additionally fulfill ``\operatorname{tr}(X) = 0``.
 
 The tolerance for the trace check of `X` can be set using `kwargs...`, which influences the `isapprox`-check.
 """
-function check_vector(M::SPDFixedDeterminant, p, X; kwargs...)
-    if !isapprox(tr(X), 0.0; kwargs...)
+function check_vector(
+    M::SPDFixedDeterminant,
+    p,
+    X::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
+    if !isapprox(tr(X), 0; atol=atol, kwargs...)
         return DomainError(
             tr(X),
             "The vector $(X) is not a tangent vector to $(p) on $(M), since it does not have a zero trace.",

src/manifolds/Spectrahedron.jl

@@ -85,10 +85,16 @@ and a $X$ has to be a symmetric matrix with trace.
 The tolerance for the base point check and zero diagonal can be set using the `kwargs...`.
 Note that symmetry of $X$ holds by construction and is not explicitly checked.
 """
-function check_vector(M::Spectrahedron, q, Y; kwargs...)
+function check_vector(
+    M::Spectrahedron,
+    q,
+    Y::T;
+    atol=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))),
+    kwargs...,
+) where {T}
     X = q * Y' + Y * q'
     n = tr(X)
-    if !isapprox(n, 0.0; kwargs...)
+    if !isapprox(n, 0; atol=atol, kwargs...)
         return DomainError(
             n,
             "The vector $(X) is not a tangent to a point on $(M) (represented py $(q) and $(Y), since its trace is nonzero.",