[Numpy-discussion] Scipy dot
Nathaniel Smith
njs at pobox.com
Fri Nov 9 03:32:52 EST 2012
On Fri, Nov 9, 2012 at 6:18 AM, Nicolas SCHEFFER
<scheffer.nicolas at gmail.com> wrote:
> Fred,
>
> Thanks for the advice.
> The code will only affect the part in _dotblas.c where gemm is called.
> There's tons of check before that make sure both matrices are of ndim 2.
> We should check though if we can do these tricks in other parts of the function.
>
> Otherwise:
> - I've built against ATLAS 3.10
> - I'm happy to add a couple more test for C and F-contiguous. I'm not
> sure how to get the third type (strided), would you have an example?
def with_memory_order(a, order):
assert order in ("C", "F", "discontig")
assert a.ndim == 2
if order in ("C", "F"):
return np.asarray(a, order=order)
else:
buf = np.empty((a.shape[0] * 2, a.shape[1] * 2), dtype=a.dtype)
buf[::2, ::2] = a
# This returns a view onto every other element of 'buf':
result = buf[::2, ::2]
assert not result.flags.c_contiguous and not result.flags.f_contiguous
return result
> The following test for instance checks integrity against
> multiarray.dot, which I believe is default when not compiled with
> BLAS.
> Dot is a hard function to test imho, so if anybody has ideas on what
> kind of test they'd like to see, please let me know.
>
> If that's ok I might now be able to:
> - Check for more bugs, I need to dig a bit more in the gemm call, make
> sure everything is ok.
> - Create an issue on github and link to this discussion
> - Make a commit in a seperate branch
> - Move forward like that.
>
> ==
> import numpy as np
> from time import time
> from numpy.testing import assert_almost_equal
>
> def test_dot_regression():
> """ Test numpy dot by comparing with multiarray dot
> """
> np.random.seed(7)
> a = np.random.randn(3, 3)
> b = np.random.randn(3, 2)
> c = np.random.randn(2, 3)
>
> _dot = np.core.multiarray.dot
>
> assert_almost_equal(np.dot(a, a), _dot(a, a))
> assert_almost_equal(np.dot(b, c), _dot(b, c))
> assert_almost_equal(np.dot(b.T, c.T), _dot(b.T, c.T))
>
> assert_almost_equal(np.dot(a.T, a), _dot(a.T, a))
> assert_almost_equal(np.dot(a, a.T), _dot(a, a.T))
> assert_almost_equal(np.dot(a.T, a.T), _dot(a.T, a.T))
You should check that the result is C-contiguous in all cases too.
for a_order in ("C", "F", "discontig"):
for b_order in ("C", "F", "discontig"):
this_a = with_memory_order(a, a_order)
this_b = with_memory_order(b, b_order)
result = np.dot(this_a, this_b)
assert_almost_equal(result, expected)
assert result.flags.c_contiguous
You could also wrap the above in yet another loop to try a few
different combinations of a and b matrices (perhaps after sticking the
code into a utility function, like run_dot_tests(a, b, expected), so
the indentation doesn't get out of hand ;-)). Then you can easily test
some of the edge cases, like Nx1 matrices.
-n
> On Thu, Nov 8, 2012 at 5:34 PM, Frédéric Bastien <nouiz at nouiz.org> wrote:
>> Hi,
>>
>> I suspect the current tests are not enought. You need to test all the
>> combination for the 3 inputs with thoses strides:
>>
>> c-contiguous
>> f-contiguous
>> something else like strided.
>>
>> Also, try with matrix with shape of 1 in each dimensions. Not all blas
>> libraries accept the strides that numpy use in that cases. Also, not all
>> blas version accept the same stuff, so if this isn't in the current version,
>> there will be probably some adjustment later on that side. What blas do you
>> use? I think ATLAS was one that was causing problem.
>>
>>
>> When we did this in Theano, it was more complicated then this diff... But
>> much of the code is boillerplate code.
>>
>> Fred
>>
>>
>>
>> On Thu, Nov 8, 2012 at 8:03 PM, Nicolas SCHEFFER
>> <scheffer.nicolas at gmail.com> wrote:
>>>
>>> Thanks Sebastien, didn't think of that.
>>>
>>> Well I went ahead and tried the change, and it's indeed straightforward.
>>>
>>> I've run some tests, among which:
>>> nosetests numpy/numpy/core/tests/test_blasdot.py
>>> and it looks ok. I'm assuming this is good news.
>>>
>>> I've copy-pasting the diff below, but I have that in my branch and can
>>> create a PR if we agree on it.
>>> I still cannot believe it's that easy (well this has been bugging me a
>>> while... ;))
>>> So I wouldn't mind waiting a day or two to see reactions on the list
>>> before moving ahead.
>>>
>>> diff --git a/numpy/core/blasdot/_dotblas.c b/numpy/core/blasdot/_dotblas.c
>>> index c73dd6a..2b4be7c 100644
>>> --- a/numpy/core/blasdot/_dotblas.c
>>> +++ b/numpy/core/blasdot/_dotblas.c
>>> @@ -770,7 +770,7 @@ dotblas_matrixproduct(PyObject *NPY_UNUSED(dummy),
>>> PyObject *args, PyObject* kwa
>>> * using appropriate values of Order, Trans1, and Trans2.
>>> */
>>>
>>> - if (!PyArray_ISCONTIGUOUS(ap2)) {
>>> + if (!PyArray_IS_C_CONTIGUOUS(ap2) &&
>>> !PyArray_IS_F_CONTIGUOUS(ap2)) {
>>> PyObject *new = PyArray_Copy(ap2);
>>>
>>> Py_DECREF(ap2);
>>> @@ -779,7 +779,7 @@ dotblas_matrixproduct(PyObject *NPY_UNUSED(dummy),
>>> PyObject *args, PyObject* kwa
>>> goto fail;
>>> }
>>> }
>>> - if (!PyArray_ISCONTIGUOUS(ap1)) {
>>> + if (!PyArray_IS_C_CONTIGUOUS(ap1) &&
>>> !PyArray_IS_F_CONTIGUOUS(ap1)) {
>>> PyObject *new = PyArray_Copy(ap1);
>>>
>>> Py_DECREF(ap1);
>>> @@ -800,6 +800,19 @@ dotblas_matrixproduct(PyObject
>>> *NPY_UNUSED(dummy), PyObject *args, PyObject* kwa
>>> lda = (PyArray_DIM(ap1, 1) > 1 ? PyArray_DIM(ap1, 1) : 1);
>>> ldb = (PyArray_DIM(ap2, 1) > 1 ? PyArray_DIM(ap2, 1) : 1);
>>> ldc = (PyArray_DIM(ret, 1) > 1 ? PyArray_DIM(ret, 1) : 1);
>>> +
>>> + /*
>>> + * Avoid temporary copies for arrays in Fortran order
>>> + */
>>> + if (PyArray_IS_F_CONTIGUOUS(ap1)) {
>>> + Trans1 = CblasTrans;
>>> + lda = (PyArray_DIM(ap1, 0) > 1 ? PyArray_DIM(ap1, 0) : 1);
>>> + }
>>> + if (PyArray_IS_F_CONTIGUOUS(ap2)) {
>>> + Trans2 = CblasTrans;
>>> + ldb = (PyArray_DIM(ap2, 0) > 1 ? PyArray_DIM(ap2, 0) : 1);
>>> + }
>>> +
>>> if (typenum == NPY_DOUBLE) {
>>> cblas_dgemm(Order, Trans1, Trans2,
>>> L, N, M,
>>>
>>> On Thu, Nov 8, 2012 at 3:58 PM, Sebastian Berg
>>> <sebastian at sipsolutions.net> wrote:
>>> > On Fri, 2012-11-09 at 00:24 +0100, Sebastian Berg wrote:
>>> >> Hey,
>>> >>
>>> >> On Thu, 2012-11-08 at 14:44 -0800, Nicolas SCHEFFER wrote:
>>> >> > Well, hinted by what Fabien said, I looked at the C level dot
>>> >> > function.
>>> >> > Quite verbose!
>>> >> >
>>> >> > But starting line 757, we can see that it shouldn't be too much work
>>> >> > to fix that bug (well there is even a comment there that states just
>>> >> > that)
>>> >> >
>>> >> > https://github.com/numpy/numpy/blob/master/numpy/core/blasdot/_dotblas.c#L757
>>> >> > I now think that should be the cleanest.
>>> >> >
>>> >> > This would only work for gemm though. I don't know what the benefit
>>> >> > is
>>> >> > for gemv for instance, but we should make that kind of changes
>>> >> > everywhere we can.
>>> >> > The evil PyArray_Copy is there twice and that's what we want to get
>>> >> > rid of.
>>> >> >
>>> >> > I'm not sure, but it looks to me that removing the copy and doing the
>>> >> > following would do the work:
>>> >> > Order = CblasRowMajor;
>>> >> > Trans1 = CblasNoTrans;
>>> >> > Trans2 = CblasNoTrans;
>>> >> > if (!PyArray_ISCONTIGUOUS(ap1)) {
>>> >> > Trans1 = CblasTrans;
>>> >> > }
>>> >> > if (!PyArray_ISCONTIGUOUS(ap2)) {
>>> >> > Trans2 = CblasTrans;
>>> >> > }
>>> >> > might be too easy to be true.
>>> >> >
>>> >>
>>> >> Sounds nice, though don't forget that the array may also be neither C-
>>> >> or F-Contiguous, in which case you need a copy in any case. So it would
>>> >> probably be more like:
>>> >>
>>> >> if (PyArray_IS_C_CONTIGUOUS(ap1)) {
>>> >> Trans1 = CblasNoTrans;
>>> >> }
>>> >> else if (PyArray_IS_F_CONTIGUOUS(ap1)) {
>>> >> Trans1 = CblasTrans;
>>> >> }
>>> >> else {
>>> >> Trans1 = CblasNoTrans;
>>> >> PyObject *new = PyArray_Copy(ap1);
>>> >> Py_DECREF(ap1);
>>> >> ap1 = (PyArrayObject *)new;
>>> >> }
>>> >>
>>> >
>>> > Well, of course I forgot error checking there, and maybe you need to set
>>> > some of the other parameters differently, but it looks like its probably
>>> > that easy, and I am sure everyone will welcome a PR with such changes.
>>> >
>>> >> Regards,
>>> >>
>>> >> Sebastian
>>> >>
>>> >> >
>>> >> >
>>> >> > On Thu, Nov 8, 2012 at 12:06 PM, Nicolas SCHEFFER
>>> >> > <scheffer.nicolas at gmail.com> wrote:
>>> >> > > I've made the necessary changes to get the proper order for the
>>> >> > > output array.
>>> >> > > Also, a pass of pep8 and some tests (fixmes are in failing tests)
>>> >> > > http://pastebin.com/M8TfbURi
>>> >> > >
>>> >> > > -n
>>> >> > >
>>> >> > > On Thu, Nov 8, 2012 at 11:38 AM, Nicolas SCHEFFER
>>> >> > > <scheffer.nicolas at gmail.com> wrote:
>>> >> > >> Thanks for all the responses folks. This is indeed a nice problem
>>> >> > >> to solve.
>>> >> > >>
>>> >> > >> Few points:
>>> >> > >> I. Change the order from 'F' to 'C': I'll look into it.
>>> >> > >> II. Integration with scipy / numpy: opinions are diverging here.
>>> >> > >> Let's wait a bit to get more responses on what people think.
>>> >> > >> One thing though: I'd need the same functionality as
>>> >> > >> get_blas_funcs in numpy.
>>> >> > >> Since numpy does not require lapack, what functions can I get?
>>> >> > >> III. Complex arrays
>>> >> > >> I unfortunately don't have enough knowledge here. If someone could
>>> >> > >> propose a fix, that'd be great.
>>> >> > >> IV. C
>>> >> > >> Writing this in C sounds like a good idea. I'm not sure I'd be the
>>> >> > >> right person to this though.
>>> >> > >> V. Patch in numpy
>>> >> > >> I'd love to do that and learn to do it as a byproduct.
>>> >> > >> Let's make sure we agree this can go in numpy first and that all
>>> >> > >> FIXME
>>> >> > >> can be fixed.
>>> >> > >> Although I guess we can resolve fixmes using git.
>>> >> > >>
>>> >> > >> Let me know how you'd like to proceed,
>>> >> > >>
>>> >> > >> Thanks!
>>> >> > >>
>>> >> > >> FIXMEs:
>>> >> > >> - Fix for ndim != 2
>>> >> > >> - Fix for dtype == np.complex*
>>> >> > >> - Fix order of output array
>>> >> > >>
>>> >> > >> On Thu, Nov 8, 2012 at 9:42 AM, Frédéric Bastien <nouiz at nouiz.org>
>>> >> > >> wrote:
>>> >> > >>> Hi,
>>> >> > >>>
>>> >> > >>> I also think it should go into numpy.dot and that the output
>>> >> > >>> order should
>>> >> > >>> not be changed.
>>> >> > >>>
>>> >> > >>> A new point, what about the additional overhead for small
>>> >> > >>> ndarray? To remove
>>> >> > >>> this, I would suggest to put this code into the C function that
>>> >> > >>> do the
>>> >> > >>> actual work (at least, from memory it is a c function, not a
>>> >> > >>> python one).
>>> >> > >>>
>>> >> > >>> HTH
>>> >> > >>>
>>> >> > >>> Fred
>>> >> > >>>
>>> >> > >>>
>>> >> > >>>
>>> >> > >>> On Thu, Nov 8, 2012 at 12:29 PM, Anthony Scopatz
>>> >> > >>> <scopatz at gmail.com> wrote:
>>> >> > >>>>
>>> >> > >>>> On Thu, Nov 8, 2012 at 7:06 AM, David Cournapeau
>>> >> > >>>> <cournape at gmail.com>
>>> >> > >>>> wrote:
>>> >> > >>>>>
>>> >> > >>>>> On Thu, Nov 8, 2012 at 12:12 PM, Dag Sverre Seljebotn
>>> >> > >>>>> <d.s.seljebotn at astro.uio.no> wrote:
>>> >> > >>>>> > On 11/08/2012 01:07 PM, Gael Varoquaux wrote:
>>> >> > >>>>> >> On Thu, Nov 08, 2012 at 11:28:21AM +0000, Nathaniel Smith
>>> >> > >>>>> >> wrote:
>>> >> > >>>>> >>> I think everyone would be very happy to see numpy.dot
>>> >> > >>>>> >>> modified to do
>>> >> > >>>>> >>> this automatically. But adding a scipy.dot IMHO would be
>>> >> > >>>>> >>> fixing
>>> >> > >>>>> >>> things
>>> >> > >>>>> >>> in the wrong place and just create extra confusion.
>>> >> > >>>>> >>
>>> >> > >>>>> >> I am not sure I agree: numpy is often compiled without
>>> >> > >>>>> >> lapack support,
>>> >> > >>>>> >> as
>>> >> > >>>>> >> it is not necessary. On the other hand scipy is always
>>> >> > >>>>> >> compiled with
>>> >> > >>>>> >> lapack. Thus this makes more sens in scipy.
>>> >> > >>>>> >
>>> >> > >>>>> > Well, numpy.dot already contains multiple fallback cases for
>>> >> > >>>>> > when it is
>>> >> > >>>>> > compiled with BLAS and not. So I'm +1 on just making this an
>>> >> > >>>>> > improvement
>>> >> > >>>>> > on numpy.dot. I don't think there's a time when you would not
>>> >> > >>>>> > want to
>>> >> > >>>>> > use this (provided the output order issue is fixed), and it
>>> >> > >>>>> > doesn't
>>> >> > >>>>> > make
>>> >> > >>>>> > sense to not have old codes take advantage of the speed
>>> >> > >>>>> > improvement.
>>> >> > >>>>>
>>> >> > >>>>> Indeed, there is no reason not to make this available in NumPy.
>>> >> > >>>>>
>>> >> > >>>>> Nicolas, can you prepare a patch for numpy ?
>>> >> > >>>>
>>> >> > >>>>
>>> >> > >>>> +1, I agree, this should be a fix in numpy, not scipy.
>>> >> > >>>>
>>> >> > >>>> Be Well
>>> >> > >>>> Anthony
>>> >> > >>>>
>>> >> > >>>>>
>>> >> > >>>>>
>>> >> > >>>>> David
>>> >> > >>>>> _______________________________________________
>>> >> > >>>>> NumPy-Discussion mailing list
>>> >> > >>>>> NumPy-Discussion at scipy.org
>>> >> > >>>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> >> > >>>>
>>> >> > >>>>
>>> >> > >>>>
>>> >> > >>>> _______________________________________________
>>> >> > >>>> NumPy-Discussion mailing list
>>> >> > >>>> NumPy-Discussion at scipy.org
>>> >> > >>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> >> > >>>>
>>> >> > >>>
>>> >> > >>>
>>> >> > >>> _______________________________________________
>>> >> > >>> NumPy-Discussion mailing list
>>> >> > >>> NumPy-Discussion at scipy.org
>>> >> > >>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> >> > >>>
>>> >> > _______________________________________________
>>> >> > NumPy-Discussion mailing list
>>> >> > NumPy-Discussion at scipy.org
>>> >> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> >> >
>>> >>
>>> >>
>>> >> _______________________________________________
>>> >> NumPy-Discussion mailing list
>>> >> NumPy-Discussion at scipy.org
>>> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> >
>>> >
>>> > _______________________________________________
>>> > NumPy-Discussion mailing list
>>> > NumPy-Discussion at scipy.org
>>> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> _______________________________________________
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
More information about the NumPy-Discussion
mailing list