Пользовательская функция, которая принимает входные данные, такие как np.dot

Question

Я хотел бы написать функцию f, которая принимает произвольную функцию g с сигнатурой g : R^n -> R^n -> int, которая "поднимает" g, так что она работает на (R^{nxm}, R^{kxm}), ведя себя как произведение точек. То есть я хочу, чтобы f имел подпись f : R^{nxm} -> R^{mxk} -> R^{nxk}, применяя g ко всем парам строк и столбцов при построении матрицы M, где M_ij = g(A[i,:], B[:,j]).

Возможно ли это?

Например, scipy.spatial.distance.cosine ожидает два вектора. Теперь я бы поднял cosine с f:

from scipy.spatial.distance import cosine

A = np.random.randint(0, 3, (3,4))
B = np.random.randint(0, 3, (5,4))

cosine_lifted = f(cosine)
cosine_lifted(A, B)

. Затем получился бы тот же результат, что и

def sim(A, B):
    ignored_states = np.seterr(divide='raise')
    return 1 - np.divide(np.dot(A, B.T), np.outer(np.linalg.norm(A, axis=1), np.linalg.norm(B, axis=1)))

, который равен sklearn.metrics.pairwise.cosine_similarity плюс * 1025. * part.

Но если бы не было sklearn.metrics.pairwise.cosine_similarity, мне пришлось бы самому реализовать эту поднятую версию cosine (что я, конечно, сделал здесь ...). Но я не хочу делать это для всех функций, которые ведут себя в основном так же, как точечный продукт, в отношении того, как они механически обрабатывают свои аргументы. Поэтому я хотел бы иметь эту функцию f.

Answer 1

Я написал другой ответ, предполагая, что ваши

np.dot(A, B.T)

с входами (3,4) и (5,4) были основной dot функциональностью, которую вы пытались эмулировать. Другими словами, (3,4), (4,5) => (3,5) с суммированием по общему размеру 4. Мой ответ показал, как это 2d-вычисление может быть выполнено с поэлементным умножением.

Для чего стоит, np.dot получает большую часть своей скорости, передавая задачу оптимизированным библиотекам BLAS (или подобным). Они были написаны на C или Фортране и оптимизированы поколениями кодировщиков численного анализа.

Но в описании вашей подписи может быть сказано другое. Это немного сбивает с толку.

g : R^n -> R^n -> int

Значит ли это, что g(x,y) принимает два (n,) массива формы и возвращает целое число? И это не может быть обобщено для работы с 2d массивами?

f : R^{nxm} -> R^{kxm} -> R^{nxm}

Означает ли это, что f(A, B) принимает форму (n, m) и (k, m) форму и возвращает ( н, м) форма? Что случилось с формой k? Это k опечатка?

В качестве альтернативы вы говорите о выполнении (я считаю)

M = np.zeros((N,N))    # (N,M) ok?
for i in range(N):
    for j in range(N):    
        x = A[i,:]; y = B[:,j]
        M[i,j] = g(x, y)

в качестве альтернативы:

M = np.array([[g(x,y) for y in B.T] for x in A])

Предполагая, что g является * Функция 1044 *, которая может работать только с 2-мя 1d-массивами (совпадающей длины) и не может быть обобщена на 2d-массивы, в numpy до compile нет механизма, описанного выше double l oop. g должно быть оценено N ** 2 раза. И если допустить, что g не является тривиальным, эти оценки N * 2 будут доминировать в общем времени оценки, а не в механизме итерации.

np.vectorize обычно принимает функцию, которая принимает скалярные входные данные, но с signature Параметр может работать с вашим g:

 f = np.vectorize(g, signature='(n),(n)')  # signature syntax may be wrong
 M = f(A, B.T)

, но в моем тестировании vectorize всегда был медленнее, чем явная итерация. С подписью еще медленнее. Так что я даже не решаюсь даже упомянуть это.

Answer 2

Вы запрашиваете функцию с такой простой сигнатурой (как умножить две матрицы) или хотите эмулировать всю np.dot поверхность API?

def lift(f):
  def dot(A, B):
    return np.array([[f(v,w) for w in zip(*B)] for v in A])
  return dot

Основное Источником неэффективности в приведенном выше коде являются выделения для всех промежуточных списков. Поскольку мы знаем окончательное возвращаемое значение, этого легко избежать:

def lift(f):
  def dot(A, B):
    result = np.empty((A.shape[0], B.shape[1]))
    for i,v in enumerate(A):
      for j,w in enumerate(zip(*B)):
        result[i,j] = f(v,w)
    return result
  return dot

Циклы довольно дорогие в Python, но, поскольку f работает с k элементами, кажется разумным предположить, что это накладных расходов мало. Вы можете уменьшить его еще больше, скомпилировав с pypy или cython.

Answer 3

matmul был разыгран как ufunc и официально имеет подпись. np.dot является более ранней версией и не имеет подписи.

Но, учитывая 2d массивы, np.dot - это фактически переданная форма умножения с последующим суммированием или «суммой произведений»:

In [587]: A = np.arange(12).reshape(3,4)                                                               
In [588]: B = np.arange(8).reshape(2,4)                                                                

In [589]: np.dot(A, B.T)                                                                               
Out[589]: 
array([[ 14,  38],
       [ 38, 126],
       [ 62, 214]])

эквивалент:

In [591]: (A[:,None,:]*B[None,:,:]).sum(axis=2)                                                        
Out[591]: 
array([[ 14,  38],
       [ 38, 126],
       [ 62, 214]])

Некоторые считают, что подпись в стиле einsum легче следовать:

In [594]: np.einsum('ij,kj->ik', A, B)                                                                 
Out[594]: 
array([[ 14,  38],
       [ 38, 126],
       [ 62, 214]])

, где повторяющиеся j сигналы dot похожи суммирование.

===

Иллюстрирование итерации в моем другом ответе:

In [601]: def g(x,y): 
     ...:     return (x*y).sum() 
     ...:                                                                                              
In [602]: A.shape, B.shape                                                                             
Out[602]: ((3, 4), (2, 4))
In [603]: np.array([[g(x,y) for y in B] for x in A])                                                   
Out[603]: 
array([[ 14,  38],
       [ 38, 126],
       [ 62, 214]])

и версия vectorize:

In [614]: f = np.vectorize(g, signature='(n),(n)->()')                                                 
In [615]: f(A[:,None,:], B[None,:,:])                                                                  
Out[615]: 
array([[ 14,  38],
       [ 38, 126],
       [ 62, 214]])

сравнительные времена:

In [616]: timeit f(A[:,None,:], B[None,:,:])                                                           
255 µs ± 6.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [617]: timeit np.array([[g(x,y) for y in B] for x in A])                                            
69.4 µs ± 116 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [618]: timeit np.dot(A, B.T)                                                                        
3.15 µs ± 128 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

и с использованием @hans' 2nd lift:

In [623]: h = lift(g)                                                                                  
In [624]: h(A,B.T)                                                                                     
Out[624]: 
array([[ 14.,  38.],
       [ 38., 126.],
       [ 62., 214.]])
In [625]: timeit h(A,B.T)                                                                              
102 µs ± 56.5 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)