Introduce defaults for atol and rtol for is_point and is_vector

[kellertuer]

This PR introduces a default atol (=0) and rtol (sqrt(eps)) for the high-level functions is_point and is_vector.

The lower level check_ functions do not have this detail, in case someone wants to be sure to not have it. Formerly I fear this is changing a default and hence breaking change – also integer vectors do not cope well with my idea of the eps-default so if there is a better idea, feel free to mention it here.

The PR yields that 1) high-level usage is easier to use and tolerant to some extend 2) low level use can still be as exact as before.

[codecov[bot]]

Codecov Report

Merging #94 (9570af5) into master (ded0ea4) will not change coverage. The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master      #94   +/-   ##
=======================================
  Coverage   99.80%   99.80%           
=======================================
  Files          15       15           
  Lines        1561     1561           
=======================================
  Hits         1558     1558           
  Misses          3        3           

Impacted Files	Coverage Δ
src/ManifoldsBase.jl	99.20% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update ded0ea4...9570af5. Read the comment docs.

[kellertuer]

I'm not sure what problem does this change actually solve? I've noticed that these defaults are exactly the same as in Julia's isapprox so all it seems to do is overriding potential manifold-specific tolerances introduce in check_vector and check_point?

Well for now they usually endet up being atol=rtol=0 which is often too strong. This introduces the same default as in Julia for our (M-prefixed) isapproxes for points and vectors.

[mateuszbaran]

Which code causes rtol to be 0? I've checked for Euclidean and Sphere and rtols are greater than 0:

julia> using Manifolds

julia> M = Euclidean(2)
Euclidean(2; field = ℝ)

julia> isapprox(M, [1, 1], [1, 1+1e-9])
true

julia> M = Sphere(1)
Sphere(1, ℝ)

julia> isapprox(M, [1, 1], [1, 1+1e-9])
true

[kellertuer]

Interesting, I am not sure, but it was the reason Grassmann failed in the slack example.

[mateuszbaran]

In this case I'd assume that just Grassmann should have less strict default tolerance.

[kellertuer]

Since I just had to set atop to 2e-14 I assume it has an rtol equal to zero, but I am not yet sure where this comes from. I also do not know why your examples do have a relative tolerance, since I do not see where our isapproxies do get that from.

[mateuszbaran]

Grassman example has a positive rtol but it doesn't matter that much when the two compared matrices are close to zero, for example:

julia> a = p'*vv
4×4 Matrix{Float64}:
 -2.45927e-16   2.48871e-17   9.14629e-17  -5.65537e-17
  1.89489e-16  -3.88603e-17  -2.96737e-16  -1.15517e-17
  2.34291e-16  -2.13559e-17  -9.71807e-17   4.89858e-17
  1.14101e-16  -1.48652e-17  -9.15055e-17   1.38072e-17

julia> b = vv'*p
4×4 Matrix{Float64}:
 -2.45927e-16   1.89489e-16   2.34291e-16   1.14101e-16
  2.48871e-17  -3.88603e-17  -2.13559e-17  -1.48652e-17
  9.14629e-17  -2.96737e-16  -9.71807e-17  -9.15055e-17
 -5.65537e-17  -1.15517e-17   4.89858e-17   1.38072e-17

See here: https://github.com/JuliaManifolds/Manifolds.jl/blob/f3e1ff5fe567d1df9dac327de22ae913dc8a20f1/src/manifolds/Grassmann.jl#L104 .

[kellertuer]

ah, so troll is fine but atol being zero is the only problem there then. Hm. But I do not want to set a nonzero atol actually. So should we leave it as is and close this PR?

[mateuszbaran]

Yes, I think this PR doesn't solve any problem. But maybe we need a nonzero default atol for that one check in Grassmann?

[kellertuer]

you mean as a default for check_point? Hm, I do not feel well with that, it breaks with the global default then? Let's keep it as is I would say.

[mateuszbaran]

Right, we probably shouldn't touch the default. One thing we can do is to write in documentation of is_vector or check_vector that for vectors close to the zero vector it may be a good idea to manually specify atol.

[kellertuer]

That sounds reasonable.

[kellertuer]

Hm, seems I am not clever enough for productRepr and such. I adapted vector_transport_direction to parallel_transport_to (and direction and both also in mutating). the non-mutating works, the mutating works if I call in the debugger from within the function the PT on a single component, but the map ends in the parallel_transport_to! method not to be found. NO clue why (since calling it actually separately in the debuger works).

My Power Manifold update might have broken essential. Hmpf. Maybe this still gets an infinite story.

[mateuszbaran]

I'll check that tomorrow or on Friday.

[kellertuer]

I at least fixed essential again – and improved error messages when checking sizes on PowerManifolds.

[JeffreySarnoff]

revised comment

This topic is one with which I have experience.

Geometric primitives that are defined/specified with floating-point values (particularly world coordinate values, which are not so well-behaved as e.g. barycentric coordinates with object relative displacements & orientations). When these primitives are used in low-level geometry constructs, the non-robustness may build.

