Вопрос по sparse-matrix, numpy, scipy, python – Замена для бесшумного вещания с использованием scipy.sparse.csc_matrix

8

В моем коде есть следующее выражение:

a = (b / x[:, np.newaxis]).sum(axis=1)

гдеb ndarray формы(M, N), а такжеx ndarray формы(M,), Сейчас,b на самом деле редко, поэтому для эффективности памяти я хотел бы заменить вscipy.sparse.csc_matrix или жеcsr_matrix, Однако вещание таким способом не реализовано (даже если деление или умножение гарантированно сохраняют разреженность) (записиx ненулевые), и поднимаетNotImplementedError, Есть лиsparse функция яЯ не знаю, что будет делать то, что я хочу? (dot() суммирует по неправильной оси.)

Чтобы было понятно, вы хотите поэлементное деление по оси 1? то есть всеN элементыb[i,:] делятся на?x[i] askewchan
Ага. "Быть понятным Вот почему я включил код. ;) Juan

Ваш Ответ

2   ответа
6

Реализоватьa = (b / x[:, np.newaxis]).sum(axis=1), ты можешь использоватьa = b.sum(axis=1).A1 / x,A1 Атрибут возвращает 1D ndarray, поэтому результатом является 1D ndarray, а неmatrix, Это краткое выражение работает, потому что вы оба масштабируетеx а также суммирование по оси 1. Например:

In [190]: b
Out[190]: 
<3x3 sparse matrix of type '<type 'numpy.float64'="">'
        with 5 stored elements in Compressed Sparse Row format>

In [191]: b.A
Out[191]: 
array([[ 1.,  0.,  2.],
       [ 0.,  3.,  0.],
       [ 4.,  0.,  5.]])

In [192]: x
Out[192]: array([ 2.,  3.,  4.])

In [193]: b.sum(axis=1).A1 / x
Out[193]: array([ 1.5 ,  1.  ,  2.25])
</type>

В целом, если вы хотите масштабировать строки разреженной матрицы с векторомxВы могли бы умножитьb слева с разреженной матрицей, содержащей1.0/x по диагонали. Функцияscipy.sparse.spdiags может быть использован для создания такой матрицы. Например:

In [71]: from scipy.sparse import csc_matrix, spdiags

In [72]: b = csc_matrix([[1,0,2],[0,3,0],[4,0,5]], dtype=np.float64)

In [73]: b.A
Out[73]: 
array([[ 1.,  0.,  2.],
       [ 0.,  3.,  0.],
       [ 4.,  0.,  5.]])

In [74]: x = array([2., 3., 4.])

In [75]: d = spdiags(1.0/x, 0, len(x), len(x))

In [76]: d.A
Out[76]: 
array([[ 0.5       ,  0.        ,  0.        ],
       [ 0.        ,  0.33333333,  0.        ],
       [ 0.        ,  0.        ,  0.25      ]])

In [77]: p = d * b

In [78]: p.A
Out[78]: 
array([[ 0.5 ,  0.  ,  1.  ],
       [ 0.  ,  1.  ,  0.  ],
       [ 1.  ,  0.  ,  1.25]])

In [79]: a = p.sum(axis=1)

In [80]: a
Out[80]: 
matrix([[ 1.5 ],
        [ 1.  ],
        [ 2.25]])
Спасибо, Уоррен! Извини, что выбрал ХаймеБолее быстрый метод ... Я действительно разрывался между скоростью и элегантностью! Оба метода великолепны и точно решают мою проблему. Заметьте также, что я несколько исказил вопрос, и мне также нужно подать заявкуxlogx() вb перед суммированием по оси (0 log (0) равно 0), так что мне все равно придется работать с b.data! Juan
Это работает даже дляM != N до тех пор, пока диагональная матрица дляx имеет форму.(M, M) askewchan
+1 Очень элегантный и чистый способ сделать это. Ницца! Jaime
Не беспокойся, @Juan. Jaime»Ответ отличный. Warren Weckesser
7

Еслиb в формате CSC, тоb.data имеет ненулевые записиb, а такжеb.indices имеет индекс строки каждой из ненулевых записей, так что вы можете сделать свое деление следующим образом:

b.data /= np.take(x, b.indices)

Это'Хакер УорренЭто элегантное решение, но, вероятно, оно также будет быстрее в большинстве случаев:

b = sps.rand(1000, 1000, density=0.01, format='csc')
x = np.random.rand(1000)

def row_divide_col_reduce(b, x):
    data = b.data.copy() / np.take(x, b.indices)
    ret = sps.csc_matrix((data, b.indices.copy(), b.indptr.copy()),
                         shape=b.shape)
    return ret.sum(axis=1)

def row_divide_col_reduce_bis(b, x):
    d = sps.spdiags(1.0/x, 0, len(x), len(x))
    return (d * b).sum(axis=1)

In [2]: %timeit row_divide_col_reduce(b, x)
1000 loops, best of 3: 210 us per loop

In [3]: %timeit row_divide_col_reduce_bis(b, x)
1000 loops, best of 3: 697 us per loop

In [4]: np.allclose(row_divide_col_reduce(b, x),
   ...:             row_divide_col_reduce_bis(b, x))
Out[4]: True

В приведенном выше примере вы можете сократить время почти вдвое, если выполните деление на месте, т.е.

def row_divide_col_reduce(b, x):
    b.data /= np.take(x, b.indices)
    return b.sum(axis=1)

In [2]: %timeit row_divide_col_reduce(b, x)
10000 loops, best of 3: 131 us per loop
Спасибо, Хайме! Я знал, что могу оперироватьb.data но я отсутствовал концептуальноnp.take вызов! Ницца! Juan
Почему вы выбралиnp.take(x, b.indices) вместо ?x[b.indices] askewchan
@askewchan Это часто быстрее, и я пытался заставить его работать как можно быстрее. Jaime

Похожие вопросы