检查一个大矩阵是否是python中的对角矩阵

Question

我计算了一个非常大的矩阵 M，其中包含许多退化特征向量（具有相同特征值的不同特征向量）。我使用 QR 分解来确保这些特征向量是正交的，所以 Q 是 M 的正交特征向量，并且 Q^{-1}MQ = D，其中 D 是对角矩阵。现在我想检查D是否真的是对角矩阵，但是当我打印D时，矩阵太大而无法显示所有它们，那么我怎么知道它是否真的是对角矩阵呢？

Answer 1

去除对角线并计算非零元素：

np.count_nonzero(x - np.diag(np.diagonal(x)))

Answer 2

不确定这与其他相比有多快，但是：

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

这样好一点。

EDIT2我想出了一个更快的方法：

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 在较大的集合上会出现内存错误。但是我的内存似乎没有他们那么多。

Answer 3

方法一：使用NumPy strides/np.lib.stride_tricks.as_strided

我们可以利用NumPy strides 为我们提供off-diag 元素作为视图。因此，那里没有内存开销和几乎免费的运行时间！这个想法之前已经在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：Hacky 方式

为了完整起见，一种巧妙的方法是临时保存诊断元素，在那里分配0s，检查所有元素是否为0。如果是这样，则发出对角矩阵的信号。最后分配回诊断元素。

实现将是 -

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's post 的其他方法捕获了其他方法 -

# @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)

Answer 4

我们实际上可以比 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

一个技巧是np.diagonal(a)实际上使用a.diagonal()，所以我们直接使用那个。但最重要的是b 的快速构建，以及b 上的就地操作。

Answer 5

我担心这是否是最有效的方法，但想法是屏蔽对角线元素并检查所有其他元素是否为零。我想这足以将矩阵标记为对角矩阵。

所以我们创建了一个与输入矩阵大小相同的虚拟数组，初始化为 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

Answer 6

我认为这是最简洁的方式：

np.allclose(np.diag(np.diag(a)), a)

Answer 7

获取真相的快速而肮脏的方式。在合理的时间内工作

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

Answer 8

import numpy as np

is_diagonal = (np.trace(mat) == np.sum(mat))