亚洲在线久爱草,狠狠天天香蕉网,天天搞日日干久草,伊人亚洲日本欧美

為了賬號安全,請及時綁定郵箱和手機立即綁定
已解決430363個問題,去搜搜看,總會有你想問的

Scipy 稀疏矩陣循環永遠持續 - 需要提高效率

Scipy 稀疏矩陣循環永遠持續 - 需要提高效率

HUX布斯 2023-02-22 15:54:02
用稀疏矩陣編寫這個循環的最有效的時間和內存方式是什么(目前使用 csc_matrix)for j in range(0, reducedsize):     xs = sum(X[:, j])     X[:, j] = X[:, j] / xs.data[0]例子:縮小尺寸 (int) - 2500X (csc_matrix) - 908x2500循環確實會迭代,但與僅使用 numpy 相比需要很長時間。
查看完整描述

2 回答

?
holdtom

TA貢獻1805條經驗 獲得超10個贊

In [388]: from scipy import sparse                                                      

制作樣本矩陣:


In [390]: M = sparse.random(10,8,.2, 'csc')                                             

我是一個矩陣:


In [393]: M.sum(axis=0)                                                                 

Out[393]: 

matrix([[1.95018736, 0.90924629, 1.93427113, 2.38816133, 1.08713479,

         0.        , 2.45435481, 0.        ]])

那些 0 在劃分時會產生警告 -nan在結果中:


In [394]: M/_                                                                           

/usr/local/lib/python3.6/dist-packages/scipy/sparse/base.py:599: RuntimeWarning: invalid value encountered in true_divide

  return np.true_divide(self.todense(), other)

Out[394]: 

matrix([[0.        , 0.        , 0.        , 0.        , 0.27079623,

                nan, 0.13752665,        nan],

        [0.        , 0.        , 0.        , 0.        , 0.        ,

                nan, 0.32825122,        nan],

        [0.        , 0.        , 0.        , 0.        , 0.        ,

                nan, 0.        ,        nan],

 ...

                nan, 0.        ,        nan]])

0 也會給您的方法帶來問題:


In [395]: for i in range(8): 

     ...:     xs = sum(M[:,i]) 

     ...:     M[:,i] = M[:,i]/xs.data[0] 

     ...:                                                                               

---------------------------------------------------------------------------

IndexError                                Traceback (most recent call last)

<ipython-input-395-0195298ead19> in <module>

      1 for i in range(8):

      2     xs = sum(M[:,i])

----> 3     M[:,i] = M[:,i]/xs.data[0]

      4 


IndexError: index 0 is out of bounds for axis 0 with size 0

但是,如果我們比較沒有 0 和的列,則值匹配:


In [401]: Out[394][:,:5]                                                                

Out[401]: 

matrix([[0.        , 0.        , 0.        , 0.        , 0.27079623],

        [0.        , 0.        , 0.        , 0.        , 0.        ],

        [0.        , 0.        , 0.        , 0.        , 0.        ],

        [0.        , 0.        , 0.        , 0.        , 0.        ],

        [0.49648886, 0.25626608, 0.        , 0.19162678, 0.72920377],

        [0.        , 0.        , 0.30200765, 0.        , 0.        ],

        [0.50351114, 0.        , 0.30445113, 0.41129367, 0.        ],

        [0.        , 0.74373392, 0.        , 0.        , 0.        ],

        [0.        , 0.        , 0.39354122, 0.        , 0.        ],

        [0.        , 0.        , 0.        , 0.39707955, 0.        ]])

In [402]: M.A[:,:5]                                                                     

Out[402]: 

array([[0.        , 0.        , 0.        , 0.        , 0.27079623],

       [0.        , 0.        , 0.        , 0.        , 0.        ],

       [0.        , 0.        , 0.        , 0.        , 0.        ],

       [0.        , 0.        , 0.        , 0.        , 0.        ],

       [0.49648886, 0.25626608, 0.        , 0.19162678, 0.72920377],

       [0.        , 0.        , 0.30200765, 0.        , 0.        ],

       [0.50351114, 0.        , 0.30445113, 0.41129367, 0.        ],

       [0.        , 0.74373392, 0.        , 0.        , 0.        ],

       [0.        , 0.        , 0.39354122, 0.        , 0.        ],

       [0.        , 0.        , 0.        , 0.39707955, 0.        ]])

