CODEWITHSUNDEEP
pythonpython-3.xnumpyarray-broadcasting
Brad Solomon
Brad Solomon

Reputation: 40878

Replace looping-over-axes with broadcasting, pt 2

Earlier I asked a similar question where the answer used np.dot, taking advantage of the fact that a dot product involves a sum of products. (To my understanding.)

Now I have a similar issue where I don't think dot will apply, because in place of a sum I want to take an element-wise diagonal. If it does, I haven't been able to apply it correctly.

Given a matrix x and array err:

x = np.matrix([[ 0.02984406, -0.00257266],
               [-0.00257266,  0.00320312]])

err = np.array([  7.6363226 ,  13.16548267])

My current implementation with loop is:

res = np.array([np.sqrt(np.diagonal(x * err[i])) for i in range(err.shape[0])])

print(res)
[[ 0.47738755  0.15639712]
 [ 0.62682649  0.20535487]]

which takes the diagonal of x.dot(i) for each i in err. Could this be vectorized? In other words, can the output of x * err be 3-dimensional, with np.diagonal then yielding a 2d array, with one element for each diagonal?

Upvotes: 2

Views: 77

Answers (2)

Divakar
Divakar

Reputation: 221574

Here's an approach making use of diagonal-view with ndarray.flat into x and then use broadcasting for element-wise multiplication, like so -

np.sqrt(x.flat[::x.shape[1]+1].A1 * err[:,None])

Sample run -

In [108]: x = np.matrix([[ 0.02984406, -0.00257266],
     ...:                [-0.00257266,  0.00320312]])
     ...: 
     ...: err = np.array([  7.6363226 ,  13.16548267])
     ...: 

In [109]: np.sqrt(x.flat[::x.shape[1]+1].A1 * err[:,None])
Out[109]: 
array([[ 0.47738755,  0.15639712],
       [ 0.62682649,  0.20535487]])

Runtime test to see how a view helps over np.diagonal that creates a copy -

In [104]: x = np.matrix(np.random.rand(5000,5000))

In [105]: err = np.random.rand(5000)

In [106]: %timeit np.diagonal(x)*err[:,np.newaxis]
10 loops, best of 3: 66.8 ms per loop

In [107]: %timeit x.flat[::x.shape[1]+1].A1 * err[:,None]
10 loops, best of 3: 37.7 ms per loop

Upvotes: 2

hamster on wheels
hamster on wheels

Reputation: 2893

Program:

import numpy as np
x = np.matrix([[ 0.02984406, -0.00257266],
               [-0.00257266,  0.00320312]])

err = np.array([  7.6363226 ,  13.16548267])
diag = np.diagonal(x)
ans = np.sqrt(diag*err[:,np.newaxis])  # sqrt of outer product
print(ans)

# use out keyword to avoid making new numpy array for many times.
ans = np.empty(x.shape, dtype=x.dtype)
for i in range(100):
    ans = np.multiply(diag, err, out=ans)
    ans = np.sqrt(ans, out=ans)

Result:

[[ 0.47738755  0.15639712]
 [ 0.62682649  0.20535487]]

Upvotes: 2

Related Questions