ENH: add matvec and vecmat gufuncs by mhvk · Pull Request #25675 · numpy/numpy

mhvk · 2024-01-23T21:57:19Z

This PR adds new matvec and vecmat gufunc to calculate the matrix-vector and vector-matrix product, to add to plain matrix multiplication with matmul and the inner vector product with vecdot.

Fixes #12348

Note that for complex numbers, vecmat is defined as x†A, i.e., the complex conjugate of the vector is taken. This seems to be the standard and is what we used for vecdot too (x†x). However, it is not what matmul does for vector-matrix or indeed vector-vector products. I think this is a bug in matmul, which I'm happy to fix here. But I'll post to the mailing list to get feedback.

Separately, with these functions available, in principle these could be used in __matmul__ and the specializations in np.matmul removed. But that can be a separate PR (if it is wanted at all).

Summary of mailing list discussion

Behaviour of np.matmul for vector-matrix should not be adjusted. Some surprise that "matrix multiplcation" deals with vectors at all.
Less obvious consensus on what @ should do.
Moderate support for having special matvec and vecmat functions.

mhvk · 2024-01-30T14:57:53Z

Refactored to avoid duplication of code: matvec now uses the matmul_noblas loops internally (as does the matrix-vector implementation in matmul itself), and vecmat creates a loop that uses the loops in vecdot, letting their outer loop deal with the matrix dimension.

mhvk · 2024-01-30T15:01:12Z

@charris - since you were the one originally asking for it, what do you think of this implementation? It is now quite lightweight, since it reuses inner loops.

mhvk · 2024-01-30T15:50:23Z

macos failure unrelated

seberg · 2024-01-31T14:41:04Z

Moving milestone to 2.1, if we get this in we should change it again, but at least I don't think it needs to be prioritized.

mhvk · 2024-01-31T16:41:12Z

Yes, not critical for 2.0 for sure, but maybe more logical to clear the milestone rather than set it to 2.1 if there still is a chance. This PR has become fairly simple after my refactor.

ev-br · 2024-08-03T14:17:52Z

Out of curiousity, what is the status of this PR?

mhvk · 2024-08-03T16:26:02Z

I think it would still be nice to get this in! I just rebased to make sure things still pass. @charris - as the person who originally suggested to have matvec and vecmat, what do you think?

mhvk · 2024-08-03T18:30:27Z

The circle-ci failure does not seem related:

/home/circleci/repo/doc/source/reference/arrays.promotion.rst:150: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:154: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:189: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:191: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:210: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:212: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:31: ERROR: Too many autonumbered footnote references: only 0 corresponding footnotes available.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:31: ERROR: Unknown target name: "hist-reasons".
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:68: ERROR: Unknown target name: "nep50".
source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.
/home/circleci/repo/doc/source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.
source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.

EDIT: also seen in #27091 so indeed unrelated.

seberg

FWIW, I am still happy with the addition. A couple of doc comments that would be nice to address. I think the main point is actually that it would be good to expand matmul docs to explain that these exist and fill a gap for stacks of vectors.

seberg · 2024-08-04T10:10:30Z

+           [ 0.,  0.,  1.],
+           [ 6.,  0.,  8.]])
+
+    .. versionadded:: 2.0.0


Needs to be updated, and I think we should move it to the end of the intro.

OK, updated. I agree the location is not very logical, but kept there anyway, since matmul and vecdot have it there too.

seberg · 2024-08-04T10:13:01Z

+    Matrix-vector dot product of two arrays.
+
+    Let :math:`\\mathbf{A}` be a maxtrix in ``x1`` and :math:`\\mathbf{v}` be
+    a vector in ``x2``. The matrix-vector product is defined as:


It's good, but if we have an idea, I don't think "a matrix in x1" is very clear about referring to stacks of matrices/vectors.

In fact it would be nice to:

Add matvec and vecmat to matmul "See Also"

Probably also to its notes on stacks (IIRC we have those).

OK, I tried to reword (will push up in a bit).

seberg · 2024-08-04T10:16:12Z

+    Rotate a set of vectors from Y to X along Z.
+
+    >>> a = np.array([[0., 1., 0.], [-1., 0., 0.], [0., 0., 1.]])
+    >>> v = np.array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [0., 6., 8.]])


Maybe we can name these Y_vectors, name a rotation_matrix, and format all as with the rows below one another (like the result).
Could also call it X_vectors = ... and then print that repr to document it.

I reformatted the arrays, but kept the names (since A and v are used in the description).

seberg · 2024-08-04T15:53:27Z

One thing that occurred to me, it might also make sense to put these into np.linalg. OTOH, since matmul and vecdot are in the main namespace, it is probably clearer there.

mhvk · 2024-08-04T16:40:29Z

One thing that occurred to me, it might also make sense to put these into np.linalg.

I think the main namespace is good enough, for now at least. It might change if the array API adopts them.

