scikit-cuda
scikit-cuda copied to clipboard
lebedov
misc.*_matvec broadcasting slightly different from numpy's
The matvec functions can align the vector with axis 0 of the matrix, but numpy always aligns from the right. Moreover, the matvec functions prefer axis 0, which can lead to discrepancies:
x22 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,2)))
a2 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,)))
# These produce different results
print skcuda.misc.add_matvec(x22, a2)
print x22.get() + a2.get()
x23 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,3)))
a2 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,)))
a3 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (3,)))
# This is fine
print skcuda.misc.add_matvec(x23, a3)
print x23.get() + a3.get()
# This causes an exception in numpy but not in add_matvec
print skcuda.misc.add_matvec(x23, a2)
print x23.get() + a2.get()
@davidweichiang I'm a bit weirded out by your last example:
# This causes an exception in numpy but not in add_matvec
print skcuda.misc.add_matvec(x23, a2)
print x23.get() + a2.get()
Do you know why numpy doesn't support this type of broadcasting? I'm guessing it's because numpy matches from the right and thus immediately throws an exception. However, since skcuda only supports 2D matrices right now, we could safely support it. But it might surprise people who are used to numpy's stricter rules (and might bit us if skcuda ever goes on to support more-dimensional tensors)
What do you guys think we should do here?
I don't know (but could take a couple of guesses) why Numpy matches from the right, but I definitely think it's good that Numpy chooses a single direction. As for the matvec functions, I don't know if I get a vote, but my suggestion would actually be to have a function add(x, y) that follows the Numpy rules exactly for any combination of 0, 1, or 2 dimensions; otherwise, delegate to GPUArray which handles 0 with 3+ and 3+ with 0.
which handles 0 with 3+ and 3+ with 0.
Sorry, but I do not understand what you mean by that?
Sorry, I mean a 0d array (or a scalar) with an array with 3 or more dimensions.
Ah, I see! The question is if there's a need for it (if someone has 3+ dimensions, they likely won't use skcuda anyhow)?
In any case, I also think that add_matvec should closely mirror numpy's broadcasting rules. If the need for a more general add (that handles 0d/3d+ arrays) also, it can probably be added later anyhow.
Agreed - mirroring numpy's rules seems desirable, and we shouldn't worry about arrays with more than 2 dimensions until we have to.
Thanks! The last case still has two different behaviors (skcuda likes it, numpy doesn't). I'll try a version and you can decide whether it's appropriate?
Are you sure you aren't working on an outdated version? Using the current master branch, both of the following lines throw a ValueError on my machine:
print(skcuda.misc.add_matvec(x23, a2))
print(x23.get() + a2.get())
Sorry, you're right! I got confused by lines 885-888, which cover the case len(a_gpu.shape) != 1. These lines may cause surprising behavior in cases like these:
x23 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,3)))
a43 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (4,3)))
a24 = pycuda.gpuarray.to_gpu(numpy.random.uniform(0., 1., (2,4)))
#print x23.get() + a43.get() # causes error
print skcuda.misc.add_matvec(x23, a43)
#print x23.get() + a24.get() # causes error
print skcuda.misc.add_matvec(x23, a24)
But since the function has "matvec" in its name, I'd suggest just throwing an exception if a is not a matrix or b is not a vector. Then lines 885-892 should not be needed unless I've missed something. Again, sorry for the false alarm.
@lebedov, please see pull request #125. I think all three of the shortcomings mentioned there could be fixed, if this kind of generality is something you want.