If the primitives use more precise floats than does a user less well-versed in compsci or float issues, things tend to work out well. An example I return to often: let the api expect Float32s while the internal work is done using Float64s. Those Float32s can be as Float32(::Float64). Of course, this may not be appropriate (I do not know enough about manifolds.jl to say).

When geometrical entities are essential to the development and utilization of a package, less expert users lean heavily on expectations of geometric "crispness," connectivity and adjacency. An implementation that works well for them ought work well for all.

---

initial comment

This topic is one with which I have experience.

"never use floating-point tests for equality" does not go far enough, often. Comparing floats that obtain as vector norms, which may appear when determining the angle between two vectors, is unsafe as it is usually done AngleBetweenVectors.jl. Your earlier approach to the determination of point and line seems to satsifice rather than resolve or, better, to extend usefulness. Whenever geometrical entities are woven into a software design, the users tend to lean heavily on their own expectations of how these entities are to behave and interwork. Almost never is the user that is expert in some other domain going to consider tolerances beyond are there? is the default good for me .. hope so.

If you are going to use floaty geometricification, which is entirely rational from a performance perspective, maybe your work would often work well for others by providing an inner accuracy that exceeds the outer accuracy to be specified. If it is possible to off-load abs and rel accuracy examination from most [not all] casual uses of your work I have not examined manifolds, so dunno, you should find describing how to best work with these pkgs much simpler.

[kellertuer]

What is floaty geometrification? I somehow seem to miss what you want to tell us here.

[mateuszbaran]

This topic is one with which I have experience.

"never use floating-point tests for equality" does not go far enough, often. Comparing floats that obtain as vector norms, which may appear when determining the angle between two vectors, is unsafe as it is usually done AngleBetweenVectors.jl. Your earlier approach to the determination of point and line seems to satsifice rather than resolve or, better, to extend usefulness. Whenever geometrical entities are woven into a software design, the users tend to lean heavily on their own expectations of how these entities are to behave and interwork. Almost never is the user that is expert in some other domain going to consider tolerances beyond are there? is the default good for me .. hope so.

If you are going to use floaty geometricification, which is entirely rational from a performance perspective, maybe your work would often work well for others by providing an inner accuracy that exceeds the outer accuracy to be specified. If it is possible to off-load abs and rel accuracy examination from most [not all] casual uses of your work I have not examined manifolds, so dunno, you should find describing how to best work with these pkgs much simpler.

Thanks for the tip, we should actually use something like AngleBetweenVectors.jl here. The main problem that started this issue was checking whether isapprox(p' * X, -conj(X' * p)) where p and X are an n by k matrices, and while p is nice (its columns are orthonormal), X can have anything in it and we don't know what a good default check for this approximate equality could look like. In this comment: https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/94#issuecomment-1031592316 I have two matrices that are expected to be approximately equal.

[JeffreySarnoff]

What is floaty geometrification?

I revised the comment.
@mateuszbaran thanks for copying my first comment. 🎩

[JeffreySarnoff]

what do those two very similar matricies each mean, what do they each signify or encode and how do they differ in a non-numerical sense?

julia> a = p'*vv
julia> b = vv'*p

oic -- each is the transpose of the other, and they are populated with values that are ~eps(Float64) or a fraction thereof

[mateuszbaran]

I've replicated a bit different but similar in style example. We have a matrix p such that p'p is approximately a unit matrix:

julia> p = [-0.04705910362254517 -0.12156937692933076; 0.24137167076604518 -0.34985331139364395; 0.6096415977032226 0.7704681918979749; 0.7535663737639254 -0.5188470992242143]
4×2 Matrix{Float64}:
 -0.0470591  -0.121569
  0.241372   -0.349853
  0.609642    0.770468
  0.753566   -0.518847

julia> p'p
2×2 Matrix{Float64}:
  1.0          -4.63875e-18
 -4.63875e-18   1.0

Now we'd like to have a random matrix X such that p'*X = -X'*p (let's take the real numbers case here), so we take a random matrix and modify it to satisfy the equation:

julia> X = randn(4, 2)
4×2 Matrix{Float64}:
 -1.46691    -0.799284
 -0.301611   -0.109984
  0.23441     0.544457
  0.0145982  -1.39093

julia> X = X - p * p' * X
4×2 Matrix{Float64}:
 -1.4043    -0.677247
 -0.178008   0.506923
 -0.209134  -0.00938728
  0.138511  -0.197069

But numerically it doesn't exactly satisfy it:

julia> p'X+X'p
2×2 Matrix{Float64}:
 -3.11344e-17  -1.74839e-17
 -1.74839e-17   3.36018e-18

julia> isapprox(p'X, -X'p)
false

[JeffreySarnoff]

You should not be treating matrices the same way that scalars test for equality [approximate or determinate].

[kellertuer]

Thanks for al l your input. I still have not yet understood what the recommendation here now would be (besides not treating them the same – how would be the right way? Up to a little too small atol the current state on master worked well and the change I proposed here is already the default actually).

[kellertuer]

Maybe we keep the current state for now - at some point it might be nice to check for a few specific manifold and update our high.level functions and their tolerance.s