WIP: ENH: Matmul as ufunc by mattip · Pull Request #11133 · numpy/numpy

mattip · 2018-05-22T00:28:04Z

The discussion of issue #9028 yielded two different possible directions for allowing ndarray subclasses to override __matmul__ via __array_ufunc__. This is option 1 - allow a single ufunc matmul which will implement vector-vector, vector-matrix, matrix-vector, and matrix-matrix multiplication. In this PR, which builds on ~~PR #11132~~ PR #11175, the implementation of matmul is converted to a ufunc.

Note the pull request incorporates the changes from pr #11132 so the changesets are repeated

What changed?

matmul is now implemented as a true ufunc, it appears in the umath's n_ops, can be set and gotten using PyArray_{G,S}etNumericOps, and its loops are defined in loops.c.src
tests were expanded to reflect this change

What is not finished?

~~since I simply copied the partial loop implementation from _umath_tests, matmul cannot accept complex or half-float dtypes~~
~~I wrote the code together with PR WIP: ENH: Expand gufunc signature to allow flexible dimension specs #11132 I probably have chages there that belong here and visa-versa~~
reuse the logic in dot to make this faster

What is next?

Get feedback on the relative merits of this approach vs. the one taken in PR WIP: ENH: allow __array_ufunc__ to override __matmul__ #11061, where I produced a demi-ufunc with a wrapper function that can be used to convert any function into a ufunc. The choice has been made, this approach wins.

Edit: update status

mattip · 2018-06-05T00:26:12Z

Rebased off PR #11175. We have loops in cblas_matrixproduct, and in vdot.c, there must be a way to reuse them nicely without having to rewrite all this

mhvk · 2018-06-06T01:13:24Z

Setup looks great; the trick will indeed be to have a nice way to choose BLAS etc for the inner loop.

mattip · 2018-06-06T20:56:54Z

Since matmul is such a core function, I ran a little comparision multiplying two 3x3 matrices via dot, einsum and matmul-as-a-ufunc. It turns out (no suprise) there is a cost to ufuncs. For this simple case, and without diving too deeply into cblas's gemm vs. a naive looping implementation, the ufunc case is roughly 2 times slower. Of course we are talking about very small numbers, around 1.2 microseconds vs. half of that, so I am going to push forward with making large arrays faster by reusing the dot loops here.

mhvk · 2018-06-06T21:11:28Z

@mattip - I think the gufunc loops need some special-case speed-ups that the ufuncs already have, in particular for the case that there is no outer looping to be done, in which case the elementary function should just be called directly. But that is independent of matmul...

One thing to possibly consider for matmul - the loop selection could include whether or not to use the BLAS routines, or that choice could be made inside the loop function. The latter may be easier to programme, but the former might be better.

mattip · 2018-06-07T00:33:57Z

One of the trickier parts of this is what to do for non-contiguous memory. Just changing the benchmark code from func(a, a) to func(b[::2], [b::2]) where a is 3x3 and b is 6x6 the times for dot and einsum doubled but the time for "dumb" matmul did not change. Increasing the size to 50x50, still with the [::2] indexing, now dot is 3 times faster than matmul, and dot is over 4 times faster than einsum

mhvk · 2018-06-07T13:36:12Z

Indeed, I think that when the strides are not contiguous, one probably should copy data over, but there is a question for what number of points this becomes useful. My guess would be that for small arrays the copy is so fast that it might not matter.

Aside: I've been implementing this in astropy - this includes a very basic 3x3 matrix multiplication routine, and in my implementation I check in the inner loop whether copies are necessary [1]. As I wrote above, I think it would be more elegant to do this outside, and select a copying or non-copying loop -- and pass on some temporary storage to the inner loop if needed -- though it may not matter that much speed-wise. Also, I think almost anything will be faster than the current call to einsum...

[1] https://gist.github.com/mhvk/a3c17f0db77b91151d50a0d278621d97#file-ufunc-c-L363

mhvk · 2018-06-08T23:41:53Z

I like the last commit - but isn't that independent of matmul? I.e., should it be build on top of #11175? It would go with the version change...

mattip · 2018-07-03T23:48:42Z

it seems the matmul vector-vector case for column-row -> 2d array was not tested. I put a marker in each case of all the matmul loops, and that one came up missing.

mattip · 2018-07-03T23:59:03Z

note this is a wip, builds on top of master and pr #11175. Once that is merged I will clean up this branch

mattip · 2018-07-05T02:10:28Z

level 1 blas_{s,d,c,z}axy can work with non-contiguous data but require zeroing out y before the call, does it give a speedup vs. writing the operation y = a * x in naive C? I chose to go with the C loop, benchmarking needed.

