so as `isequal(+0.0im,-0.0im)`

.

I discovered this during a painful debug process.

I really want to understand the philosophy behind this designing, and I wish to get some advice on numerical stability. Thank you @everyone !

so as `isequal(+0.0im,-0.0im)`

.

I discovered this during a painful debug process.

I really want to understand the philosophy behind this designing, and I wish to get some advice on numerical stability. Thank you @everyone !

If you want `+0.0`

and `-0.0`

to be equal, use `==`

, which implements IEEE floating-point equality rules. `isequal`

is used more for hash tables and other purposes that require finer-grained distinctions.

(In Julia, there are three different equality predicates: `==`

, `===`

, and `isequal`

, to express slightly different notions of equality.)

7 Likes

As an aside, I learned recently from Valentin that `-0.0`

is the additive identity, while `0.0`

is not:

```
julia> mymulv1(a,b) = muladd(a,b,0.0)
mymulv1 (generic function with 1 method)
julia> mymulv2(a,b) = muladd(a,b,-0.0)
mymulv2 (generic function with 1 method)
julia> @code_native debuginfo=:none mymulv1(1.2, 3.4)
.text
vxorpd %xmm2, %xmm2, %xmm2
vfmadd213sd %xmm2, %xmm1, %xmm0 # xmm0 = (xmm1 * xmm0) + xmm2
retq
nopw (%rax,%rax)
julia> @code_native debuginfo=:none mymulv2(1.2, 3.4)
.text
vmulsd %xmm1, %xmm0, %xmm0
retq
nopw %cs:(%rax,%rax)
nop
```

So the compiler eliminates `-0.0`

, but not `0.0`

.

It would be nicer if it were the other way around, because `0.0`

is faster to create (via `xor`

instructions), making it the better choice if you canâ€™t trust it to be eliminated. Itâ€™s also what you get from `zero(::Float64)`

or `zero(::Type{Float64})`

.

Although the people came up with the IEEE rules thought about this a lot more than I have, so thereâ€™s likely a good reason I donâ€™t know about.

9 Likes

Thanks Steven!

In fact I was struggling with a version of `unique()`

with certain tolerance like `1e-16`

. I tracked the Julia implementation of `unique()`

in `base.jl`

line 1403:

`unique(A::AbstractArray; dims::Union{Colon,Integer} = :) = _unique_dims(A, dims)`

In the lines below that, I saw the actual code of `_unique_dims()`

. It seems to me that the implementation of `_unique_dims()`

is more â€śhashâ€ť.

My version of `unique()`

with tolerance is

`@inline uniquen(v,n) = unique( v .|> (x->round(x,digits=n)) )`

I simply round the digits beyond certain precision.

Do you think this would be stable enough?

PS:

I just realized that

`unique([-0.0,+0.0])`

is a `Vector{Float64}`

with TWO elements `[0.0,0.0]`

1 Like

are you from Mars?

1 Like

See the discussion here: https://github.com/JuliaLang/julia/issues/19147

A fundamental issue is that approximate equality is not an equivalence relation, so a `unique`

function with a tolerance would have to give implementation-dependent results, or rather its exact behavior would be so complicated to describe that it would effectively be defined by its implementation.

1 Like

thank you, Iâ€™ve learnt this way of thinking now!

What you might do is something like:

```
unique(x -> round(x, sigdigits=13), myarray)
```

which would perform the `unique`

operation based on the first 13 significant digits. (You could use `digits`

instead of `sigdigits`

if you want an absolute tolerance).

But this might not necessarily be what you want, e.g. `unique(x -> round(x, sigdigits=13), [1.000000000000501, 1.0000000000005])`

gives `[1.000000000000501, 1.0000000000005]`

because those two values round differently in the 13th decimal digit. (Again, the fact that approximate equality is not an equivalence relation â€” itâ€™s not transitive â€” means that *any* implementation of `unique`

with a tolerance will have some odd behaviors.)

4 Likes

I see â€¦ I will investigate further with actual examples. I am dealing with degenerate eigenvalues, whose behaviour depends largely on the input matrix.

To get rid of the `unique([-0.0,+0.0])`

trouble, I am currently using the following implementation:

`@inline correct_minus_00(x) = (real(x)==0.0 ? +0.0+imag(x)*im : x)`

`@inline uniquen(myarray,n) = unique(x -> correct_minus_00.(round(x, digits=n)), myarray)`

which adds a sign correction to your suggestion

Hi Elrod, do you have more elegant solution of the `-0.0`

problem? I mean, `-0.0`

is practically equal to `0.0`

in my code, so I use this

`@inline correct_minus_00(x) = (real(x)==0.0 ? +0.0+imag(x)*im : x)`

to correct the sign of `0.0`

. But it looks silly â€¦

If the sign of zero doesnâ€™t matter in your code, why do you need to â€ścorrectâ€ť it?

I am dealing with degenerate eigenvalues, whose behaviour depends largely on the input matrix.

I would very strongly recommend that you re-think your algorithms to not rely on detecting exactly degenerate eigenvalues â€¦ sounds numerically unstable. What are you actually trying to compute?

1 Like

I am trying to identify degenerate eigenvectors within a certain tolerance. These degenerate eigenvectors means quantum states of same energy, they are physically very interesting. I just need to identify and collect them, but my current algorithm is not very stable and generates wrong collections. I will flatten all our discussions in my mind â€¦

I have the following list of vectors

```
testv = [ Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[1.0+1.55431e-15im, 1.0-1.30451e-15im],
Complex{Float64}[0.0+1.0im, -0.0+1.0im, 0.0-1.0im, -0.0-1.0im],
Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im, -1.0-4.79e-16im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im],
Complex{Float64}[-1.0+1.22125e-15im, -1.0-1.35308e-15im],
Complex{Float64}[-1.0+0.0im, -1.0-0.0im], ]
```

within a tolerance of `1e-14`

, I hope to reduce them to distinct values. Because `unique([-0.0,0.0])`

gives two elements, I couldnâ€™t use it in an intuitive manner.

But Iâ€™ve got my solution now:

```
@inline uniquen(myarray,n) = unique(myarray .|> (x->(round(x,digits=n)+(0.0+0.0im))))
uniquen.(testv,12)
```

yields desired results.

BTW,

```
f1(myarray) = unique(x->round(x,digits=12), myarray)
```

and

```
f2(myarray) = unique(myarray .|> (x->round(x,digits=12)))
```

behaves differently. For example:

```
julia> tv = Complex{Float64}[-1.0+1.11022e-15im, -1.0-1.04777e-15im, -1.0-4.79e-16im] ;
julia> f1(tv)
2-element Array{Complex{Float64},1}:
-1.0 + 1.11022e-15im
-1.0 - 1.04777e-15im
julia> f2(tv)
2-element Array{Complex{Float64},1}:
-1.0 + 0.0im
-1.0 - 0.0im
```

I wished that `f1`

could give me the result of `f2`

(shown above), but it didnâ€™t.

If you are looking for non-accidental degeneracies (for a Hamiltonian that commutes with some symmetry group), then a more reliable way is to project the eigenvector onto irreps of the symmetry group.

3 Likes

yes I am doing that in the afternoon

Hereâ€™s my obligatory mention of one of my favorite-titled papers of all time, Branch cuts for complex elementary functions, or much ado about nothingâ€™s sign bit, by the estimable William Kahan.

6 Likes

Thanks for the paper, I realized that handling zeros is a highly non-trivial task indeed.

1 Like