mhvk · 2024-09-10T15:51:21Z

@seberg - shall we get this in? I fixed the error I had made in the change-log entry so now also circle-ci passes.

seberg · 2024-09-10T16:14:12Z

I am happy to get it in, someone should maybe confirm they like it as well. Am travelling right now, so may just punt for another while before doing another pass and merging, sorry.

mattip · 2024-09-10T17:17:46Z

Did you compare performance of calls to matmul and vecmat/matvec? I would be very surprised if they were different.

mhvk · 2024-09-10T19:37:03Z

There should be no performance difference - it uses the same code. The advantage of matvec and vecmat is that one can pass in stacks of vectors and stacks of matrices instead of being stuck to a single vector.

ping also @charris who originally suggested we have these functions...

mattip · 2024-11-13T18:58:35Z

We discussed this at a recent triage meeting and would like to put it in. It needs a final review.

seberg

I looked through, and the only thing that stood out a bit is that vecmat and matvec may have subtly different behavior when it comes to accuracy for float32.

Overall, still happy to put this in.

seberg · 2024-11-21T10:32:22Z

+    npy_intp vecdot_steps[5] = {0, is2_m, os_m, is1_n, is2_n};
+    for (npy_intp i = 0; i < n_outer; i++, args[0] += s0, args[1] += s1, args[2] += s2) {
+        char *vecdot_args[3] = {args[0], args[1], args[2]};
+        @TYPE@_vecdot(vecdot_args, vecdot_dims, vecdot_steps, NULL);


OK, I have two comments, that I think I can glance over if we want this in the 2.2 release, but I also think should get an issue for follow up:

The vecmat code uses our vector multiplication function. This means, it ensures a float64 accumulation (at least partially), and correctly doesn't need to fall back to our implementation for huge arrays. (but I think/hope that 64bit blas is starting to be common).

We probably should be able to remove the outer loop here using C transpose parameter to some blas functions.

mhvk · 2024-11-21T17:28:00Z

@seberg - I hadn't thought about implementation differences for float32 at all...

I'm fairly sure I initially had vecmat quite similar to matvec here, for the blas call just ensuring the conjugate, but with a separate no-blas implementation that ensured conjugation of the vector (since one couldn't use the matmat one). Then, realizing that vecdot had more or less what was needed already, it just seemed simpler to use it, reducing the amount of extra code...

I do think it is better if the two functions behave the same way. My tendency would be to use the better precision sums for both, which would mean repeating the complex vecdot implementation without conjugation for matvec. That said, it is easier to imaging vector multiplication where the precision matters than matrix-vector, since the matrices would generally become quite huge.

I'll look a bit at whether one can do this fairly easily, but let me know if you have suggestions...

mhvk · 2024-11-21T22:22:49Z

@seberg - after spending most of the day on this (grr... still do not quite get the blas dimension stuff - too much switching between fortran and C...), I ended up with a bit of a frankencode: now the blas matrix-vector gemv will be called for float arguments, but not for complex. In principle, it works fine for complex for matvec, but it cannot handle the conjugation of the vector in vecmat; for symmetry, I thought it might be best to just fall through to my default implementation in that case -- which uses the various *_dot, so can still be blas-optimized in some cases.

What do you think?

seberg · 2024-11-22T11:45:58Z

Sorry, I didn't mean to push you into a rabbit hole :/, just realized that there likely is a difference due to the discussion about sum and float32 and was wondering if we should worry.
I think you are right, gemm could be used to achieve this probably, but gemv doesn't seem to have the transposes we would need.

I suspect your frankencode is fine (and the original would probably also have been fine). @mattip should we just go ahead with this?

Pushed a commit to fix the __module__ thing that failed CI in main. Points at one other thing: I would be totally fine with only adding these to functions to numpy.linalg (no strong opinion, though).

EDIT: Of course feel free to force push my commits away.

mhvk · 2024-11-22T14:47:50Z

Thanks for pushing the fixes - I did squash them together so it is one commit again.

Given that vecdot and matmul are in the main namespace, I think these might as well too...

seberg · 2024-11-22T14:55:16Z

