Diagonal Matrix Exponential in Python

Question

I'm writing a numerical algorithm with speed in mind. I've come across the two matrix exponential functions in scipy/numpy (scipy.linalg.expm2, scipy.linalg.expm). However I have a matrix that I know to be diagonal beforehand. Do these scipy functions check if the matrix is diagonal before they run? Obviously the exponentiation algorithm can be much faster for a diagonal matrix, and I just want to make sure that these are doing something smart with that - if they aren't, is there an easy way to do it?

Answer 1

If a matrix is diagonal, then its exponential can be obtained by just exponentiating every entry on the main diagonal, so you can calculate it by:

np.diag(np.exp(np.diag(a)))

Answer 2

If you know A is diagonal and you want the k-th power:

def dpow(a, k):
    return np.diag(np.diag(a) ** k)

Check if a matrix is diagonal:

def isdiag(a):
    return np.all(a == np.diag(np.diag(a)))

so :

def pow(a, k):
    if isdiag(a):
        return dpow(a, k)
    else:
        return np.asmatrix(a) ** k

Similarly for exponential (which you can get mathematically from the expansion of a suite of pow) you can do:

def dexp(a, k):
    return np.diag(np.exp(np.diag(a)))

def exp(a, k):
    if isdiag(a):
        return dexp(a, k)
    else:
        #use scipy.linalg.expm2 or whatever

Answer 3

I've developed a tool that can help being faster doing the same as HYRY but by doing it in-place:

def diagonal(array):
    """ Return a **view** of the diagonal elements of 'array' """
    from numpy.lib.stride_tricks import as_strided
    return as_strided(array,shape=(min(array.shape),),strides=(sum(array.strides),))

# generate a random diagonal array
d  = np.diag(np.random.random(4000))

# in-place exponent of the diagonal elements
ddiag = diagonal(d)
ddiag[:] = np.exp(ddiag)

# timeit comparison with HYRY's method
%timeit -n10 np.diag(np.exp(np.diag(d)))   
    # out> 10 loops, best of 3: 52.1 ms per loop
%timeit -n10 ddiag = diagonal(d); ddiag[:] = np.exp(ddiag)
    # out> 10 loops, best of 3: 108 µs per loop

Now,

• HYRY's method is quadratic w.r.t the diagonal length (probably because of the new array memory allocation), and so if your matrices are of little dimension, the difference might not be as big.

• you need to be alright with in-place computation

• Finally, the off-diagonal elements are 0, so their exponential should be 1, shouldn't it ? In both our method the off-diagonal are 0.

For that last part, if you want all off-diagonal elements to be 1, then you can do:

d2 = np.ones_like(d); 
diagonal(d2)[:] = np.exp(np.diag(d))

print (d2==np.exp(d)).all()  # True

But this is linear w.r.t to array size, so quadratic w.r.t to diagonal length. The timeit gives ~90ms for a 4000x4000 array and 22.3ms for a 2000x2000.

Finally, you can also do it in-place to get a little speed up:

diag = np.diag(d)
d[:]=1
diagonal(d)[:] = np.exp(diag)

Timeit gives 66.1ms for 4000^2 array, and 16.8ms for 2000^2