Checking if a list of vectors are all identical

General Usage Performance
1

(Looking for the most idiomatic way to solve this problem. It’s inside several nested loops.)

I have a situation where I have two lists of 3-vectors. Often these two lists are identical except for displacement. (One set of points is merely a shifted copy of the other.) To check for identity, I just subtract them and then check to see that the component-wise difference is the same for all vectors. So the difference between two sets of vectors might be something like this (just a MWE, not realistic)

julia> t
4-element Vector{Vector{Float64}}:
 [0.1, 0.2, 3.3]
 [0.2, -0.2, 1.0]
 [0.3, 0.25, -1.0]
 [-0.1, -0.2, 0.3]

julia> s
4-element Vector{Vector{Float64}}:
 [-0.1, 0.5, 4.6]
 [0.0, 0.1, 2.3]
 [0.1, 0.55, 0.3]
 [-0.3, 0.1, 1.6]

julia> diffs=t-s
4-element Vector{Vector{Float64}}:
 [0.2, -0.3, -1.2999999999999998]
 [0.2, -0.30000000000000004, -1.2999999999999998]
 [0.19999999999999998, -0.30000000000000004, -1.3]
 [0.19999999999999998, -0.30000000000000004, -1.3]

So then I check to see that each difference vector matches the others by:

julia> all([diffs[1]≈i for i in diffs[2:end]])

Not sure if this is the best way to do this. Is there a more idiomatic way? It seems verbose with the list comprehension…

(I could perhaps store the points as columns in a matrix instead of as a vector of vectors…)

2

If allocations/performance don’t matter, I propose

diffs = t-s
sum(sum.(diff(diffs))) ≈ 0

but not sure whether this achieves the maximal clarity :laughing:

3

only slightly better:

julia> diffs = t-s;

julia> all(≈(first(diffs)), diffs)
true
3 Likes
4

First of all, remove the brackets:

all(diffs[1]≈i for i in diffs[2:end])

to avoid (one of) the redundant allocation(s).

In general, though I would try to avoid subtracting all those vectors, don’t create diffs in the first place, instead iterate, comparing subtractions on the fly. Ideally, there should be zero allocations for this.

3 Likes
5

For example:

all(tᵢ - sᵢ ≈ t[1] - s[1] for (tᵢ, sᵢ) in zip(t,s))
6 
t .- [t[1]] ≈ s .- [s[1]]
1 Like
7

@rocco_sprmnt21, are you sure it is computing the same thing?
From the docs, it looks more like it uses the norm:
gif

1 Like
8

Using ≈ elementwise like this may be unreliable. Approximate comparison ≈ is not something you can use blindly — you are asking that the difference should be small, but you should always think “small compared to what”?

The problem with doing approximate comparison elementwise is that a ≈ b checks by default that |a - b| is small compared to |a| or |b|, and by default to about 8 digits. But what if t[1] - s[1] happens to be zero (or nearly zero)? Then tᵢ - sᵢ ≈ t[1] - s[1] will almost always be false, because nothing is small compared to zero (except for zero itself).

In your case, it’s a tricky question — what is the right “yardstick” to compare to? If you know something from your problem (e.g. some overall scale of your coordinate system), then you can use it to pass an absolute tolerance atol to isapprox (the same function as ≈). Otherwise, you would want to compute some relevant scale, e.g. the diameter of your point cloud.

(There’s also the delicate question of choosing an appropriate tolerance. atol = lengthscale * small, but what is small? √eps() ≈ 1.5e-8 is the default relative tolerance rtol for ≈, corresponding to half the significant digits, because this is a reasonable rule of thumb for unit tests: small enough to catch most bugs, but large enough to be resilient to roundoff errors. But if your application is not a unit test, your needed tolerance may be very different.)

I have a situation where I have two lists of 3-vectors

Working with this data structure will probably be much, much faster if you use arrays of SVector{3,Float64} from StaticArrays.jl. Then you can use natural vector operations, like subtracting two vectors, and it will be (a) allocation-free and (b) unrolled by the compiler. See the performance tips.

5 Likes
9