mattip · 2018-07-05T15:34:50Z

Just noticed there is an optimization for a @ a.T in #6932 and a request for a@ a.H in #6951. Make sure these are not ruined by this refactoring

codecov-io · 2018-07-23T15:14:17Z

Codecov Report

Merging #11133 into master will increase coverage by 0.02%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master   #11133      +/-   ##
==========================================
+ Coverage    85.7%   85.73%   +0.02%     
==========================================
  Files         327      327              
  Lines       81995    82139     +144     
==========================================
+ Hits        70274    70418     +144     
  Misses      11721    11721

Impacted Files	Coverage Δ
numpy/core/numeric.py	`90.96% <ø> (ø)`	⬆️
numpy/core/_add_newdocs.py	`100% <ø> (ø)`	⬆️
numpy/core/setup.py	`0% <ø> (ø)`	⬆️
numpy/core/tests/test_ufunc.py	`98.95% <100%> (+0.1%)`	⬆️
numpy/core/tests/test_multiarray.py	`98.27% <100%> (+0.01%)`	⬆️
numpy/core/tests/test_deprecations.py	`91.73% <0%> (-0.11%)`	⬇️
numpy/linalg/linalg.py	`91.07% <0%> (ø)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0124839...4254dc0. Read the comment docs.

tylerjereddy

Please rebase on top of the latest master branch so that `.codecov.yml' setting disabling auto-comments can be detected in the PR.

jakirkham · 2018-07-28T03:42:03Z

Just noticed there is an optimization for a @ a.T in #6932 and a request for a@ a.H in #6951. Make sure these are not ruined by this refactoring

Thanks for taking these into account. :)

We got a @ a.T to work before by checking whether the data was the same and only striding changed. a @ a.H was a bit hard for us to figure out before given how conjugating will generate a copy.

tylerjereddy · 2018-08-04T03:54:18Z

+ *        out[n, 2] = in1[n, 0]*in2[n, 1] - in1[n, 1]*in2[n, 0]
+ */
+static void
+@TYPE@_cross1d(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))


I may as well just ask--do we have some docs that describe the mechanics of the @ symbol usage in NumPy C source? It appears to be some kind of mechanism for generalization / substitution, but unless I missed a chapter in the C text book, not standard?

Numpy uses a custom templating system. The parser is defined in conv_template.py.

Yup, I know that now, but it is not well exposed documentation -- I've been working on #11689 & to be fair the distutils docs template section is nothing more than a placeholder.

tylerjereddy · 2018-08-04T03:58:44Z

