我有一個很大的N = 100003d 向量矩陣。為簡化起見,我將在此處使用 10 x 3 矩陣作為示例:
import numpy as np
A = np.array([[1.2, 2.3, 0.8],
[3.2, 2.1, 0.5],
[0.8, 4.4, 4.4],
[-0.2, -1.1, -1.1],
[2.4, 4.6, 1.6],
[0.5, 0.96, 0.33],
[1.1, 2.2, 3.3],
[-2.2, -4.41, -6.62],
[3.4, 5.5, 3.8],
[-5.1, -28., -28.1]])
我想找到所有幾乎彼此平行的唯一向量對。需要使用容差測量,我想獲得所有唯一的行索引對(無論順序如何)。我設法撰寫了以下代碼:
def all_parallel_pairs(A, tol=0.1):
res = set()
for i, v1 in enumerate(A):
for j, v2 in enumerate(A):
if i == j:
continue
norm = np.linalg.norm(np.cross(v1, v2))
if np.isclose(norm, 0., rtol=0, atol=tol):
res.add(tuple(sorted([i, j])))
return np.array(list(res))
print(all_parallel_pairs(A, tol=0.1))
out[1]: [[0 4]
[2 3]
[6 7]
[4 5]
[0 5]]
但是,因為我使用了兩個 for 回圈,所以當N它很大時會變慢。我覺得應該有更有效和更 Numpyic 的方法來做到這一點。有什么建議?
uj5u.com熱心網友回復:
請注意,該函式np.cross從檔案中接收一個向量陣列:
回傳兩個(陣列)向量的叉積。
因此,一種方法是使用 numpy 高級索引來找到需要計算叉積的正確向量:
# generate the i, j indices (note that only the upper triangular matrices of indices is needed)
rows, cols = np.triu_indices(A.shape[0], 1)
# find the cross products using numpy indexing on A, and the np.cross can take array of vectors
cross = np.cross(A[rows], A[cols])
# find the values that are close to 0
arg = np.argwhere(np.isclose(0, (cross * cross).sum(axis=1) ** 0.5, rtol=0, atol=0.1))
# get the i, j indices where is 0
res = np.hstack([rows[arg], cols[arg]])
print(res)
輸出
[[0 4]
[0 5]
[2 3]
[4 5]
[6 7]]
表達方式:
(cross * cross).sum(axis=1) ** 0.5
是應用于np.linalg.norm向量陣列的更快替換。
uj5u.com熱心網友回復:
作為對Dani Masejo answer的改進更新,您可以使用 GPU_aided 或 TPU_aided 庫,例如JAX:
from jax import jit
@jit
def test_jit():
rows, cols = np.triu_indices(A.shape[0], 1)
cross = np.cross(A[rows], A[cols])
arg = np.argwhere(np.isclose(0, (cross * cross).sum(axis=1) ** 0.5, rtol=0, atol=0.1))
res = np.hstack([rows[arg], cols[arg]])
return res
print(test_jit())
使用 google colab TPU 運行時的結果如下:
100 loops, best of 5: 12.2 ms per loop # the question code
100 loops, best of 5: 152 μs per loop # Dani Masejo code
100 loops, best of 5: 81.5 μs per loop # using jax library
當資料量增加時,差異將顯著。
轉載請註明出處,本文鏈接:https://www.uj5u.com/yidong/366735.html
上一篇:來自檔案的數學.txt
