CODEWITHSUNDEEP
pythonnumpy
ale
ale

Reputation: 11820

Determine whether two complex numbers are equal

The following code causes the print statements to be executed:

import numpy as np
import math

foo = np.array([1/math.sqrt(2), 1/math.sqrt(2)], dtype=np.complex_)

total = complex(0, 0)
one = complex(1, 0)
for f in foo:
   total = total + pow(np.abs(f), 2)
   if(total != one):
      print str(total) + " vs " + str(one)
      print "NOT EQUAL"

However, my input of [1/math.sqrt(2), 1/math.sqrt(2)] results in the total being one:

(1+0j) vs (1+0j) NOT EQUAL

Is it something to do with mixing NumPy with Python's complex type?

Upvotes: 1

Views: 2036

Answers (3)

Antony Hatchkins
Antony Hatchkins

Reputation: 33994

TLDR:

The proper way of comparing either floats or complex numbers is:

def isclose(a, b, rtol=1e-5, atol=1e-8):
    return abs(a-b) < atol + rtol * abs(b)

which is essentially what np.isclose() does under the hood, with np.isclose() being preferable as it takes care about infinities, not-a-numbers, etc.

Details:

The particular case from the question is not specific to complex numbers. If you replace

total = complex(0, 0)
one = complex(1, 0)

with its floating point equivalent

total = 0
one = 1

you'll get exactly the same result:

0.4999999999999999 vs 1
NOT EQUAL
0.9999999999999998 vs 1
NOT EQUAL

The criterion abs(a-b) < epsilon suggested by @bereal works in some cases, but if you look at the error in different scales:

>>> np.sqrt(1.234)**2-1.234
2.220446049250313e-16
>>> np.sqrt(12.34)**2-12.34
1.7763568394002505e-15
>>> np.sqrt(123.4)**2-123.4
-1.4210854715202004e-14

you'll see that it linearly increases for large numbers (which is not really surprising as floating point only keeps 15 decimal digits for float64), so it would be more reasonable to use a relative difference rather than the absolute one:

abs(a-b)/a < epsilon

It resolves an issue of the increasing error, but even a quick glance reveals problems with a==0. And no matter which flavor you prefer:
  • abs(a-b)/max(a,b) (used in math.isclose()),
  • abs(a-b)/(a+b)
or what not, they all have the same problem: they fail when a and/or b is zero.

To address this issue, it is a common practice to have two 'epsilons': an absolute epsilon (aka atol, absolute tolerance) for comparing with zero and a relative one (aka rtol, relative tolerance) for comparing with anything else.

As a bottomline:
  • np.isclose() does abs(a-b) <= atol + rtol * abs(b)
  • math.isclose() does abs(a-b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)

The math version is more symmetric, but the numpy version is faster.

Upvotes: 0

sebix
sebix

Reputation: 3240

When using floating point numbers it is important to keep in mind that working with these numbers is never accurate and thus computations are every time subject to rounding errors. This is caused by the design of floating point arithmetic and currently the most practicable way to do high arbitrary precision mathematics on computers with limited resources. You can't compute exactly using floats (means you have practically no alternative), as your numbers have to be cut off somewhere to fit in a reasonable amount of memory (in most cases at maximum 64 bits), this cut-off is done by rounding it (see below for an example).

To deal correctly with these shortcomings you should never compare to floats for equality, but for closeness. Numpy provides 2 functions for that: np.isclose for comparison of single values (or a item-wise comparison for arrays) and np.allclose for whole arrays. The latter is a np.all(np.isclose(a, b)), so you get a single value for an array.

>>> np.isclose(np.float32('1.000001'), np.float32('0.999999'))
True

But sometimes the rounding is very practicable and matches with our analytical expectation, see for example:

>>> np.float(1) == np.square(np.sqrt(1))
True

After squaring the value will be reduced in size to fit in the given memory, so in this case it's rounded to what we would expect.

These two functions have built-in absolute and relative tolerances (you can also give then as parameter) that are use to compare two values. By default they are rtol=1e-05 and atol=1e-08.

Also, don't mix different packages with their types. If you use Numpy, use Numpy-Types and Numpy-Functions. This will also reduce your rounding errors.

Btw: Rounding errors have even more impact when working with numbers which differ in their exponent widely.

Upvotes: 3

bereal
bereal

Reputation: 34282

I guess, the same considerations as for real numbers are applicable: never assume they can be equal, but rather close enough:

eps = 0.000001
if abs(a - b) < eps:
    print "Equal"

Upvotes: 2

Related Questions