我计算了一个非常大的矩阵M，其中包含许多退化的特征向量（不同特征向量具有相同的特征值）。 我使用QR分解来确保这些特征向量是正交的，因此Q是M的正交特征向量，Q^{-1}MQ=D，其中D是对角矩阵。 现在我想检查D是否真的是对角矩阵，但当我打印D时，该矩阵太大无法显示全部内容，那我如何知道它是否真的是对角矩阵？
9个回答
17
去掉对角线并计算非零元素的个数：
np.count_nonzero(x - np.diag(np.diagonal(x)))
- donkopotamus
6
我不确定这与其他技术相比有多快，但：
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。
- Daniel F
5
方法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)
- Divakar
2
2
我们可以做得比Daniel F建议的更好：
最初的回答
最初的回答
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的原地操作。
注：原地操作是指在不创建新变量的情况下直接修改原始变量。
- Make42
3
我在我的系统上看不到加速，但是后来我找到了一种超过两者的方法。请检查我的编辑。 - Daniel F
我没有使用Jupyter - 我如何测试这个（%timeit在PyCharm控制台中无效）。 - Make42
我只是在Spyder中进行测试。 - Daniel F
0
由@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
- Sun Bear
0
快速而简单的方法获取真相。在合理的时间内运作良好。
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
- Michael Alemu
0
我相信这是最简洁的方式：
np.allclose(np.diag(np.diag(a)), a)
- liwt31
3
比起@DanielF建议的方法要慢得多。 - Make42
@make42 是的，既然我们只是复制/粘贴值而不进行计算，我们实际上并不需要
np.allclose 的强大功能。但我可以理解为什么有些人会被浮点数不等式所困扰。 - Daniel F@DanielF：我明白了 :-D，好观点。可以先使用我的方法（看看我能获得的速度增益:-)），如果答案为False，则只需将np.all()替换为np.allclose()。 - Make42
0
我担心这是否是最有效的方法，但想法是掩盖对角线元素并检查所有其他元素是否为零。我猜这足以检查将矩阵标记为对角矩阵。
因此，我们创建一个与输入矩阵大小相同的虚拟数组，并用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
- ZdaR
-3
import numpy as np
is_diagonal = (np.trace(mat) == np.sum(mat))
- user3166559
2
3虽然这段代码可能解决了问题，但包括解释它是如何解决问题的，以及为什么这样做可以提高你的帖子质量，并且可能会得到更多的赞。请记住，你正在回答未来的读者，而不仅仅是现在提问的人。请编辑你的答案以添加解释，并指出适用的限制和假设。 - fcdt
1一个简单的反例：
mat = np.array([[1,-1,0],[0,1,1],[0,0,1]]) - anon01网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
a.shape[0]>a.shape[1]，它将跑出缓冲区的末尾。我想赏金只适用于方形输入，但是有些情况下非方形对角张量很有用。我想我的安全阀“assert”是我时间滞后的一部分。 - Daniel Fassert的时间记录。 - Divakar