ENH: make matmul into a ufunc by mattip · Pull Request #12219 · numpy/numpy

mattip · 2018-10-19T12:49:06Z

Replaces #11133.
Uses the extended generalized ufunc signature (n?,k),(k,m?)->(n?,m?), a custom PyUFunc_MatmulTypeResolver function, and new inner loop functions to implement matmul as a true ufunc.

Still known todo:

PR MAINT: refactor ufunc iter operand flags handling #11580 is needed to pass the failing test. The TypeResolver checks overlap. If the output and input overlap (via an out kwarg), the ufunc->op_flag will be set to use writeback semantics. On master, this currently does not work since all the flags are overridden.
The docstring processing needs refactoring, but I wanted to change it in two commits: the first to modify it for the new behavior, then to move it from core/_add_newdocs.py to ufunc_docstrings
~~The TypeResolver should already choose the proper vector-vector, vector-matrix, or matrix-matrix inner loop based on shape (maybe, needs benchmarking to see which is faster)~~ Leave this for a future PR
Make sure the code handles non-contiguous inputs and output
Make sure the conversions from npy_intp to blas int do not overflow
General performance improvements to make sure the new code is no slower than the old code, which went through einsum.

Edit: update after commit 207fd8b3

charris · 2018-11-06T21:19:41Z

It would be good to be able to deal with stacked vectors as well as stacked matrices. I started with four functions: vecvec, matvec, vecmat, matmat. Now that doesn't give matmul as a ufunc, but it does allow defining it using all ufuncs, which should be just a good. One could also just implement matmat, and do everything else by adding, removing a dimension, so that would also go through a ufunc. In all those versions, matmul should not be covered by __array_function__, it is just to supply a function for the @ operator for Python < 3.5.

@mhvk, @shoyer Thoughts?

shoyer · 2018-11-06T21:43:18Z

If I recall correctly, the agreement on the mailing list was to make np.matmul a real ufunc (with optional dimensions), and use it for implementing @. This will let us complete the table of operators and ufuncs from NEP 13.

We should think about how single gufunc interacts with axis/axes arguments -- do the selected axes determine which function gets dispatched? If so, np.matmul(x, y, axes=[(-1,), (-1,)]) might be enough for specifying vector-vector multiplication.

For the sake of clarity it still might be better to have alternative gufuncs like vecvec and matvec. But this should not change the np.matmul gufunc signature with (m?,n),(n,p?)->(m?,p?), which we really want so we can implement __matmul__/@ consistently with other arithmetic operators.

mattip · 2018-11-07T00:25:56Z

Here is the asv benchmark comparison run from commit 3703c095 to 97df928 (master just before this commit)

