在我的代碼中，我有一個點位置陣列和一個稀疏關聯矩陣。對于矩陣的每個非零 ij 元素，我想計算 i 點和 j 點之間的距離。
目前我正在使用“nonzero()”提取索引，但是這種方法對輸出索引進行排序，并且這種排序占用了我整個應用程式的大部分執行時間。
import numpy as np
import scipy as sc
import scipy.sparse
#number of points
n = 100000
#point positions
pos = np.random.rand(n,3)
#building incidence matrix, defining the pair wise calculation
density = 1e-3
n_entries = int(n**2*density)
indx = np.random.randint(0,n, n_entries)
indy = np.random.randint(0,n, n_entries)
Pairs = scipy.sparse.csc_matrix((np.ones(n_entries) ,(indx,indy)) , (n,n))
#current method
ii,jj = Pairs.nonzero() #this part is slow
d = np.linalg.norm(pos[ii] - pos[jj], axis = 1)
distance_matrix = scipy.sparse.csr_matrix((d , (ii,jj)) , (n,n))
uj5u.com熱心網友回復：
修改答案
主要問題是indx, indy. 出于這個原因，我們不能簡單地做：
d2 = np.linalg.norm(pos[indx] - pos[indy], axis=1)
distance_matrix2 = scipy.sparse.csr_matrix((d2, (indx, indy)), (n, n))
...因為重復索引處的值將乘以重復項的數量。
下面的原始答案是首先對索引進行重復資料洗掉。但即使速度更快pd_unique2_idx，整個操作最終也比原來花費更多的時間。
相反，由于您已經計算過Pairs（如果您注意到，它也存在重復索引問題：Pairs.max()通常給出超過 1.0），讓我們嘗試獲取非零值的索引而不進行任何排序，僅使用結構中的內容。這是一種方法，遵循indptr和indices 在檔案中的定義：
ii = Pairs.indices
jj = np.repeat(np.arange(len(Pairs.indptr) - 1), np.diff(Pairs.indptr))
您得到的索引與您自己的索引相同ii, jj，但順序不同（列優先而不是行優先）。
定時
只看得到的部分ii, jj：
%timeit Pairs.nonzero()
# 1.27 s ± 5.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit np.repeat(np.arange(len(Pairs.indptr) - 1), np.diff(Pairs.indptr))
# 41.3 ms ± 49.7 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
在背景關系中：
def orig(pos, Pairs, n):
#current method
ii, jj = Pairs.nonzero() #this part is slow
d = np.linalg.norm(pos[ii] - pos[jj], axis=1)
distance_matrix = scipy.sparse.csr_matrix((d, (ii, jj)), (n, n))
return distance_matrix
def modified(pos, Pairs, n):
# only change: these two lines:
ii = Pairs.indices
jj = np.repeat(np.arange(len(Pairs.indptr) - 1), np.diff(Pairs.indptr))
# unchanged
d = np.linalg.norm(pos[ii] - pos[jj], axis=1)
distance_matrix = scipy.sparse.csr_matrix((d, (ii, jj)), (n, n))
return distance_matrix
在您的資料上（np.random.seed(0)以開始，以獲得可重復的結果）：
%timeit orig(pos, Pairs, n)
# 2.16 s ± 15.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit modified(pos, Pairs, n)
# 1.11 s ± 4.17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
正確性
dm = orig(pos, Pairs, n)
dmod = modified(pos, Pairs, n)
np.all(np.c_[dm.nonzero()] == np.c_[dmod.nonzero()])
# True
i, j = dm.nonzero()
np.all(dm[i, j] == dmod[i, j])
# True
原答案
您可以按indx, indy順序計算另一個距離串列，然后制作一個稀疏矩陣：
d2 = np.linalg.norm(pos[indx] - pos[indy], axis=1)
distance_matrix2 = scipy.sparse.csr_matrix((d2, (indx, indy)), (n, n))
但是，與您自己的距離矩陣的比較顯示了一些罕見但可能存在的差異：它們來自您有重復位置的情況，indx, indy對應的Pairs.toarray()[indx, indy]是重復次數。這源于的一個眾所周知的“功能”csc和csr稀疏矩陣重復的元素的總和。例如這里或這里。
For example:
>>> indx = np.array([0,1,0])
>>> indy = np.array([2,0,2])
>>> scipy.sparse.csc_matrix((np.ones(3), (indx, indy))).toarray()
array([[0., 0., 2.],
[1., 0., 0.]])
Unfortunately, the solution proposed by @hpaulj in the posts linked above using todok() do not work anymore with recent versions of scipy.
Using np.unique() is slow, as it sorts the values. pd.unique() is faster (no sort), but accepts only 1D arrays. Here is a way to do a 2D unique that is fast. Note that we need the indices of the unique rows, so that we can index indx and indy, but also the distance matrix if that was already calculated.
def pd_unique2_idx(indx, indy):
df = pd.DataFrame(np.c_[indx, indy])
return np.nonzero(~df.duplicated().values)
# by contrast (slower on long sequences, because it sorts)
def np_unique2_idx(indx, indy):
rc = np.vstack([indx, indy]).T.copy()
dt = rc.dtype.descr * 2
i = np.unique(rc.view(dt), return_index=True)[1]
return i
On 1M elements for both indx and indy, I see:
pd_unique2_idx: 52.1 ms ± 529 μs per loopnp_unique2_idx: 972 ms ± 19.5 ms per loop
How to use this? Following your setup:
i = pd_unique2_idx(indx, indy)
indx = indx[i]
indy = indy[i]
d2 = np.linalg.norm(pos[indx] - pos[indy], axis=1)
distance_matrix2 = scipy.sparse.csr_matrix((d2, (indx, indy)), (n, n))
uj5u.com熱心網友回復：
使用您的樣本陣列：
In [104]: Pairs = sparse.csc_matrix((np.ones(n_entries) ,(indx,indy)) , (n,n))
創建csc本身的時間并非微不足道：
In [105]: timeit Pairs = sparse.csc_matrix((np.ones(n_entries) ,(indx,indy)) , (
...: n,n))
1.08 s ± 44.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
nonzero將矩陣轉換為coo，然后對 0 進行額外檢查：
def nonzero(self):
A = self.tocoo()
nz_mask = A.data != 0
return (A.row[nz_mask], A.col[nz_mask])
該coo轉換是比較快的。我認為不涉及任何排序。
In [106]: timeit Pairs.tocoo()
221 ms ± 17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
額外的掩蔽確實增加了時間。
In [107]: timeit Pairs.nonzero()
1.86 s ± 37.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
由于新制作的csc不會有任何明確的 0，您可以執行以下操作：
M = Pairs.tocoo()
ii,jj = M.row, M.col
coo直接創建矩陣更快：
In [108]: timeit sparse.coo_matrix((np.ones(n_entries) ,(indx,indy)) , (n,n))
158 ms ± 15.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
此套row和col至indx和indy直接（或至多與D型的變化）。它是對csc索引進行排序（按行然后按列）的創建，并進行任何重復的求和。
轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/396610.html