Scipy разреженный ... массивы?

Question

Итак, я делаю некоторую классификацию Kmeans, используя довольно редкие числовые массивы - множество нулей. Я решил, что воспользуюсь «разреженным» пакетом scipy, чтобы уменьшить накладные расходы на хранение, но меня немного смущает вопрос о том, как создавать массивы, а не матрицы.

Я прошел этот урок о том, как создавать разреженные матрицы: http://www.scipy.org/SciPy_Tutorial#head-c60163f2fd2bab79edd94be43682414f18b90df7

Чтобы имитировать массив, я просто создаю матрицу 1xN, но, как вы можете догадаться, Asp.dot (Bsp) не совсем работает, потому что вы не можете умножить две матрицы 1xN. Я должен был бы транспонировать каждый массив в Nx1, и это довольно неубедительно, так как я делал бы это для каждого вычисления точечного произведения.

Затем я попытался создать матрицу NxN, где столбец 1 == строка 1 (чтобы вы могли умножить две матрицы и просто взять верхний левый угол в качестве точечного произведения), но это оказалось действительно неэффективным .

Я бы хотел использовать редкий пакет scipy в качестве волшебной замены для array () numpy, но пока я не совсем уверен, что делать.

Любой совет?

Answer 1

Используйте формат scipy.sparse, основанный на строках или столбцах: csc_matrix и csr_matrix.

В них используются эффективные реализации C под капотом (включая умножение), а транспонирование не используется (особенно если вы вызываете transpose(copy=False)), как и в случае с массивными массивами.

РЕДАКТИРОВАТЬ: некоторые тайминги через ipython :

import numpy, scipy.sparse
n = 100000
x = (numpy.random.rand(n) * 2).astype(int).astype(float) # 50% sparse vector
x_csr = scipy.sparse.csr_matrix(x)
x_dok = scipy.sparse.dok_matrix(x.reshape(x_csr.shape))

Теперь x_csr и x_dok на 50% разрежены:

print repr(x_csr)
<1x100000 sparse matrix of type '<type 'numpy.float64'>'
        with 49757 stored elements in Compressed Sparse Row format>

И время:

timeit numpy.dot(x, x)
10000 loops, best of 3: 123 us per loop

timeit x_dok * x_dok.T
1 loops, best of 3: 1.73 s per loop

timeit x_csr.multiply(x_csr).sum()
1000 loops, best of 3: 1.64 ms per loop

timeit x_csr * x_csr.T
100 loops, best of 3: 3.62 ms per loop

Похоже, я солгал. Транспозиция очень дешева , но эффективной реализации Csr * csc не существует (в последней версии scipy 0.9.0). Новый объект CSR создается в каждом вызове: - (

Как хак (хотя scipy в наши дни относительно стабилен), вы можете использовать скалярное произведение непосредственно на разреженных данных:

timeit numpy.dot(x_csr.data, x_csr.data)
10000 loops, best of 3: 62.9 us per loop

Обратите внимание, что этот последний подход снова делает умноженное плотное умножение. Разреженность составляет 50%, поэтому она на самом деле быстрее, чем dot(x, x) в 2 раза.

Answer 2

Я не уверен, что это действительно намного лучше или быстрее, но вы могли бы сделать это, чтобы избежать использования транспонирования:

Asp.multiply(Bsp).sum()

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

Однако, вероятно, их проще просто транспонировать:

Asp*Bsp.T

Это не так уж и много, но вы также можете создать подкласс и изменить метод mul ().

Answer 3

Вы можете создать подкласс одного из существующих 2d разреженных массивов

from scipy.sparse import dok_matrix

class sparse1d(dok_matrix):
    def __init__(self, v):
        dok_matrix.__init__(self, (v,))
    def dot(self, other):
        return dok_matrix.dot(self, other.transpose())[0,0]

a=sparse1d((1,2,3))
b=sparse1d((4,5,6))
print a.dot(b)