Joe*_*mes 5 python numpy matrix
我已经计算了一个非常大的矩阵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
这样好一些。
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 Let's check times来说更胖：
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 在较大的集合上给出了内存错误。但后来我似乎没有他们那么多的内存。
方法#1：使用 NumPy strides/np.lib.stride_tricks.as_strided
我们可以利用NumPy strides提供非诊断元素作为视图。因此，没有内存开销，并且几乎免费的运行时间！这个想法之前已经被探索过this post。
因此，我们有——
\n# /sf/answers/3063335901/ @Divakar\ndef nodiag_view(a):\n m = a.shape[0]\n p,q = a.strides\n return np.lib.stride_tricks.as_strided(a[:,1:], (m-1,m), (p+q,q))\n示例运行以展示其用法 -
\nIn [175]: a\nOut[175]: \narray([[ 0, 1, 2, 3],\n [ 4, 5, 6, 7],\n [ 8, 9, 10, 11],\n [12, 13, 14, 15]])\n\nIn [176]: nodiag_view(a)\nOut[176]: \narray([[ 1, 2, 3, 4],\n [ 6, 7, 8, 9],\n [11, 12, 13, 14]])\n让我们通过在大型数组上使用它来验证免费运行时和无内存开销声明 -
\nIn [182]: a = np.zeros((10000,10000), dtype=int)\n ...: np.fill_diagonal(a,np.arange(len(a)))\n\nIn [183]: %timeit nodiag_view(a)\n6.42 \xc2\xb5s \xc2\xb1 48.2 ns per loop (mean \xc2\xb1 std. dev. of 7 runs, 100000 loops each)\n\nIn [184]: np.shares_memory(a, nodiag_view(a))\nOut[184]: True\n现在，我们在这里如何使用它？只需检查所有nodiag_view元素是否均为 0，表示这是一个对角矩阵！
因此，为了解决我们这里的情况，对于输入数组a，它将是 -
isdiag = (nodiag_view(a)==0).all()\n方法#2：Hacky 方式
\n为了完整起见，一种巧妙的方法是临时保存 diag 元素，0s在那里分配，检查所有元素是否为 0。如果是，则发出对角矩阵信号。最后分配回 diag 元素。
实施将是 -
\ndef hacky_way(a):\n diag_elem = np.diag(a).copy()\n np.fill_diagonal(a,0)\n out = (a==0).all()\n np.fill_diagonal(a,diag_elem)\n return out\n让我们花时间在一个大型数组上，看看它们的性能比较如何 -
\nIn [3]: a = np.zeros((10000,10000), dtype=int)\n ...: np.fill_diagonal(a,np.arange(len(a)))\n\nIn [4]: %timeit (nodiag_view(a)==0).all()\n52.3 ms \xc2\xb1 393 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n\nIn [5]: %timeit hacky_way(a)\n51.8 ms \xc2\xb1 250 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n@Daniel F\'s 帖子中的其他方法捕获了其他方法 -
\n# @donkopotamus solution improved by @Daniel F\ndef b(M):\n return np.all(M == np.diag(np.diagonal(M)))\n\n# @Daniel F\'s soln without assert check\ndef isDiag2(M):\n i, j = M.shape\n test = M.reshape(-1)[:-1].reshape(i-1, j+1)\n return ~np.any(test[:, 1:])\n\n# @Make42\'s soln\ndef Make42(m):\n b = np.zeros(m.shape)\n np.fill_diagonal(b, m.diagonal())\n return np.all(m == b)\n计时设置与之前相同 -
\nIn [6]: %timeit b(a)\n ...: %timeit Make42(a)\n ...: %timeit isDiag2(a)\n218 ms \xc2\xb1 1.68 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n302 ms \xc2\xb1 1.25 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n67.1 ms \xc2\xb1 1.35 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n