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