Unable to broadcast argument 1 to output array

[andrew-pomerance]

This code

@jit(nopython=True)
def _sgd_batch(learning_rate, w, xtrain, ytrain):
    for i in range(xtrain.shape[1]):
        wx = w*xtrain[:,i]
        error = wx - ytrain[:,i]
        for j in range(w.shape[0]):
            w[j,:] = w[j,:]-learning_rate*error[j]*np.transpose(xtrain[:,i])
    return w

throws a ValueError: unable to broadcast argument 1 to output array exception when called with the following parameters: learning_rate is a scalar, w is a 3x1000 matrix, xtrain is a 1000 x 10k matrix, and ytrain is a 3 x 100k matrix. When I run it without @jit it works just fine. It doesn't appear to have anything to do with w because reordering the parameters doesn't change anything.

There doesn't seem to be any google-discoverable documentation on this error message. I'd be happy to write up a couple paragraphs if someone can explain to me what's going on and what causes this error message. Again, without @jit it is fine so it's not an underlying Numpy broadcast problem.

Thanks,
Andrew

[rjenc29]

Hi, what happens if you execute the following?

@jit(nopython=True)
def _sgd_batch(learning_rate, w, xtrain, ytrain):
    for i in range(xtrain.shape[1]):
        wx = w*xtrain[:,i]
        error = wx - np.atleast_2d(ytrain[:,i]).T
        for j in range(w.shape[0]):
            w[j,:] = w[j,:]-learning_rate*error[j]*np.transpose(xtrain[:,i])
    return w

I'm not sure but it looks like the error term might be the problem.

[andrew-pomerance]

Thanks! It does something... weird... The code you wrote compiles with Numba just fine but gives the wrong answer (order ~1e+18 instead of ~100). If I take out the @jit, I get ValueError: shapes (1,3) and (1,1000) not aligned: 3 (dim 1) != 1 (dim 0) from NumPy.

I did a bit of digging and according to some debug statements placed in numpy_ufunc_kernel() in npyimpl.py, the expected type of the line error = (w*xtrain[:,i]) - ytrain[:,i] (I refactored) is (array(float64, 1d, F), array(float64, 1d, F), array(float64, 2d, C)) -> array(float64, 2d, C) for some reason. That suggested that I should try this function:

@jit(nopython=True)
def _sgd_batch(w, learning_rate,  xtrain, ytrain):
    for i in range(xtrain.shape[1]):
        error = (w*xtrain[:,i])[:,0] - ytrain[:,i]  # note forcing the implied 2D array down to 1D
        for j in range(w.shape[0]):
            w[j,:] -= learning_rate*error[j]*np.transpose(xtrain[:,i])
    return w

This function compiles but gives the wrong answer when @jitted, and works just fine when not. Then I tried this:

@jit(nopython=True)
def _sgd_batch(w, learning_rate,  xtrain, ytrain):
    for i in range(xtrain.shape[1]):
        for j in range(w.shape[0]):
            error = np.dot(w[j,:], xtrain[:,i]) - ytrain[j,i]
            w[j,:] -= learning_rate*error*np.transpose(xtrain[:,i])
    return w

It works! However, is this a bug? If not, what is the logic that suggests that my original function needs to be written like this to work? Also, the fact that one function compiles in numba but fails in NumPy and another gives different answers between the two makes me think this is a bug.

Thanks,
Andrew

[esc]

@andrew-pomerance thank you very much for your bug report and taking the time to perform such a detailed analysis. I expect this will take some digging, so I will label it as 'needtriage' for the time being.

[stuartarchibald]

This doesn't seem to work in plain NumPy:

from numba import njit
import numpy as np

##### NOT USING JIT! @njit
def _sgd_batch(learning_rate, w, xtrain, ytrain):
    for i in range(xtrain.shape[1]):
        wx = w*xtrain[:,i]
        error = wx - ytrain[:,i]
        for j in range(w.shape[0]):
            w[j,:] = w[j,:]-learning_rate*error[j]*np.transpose(xtrain[:,i])
    return w

learning_rate = 1.2
factor = 100
a = 3
b = 1000
c  = 10000
d  = 100000

w = np.zeros((a, b))
xtrain = np.zeros((b, c))
ytrain = np.zeros((a, d))

def print_shp(name, thing):
    print("%s: %s" % (name, thing.shape))


print_shp('w', w)
print_shp('xtrain', xtrain)
print_shp('ytrain', ytrain)

_sgd_batch(learning_rate, w, xtrain, ytrain)

gives:

$ python issue4632.py 
w: (3, 1000)
xtrain: (1000, 10000)
ytrain: (3, 100000)
Traceback (most recent call last):
  File "issue4632.py", line 32, in <module>
    _sgd_batch(learning_rate, w, xtrain, ytrain)
  File "issue4632.py", line 8, in _sgd_batch
    error = wx - ytrain[:,i]
ValueError: operands could not be broadcast together with shapes (3,1000) (3,) 

[andrew-pomerance]

I call it with xtrain/ytrain as matrices. I have no idea if that is what I should be doing, but it seems that NumPy broadcasting is doing some thing that I can't seem to understand which makes the shapes of slices different than what I expect. But np.matrices work as expect so I call with:

_sgd_batch(learning_rate, np.matrix(w), np.matrix(xtrain), np.matrix(ytrain))

and it works fine. Maybe that's part of the problem?

[stuartarchibald]

I suspect there's some nuance to do with the slight difference in semantics that arises from the use of np.matrix that Numba is not capturing.

[lanougue]

Hello Everyone,
Is there any update about this issue ?
I tried the latest version of numba (0.55.1) but I still have the same problem.
Thanks for the help

[qris]

Thanks, this gave me the hint I needed to work around the problem - pass the scalar parameters as a 1x1 Numpy array instead of a Python variable.