Replace division by zero numpy

Question

I'm doing a matrix by matrix pointwise division however there are some zeros in the divisor matrix. This results in a warning and in some NaNs. I want these to map to 0, which I can do like this:

edge_map = (xy/(x_norm*y_norm))
edge_map[np.isnan(edge_map)] = 0

However there are two issues with this, first of all it still gives a warning (I don't like warnings) and second of all this requires a second pass over the matrix (not sure if this is unavoidable) and efficiency is very important for this part of the code. Ideas?

Yes so this would mean I do need the second pass, any clue about the warning? — Jan van der Vegt, May 30 '16 at 20:16
http://stackoverflow.com/a/26248892/382936 suggests using the `numpy.errstate(divide='ignore')` context to suppress the warning. Also see http://docs.scipy.org/doc/numpy/reference/generated/numpy.errstate.html#numpy.errstate — Seth Difley, May 30 '16 at 20:20
This is a good question. A `where` function also triggers the warning as the result is calculated. Masks are probably your best option. — Chiel, May 30 '16 at 20:23

Chiel · Accepted Answer · 2016-05-30T20:30:36.070

5

This is probably the fastest solution, but the where function does trigger the error as it precalculates the solutions:

import numpy as np

n = 4

xy = np.random.randint(4, size=(n,n)).astype(float)
x_norm = np.random.randint(4, size=(n,n)).astype(float)
y_norm = np.random.randint(4, size=(n,n)).astype(float)

xy_norm = x_norm*y_norm

edge_map = np.where(xy_norm == 0, xy_norm, xy/xy_norm)

print(xy)
print(xy_norm)
print(edge_map)

edited May 30 '16 at 20:30

answered May 30 '16 at 20:24

Chiel

6,006
2
32
57

3

This combined with "with numpy.errstate(divide='ignore'):" solves the problem and is faster than my own solution, thanks – Jan van der Vegt May 30 '16 at 20:29
I have added type conversions, otherwise the divisions all give zero. – Chiel May 30 '16 at 20:30

score 1 · Answer 2 · answered Jan 13 '20 at 11:06

While Chiel's solution works for this example, it might get expensive as np.where is not lazy. You have to calculate both possible outputs xy_norm and xy/xy_norm for all arguments. If this is not desired, you can just calculate the output, which typically requires a few more lines:

import numpy as np

n = 4

xy = np.random.randint(4, size=(n,n)).astype(float)
x_norm = np.random.randint(4, size=(n,n)).astype(float)
y_norm = np.random.randint(4, size=(n,n)).astype(float)

xy_norm = x_norm*y_norm

edge_map = np.zeros_like(xy_norm)
non_zero = xy_norm != 0
edge_map[non_zero] = xy[non_zero]/xy_norm[non_zero]

print(xy)
print(xy_norm)
print(edge_map)

Here we avoid the division for all elements where xy_norm is zero. As there typically are few zeros and division is not terribly expensive, it's probably not worth from a performance point of view. Alternatively you can change

edge_map[non_zero] = xy[non_zero]/xy_norm[non_zero]

to

edge_map = np.divide(xy, xy_norm, out=edge_map, where=non_zero)

in case you find that more readable (see this).

Replace division by zero numpy

2 Answers2