假設我有一個num_indices * n索引矩陣 (in range(m)) 和一個num_indices * n值矩陣，即
m, n = 100, 50
num_indices = 100000
indices = np.random.randint(0, m, (num_indices, n))
values = np.random.uniform(0, 1, (num_indices, n))
我想得到一個 shape 的結果矩陣m * n，這樣每一列都被索引并在索引和值矩陣中的相應列上求和。以下是一些實作：
- 基本的for回圈
tic = time.time()
res1 = np.zeros((m, n))
for i in range(num_indices):
for j in range(n):
res1[indices[i, j], j] = values[i, j]
print(f'element-wise for loop used {time.time() - tic:.5f}s')
- np.add.at，多維
tic = time.time()
res2 = np.zeros((m, n))
np.add.at(res2, (indices, np.arange(n)[None, :]), values)
print(f'np.add.at used {time.time() - tic:.5f}s')
- np.add.at，帶有 for 回圈
tic = time.time()
reslst = []
for colid in range(n):
rescol = np.zeros((m, 1))
np.add.at(rescol, (indices[:, colid], np.zeros(num_indices, dtype=int)), values[:, colid])
reslst.append(rescol)
res3 = np.hstack(reslst)
print(f'np.add.at with loop used {time.time() - tic:.5f}s')
- scipy.sparse，多維
from scipy import sparse
tic = time.time()
res4 = sparse.csr_matrix((values.ravel(), (indices.ravel(), np.tile(np.arange(n), num_indices))), shape=(m, n)).toarray()
print(f'sparse.csr used {time.time() - tic:.5f}s')
- scipy.sparse，帶有 for 回圈
tic = time.time()
reslst = []
for colid in range(n):
rescol = sparse.csr_matrix((values[:, colid], (indices[:, colid], np.zeros(num_indices, dtype=int))), shape=(m, 1)).toarray()
reslst.append(rescol)
res5 = np.hstack(reslst)
print(f'sparse.csr with loop used {time.time() - tic:.5f}s')
結果都是一樣的：
assert np.isclose(res2, res1).all()
assert np.isclose(res3, res1).all()
assert np.isclose(res4, res1).all()
assert np.isclose(res5, res1).all()
但是，速度很奇怪，并不令人滿意：
for loop used 1.69889s
np.add.at used 0.17287s
np.add.at with loop used 0.47767s
sparse.csr used 0.16847s
sparse.csr with loop used 0.09845s
我的基本問題是：
- 為什么
np.add.at這么慢——甚至還慢scipy.sparse？我雖然 numpy 會受益很多，尤其是在多維情況下。 - 對于
scipy.sparse，為什么多維比 for 回圈還要慢？不使用并發嗎？（為什么對于 numpy，多維更快？） - 總的來說，有沒有更快的方法來實作我想要的？
uj5u.com熱心網友回復：
令人驚訝的是，事實證明，np.add.at在 Numpy 中實作的效率非常低。
關于scipy.sparse，我無法在我的機器上使用 Scipy 1.7.1 重現相同的性能改進：它幾乎沒有更快。像 Numpy's 一樣np.add.at，它遠非快速。
您可以使用Numba有效地實作此操作：
import numba as nb
@nb.njit('(int_[:,::1], float64[:,::1], int_, int_)')
def compute(indices, values, n, m):
res = np.zeros((m, n), dtype=np.float64)
for i in range(num_indices):
for j in range(n):
#assert 0 <= indices[i, j] < m
res[indices[i, j], j] = values[i, j]
return res
tic = time.time()
result = compute(indices, values, n, m)
print(f'Numba used {time.time() - tic:.5f}s')
請注意，該函式假定indices包含有效值（即沒有越界）。否則，該函式可能會崩潰或默默地破壞程式。如果您不確定以較慢的計算為代價（在我的機器上慢兩倍），您可以啟用斷言。
請注意，此實作速度很快，只要8 * n * m足夠小以使輸出陣列適合 L1 快取（通常從 16 KB 到 64 KB）。否則，由于低效的隨機訪問模式，它可能會慢一些。如果輸出陣列不適合 L2 快取（通常為幾百 KB），那么這將顯著變慢。
結果
element-wise for loop used 2.45158s
np.add.at used 0.28600s
sparse.csr used 0.19600s
sparse.csr with loop used 0.18900s
Numba used 0.00500s
因此，Numba比最快的實作快約 40 倍，比最慢的實作快約 500 倍。numba 代碼要快得多，因為與 Numpy 和 Scipy 相比，代碼被編譯為優化的本機二進制檔案，沒有額外的開銷（即邊界檢查、臨時陣列、函式呼叫）。
uj5u.com熱心網友回復：
復制它以更好地了解正在發生的事情：
In [48]: m, n = 100, 50
...: num_indices = 100000
...: indices = np.random.randint(0, m, (num_indices, n))
...: values = np.random.uniform(0, 1, (num_indices, n))
add.at案例_
In [49]: res2 = np.zeros((m, n))
...: np.add.at(res2, (indices, np.arange(n)[None, :]), values)
In [50]: res2
Out[50]:
array([[472.25600004, 526.36528387, 493.10320355, ..., 548.43850955,
497.3412013 , 508.51132267],
...])
In [51]: %%timeit
...: res2 = np.zeros((m, n))
...: np.add.at(res2, (indices, np.arange(n)[None, :]), values)
613 ms ± 11.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
稀疏方法
In [52]: from scipy import sparse
In [53]: res4 = sparse.csr_matrix((values.ravel(), (indices.ravel(), np.tile(np.
...: arange(n), num_indices))), shape=(m, n)).toarray()
In [54]: res4
Out[54]:
array([[472.25600004, 526.36528387, 493.10320355, ..., 548.43850955,
497.3412013 , 508.51132267],
[....)
值是相同的，它們沒有什么病態的（即并非全為 0。只是許多浮點數的總和。）
確實，稀疏時間要好一些：
In [55]: timeit res4 = sparse.csr_matrix((values.ravel(), (indices.ravel(), np.t
...: ile(np.arange(n), num_indices))), shape=(m, n)).toarray()
402 ms ± 22.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
通常我發現創建稀疏矩陣比創建密集矩陣需要更長的時間，尤其是當大小僅為 (100,50) 時。從某種意義上說，它確實需要很長時間，半秒。但與add.at.
用 1替換values，我確認平均每個陣列元素有 1000 個重復項（我可以從res2值及其 (0,1) 隨機分布中推斷出）。
In [57]: res5 = sparse.csr_matrix((np.ones_like(indices.ravel()), (indices.ravel
...: (), np.tile(np.arange(n), num_indices))), shape=(m, n)).toarray()
In [58]: res5
Out[58]:
array([[ 971, 1038, 1004, ..., 1056, 988, 1030],
[1016, 992, 949, ..., 994, 982, 1011],
[ 984, 1014, 980, ..., 1001, 1053, 1057],
...,])
如果我正確回憶過去的稀疏代碼探索，則使用這些輸入，csr_matrix首先創建一個coo矩陣，然后按詞法對索引進行排序——即對行進行排序，并在列內進行排序。這將重復的索引放在一起，它可以很容易地求和。鑒于大量重復，這相對有效也就不足為奇了。
隨機性只是在行索引中。sparse_with_loop 利用了這一點，通過迭代n, 50。因此創建每個 (100,1) 稀疏矩陣更容易；再次只是對這些行索引進行排序并求和。
In [61]: %%timeit
...: reslst = []
...: for colid in range(n):
...: rescol = sparse.csr_matrix((values[:, colid], (indices[:, colid], n
...: p.zeros(num_indices, dtype=int))), shape=(m, 1)).toarray()
...: reslst.append(rescol)
...: res5 = np.hstack(reslst)
258 ms ± 452 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
我可能會補充一點np.add.at，雖然比完整的 python 回圈快，但通常比 =它糾正的緩沖慢。
In [72]: %%timeit
...: res10 = np.zeros((m,n))
...: res10[indices, np.arange(n)] = values
103 ms ± 55.1 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
為了正確看待這些時間，生成values陣列的時間是：
In [75]: timeit values = np.random.uniform(0, 1, (num_indices, n))
97.5 ms ± 90.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
大約 1/6 的add.at時間。以及從中創建稀疏矩陣的時間：
In [76]: timeit sparse.csr_matrix(np.random.uniform(0, 1, (num_indices, n)))
379 ms ± 4.51 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
因此，執行索引加法的時間大約為 300 毫秒并不少見。無論如何，處理 (m,n) 形狀的陣列都需要時間。
res5 情況下，計算重復項可以通過以下方式完成bincount：
In [88]: np.array([np.bincount(indices[:,i]) for i in range(n)]).T
Out[88]:
array([[ 971, 1038, 1004, ..., 1056, 988, 1030],
[1016, 992, 949, ..., 994, 982, 1011],
[ 984, 1014, 980, ..., 1001, 1053, 1057],
...,])
In [89]: timeit np.array([np.bincount(indices[:,i]) for i in range(n)]).T
64.7 ms ± 100 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
編輯
有趣的是，如果我做的coo_matrix時間會快很多：
In [95]: timeit res4 = sparse.coo_matrix((values.ravel(), (indices.ravel(), np.t
...: ile(np.arange(n), num_indices))), shape=(m, n)).toarray()
78.3 ms ± 46.3 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
它仍在執行重復項的總和步驟，但作為該to_array()步驟的一部分。但它不必做創建csr矩陣的額外作業。
編輯
我只記得以前的add.atSO 使用bincount過weights：
In [105]: res0 = np.array([np.bincount(indices[:,i],weights=values[:,i]) for i i
...: n range(n)]).T
In [106]: np.allclose(res2,res0)
Out[106]: True
In [107]: timeit res0 = np.array([np.bincount(indices[:,i],weights=values[:,i])
...: for i in range(n)]).T
142 ms ± 133 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/467437.html