How to speed up my numpy loop using numpy.where ()
I recently wrote a function about the ordered logit model.
But I need a lot of time when dealing with big data.
So I want to rewrite the code and replace the numpy.where function with an if .
I have some problems with my new code, I don't know how to do it.
If you know please help me. Many thanks!
This is my original function.
import numpy as np
from scipy.stats import logistic
def func(y, X, thresholds):
ll = 0.0
for row in zip(y, X):
if row[0] == 0:
ll += logistic.logcdf(thresholds[0] - row[1])
elif row[0] == len(thresholds):
ll += logistic.logcdf(row[1] - thresholds[-1])
else:
for i in xrange(1, len(thresholds)):
if row[0] == i:
diff_prob = logistic.cdf(thresholds[i] - row[1]) - logistic.cdf(thresholds[i - 1] - row[1])
if diff_prob <= 10 ** -5:
ll += np.log(10 ** -5)
else:
ll += np.log(diff_prob)
return ll
y = np.array([0, 1, 2])
X = [2, 2, 2]
thresholds = np.array([2, 3])
print func(y, X, thresholds)
This is new, but not perfect code.
y = np.array([0, 1, 2]) X = [2, 2, 2] thresholds = np.array([2, 3]) ll = np.where(y == 0, logistic.logcdf(thresholds[0] - X), np.where(y == len(thresholds), logistic.logcdf(X - thresholds[-1]), np.log(logistic.cdf(thresholds[1] - X) - logistic.cdf(thresholds[0] - X)))) print ll.sum()
The problem is that I don't know how to rewrite the subcycle ( for I am in xrange (1, len (thresholds)):) .
source to share
I think asking how to implement it using np.where
is a bit of an X / Y problem .
So, I will try to explain how to approach optimization of this function.
My first instinct is to get rid of the cycle for
that was causing the pain anyway:
import numpy as np
from scipy.stats import logistic
def func1(y, X, thresholds):
ll = 0.0
for row in zip(y, X):
if row[0] == 0:
ll += logistic.logcdf(thresholds[0] - row[1])
elif row[0] == len(thresholds):
ll += logistic.logcdf(row[1] - thresholds[-1])
else:
diff_prob = logistic.cdf(thresholds[row[0]] - row[1]) - \
logistic.cdf(thresholds[row[0] - 1] - row[1])
diff_prob = 10 ** -5 if diff_prob < 10 ** -5 else diff_prob
ll += np.log(diff_prob)
return ll
y = np.array([0, 1, 2])
X = [2, 2, 2]
thresholds = np.array([2, 3])
print(func1(y, X, thresholds))
I just replaced i
with row[0]
without changing the semantics of the loop. So one for a loop less.
Now I would like the form of statements in different branches to if-else
be the same. To this end:
import numpy as np
from scipy.stats import logistic
def func2(y, X, thresholds):
ll = 0.0
for row in zip(y, X):
if row[0] == 0:
ll += logistic.logcdf(thresholds[0] - row[1])
elif row[0] == len(thresholds):
ll += logistic.logcdf(row[1] - thresholds[-1])
else:
ll += np.log(
np.maximum(
10 ** -5,
logistic.cdf(thresholds[row[0]] - row[1]) -
logistic.cdf(thresholds[row[0] - 1] - row[1])
)
)
return ll
y = np.array([0, 1, 2])
X = [2, 2, 2]
thresholds = np.array([2, 3])
print(func2(y, X, thresholds))
Now the expression in each branch looks like ll += expr
.
There are several different paths in this piont that might need optimization. You can try to optimize the loop by writing it as an understanding, but I suspect it won't give you much speed boost.
An alternative way is to pull the conditions if
out of the loop. This was your intention with np.where
:
import numpy as np
from scipy.stats import logistic
def func3(y, X, thresholds):
y_0 = y == 0
y_end = y == len(thresholds)
y_rest = ~(y_0 | y_end)
ll_1 = logistic.logcdf(thresholds[0] - X[ y_0 ])
ll_2 = logistic.logcdf(X[ y_end ] - thresholds[-1])
ll_3 = np.log(
np.maximum(
10 ** -5,
logistic.cdf(thresholds[y[ y_rest ]] - X[ y_rest ]) -
logistic.cdf(thresholds[ y[y_rest] - 1 ] - X[ y_rest])
)
)
return np.sum(ll_1) + np.sum(ll_2) + np.sum(ll_3)
y = np.array([0, 1, 2])
X = np.array([2, 2, 2])
thresholds = np.array([2, 3])
print(func3(y, X, thresholds))
Note that I turned X
in np.array
to be able to use immanent indexing.
At this point, I would say that it is fast enough for my purposes. However, you can stop earlier or above this point, depending on your requirements.
On my computer, I get the following results:
y = np.random.random_integers(0, 10, size=(10000,))
X = np.random.random_integers(0, 10, size=(10000,))
thresholds = np.cumsum(np.random.rand(10))
%timeit func(y, X, thresholds) # Original
1 loops, best of 3: 1.51 s per loop
%timeit func1(y, X, thresholds) # Removed for-loop
1 loops, best of 3: 1.46 s per loop
%timeit func2(y, X, thresholds) # Standardized if statements
1 loops, best of 3: 1.5 s per loop
%timeit func3(y, X, thresholds) # Vectorized ~ 500x improvement
100 loops, best of 3: 2.74 ms per loop
source to share