我想對數組的每個對角線應用一個函數（即 np.var，但一般方法會很好）。我的陣列是方形的。我可以這樣做：offset_list = np.arange(-1 * len(arr) + 1, len(arr))diag_var_list = [np.var(np.diagonal(arr, k)) for k in offset_list]如果我想要對角線的一個子集，我可以更改 offset_list。但是使用列表理解似乎效率低下，因為我將在許多大型數組上執行此操作。有沒有更好的辦法？
2 回答
慕的地8271018
TA貢獻1796條經驗 獲得超4個贊
您可以使用以下想法。如果major, minor = A.strides然后將步幅設置A為major + minor, minor（應小心操作以避免超出數組邊界），您將獲得每列為對角線的數組。通過這種方式，A.sum(axis=0)您可以計算對角線之和。對于意味著您可以使用相同但乘以某些值，例如A.shape[0] / [1, 2, ... A.shape[0], ... 2, 1]修復對角線長度的變化。對于方差，您可以使用它variance = <(x - <x>)**2> = <x**2> - <x>**2。
import numpy as np
def rot45(A):
"""
>>> A = np.triu(np.arange(25).reshape(5, 5), 1)
>>> print(A)
[[ 0 1 2 3 4]
[ 0 0 7 8 9]
[ 0 0 0 13 14]
[ 0 0 0 0 19]
[ 0 0 0 0 0]]
>>> print(rot45(A))
[[ 1 2 3 4]
[ 7 8 9 0]
[13 14 0 0]
[19 0 0 0]]
"""
major, minor = A.strides
strides = major + minor, minor
shape = A.shape[0] - 1, A.shape[1]
return np.lib.stride_tricks.as_strided(A, shape, strides)[:, 1:]
def apply_diag(A, func):
"""
>>> A = np.arange(25).reshape(5, 5)
>>> print(A)
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]
[20 21 22 23 24]]
>>> offset_list = np.arange(-1 * len(A) + 1, len(A))
>>> diag_var_list = [np.sum(np.diagonal(A, k)) for k in offset_list]
>>> diag_var_list
[20, 36, 48, 56, 60, 40, 24, 12, 4]
>>> print(apply_diag(A, np.sum))
[20, 36, 48, 56, 60, 40, 24, 12, 4]
"""
U = np.triu(A, 1)
U = rot45(U)
D = np.tril(A, -1).T.copy()
D = rot45(D)
return func(D, axis=0)[::-1].tolist() + [func(np.diag(A))] + func(U, axis=0)[::-1].tolist()[::-1]
def using_numpy(A):
"""
>>> A = np.arange(25).reshape(5, 5)
>>> print(A)
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]
[20 21 22 23 24]]
>>> offset_list = np.arange(-1 * len(A) + 1, len(A))
>>> diag_var_list = [np.var(np.diagonal(A, k)) for k in offset_list]
>>> diag_var_list
[0.0, 9.0, 24.0, 45.0, 72.0, 45.0, 24.0, 9.0, 0.0]
>>> print(using_numpy(A))
[ 0. 9. 24. 45. 72. 45. 24. 9. 0.]
"""
multiply = (A.shape[0] - 1) / np.r_[1:A.shape[0], A.shape[0] - 1, A.shape[0] - 1:0:-1]
return multiply * apply_diag(A ** 2, np.mean) - (multiply * apply_diag(A, np.mean))**2
def list_comp(A, func):
"""
>>> A = np.arange(25).reshape(5, 5)
>>> list_comp(A, np.sum)
[20, 36, 48, 56, 60, 40, 24, 12, 4]
"""
offset_list = np.arange(-1 * len(A) + 1, len(A))
return [func(np.diagonal(A, k)) for k in offset_list]
對于 (100, 100) 大小的矩陣，速度似乎提高了 10 倍，但對于更大的矩陣，速度差異下降得更低。
In [9]: A = np.random.randn(100, 100)
In [10]: %timeit using_numpy(A)
761 μs ± 3.35 μs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [11]: %timeit list_comp(A, np.var)
9.57 ms ± 19.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [12]: A = np.random.randn(1000, 1000)
In [13]: %timeit using_numpy(A)
37.4 ms ± 125 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [14]: %timeit list_comp(A, np.var)
112 ms ± 927 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
呼如林
TA貢獻1798條經驗 獲得超3個贊
線性代數中使用的稀疏矩陣通常具有對角線結構（盡管并非所有對角線）。scipy.sparse 有兩種指定對角線的方法 - 作為數組列表，每個數組的長度不同，以及作為二維填充數組。
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.diags.html
對于numpy（密集）數組，對角線并不是那么特別。所有函數一次只處理一個對角線（可以偏移）。 np.diag*
添加回答
舉報
0/150
提交
取消