Check if a large matrix is diagonal matrix in python

We can actually do quite a bit better than what Daniel F suggested:

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)

results in

2.5737204551696777
0.6501829624176025

One tricks is that np.diagonal(a) actually uses a.diagonal(), so we use that one directly. But what takes the cake the the fast build of b, combined with the in-place operation on b.

Not sure how fast this is compared to the others, but:

def isDiag(M):
    i, j = np.nonzero(M)
    return np.all(i == j)

EDIT Let's time things:

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

Hmm, that's slower, probably from constructing i and j

Let's try to improve the @donkopotamus solution by removing the subtraction step:

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

That's a bit better.

EDIT2 I came up with an even faster method:

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:])

This isn't doing any calculations, just reshaping. Turns out reshaping to +1 rows on a diagonal matrix puts all the data in the first column. You can then check a contiguous block for any nonzeros which is much fatser for 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

Seems my original is faster than @Make42 for smaller sets

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

And @Make42 gives memory error on the larger set. But then I don't seem to have as much RAM as they do.

Remove the diagonal and count the non zero elements:

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

Appproach #1 : Using NumPy strides/np.lib.stride_tricks.as_strided

We can leverage NumPy strides to give us the off-diag elements as a view. So, no memory overhead there and virtually free runtime! This idea has been explored before in this post.

Thus, we have -

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

Sample run to showcase its usage -

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]])

Let's verify the free runtime and no memory overhead claims, by using it on a large array -

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

Now, how do we use it here? Simply check if all nodiag_view elements are 0s, signalling a diagonal matrix!

Hence, to solve our case here, for an input array a, it would be -

isdiag = (nodiag_view(a)==0).all()

Appproach #2 : Hacky way

For completeness, one hacky way would be to temporarily save diag elements, assign 0s there, check for all elements to be 0s. If so, signalling a diagonal matrix. Finally assign back the diag elements.

The implementation would be -

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

Benchmarking

Let's time on a large array and see how these compare on performance -

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)

Other approaches from @Daniel F's post that captured other approaches -

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

Timings with same setup as earlier -

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)