Как я могу вызвать список функций numpy без цикла for?

Question

Я делаю анализ данных, который включает в себя минимизацию ошибки наименьших квадратов между набором точек и соответствующим набором ортогональных функций.Другими словами, я беру набор значений y и набор функций, и пытаюсь сосредоточиться на значении x, которое получает все функции, наиболее близкие к их соответствующему значению y.Все делается в классе data_set.Все функции, с которыми я сравниваю, хранятся в одном списке, и я использую метод класса для вычисления общей ошибки lsq для всех них:

self.fits = [np.poly1d(np.polyfit(self.x_data, self.y_data[n],10)) for n in range(self.num_points)]

def error(self, x, y_set):
    arr = [(y_set[n] - self.fits[n](x))**2 for n in range(self.num_points)]
    return np.sum(arr)

Это было нормально, когда я имелзначительно больше времени, чем данных, но теперь я беру тысячи значений x, каждое из которых имеет тысячи значений y, и этот цикл for является недопустимо медленным.Я пытался использовать np.vectorize:

#global scope
def func(f,x):
    return f(x)
vfunc = np.vectorize(func, excluded=['x'])
…
…
#within data_set class
    def error(self, x, y_set):
        arr = (y_set - vfunc(self.fits, x))**2
        return np.sum(arr)

func(self.fits[n], x) работает нормально, пока действует n, и, насколько я могу судить из документов , vfunc(self.fits, x) должно быть эквивалентно

[self.fits[n](x) for n in range(self.num_points)]

, но вместо этого выдается:

ValueError: cannot copy sequence with size 10 to array axis with dimension 11

10 - степень соответствия полинома, а 11 - (по определению) числоусловия в нем, но я понятия не имею, почему они появляются здесь.Если я изменяю порядок подгонки, сообщение об ошибке отражает изменение.Кажется, что np.vectorize принимает каждый элемент self.fits в виде списка, а не np.poly1d функции.

В любом случае, если кто-то может помочь мне лучше понять np.vectorize или предложить другой способустранить эту петлю, это было бы здорово.

Answer 1

Поскольку все рассматриваемые функции имеют очень похожую структуру, мы можем «векторизовать» их вручную, как только мы извлечем поли-коэффициенты.Фактически, тогда функция представляет собой довольно простой однострочный текст, eval_many ниже:

import numpy as np

def poly_vec(list_of_polys):
    O = max(p.order for p in list_of_polys)+1
    C = np.zeros((len(list_of_polys), O))
    for p, c in zip(list_of_polys, C):
        c[len(c)-p.order-1:] = p.coeffs
    return C

def eval_many(x,C):
    return C@np.vander(x,11).T

# make example
list_of_polys = [np.poly1d(v) for v in np.random.random((1000,11))]
x = np.random.random((2000,))

# put all coeffs in one master matrix
C = poly_vec(list_of_polys)

# test
assert np.allclose(eval_many(x,C), [p(x) for p in list_of_polys])

from timeit import timeit

print('vectorized', timeit(lambda: eval_many(x,C), number=100)*10)
print('loopy     ', timeit(lambda: [p(x) for p in list_of_polys], number=10)*100)

Пример выполнения:

vectorized 6.817315469961613
loopy      56.35076989419758