CODEWITHSUNDEEP
numpymatrix-multiplication
LMB
LMB

Reputation: 464

Python numpy matrix multiplication with one diagonal matrix

I have two arrays A (4000,4000) of which only the diagonal is filled with data, and B (4000,5), filled with data. Is there a way to multiply (dot) these arrays that is faster than the numpy.dot(a,b) function?

So far I found that (A * B.T).T should be faster (where A is one dimensional (4000,), filled with the diagonal elements), but it turned out to be roughly twice as slow.

is there a faster way to calculate B.dot(A) in the case where A is a diagnal array?

Upvotes: 2

Views: 11935

Answers (2)

LMB
LMB

Reputation: 464

Appears that my initial claim of (A * B.T).T being slower is incorrect.

from timeit import default_timer as timer
import numpy as np

##### Case 1
a = np.zeros((4000,4000))
np.fill_diagonal(a, 10)
b = np.ones((4000,5))

dot_list = []

def time_dot(a,b):
    start = timer()
    c = np.dot(a,b)
    end = timer()
    return end - start

for i in range(100):
    dot_list.append(time_dot(a,b))

print np.mean(np.asarray(dot_list))

##### Case 2
a = np.ones((4000,))
a = a * 10
b = np.ones((4000,5))

shortcut_list = []

def time_quicker(a,b):
    start = timer()
    c = (a*b.T).T
    end = timer()
    return end - start

for i in range(100):
    shortcut_list.append(time_quicker(a,b))

print np.mean(np.asarray(shortcut_list))


##### Case 3
a = np.zeros((4000,4000)) #diagonal matrix
np.fill_diagonal(a, 10)
b = np.ones((4000,5))

case3_list = []

def function(a,b):
    start = timer()
    np.multiply(b.T,np.diag(a))
    end = timer()
    return end - start

for i in range(100):
    case3_list.append(function(a,b))

print np.mean(np.asarray(case3_list))

results in :

0.119120892431

0.00010633951868

0.00214490709662

so the second method is fastest

Upvotes: 0

Divakar
Divakar

Reputation: 221504

You could simply extract the diagonal elements and then perform broadcasted elementwise multiplication.

Thus, a replacement for B*A would be -

np.multiply(np.diag(B)[:,None], A)

and for A.T*B -

np.multiply(A.T,np.diag(B))

Runtime test -

In [273]: # Setup
     ...: M,N = 4000,5
     ...: A = np.random.randint(0,9,(M,N)).astype(float)
     ...: B = np.zeros((M,M),dtype=float)
     ...: np.fill_diagonal(B, np.random.randint(11,99,(M)))
     ...: A = np.matrix(A)
     ...: B = np.matrix(B)
     ...: 

In [274]: np.allclose(B*A, np.multiply(np.diag(B)[:,None], A))
Out[274]: True

In [275]: %timeit B*A
10 loops, best of 3: 32.1 ms per loop

In [276]: %timeit np.multiply(np.diag(B)[:,None], A)
10000 loops, best of 3: 33 µs per loop

In [282]: np.allclose(A.T*B, np.multiply(A.T,np.diag(B)))
Out[282]: True

In [283]: %timeit A.T*B
10 loops, best of 3: 24.1 ms per loop

In [284]: %timeit np.multiply(A.T,np.diag(B))
10000 loops, best of 3: 36.2 µs per loop

Upvotes: 8

Related Questions