我计算了一个非常大的矩阵M,其中包含许多退化的特征向量(不同特征向量具有相同的特征值)。 我使用QR分解来确保这些特征向量是正交的,因此Q是M的正交特征向量,Q^{-1}MQ=D,其中D是对角矩阵。 现在我想检查D是否真的是对角矩阵,但当我打印D时,该矩阵太大无法显示全部内容,那我如何知道它是否真的是对角矩阵?
去掉对角线并计算非零元素的个数:
np.count_nonzero(x - np.diag(np.diagonal(x)))
我不确定这与其他技术相比有多快,但:
def isDiag(M):
i, j = np.nonzero(M)
return np.all(i == j)
编辑 让我们计时:
M = np.random.randint(0, 10, 1000) * np.eye(1000)
def a(M): #donkopotamus solution
return np.count_nonzero(M - np.diag(np.diagonal(M)))
%timeit a(M)
100 loops, best of 3: 11.5 ms per loop
%timeit is_diagonal(M)
100 loops, best of 3: 10.4 ms per loop
%timeit isDiag(M)
100 loops, best of 3: 12.5 ms per loop
嗯,这个比较慢,可能是由于构建 i
和 j
导致的。
让我们试着通过删除减法步骤来改进 @donkopotamus 的解决方案:
def b(M):
return np.all(M == np.diag(np.diagonal(M)))
%timeit b(M)
100 loops, best of 3: 4.48 ms per loop
这稍微好一些。
编辑2 我想出了一种更快的方法:
def isDiag2(M):
i, j = M.shape
assert i == j
test = M.reshape(-1)[:-1].reshape(i-1, j+1)
return ~np.any(test[:, 1:])
这并没有进行任何计算,只是进行了重塑。事实证明,将一个对角矩阵重塑为+1行会将所有数据放在第一列中。然后,您可以检查任何非零的连续块,这对于numpy
来说速度要快得多。让我们看一下时间:
def Make42(m):
b = np.zeros(m.shape)
np.fill_diagonal(b, m.diagonal())
return np.all(m == b)
%timeit b(M)
%timeit Make42(M)
%timeit isDiag2(M)
100 loops, best of 3: 4.88 ms per loop
100 loops, best of 3: 5.73 ms per loop
1000 loops, best of 3: 1.84 ms per loop
似乎在小数据集上,我的原始代码比 @Make42 的更快。
M = np.diag(np.random.randint(0,10,10000))
%timeit b(M)
%timeit Make42(M)
%timeit isDiag2(M)
The slowest run took 35.58 times longer than the fastest. This could mean that an intermediate result is being cached.
1 loop, best of 3: 335 ms per loop
<MemoryError trace removed>
10 loops, best of 3: 76.5 ms per loop
在更大的数据集上,@Make42会出现内存错误。但我似乎没有和他们一样多的RAM。
方法1:使用NumPy的strides/np.lib.stride_tricks.as_strided
我们可以利用NumPy strides
来将非对角线元素作为视图呈现。因此,这里没有内存开销,而且几乎是免费运行!这个想法在this post
中已经被探索过。
因此,我们有 -
# https://stackoverflow.com/a/43761941/ @Divakar
def nodiag_view(a):
m = a.shape[0]
p,q = a.strides
return np.lib.stride_tricks.as_strided(a[:,1:], (m-1,m), (p+q,q))
展示使用方法的示例运行 -
In [175]: a
Out[175]:
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
In [176]: nodiag_view(a)
Out[176]:
array([[ 1, 2, 3, 4],
[ 6, 7, 8, 9],
[11, 12, 13, 14]])
我们通过在一个大数组上使用它来验证免费运行时和无内存开销的声明 -
In [182]: a = np.zeros((10000,10000), dtype=int)
...: np.fill_diagonal(a,np.arange(len(a)))
In [183]: %timeit nodiag_view(a)
6.42 µs ± 48.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [184]: np.shares_memory(a, nodiag_view(a))
Out[184]: True
那么,我们在这里怎么使用它呢?只需检查所有的nodiag_view
元素是否为0,表示对角矩阵!
因此,要解决我们这里的问题,对于输入数组a
,可以这样写:
isdiag = (nodiag_view(a)==0).all()
方法 #2:诡异的方式
为了完整性,一种诡异的方式是临时保存对角线元素,将其赋值为0s
,检查所有元素是否都为0s
。如果是,表示这是一个对角矩阵。最后将对角元素重新赋值回去。
实现方法如下:
def hacky_way(a):
diag_elem = np.diag(a).copy()
np.fill_diagonal(a,0)
out = (a==0).all()
np.fill_diagonal(a,diag_elem)
return out
让我们在一个大数组上进行时间测试,并比较它们的性能-
In [3]: a = np.zeros((10000,10000), dtype=int)
...: np.fill_diagonal(a,np.arange(len(a)))
In [4]: %timeit (nodiag_view(a)==0).all()
52.3 ms ± 393 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [5]: %timeit hacky_way(a)
51.8 ms ± 250 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
其他方法来自于@Daniel F的帖子,涉及到其他方法 -
# @donkopotamus solution improved by @Daniel F
def b(M):
return np.all(M == np.diag(np.diagonal(M)))
# @Daniel F's soln without assert check
def isDiag2(M):
i, j = M.shape
test = M.reshape(-1)[:-1].reshape(i-1, j+1)
return ~np.any(test[:, 1:])
# @Make42's soln
def Make42(m):
b = np.zeros(m.shape)
np.fill_diagonal(b, m.diagonal())
return np.all(m == b)
使用与之前相同的设置进行计时 -
In [6]: %timeit b(a)
...: %timeit Make42(a)
...: %timeit isDiag2(a)
218 ms ± 1.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
302 ms ± 1.25 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
67.1 ms ± 1.35 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
最初的回答
import numpy as np
import time
a = np.diag(np.random.random(19999))
t1 = time.time()
np.all(a == np.diag(np.diagonal(a)))
print(time.time()-t1)
t1 = time.time()
b = np.zeros(a.shape)
np.fill_diagonal(b, a.diagonal())
np.all(a == b)
print(time.time()-t1)
导致
2.5737204551696777
0.6501829624176025
其中一个技巧是使用a.diagonal()
而不是np.diagonal(a)
,因为后者实际上调用了前者。但最厉害的是快速构建b
并结合对b
的原地操作。
注:原地操作是指在不创建新变量的情况下直接修改原始变量。
由@Divakar提供的解决方案nodiag_view
和hacky_way
以及@Daniel F提供的解决方案isDiag2
都非常高效。其他解决方案则比较慢。
被接受的解决方案是@donkopotamus提供的,虽然易于实现,但是速度仅次于最快答案7.55倍之多。
from timeit import timeit
setup = '''
import numpy as np
d = np.zeros((10000,10000), dtype=int)
np.fill_diagonal(d,np.arange(len(d)))
def nodiag_view(a):
m = a.shape[0]
p,q = a.strides
return (np.lib.stride_tricks.as_strided(a[:,1:], (m-1,m), (p+q,q)) == 0).all()
def hacky_way(a):
diag_elem = np.diag(a).copy()
np.fill_diagonal(a,0)
out = (a==0).all()
np.fill_diagonal(a,diag_elem)
return out
def isDiag(M):
i, j = np.nonzero(M)
return np.all(i == j)
def isDiag2(M):
i, j = M.shape
assert i == j
test = M.reshape(-1)[:-1].reshape(i-1, j+1)
return ~np.any(test[:, 1:])
def Make42(m):
b = np.zeros(m.shape)
np.fill_diagonal(b, m.diagonal())
return np.all(m == b)
def by_donkopotamus(a):
return np.count_nonzero(a - np.diag(np.diag(a))) == 0
def by_liwt31(a):
np.allclose(np.diag(np.diag(a)), a)
'''
test0 = '''nodiag_view(d)'''
test1 = '''hacky_way(d)'''
test2 = '''isDiag(d)'''
test3 = '''isDiag2(d)'''
test4 = '''Make42(d)'''
test5 = '''by_donkopotamus(d)'''
test6 = '''by_liwt31(d)'''
print('test0:', timeit(test0, setup, number=100))
print('test1:', timeit(test1, setup, number=100))
print('test2:', timeit(test2, setup, number=100))
print('test3:', timeit(test3, setup, number=100))
print('test4:', timeit(test4, setup, number=100))
print('test5:', timeit(test5, setup, number=100))
print('test6:', timeit(test6, setup, number=100))
结果:
test0: 4.194842008990236
test1: 4.11843847198179
test2: 28.11888137299684
test3: 5.095675196003867
test4: 22.56097131301067
test5: 31.05823188900831
test6: 106.19386338599725
for i in range(0, len(matrix[0])):
for j in range(0, len(matrix[0])):
if ((i != j) and
(matrix[i][j] != 0)) :
return False
return True
我相信这是最简洁的方式:
np.allclose(np.diag(np.diag(a)), a)
np.allclose
的强大功能。但我可以理解为什么有些人会被浮点数不等式所困扰。 - Daniel F我担心这是否是最有效的方法,但想法是掩盖对角线元素并检查所有其他元素是否为零。我猜这足以检查将矩阵标记为对角矩阵。
因此,我们创建一个与输入矩阵大小相同的虚拟数组,并用1初始化。然后将对角线元素替换为零。现在我们执行输入矩阵和虚拟矩阵的逐元素乘法。因此,在这里,我们将输入矩阵的对角线元素替换为零,并将其他元素保持不变。
现在最后,我们检查是否有任何非零元素。
def is_diagonal(matrix):
#create a dummy matrix
dummy_matrix = np.ones(matrix.shape, dtype=np.uint8)
# Fill the diagonal of dummy matrix with 0.
np.fill_diagonal(dummy_matrix, 0)
return np.count_nonzero(np.multiply(dummy_matrix, matrix)) == 0
diagonal_matrix = np.array([[3, 0, 0],
[0, 7, 0],
[0, 0, 4]])
print is_diagonal(diagonal_matrix)
>>> True
random_matrix = np.array([[3, 8, 0],
[1, 7, 8],
[5, 0, 4]])
print is_diagonal(random_matrix)
>>> False
import numpy as np
is_diagonal = (np.trace(mat) == np.sum(mat))
mat = np.array([[1,-1,0],[0,1,1],[0,0,1]])
- anon01
a.shape[0]>a.shape[1]
,它将跑出缓冲区的末尾。我想赏金只适用于方形输入,但是有些情况下非方形对角张量很有用。我想我的安全阀“assert”是我时间滞后的一部分。 - Daniel Fassert
的时间记录。 - Divakar