Replace looping-over-axes with broadcasting, pt 2

Question

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?

hamster on wheels · Accepted Answer · 2017-07-19T22:41:00.407

2

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

edited Jul 19 '17 at 22:41

answered Jul 19 '17 at 22:19

hamster on wheels

2,771
17
50

could use out keyword to avoid allocating new numpy array, which is slow. – hamster on wheels Jul 19 '17 at 22:20
sqrt faster than **0.5 in old version of python. not sure about new. https://stackoverflow.com/questions/327002/which-is-faster-in-python-x-5-or-math-sqrtx I am usually too lazy to type sqrt – hamster on wheels Jul 19 '17 at 22:22
don't worry about allocating new array if you only do this for once, instead of many times. – hamster on wheels Jul 19 '17 at 22:37
are you writing some kind of least square fitting? – hamster on wheels Jul 19 '17 at 22:41
Yes, with the intention of generalizing to cases with more than one `y` vector. This finds the standard errors of the coefficients. Why do you ask? – Brad Solomon Jul 20 '17 at 00:41
Have you searched if anyone has written code/library for that? writing new code is the last resort. – hamster on wheels Jul 20 '17 at 01:04
Yes, there are `statsmodels` and `sklearn` but I'm unsatisfied with both and in the process of tweaking them – Brad Solomon Jul 20 '17 at 01:12

score 2 · Answer 2 · answered Jul 20 '17 at 09:55

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

Replace looping-over-axes with broadcasting, pt 2

2 Answers2