2D разреженная матрица, умноженная на 3D массив

Question

Можно ли как-нибудь умножить разреженную 2D-матрицу на 3D-массив?Например, у меня есть эта функция

def myFun(x, p):
    r = 2
    out = x * np.log(p) + r * np.log(1-p)
    return out

, где x - массив измерений 3500, 90, а p - другой массив с измерениями 3500, 90, 70.На данный момент и x, и p являются плотными массивами, и я просто транслирую, когда вызываю функцию:

out = myFun(x[..., None], p)

Однако массив x довольно редок, только 7% его элементовненулевая.С другой стороны, p не имеет нулей, а только плавает между нулем и единицей.Хотя я надеюсь, что с разреженной матрицей (от scipy.sparse, вероятно) я увижу улучшение скорости.Тем не менее, я не знаю, как сделать эту операцию или, если это более эффективно, пожалуйста.

Я использую Python 3.

Большое спасибо

Answer 1

Вы можете попробовать следующую реализацию.Для этой простой функции это выглядит немного преувеличенным, но у меня также были проблемы с получением numexpr для работы с Intel SVML (в противном случае я бы предпочел цифру pr).Это решение должно давать 0,07 с на вызов на Quadcore i7 и должно довольно хорошо масштабироваться на большем количестве ядер.Также обратите внимание, что при первом вызове затраты на компиляцию составляют около 0,5 с.

Установка Intel SVML

import numpy as np
import numba as nb

x = np.random.uniform(-13, 1, (3500, 90, 1)).clip(0, None)
p = np.random.random((3500, 90, 70))

@nb.njit(parallel=True,fastmath=True)
def nb_myFun_sp(x, p):
    out=np.empty(p.shape,p.dtype)
    r = 2.

    for i in nb.prange(p.shape[0]):
      for j in range(p.shape[1]):
        if x[i,j,0]!=0.:
          x_=x[i,j,0]
          for k in range(p.shape[2]):
            out[i,j,k] = x_ * np.log(p[i,j,k]) + r * np.log(1.-p[i,j,k])
        else:
          for k in range(p.shape[2]):
            out[i,j,k] = r * np.log(1.-p[i,j,k])

    return out

@nb.njit(parallel=True,fastmath=True)
def nb_myFun(x, p):
    out=np.empty(p.shape,p.dtype)
    r = 2.

    for i in nb.prange(p.shape[0]):
      for j in range(p.shape[1]):
        x_=x[i,j,0]
        for k in range(p.shape[2]):
          out[i,j,k] = x_ * np.log(p[i,j,k]) + r * np.log(1.-p[i,j,k])
    return out

Answer 2

Вы можете использовать разреженность x, используя ключевое слово where.

def sprse(x, p):
    r = 2
    out = x * np.log(p, where=x.astype(bool)) + r * np.log(1-p)
    return out

from timeit import timeit
x = np.random.uniform(-13, 1, (3500, 90, 1)).clip(0, None)
p = np.random.random((3500, 90, 70))
assert np.all(sprse(x, p)==myFun(x, p))
def f():
    return myFun(x, p)
print(timeit(f, number=3))
def f():
    return sprse(x, p)
print(timeit(f, number=3))

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

5.171174691990018
3.2122434769989923