A sum of numpy elements with a given index

Question

This code implements a single variable regression without an intercept. It works, but I can't figure out a way to do it without resorting to using slow python iteration. Any ideas?

# y: numpy array of n values, for large n
# coeff: numpy array of n values, for large n
# L : size of result
# l_indices : numpy array of indices from 0 to L-1

def simple_regression(y, coeff, L, l_index):
    numerator = y*coeff
    denominator = np.square(coeff)
    numsum = np.zeros(L)
    denomsum = np.zeros(L)
    for (n,d,l) in zip(numerator,denominator,l_index):
        numsum[l] += n
        denomsum[l] += d
    return numsum / denomsum

fundamentally an operation like the following, that doesn't do a bunch of memory allocation:

 numsum[l] = np.sum(numerator[l_index == l])

(Doing it that way is much lower then my first code)

Paul Panzer · Accepted Answer · 2018-06-11 16:42:41Z

You can use numpy.bincount:

import numpy as np

def simple_regression(y, coeff, L, l_index):
    numerator = y*coeff
    denominator = np.square(coeff)
    numsum = np.zeros(L)
    denomsum = np.zeros(L)
    for (n,d,l) in zip(numerator,denominator,l_index):
        numsum[l] += n
        denomsum[l] += d
    return numsum / denomsum

def simple_regression_pp(y, coeff, L, l_index):
    numerator = y*coeff
    denominator = np.square(coeff)
    numsum = np.bincount(l_index, numerator, L)
    denomsum = np.bincount(l_index, denominator, L)
    return numsum / denomsum

def simple_regression_br(y, coeff, L, l_index):
    numerator = y*coeff
    denominator = np.square(coeff)
    numsum = np.zeros(L)
    denomsum = np.zeros(L)
    np.add.at(numsum, l_index, numerator)
    np.add.at(denomsum, l_index, denominator)
    return numsum / denomsum

L, N = 1_000, 1_000_000
y, coeff = np.random.random((2, N))
l_index = np.random.randint(0, L, (N,))

from timeit import timeit

print('OP', timeit("simple_regression(y, coeff, L, l_index)", globals=globals(),
                   number=10), 'sec')
print('pp', timeit("simple_regression_pp(y, coeff, L, l_index)",
                   globals=globals(), number=10), 'sec')
print('br', timeit("simple_regression_br(y, coeff, L, l_index)",
                   globals=globals(), number=10), 'sec')

Sample run:

OP 6.602819449035451 sec
pp 0.12009818502701819 sec
br 1.5504542298149318 sec

Bi Rico · Accepted Answer · 2018-06-11 16:31:57Z

1

if you know your index l_index only has unique values, you can do:

numsum[l_index] += numerator
denomsum[l_index] += denominator

if your index is not known to be unique you can do the same thing using numpy.add.at:

numpy.add.at(numsum, l_index, numerator)
numpy.add.at(denomsum, l_index, denominator)

answered Jun 11, 2018 at 16:31

Bi Rico

25.9k3 gold badges57 silver badges75 bronze badges

Collectives™ on Stack Overflow

A sum of numpy elements with a given index

2 Answers 2

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related