Note that one of the benchmarks does not run. This is due to the following code segfaulting. Now debugging this (only fails for matmul ufunc, add for instance passes

a = np.arange(240000.0).reshape(400, 600);
b = np.arange(600); 
np.matmul(a, b)

mhvk · 2018-11-07T01:20:09Z

@shoyer - that is an intriguing suggestion, to use axes to indicate that one wants the 1-D version (as it, that doesn't work, but it might be possible to make it work...).

Do agree that having matmul override is important; we put the whole ? in the signature for that to become possible...

mattip · 2018-11-07T02:58:54Z

valgrind says it is cblas_dgemv (dscal_k_HASWELL) writing past the data in libopenblasp-r0.2.20.so. Maybe the loop code is wrong? This is in matvec for double dtype.

mhvk

Tried finding your bug, and found one mistake, but I think it is unrelated. Could it be that you should use is2_n? (see in-line comment) Or is perhaps the assumption of contiguous data broken?

mattip · 2018-11-07T15:55:18Z

I found the segfault, but in doing so discovered that my cblas routines do not work with non-contiguous data.

mattip · 2018-11-07T19:53:31Z

All the benchmarks run to completion, here are the latest results. There are a few regressions, and a few slight improvements

mhvk · 2018-11-07T20:34:02Z

Only the first seems worth checking in a bit more detail. That's a large regression for a fairly common case (A * AT)

mattip · 2018-11-07T21:34:57Z

Commit cafcfa4dc improved the A @ A.T benchmark, ~~but einsum is still faster~~, there is no significant difference.

Edit: misinterpreted the results

charris · 2018-11-07T22:41:32Z

If I recall correctly, the agreement on the mailing list was to make np.matmul a real ufunc (with optional dimensions), and use it for implementing @.

One of the problems is that there are a number of optimizations that can be made that do not fit nicely into the ufunc scheme, ufuncs are too low level. For instance, a stack_of_matrices X single_matrix can be optimized if the stacked matrices are reshaped so that it is simply matrix X matrix, same with stack of vectors. My own preference would be to implement a single matmat ufunc and implement @ using it by jiggling the dimensions.

mattip · 2018-11-08T17:48:45Z

a stack_of_matrices X single_matrix can be optimized if the stacked matrices are reshaped so that it is simply matrix X matrix

Isn't that essentially the difference between dot and matmul?

mattip · 2018-11-08T17:52:54Z

@charris I opened #12348 to discuss further enhancements. Or are you suggesting to close this PR?

charris · 2018-11-08T18:53:16Z

@mattip I still need to make a careful review, before writing out my thoughts.

mattip · 2018-11-09T20:05:17Z

I think matmul as a ufunc is ready for review.

The loops cast M, N, and P from npy_intp (8 bytes) to int (4 bytes) in the cblas matrixmatrix routines, but did so in the past as well, and dot also suffers from this casting, so I am leaving that for a future issue/PR.

Speed improvements could be made, but the code seems to be as fast as the versions on master, so this too should be dealt with in subsequent work.

jakirkham · 2018-11-10T05:00:02Z

Sign in to view

Should there be some sort of error handling here or should these be dropped? Just asking as the comment says this shouldn't be run yet the branch is here.

Also same question for the branch below.

I'd use an assert here too, and also move the default definitions of trans1 and lda here (from lines 111 and 118).

assert in C will abort in debug and be ignored in release builds. There is no error handling in inner loops. I added a check before this call to ensure the input is contiguous.

Also moved the definitions.

eric-wieser · 2018-11-10T05:23:28Z

Is this worth it vs just letting python emit the error?

I tend to agree to just let python error.

The error for np.arange(10) @ 3 then becomes "ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)"

Will the typical user understand what went wrong?

3 @ 3 errors with "TypeError: unsupported operand type(s) for @: 'int' and 'int'"

The TypeError would definitely be fine. For the other, it does support keeping the error; on the other hand, really, when something goes wrong, how hard is it to do np.matmul?...

np.matmul(np.arange(10), 3) will also reach this call, and will error.

In commit 98c8ac4cd221 I removed the call to PyErr_SetString, we can revert it if desired.

mhvk

@mattip - I went through with a reasonably fine-toothed comb, and found some little bits, but nothing too major. The main question I have is whether if cblas is used, it is reasonable to copy the whole input array (a possibly very large stack of matrices) rather than copy matrix by matrix. I must say my tendency would be to just not use cblas for non-contiguous data. Anyway, that perhaps can be decided after - since I think you currently just reproduce what was done before.

mhvk · 2018-11-10T18:20:34Z

dot(a,b) -> matmul(x1, x2)? Though I think it is also fine to rewrite the sentence to

it must be C-contiguous and have both the right type and the right dtype```

mhvk · 2018-11-10T18:23:01Z

maybe "but a subsequent call to matmul will..." (I'd parenthesize that, since it is somewhat of an aside).

mhvk · 2018-11-10T19:00:30Z

Why the strange indentation?

mattip · 2018-11-13T23:26:53Z

@mhvk I would like to rebase and squash this since it has conflicts. Since that will make review more difficult, should I wait for more comments?

shoyer · 2018-11-13T23:40:04Z

Consider doing a merge instead of a rebase -- github preserves history for reviewing better that way.

…

On Tue, Nov 13, 2018 at 3:26 PM Matti Picus ***@***.***> wrote: @mhvk <https://github.com/mhvk> I would like to rebase and squash this since it has conflicts. Since that will make review more difficult, should I wait for more comments? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#12219 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABKS1rpt6_N9b2lcbqYfQBs1K9tsUMAaks5uu1U_gaJpZM4XwSxP> .

eric-wieser · 2018-11-27T23:03:52Z

Are you sure these ns and ms aren't swapped?

the assert was failing in the debug CI test, fixing.

I'm wrong, these ms and ns are correct.

What's wrong is that the asserts on i1 should be within the if or if-else

I rewrote the macro instead

eric-wieser · 2018-11-27T23:50:00Z

Perhaps rename to blasable_matrix_vector, and vector_blasable_matrix above?

will do if there are more fixes

eric-wieser · 2018-11-30T06:24:31Z

This can be removed

eric-wieser · 2018-11-30T06:24:36Z

Not strictly true, when you consider the 1d overloads. Maybe fine as is anyway.

leaving this for now

eric-wieser · 2018-11-30T06:24:40Z

I think you want double backticks here, since single backticks tries to make a reference.

eric-wieser · 2018-11-30T06:24:43Z

For a separate doc cleanup PR: This would be a lot clearer if the dimensions were all different. (2, 3, 4) and (2, 4, 5) would show the behavior nicely.

eric-wieser · 2018-11-30T06:24:44Z

For a separate doc cleanup PR: spaces around commas for PEP8 - let's not train new users into going against PEP8 if we can avoid it.

easy enough to clean up now

eric-wieser · 2018-11-30T06:25:09Z

Not a fan of this macro. I think I'd prefer a form that takes all the parameters, with a little more explanation. Also, may as well just use an inline function.

Perhaps something like:

/* -1 to be conservative, in case blas internally uses a for loop with an inclusive upper bound */ #define BLAS_MAXSIZE (NPY_MAX_INT - 1) /** * Return true if matrix with strides (s1, s2) and dimensions (n1, n2) can be used with blas */ static inline npy_bool is_corder_blasable_2d(npy_intp s1, npy_intp s2, npy_intp n1, npy_intp n2, npy_intp itemsize) { npy_intp step1; if (!(s2 == itemsize)) { /* second axis must be contiguous */ return NPY_FALSE; } if (!(s1 % itemsize == 0)) { /* strides in the first axis must count whole elements */ return NPY_FALSE; } step1 = s1 / itemsize; if (!(n1 <= step1) { /* number of elements must not cause slices of the last axis to overlap */ return NPY_FALSE; } /* todo: can n1 or n2 be 0? */ if (!(n1 <= BLAS_MAXSIZE && n2 <= BLAS_MAXSIZE && step1 <= BLAS_MAXSIZE) { /* dimensions must be representable in BLAS */ return NPY_FALSE; } return NPY_TRUE; } static inline npy_bool is_forder_blasable_2d(npy_intp s1, npy_intp s2, npy_intp n1, npy_intp n2, npy_intp itemsize) { /* just transpose the array */ return is_corder_blasable_2d(s2, s1, n2, n1, itemsize); } static inline npy_bool is_blasable_1d(npy_intp s, npy_intp n, npy_intp itemsize) { npy_intp step; if (!(s1 % itemsize == 0)) { /* strides must count whole elements */ return NPY_FALSE; } /* todo: can INCX be 0 or negative? */ /* todo: can N be 0? */ if (!(n <= BLAS_MAXSIZE && step <= BLAS_MAXSIZE) { /* dimensions must fit in an int */ return NPY_FALSE; } return NPY_TRUE; }

Yes, this is far more verbose, but I'd hope they'd be optimized to the same code, and this gives far more places for comments.

It's not clear to me whether c / f order is the right way to name these function, or whether is_forder_blasable is more useful than explicitly calling is_corder_blasable with the arguments in the other order.

adopted a modified version of this

eric-wieser · 2018-11-30T06:25:12Z

Why the 2 here?

eric-wieser · 2018-11-30T06:25:15Z

could omit the == 1 after you remove the 2 above

eric-wieser · 2018-12-01T04:14:51Z

I don't think these three lines are correct, which is why I was pushing for a distinction between "is supported by blas" and "is supported by blas if I transpose it".

ok, will extend the function

eric-wieser · 2018-12-01T23:01:44Z

This is not true - it needn't be a multiple, it just needs to exceed the length of the fast axis

eric-wieser · 2018-12-01T23:02:38Z

This should move into the if above.

eric-wieser · 2018-12-01T23:50:04Z

Should this still be here?

eric-wieser · 2018-12-01T23:54:33Z

Probably safer to use BLAS_MAX_SIZE here

eric-wieser · 2018-12-03T03:13:08Z

+            if (op[i] == NULL) {
+                op[i] = operands[i];
+                Py_INCREF(op[i]);
+            }


These lines are unrelated cleanup, right?

eric-wieser

Looks ready now - just one comment left about documenting the need for the non-default ufunc flags, which you explained well in person.

mattip · 2018-12-03T18:38:05Z

Fixed comment on non-default flags, ready for merge

eric-wieser · 2018-12-04T16:27:57Z

Thanks for the persistence, @mattip!

shoyer · 2018-12-04T19:56:30Z

Woohoo, thank you @mattip!

mhvk · 2018-12-04T19:57:06Z

Yay!! Thanks, @mattip!

jakirkham · 2018-12-05T03:29:28Z

Thanks so much for working on this, @mattip! 😄 🎉

mattip added 01 - Enhancement 56 - Needs Release Note. Needs an entry in doc/release/upcoming_changes labels Oct 19, 2018

tylerjereddy added the 25 - WIP label Oct 19, 2018

mattip added this to the 1.16.0 release milestone Nov 4, 2018

mattip force-pushed the matmul-as-ufunc2 branch from df5e819 to 3703c09 Compare November 6, 2018 01:49

mhvk reviewed Nov 7, 2018

View reviewed changes

Comment thread numpy/core/src/umath/matmul.c.src Resolved

Comment thread numpy/core/src/umath/matmul.c.src Resolved

Comment thread numpy/core/src/umath/matmul.c.src Resolved

Comment thread numpy/core/src/umath/matmul.c.src Resolved

mattip mentioned this pull request Nov 7, 2018

ENH: expose matmat, matvec, vecmat ufuncs #12348

Closed

mattip force-pushed the matmul-as-ufunc2 branch from b14424e to 7c62240 Compare November 8, 2018 21:54

mattip mentioned this pull request Nov 9, 2018

Numpy passes invalid arguments to DGEMV when stride = 0 #12156

Closed

mattip changed the title ~~WIP, ENH: make matmul into a ufunc~~ ENH: make matmul into a ufunc Nov 9, 2018

jakirkham reviewed Nov 10, 2018

View reviewed changes

eric-wieser reviewed Nov 10, 2018

View reviewed changes

mhvk reviewed Nov 10, 2018

View reviewed changes

mattip mentioned this pull request Nov 11, 2018

ENH: improve handling of non-contiguous operands in matmul #12365

Closed

eric-wieser reviewed Nov 27, 2018

View reviewed changes

Comment thread numpy/core/src/umath/matmul.c.src Resolved

eric-wieser reviewed Nov 27, 2018

View reviewed changes

Comment thread numpy/core/src/umath/matmul.c.src Resolved

eric-wieser reviewed Nov 30, 2018

View reviewed changes

eric-wieser reviewed Dec 1, 2018

View reviewed changes

Comment thread numpy/core/tests/test_multiarray.py Resolved

mattip force-pushed the matmul-as-ufunc2 branch 2 times, most recently from 962342f to ae9bb7b Compare December 1, 2018 23:42

eric-wieser reviewed Dec 1, 2018

View reviewed changes

Comment thread numpy/core/src/umath/matmul.c.src Outdated

eric-wieser Dec 1, 2018

Copy link
Copy Markdown

Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this still be here?

eric-wieser reviewed Dec 1, 2018

View reviewed changes

Comment thread numpy/core/src/umath/matmul.c.src Outdated

eric-wieser Dec 1, 2018

Copy link
Copy Markdown

Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably safer to use BLAS_MAX_SIZE here

mattip force-pushed the matmul-as-ufunc2 branch from ae9bb7b to 2bcc6e3 Compare December 1, 2018 23:58

ENH: corrections and fixes from review

72f2abf

mattip force-pushed the matmul-as-ufunc2 branch from 2bcc6e3 to 72f2abf Compare December 2, 2018 00:14

eric-wieser reviewed Dec 3, 2018

View reviewed changes

Comment thread numpy/core/src/umath/ufunc_object.c Resolved

eric-wieser approved these changes Dec 3, 2018

View reviewed changes

DOC: add a comment detailing the output ufunc flag choices

18476f5

eric-wieser merged commit 82428bb into numpy:master Dec 4, 2018

shoyer mentioned this pull request Dec 4, 2018

Update NDArrayOperatorsMixin to implement __matmul__ using np.matmul #12482

Closed

mattip deleted the matmul-as-ufunc2 branch March 4, 2019 13:19

ahaldane mentioned this pull request Mar 29, 2019

var, std memory consumption #13199

Open

mattip mentioned this pull request Aug 20, 2019

np.matmul(a.T, b) 1000x slower depending on data order, .copy() "fixes" it #14306

Closed

Conversation

mattip commented Oct 19, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

charris commented Nov 6, 2018

Uh oh!

shoyer commented Nov 6, 2018

Uh oh!

mattip commented Nov 7, 2018

Uh oh!

mhvk commented Nov 7, 2018

Uh oh!

mattip commented Nov 7, 2018

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

mattip commented Nov 7, 2018

Uh oh!

mattip commented Nov 7, 2018

Uh oh!

mhvk commented Nov 7, 2018

Uh oh!

mattip commented Nov 7, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

charris commented Nov 7, 2018

Uh oh!

mattip commented Nov 8, 2018

Uh oh!

mattip commented Nov 8, 2018

Uh oh!

charris commented Nov 8, 2018

Uh oh!

mattip commented Nov 9, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

This comment was marked as outdated.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

mattip commented Oct 19, 2018 •

edited

Loading

mattip commented Nov 7, 2018 •

edited

Loading