回到 [394] 我應該首先將矩陣和轉換為稀疏,所以結果也將是稀疏的。稀疏沒有逐元素除法,所以我必須先求密集矩陣的逆。0 仍然很麻煩。


In [409]: M.multiply(sparse.csr_matrix(1/Out[393]))                                     

...

Out[409]: 

<10x8 sparse matrix of type '<class 'numpy.float64'>'

    with 16 stored elements in Compressed Sparse Column format>


查看完整回答
反對 回復 2023-02-22
?
德瑪西亞99

TA貢獻1770條經驗 獲得超3個贊

如果您想在沒有任何內存開銷的情況下執行此操作(就地)

始終考慮數據的實際存儲方式。csc 矩陣的一個小例子。


shape=(5,5)

X=sparse.random(shape[0], shape[1], density=0.5, format='csc')

print(X.todense())


[[0.12146814 0.         0.         0.04075121 0.28749552]

 [0.         0.92208639 0.         0.44279661 0.        ]

 [0.63509196 0.42334964 0.         0.         0.99160443]

 [0.         0.         0.25941113 0.44669367 0.00389409]

 [0.         0.         0.         0.         0.83226886]]


i=0 #first column

print(X.data[X.indptr[i]:X.indptr[i+1]])

[0.12146814 0.63509196]

一個 Numpy 解決方案


所以我們在這里唯一要做的就是逐列修改非零條目。這可以使用部分矢量化的 numpy 解決方案輕松完成。data只是包含所有非零值的數組,indptr存儲每列開始和結束的信息。


def Numpy_csc_norm(data,indptr):

    for i in range(indptr.shape[0]-1):

        xs = np.sum(data[indptr[i]:indptr[i+1]])

        #Modify the view in place

        data[indptr[i]:indptr[i+1]]/=xs

關于性能,這個就地解決方案已經不錯了。如果你想進一步提高性能,你可以使用 Cython/Numba/ 或其他一些可以或多或少容易地用 Python 包裝的編譯代碼。


一個 Numba 解決方案


@nb.njit(fastmath=True,error_model="numpy",parallel=True)

def Numba_csc_norm(data,indptr):

    for i in nb.prange(indptr.shape[0]-1):

        acc=0

        for j in range(indptr[i],indptr[i+1]):

            acc+=data[j]

        for j in range(indptr[i],indptr[i+1]):

            data[j]/=acc

表現


#Create a not to small example matrix

shape=(50_000,10_000)

X=sparse.random(shape[0], shape[1], density=0.001, format='csc')


#Not in-place from hpaulj

def hpaulj(X):

    acc=X.sum(axis=0)

    return X.multiply(sparse.csr_matrix(1./acc))


%timeit X2=hpaulj(X)

#6.54 ms ± 67.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)


#All 2 variants are in-place, 

#but this shouldn't have a influence on the timings


%timeit Numpy_csc_norm(X.data,X.indptr)

#79.2 ms ± 914 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)


#parallel=False -> faster on tiny matrices

%timeit Numba_csc_norm(X.data,X.indptr)

#626 μs ± 30.6 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)


#parallel=True -> faster on larger matrices

%timeit Numba_csc_norm(X.data,X.indptr)

#185 μs ± 5.39 μs per loop (mean ± std. dev. of 7 runs, 10000 loops each)


查看完整回答
反對 回復 2023-02-22
  • 2 回答
  • 0 關注
  • 116 瀏覽
慕課專欄
更多

添加回答

舉報

0/150
提交
取消
微信客服

購課補貼
聯系客服咨詢優惠詳情

幫助反饋 APP下載

慕課網APP
您的移動學習伙伴

公眾號

掃描二維碼
關注慕課網微信公眾號