[Numpy-discussion] Scipy dot
George Nurser
gnurser at gmail.com
Fri Nov 9 04:32:53 EST 2012
Hi,
It's really good to see this work being done.
It would be great if this could somehow be put also into the np.einsum
function, which currently doesn't even use blas, and is consequently
substantially slower than current np.dot (timings on
http://mail.scipy.org/pipermail/numpy-discussion/2012-October/064259.html ).
I really do think that the Einstein summation convention is a really clear
way of doing array operations, and its use should be encouraged.
--George
On 9 November 2012 06:18, 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?
>
> 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))
>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20121109/99ca829c/attachment.html>
More information about the NumPy-Discussion
mailing list