我正在嘗試采用兩個 matices 的hilbert schmidt 內積。
對于兩個矩陣 A、B，此操作需要 A 與 B 的 Hermitian 共軛的矩陣乘積，然后是跡線（對角線求和）。由于您只需要對角線項，因此找到完整的矩陣乘積毫無意義，因為不需要非對角線項。
實際上，需要計算：
????(??? ??)=∑_{????} (??_????* ??_????)
其中 i 索引行和 j 索引列。
對于稀疏矩陣，最快的方法是什么？我找到了以下類似的文章：
在numpy中計算矩陣乘積的最佳方法是什么？
目前，我正在做：
def hilbert_schmidt_inner_product(mat1, mat2):
## find nonzero ij indices of each matrix
mat1_ij = set([tuple(x) for x in np.array(list(zip(*mat1.nonzero())))])
mat2_ij = set([tuple(x) for x in np.array(list(zip(*mat2.nonzero())))])
## find common ij indices between these that are both nonzero
common_ij = np.array(list(mat1_ij & mat2_ij))
## select common i,j indices from both (now will be 1D array)
mat1_survied = np.array(mat1[common_ij[:,0], common_ij[:,1]])[0]
mat2_survied = np.array(mat2[common_ij[:,0], common_ij[:,1]])[0]
## multiply the nonzero ij common elements (1D dot product!)
trace = np.dot(mat1_survied.conj(),mat2_survied)
return trace
然而，這比：
import numpy as np
sum((mat1.conj().T@mat2).diagonal())
它在進行跟蹤之前會執行完整的 matix 乘積，因此會執行無意義的操作來查找非對角線元素。有沒有更好的方法來做到這一點？
我正在使用以下內容進行基準測驗：
import numpy as np
from scipy.sparse import rand
Dimension = 2**12
A = rand(Dimension, Dimension, density=0.001, format='csr')
B = rand(Dimension, Dimension, density=0.001, format='csr')
運行一些測驗，我發現：
%timeit hilbert_schmidt_inner_product(A,B)
49.2 ms ± 3.13 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit sum((A.conj().T@B).diagonal())
1.48 ms ± 32 μs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit np.einsum('ij,ij->', A.conj().todense(), B.todense())
53.9 ms ± 2.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
uj5u.com熱心網友回復：
另一種選擇是(A.conj().multiply(B)).sum()。
In [111]: Dimension = 2**12
In [112]: A = rand(Dimension, Dimension, density=0.001, format='csr')
...: B = rand(Dimension, Dimension, density=0.001, format='csr')
比較sum((A.conj().T @ B).diagonal())：
In [113]: sum((A.conj().T @ B).diagonal())
Out[113]: 4.152218112255467
In [114]: (A.conj().multiply(B)).sum()
Out[114]: 4.152218112255466
In [115]: %timeit sum((A.conj().T @ B).diagonal())
2.7 ms ± 11.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [116]: %timeit (A.conj().multiply(B)).sum()
1.12 ms ± 4.39 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
當然，對于較大的 值Dimension，相對性能差異要大得多（O(Dimension**3)對于全矩陣乘法與O(Dimension**2)逐元素乘法）：
In [119]: Dimension = 2**14
In [120]: A = rand(Dimension, Dimension, density=0.001, format='csr')
...: B = rand(Dimension, Dimension, density=0.001, format='csr')
In [121]: sum((A.conj().T @ B).diagonal())
Out[121]: 69.23254213582365
In [122]: (A.conj().multiply(B)).sum()
Out[122]: 69.23254213582364
In [123]: %timeit sum((A.conj().T @ B).diagonal())
124 ms ± 1.22 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [124]: %timeit (A.conj().multiply(B)).sum()
8.67 ms ± 63.5 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
轉載請註明出處，本文鏈接：https://www.uj5u.com/yidong/519644.html
上一篇：Numpy轉換資料
下一篇：創建一個回圈以計數0-6100次