Thanks Steve. Would something like this be better? (*edited using: atol)

all(isapprox(tᵢ - sᵢ, t[1] - s[1], atol=1e-3) for (tᵢ, sᵢ) in zip(t,s))
10

How does:

all(s .+ Ref(mean(t-s)) .≈ t)

look?

11

No, because you are passing only a relative tolerance rtol. That means it will still probably fail if t[1] - s[1] happens to be nearly zero.

Yes, using atol is safer, but requires you to pick a reasonable atol, and that depends on the problem because atol is dimensionful.

For example, suppose your vectors represent positions of stars in the galaxy, with coordinates measured in meters. Then atol = 1e-3 (1mm) is probably a ridiculously small tolerance. On the other hand, if you are doing atomic physics with distances in meters, where the typical scale is angstroms, then atol = 1e-3 is a ridiculously large tolerance.

This is still likely to fail if any element of t happens to be zero or nearly zero.

Note, by the way, that ≈ works perfectly well on arrays, or arrays of arrays (anything with a norm), so using ≈ with all is generally unnecessary (and more susceptible to scaling problems). You can just do s .+ Ref(mean(t-s)) ≈ t (no dot on the ≈), which is more robust although it will still be problematic if norm(t) == 0 (maybe that is less likely … but ultimately you should still be thinking about an appropriate problem-dependent atol).

1 Like
12

I’m not sure I understand the question/observation. I definitely don’t understand the reference to the norm.
What I wrote (which is ultimately equivalent to the other proposed forms. I don’t know if from a numerical point of view it is more robust to Stevenj’s observations) is based on the fact that if the two vectors are equivalent under a translation, then translating them into so that the first element of both goes to the origin, the two vectors must “overlap” element by element

13

:person_facepalming: Yes, I didn’t fix the problem.
But:

using LinearAlgebra

all(s - t .≈ Ref(mean(s - t)))

might fix it: isapprox defaults to use norm when comparing vectors, and mean(s-t)'s norm is a good measuring yardstick for how small the dispersion of s - t should be (limited through rtol).

But I’m still not sure enough. Doing a more explicit and simpler test would avoid bugs.

14

Still seems wrong. The mean could easily be zero, no?

15

Your code is better than elementwise — as I wrote above, ≈ is generally better to do on whole arrays — the more values the better to gauge the correct overall scale. But the problem with your code is still the implicit relative tolerance on the ≈, which can make the approximate comparison fail unexpectedly.

For example, your code returns false in this case:

julia> s = [1, 1, 1]; t = [1, 1, 1 + 1e-15]
3-element Vector{Float64}:
 1.0
 1.0
 1.000000000000001

julia> t .- [t[1]] ≈ s .- [s[1]]
false

even though s ≈ t to 15 digits. (It is giving a correct answer to the question you asked! If you compare [0,0,1e-15] to [0,0,0], they don’t match in the third component to any significant digits. You have given isapprox no other yardstick by which to say that 1e-15 is small in this problem.)

To get a more reliable comparison, you really need to think about a problem-specific atol.

1 Like
16

Assuming that 1e-3 is a good yardstick for the OP’s data, then your preference would be to modify @rocco_sprmnt21’s snippet as follows?

isapprox(t .- [t[1]], s .- [s[1]], atol=1e-3)
17

Oh this looks pretty sweet and works even in some edge cases:

julia> diff(s) ≈ diff(t)
true
1 Like
18

If they’re all three vectors, then you might want to use a static length structure for them all… and then that lends itself nicely to being reinterpreted as a matrix… and then there’s Linear Algebra fun to be had with linear dependencies. :slight_smile:

19

I’m afraid not. The counterexample proposed by Stevengj also applies to this nice solution.
This leads me to advance the following conjecture: for every proposed solution, in the “numerical context” there is always a counterexample that makes it fail.

1 Like
20

Now I think I understand your previous observation. I was looking for an answer to the OP’s request for an “idiomatic solution”, but I didn’t take the numerical aspects into consideration at all (not voluntarily: just out of ignorance. I don’t have any experience with software for “real” uses).