On Tue, 23 Jan 2024 at 23:13, Marten van Kerkwijk <mhvk@astro.utoronto.ca> wrote:

...

...

I also 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 (remember that those are possible only if the vector is one-dimensional, i.e., not with a stack of vectors). I think this is a bug in matmul, which I'm happy to fix. But I'm posting here in part to get feedback on that.

Does matmul not mean "matrix multiplication"?

There is no conjugation in matrix multiplication.

Have I misunderstood what matmul is supposed to be? (If it does mean anything other than matrix multiplication then it is an extremely poor choice of name.)

Right now, there is this weird special-casing of one-1 arguments which are treated as vectors that are promoted to matrices [1]. I agree that it is illogical to have, which is why I had the postscript about removing it from np.matmul, and instead letting `@` call vecdot, vecmat, matvec, or matmul as appropriate.

-- Marten

[1] From the docstring of np.matmul """ ... The behavior depends on the arguments in the following way.

- If both arguments are 2-D they are multiplied like conventional matrices. - If either argument is N-D, N > 2, it is treated as a stack of matrices residing in the last two indexes and broadcast accordingly. - If the first argument is 1-D, it is promoted to a matrix by prepending a 1 to its dimensions. After matrix multiplication the prepended 1 is removed. - If the second argument is 1-D, it is promoted to a matrix by appending a 1 to its dimensions. After matrix multiplication the appended 1 is removed. """

Note that the description is consistent with the current behaviour, but I don't think it is correct (certainly for two 1-d arguments) and my suggestion corresponds to changing the third item to,

- If the first argument is 1-D, it is promoted to a matrix by prepending a 1 to its dimensions and taking the conjugate if complex. After matrix multiplication the prepended 1 is removed.

I can understand the desire to generalise the idea of matrix multiplication for when the arrays are not both 2-D but taking the complex conjugate makes absolutely no sense in the context of matrix multiplication. You note above that "vecmat is defined as x†A" but my interpretation of that is that vecmat(x, A) == matmul(conjugate(transpose(x)), A). If you want to define vecmat like that then maybe that makes sense in some contexts but including the conjugate as an implicit part of matmul is something that I would find very confusing: such a function should not be called matmul. -- Oscar

FWIW, +1 for matvec & vecmat to complement matmat (erm, matmul). Having a binop where one argument is a matrix and the other is a stack/batch of vectors is indeed awkward otherwise, and a dedicated function to clearly distinguish "two matrices" from "a matrix and a batch of vectors" sounds great from a user perspective. As to vecmat doing hermitian conjugation of the vector argument --- I'd be in favor (because <\psi | \hat H | \psi>) but this is a weak preference. ср, 24 янв. 2024 г., 21:28 Marten van Kerkwijk <mhvk@astro.utoronto.ca>:

...

Why do these belong in NumPy? What is the broad field of application of these functions? And, does a more general concept underpin them?

Multiplication of a matrix with a vector is about as common as matrix with matrix or vector with vector, and not currently easy to do for stacks of vectors, so I think the case for matvec is similarly strong as that for matmul and vecdot.

Arguably, vecmat is slightly less common, though completes the quad.

-- Marten

_______________________________________________ NumPy-Discussion mailing list -- numpy-discussion@python.org To unsubscribe send an email to numpy-discussion-leave@python.org https://mail.python.org/mailman3/lists/numpy-discussion.python.org/ Member address: evgeny.burovskiy@gmail.com

On Wed, 24 Jan 2024 at 19:29, Marten van Kerkwijk <mhvk@astro.utoronto.ca> wrote:

...

Why do these belong in NumPy? What is the broad field of application of these functions? And, does a more general concept underpin them?

Multiplication of a matrix with a vector is about as common as matrix with matrix or vector with vector, and not currently easy to do for stacks of vectors, so I think the case for matvec is similarly strong as that for matmul and vecdot.

If either argument is conjugated then I would not describe this simply as "multiplication". A different terminology would be better (and perhaps different function names). -- Oscar

Could you please offer some code or math notation to help communicate this? I am forced to guess at the need.

The words "matrix" and "vector" are ambiguous. After all, matrices (of given shape) are a type of vector (i.e., can be added and scaled.) So if by "matrix" you mean "2d array" and by "stack of vectors" you effectively mean "2d array", this sounds like a use for np.dot (possibly after a transpose). However I am going to guess that here by "vector" you actually mean a matrix (i.e., a 2d array) with only one row or only one column, so a "stack" of them is actually 3d. Perhaps the needless dimension is then the real problem and can either not be produced or can be squeezed away..

Stack of matrices in this context is a an ndarray in which the last two dimensions are interpreted as columns and rows of matrices (in that order), stack of vectors as an ndarray in which the last dimension is interpreted as column vectors. (Well, axes in all these functions can be chosen arbitrarily, but those are the defaults.) So a simple example for matvec would be a rotation matrix that I'd like to apply to a large number of vectors (say to points in space); this is currently not easy. Complex vectors might be Fourier components. (Note that I was rather sloppy in calling it multiplication rather than using the term vector-matrix product, etc.; it is definitely not element-wise!). The vector-matrix product comes up a bit less, but as mentioned by Evgeni in physics there is, e.g., the bra-ket stuff with Hamiltonians (<\psi | \hat H | \psi>), and in linear least squares one often gets terms which for real data are written as Xᵀ Z, but for complex would be Xᴴ Z [1] Hope this clarifies things! -- Marten [1] https://en.wikipedia.org/wiki/Linear_least_squares

Hi Alan, The problem with .dot is not that it is not possible, but more that it is not obvious exactly what will happen given the overloading of multiple use cases; indeed, this is why `np.matmul` was created. For the stacks of vectors case in particular, it is surprising that the vector dimension is the one but last rather than the last dimension (at least if the total dimension is more than 2). The idea here is to have a function that does just one obvious thing (sadly, for matmul we didn't quite stick to that...). All the best, Marten

I tend to agree with not using the complex conjugate for vecmat, but would prefer having separate functions for that that make it explicit in the name. I also note that mathematicians use sesquilinear forms, which have the vector conjugate on the other side, so there are different conventions. I prefer the Dirac convention myself, but many mathematical methods texts use the opposite. It is tricky for the teacher in introductory courses, right up there with vectors being called contravariant when they are actually covariant (the coefficients are contravariant). Anyway, I think having the convention explicit in the name will avoid confusion.

Hi Chuck, Indeed, which argument to conjugate seems not generally agreed on -- definitely found wikipedia pages using both convention! But I think numpy made the choice with vdot and vecdot, so that seems mostly a matter of documentation. I'm slightly confused whether or not you agree that the complex conjugate is useful, but can see your point of a name that makes it clear. Or perhaps one could have a wrapper function that let's one switch between conjugating the vector or not? One thing perhaps worth mentioning is that if we decide that vecmat does conjugate, the non-conjugating version can still trivially be done by a transpose using matvec(matrix.mT, vector) [1], while if vecmat does not conjugate by default, getting the conjugate version using the "obvious" vecmat(vector.conj(), matrix) implies a potentially costly copy of the whole vector array. (This can be avoided using the much less obvious np.vecdot(vec[..., np.newaxis, :], mat.mT), but I'm not sure that one really wants to suggest that...) More generally, in my own experience working with complex matrices and vectors, just about anytime one needs a transpose, a conjugate is generally needed too. All the best, Marten [1] Indeed, the implementation of vecmat uses the matvec loops with the matrix transposed for all but the complex numbers.