Evaluating function using numpy array returns inf and nan

Question

I have the following numpy array and function:

import numpy
import scipy.special

i = numpy.linspace(0, 500, 501, dtype = int)

def func(x, n):
  return scipy.special.eval_hermite(n, x)

I evaluate the function for every element in my numpy array i using two different appraoches:

Approach 1:

hermites = func(0, i)

Approach 2:

hermites = [func(0, idx) for idx in i]

The two approaches result in two different answers. They are different because after element 50, approach 1 starts to return inf and nan. Approach 2 also does not give the correct value for every element of i. However, it is able to calculate more of them. Approach 2 fails for i >= 64.

Both approaches give me an answer in approximately the same time (0.7 s for len(i) = 15000, determined using timeit). What I do not understand are the different results. This is because I learned to avoid for loops in python as much as possible. This does not seems to be the case this time.

The idea that it had to do with memory also crossed my mind. However, evaluating one single element, i.e. print func(0, 64) also returns 0. (equal to the output of approach 2).

What is going on?

Answer 1

This is a bug in scipy created by the occasionally surprising casting rules of numpy's "ufuncs". The problem is that, in scipy version 0.16 and older, when the first argument of eval_hermite is an integer array and the second is an integer scalar, the data type of the return value is single precision floating point (numpy.float32). When the second argument is a 64 bit floating point value, the return type is numpy.float64. The largest value that can be represented with float32 is much less than that of float64, so when the second argument is an integer, eval_hermite overflows to infinity much sooner.

For example, this is with scipy 0.16.0 and numpy 1.10.1:

In [26]: from scipy.special import eval_hermite

Note that the data type of the return value is float32:

In [27]: eval_hermite([20, 50, 100, 200], 0)
Out[27]: 
array([  6.70442586e+11,             -inf,              inf,
                    inf], dtype=float32)

If the second argument is floating point, the return type is float64, and the large values can be represented:

In [28]: eval_hermite([20, 50, 100, 200], 0.0)
Out[28]: 
array([  6.70442573e+011,  -1.96078147e+039,   3.06851876e+093,
         8.45055019e+216])

The work-around for your code is to always ensure that the second argument to eval_hermite is a floating point value. For example,

hermites = func(0.0, i)

This issue has been fixed in the scipy development version (see https://github.com/scipy/scipy/pull/4896), so scipy 0.17 should not have this problem when it is released.