+    else
+#endif
+    {
+        @type@ sum = (@type@)0;  /* could make this double */


The comment here means that there's no compelling advantage to having @type@ of sum be flexible, or something else?

tylerjereddy · 2018-08-05T03:32:11Z

   elementary function's implementation.
 #. White spaces are ignored.
+#. An integer as a dimension name freezes that dimension to the value,
+#. The name can be suffixed with a question mark, this make the dimension a


make -> makes I think

tylerjereddy · 2018-08-05T03:52:00Z

+                            (void*)op, dm, dn);
+            } else {
+                /* cannot reach here */
+                /* printf("inconceivable to reach this line m=%ld, n=%ld, p=%ld\n",


remove debug code and raise some kind of exception here?

We cannot error inside a loop. I could just remove it, but it is helpful for reasoning about the code and can be enabled for debugging

tylerjereddy · 2018-08-05T03:52:29Z

+        /* npy_intp ib1_m = is1_m * dm;
+        npy_intp ob_m = os_m * dp; */
+        /*
+         * TODO: refactor this out to a inner_loop_selector, instead of


Still needs to be done?

Maybe, but not in this iteration of the code

tylerjereddy · 2018-08-05T03:55:59Z

+        assert_equal(c, b)
+        a = np.arange(18).reshape(2, 3, 3)
+        c = np.matmul(a, a, out=a[::-1, ...])
+        assert(c.base is a.base)


why the switch from assert_ above to plain assert here?

tylerjereddy · 2018-08-05T04:01:56Z

+        assert_raises(ValueError, umt.cross1d, a, np.arange(4.))
+        assert_raises(ValueError, umt.cross1d, a, np.arange(3.), np.zeros((3, 4)))
+
+    def test_can_ignore_signature(self):


Any way to split this 70-line "unit test" up into smaller pieces logically?

This come from PR #11175 which must be merged before this one. @mhvk any thoughts?

Yes, my fault. Might be good to refactor, but out of scope for this PR, I think.

mattip · 2018-10-19T08:51:06Z

Closing this, will reopen as a new, cleaner PR

charris added 01 - Enhancement 23 - Wish List 15 - Discussion and removed 23 - Wish List labels May 23, 2018

mattip force-pushed the matmul-as-ufunc branch 2 times, most recently from 59ef4ea to 7ab13be Compare May 23, 2018 20:12

mattip mentioned this pull request May 29, 2018

ENH: Generalized ufunc signature expansion for frozen and flexible dimensions #11175

Merged

mattip force-pushed the matmul-as-ufunc branch from 7ab13be to e310ed8 Compare June 5, 2018 00:23

mattip added the component: numpy.ufunc label Jun 8, 2018

mattip mentioned this pull request Jun 25, 2018

ENH: implement nep 0015: merge multiarray and umath #10915

Merged

mattip force-pushed the matmul-as-ufunc branch 2 times, most recently from 7700bf9 to 7efa656 Compare June 28, 2018 19:13

mattip force-pushed the matmul-as-ufunc branch 3 times, most recently from 5cf6b60 to 2bcf710 Compare July 4, 2018 17:54

This was referenced Jul 5, 2018

np.dot deprecations #5859

Closed

ENH: allocate working buffers outside ufunc's inner loop #11510

Open

mattip force-pushed the matmul-as-ufunc branch from 2515ae1 to 48de0d4 Compare July 5, 2018 23:04

mattip force-pushed the matmul-as-ufunc branch from 5c41a57 to c80c88b Compare July 23, 2018 06:20

mattip force-pushed the matmul-as-ufunc branch from c80c88b to 089b1a2 Compare July 23, 2018 15:13

tylerjereddy reviewed Jul 23, 2018

View reviewed changes

mattip force-pushed the matmul-as-ufunc branch from d916028 to 8567f2f Compare July 26, 2018 03:25

mattip mentioned this pull request Jul 26, 2018

numpy.matmul is slow #7569

Closed

jaimefrio and others added 12 commits July 26, 2018 19:56

ENH: test, document, implement flexible signature (PR numpy#11175)

935480e

ENH: make matmul a full-fledged ufunc

a820bc8

add half and complex loops, start thinking about cblas, expand tests

3788c9c

backward compatible scalar error message

c8abc54

refactor loops to prepare for cblas integration

726a821

add many cases and use first BLAS function, add test for missing case

a022d4d

test and fix array([2]) @ array([3, 4]), decide against level 1 blas

412cb35

test, fix overlapping and non-contiguous via special type_resolver

313d898

add matrix-vector BLAS gemv calls

cf8bf84

fix non-blas compilation, cast to void* to avoid warning

37d0781

codify inability to use out=a for non-contiguous or wrong type a

9c1b5f8

add blas level 3 and remove redundant threading directives

4254dc0

mattip force-pushed the matmul-as-ufunc branch from 8567f2f to 4254dc0 Compare July 26, 2018 23:57

tylerjereddy reviewed Aug 4, 2018

View reviewed changes

tylerjereddy reviewed Aug 5, 2018

View reviewed changes

fixes from review

af51ac4

mattip mentioned this pull request Sep 8, 2018

ENH: Support object arrays in matmul #11332

Closed

This was referenced Oct 12, 2018

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

Closed

Tracking issue for implementation of NEP-18 (__array_function__) #12028

Closed

mattip closed this Oct 19, 2018

mattip mentioned this pull request Oct 19, 2018

ENH: make matmul into a ufunc #12219

Merged

6 tasks done

mattip mentioned this pull request Oct 27, 2018

MAINT: refactor ufunc iter operand flags handling #11580

Merged

mattip deleted the matmul-as-ufunc branch March 4, 2019 13:16

Conversation

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

Uh oh!

Uh oh!

mattip commented Jun 5, 2018

Uh oh!

mhvk commented Jun 6, 2018

Uh oh!

mattip commented Jun 6, 2018

Uh oh!

mhvk commented Jun 6, 2018

Uh oh!

mattip commented Jun 7, 2018

Uh oh!

mhvk commented Jun 7, 2018

Uh oh!

mhvk commented Jun 8, 2018

Uh oh!

mattip commented Jul 3, 2018

Uh oh!

mattip commented Jul 3, 2018

Uh oh!

mattip commented Jul 5, 2018

Uh oh!

mattip commented Jul 5, 2018

Uh oh!

codecov-io commented Jul 23, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

tylerjereddy left a comment

Choose a reason for hiding this comment

Uh oh!

jakirkham commented Jul 28, 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!

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!

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!

mattip commented Oct 19, 2018

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

8 participants

mattip commented May 22, 2018 •

edited

Loading

codecov-io commented Jul 23, 2018 •

edited

Loading