Yeah, I think we can just put this in for 2.2. @charris wanted to branch this weekend. I pinged the mailing list; but if anyone shouts there (or here), we can always do a backport (it'll be a few weeks before the actual release anyway).

mhvk · 2024-11-22T15:11:08Z

With a fresh look, I now can get gemm to work. Also noticed a slight omission in the blasable tests. So, wait for merging just a little!

Internally, they mostly just call the relevant blas, or vecdot routines.

mhvk · 2024-11-22T15:56:15Z

OK, now did the gemm for vecmat complex - should be ready!

charris · 2024-11-25T16:50:20Z

Thanks Marten.

mhvk added 01 - Enhancement component: numpy._core labels Jan 23, 2024

mhvk added this to the 2.0.0 release milestone Jan 23, 2024

mhvk force-pushed the matvec-vecmat-ufuncs branch 2 times, most recently from 3647a3c to a7ac211 Compare January 23, 2024 22:11

mhvk mentioned this pull request Jan 25, 2024

ENH: expose matmat, matvec, vecmat ufuncs #12348

Closed

mhvk force-pushed the matvec-vecmat-ufuncs branch from a7ac211 to 5e5c51a Compare January 30, 2024 14:54

mhvk force-pushed the matvec-vecmat-ufuncs branch from 5e5c51a to 5b5d64f Compare January 30, 2024 15:00

mhvk force-pushed the matvec-vecmat-ufuncs branch from 07487e7 to 09d4504 Compare January 30, 2024 15:12

seberg modified the milestones: 2.0.0 release, 2.1.0 release Jan 31, 2024

mhvk removed this from the 2.1.0 release milestone Jan 31, 2024

mhvk force-pushed the matvec-vecmat-ufuncs branch from 09d4504 to c373c4a Compare August 3, 2024 16:24

mhvk force-pushed the matvec-vecmat-ufuncs branch 2 times, most recently from 1b6ca20 to a1e66f3 Compare August 3, 2024 17:36

seberg approved these changes Aug 4, 2024

View reviewed changes

mhvk force-pushed the matvec-vecmat-ufuncs branch from a1e66f3 to 16023ac Compare August 4, 2024 15:39

mhvk force-pushed the matvec-vecmat-ufuncs branch from 16023ac to 27ae492 Compare September 2, 2024 14:44

mhvk force-pushed the matvec-vecmat-ufuncs branch from 27ae492 to b09440e Compare September 10, 2024 13:02

mhvk added this to the 2.2.0 release milestone Sep 30, 2024

seberg approved these changes Nov 21, 2024

View reviewed changes

mhvk force-pushed the matvec-vecmat-ufuncs branch from b09440e to 89ffb69 Compare November 21, 2024 22:19

mhvk force-pushed the matvec-vecmat-ufuncs branch from 3825a86 to 981dd6e Compare November 22, 2024 14:10

ENH: define matvec and vecmat gufuncs

8cec646

Internally, they mostly just call the relevant blas, or vecdot routines.

mhvk force-pushed the matvec-vecmat-ufuncs branch from 981dd6e to 8cec646 Compare November 22, 2024 15:55

charris merged commit cd84af2 into numpy:main Nov 25, 2024

charris mentioned this pull request Nov 25, 2024

ENH: add matvec and vecmat gufuncs #27846

Merged

mhvk deleted the matvec-vecmat-ufuncs branch November 25, 2024 23:05

neutrinoceros mentioned this pull request Nov 27, 2024

BUG: add missing ufunc helpers for np.matvec and np.vecmat astropy/astropy#17464

Merged

1 task

rgommers mentioned this pull request Dec 26, 2024

RFC: consider adding the specialized methods for missing matrix products data-apis/array-api#877

Closed

Conversation

mhvk commented Jan 23, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Jan 30, 2024

Uh oh!

mhvk commented Jan 30, 2024

Uh oh!

mhvk commented Jan 30, 2024

Uh oh!

seberg commented Jan 31, 2024

Uh oh!

mhvk commented Jan 31, 2024

Uh oh!

ev-br commented Aug 3, 2024

Uh oh!

mhvk commented Aug 3, 2024

Uh oh!

mhvk commented Aug 3, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

seberg left a comment

Choose a reason for hiding this comment

Uh oh!

seberg Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

mhvk Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

Uh oh!

seberg Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

mhvk Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

seberg Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

mhvk Aug 4, 2024

Choose a reason for hiding this comment

Uh oh!

seberg commented Aug 4, 2024

Uh oh!

mhvk commented Aug 4, 2024

Uh oh!

mhvk commented Sep 10, 2024

Uh oh!

seberg commented Sep 10, 2024

Uh oh!

mattip commented Sep 10, 2024

Uh oh!

mhvk commented Sep 10, 2024

Uh oh!

mattip commented Nov 13, 2024

Uh oh!

seberg left a comment

Choose a reason for hiding this comment

Uh oh!

seberg Nov 21, 2024

Choose a reason for hiding this comment

Uh oh!

mhvk commented Nov 21, 2024

Uh oh!

mhvk commented Nov 21, 2024

Uh oh!

seberg commented Nov 22, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented Nov 22, 2024

Uh oh!

seberg commented Nov 22, 2024

Uh oh!

mhvk commented Nov 22, 2024

Uh oh!

mhvk commented Nov 22, 2024

Uh oh!

charris commented Nov 25, 2024

Uh oh!

mhvk commented Jan 23, 2024 •

edited

Loading

mhvk commented Aug 3, 2024 •

edited

Loading

seberg commented Nov 22, 2024 •

edited

Loading