ENH: define matvec and vecmat gufuncs

Internally, they mostly just call the relevant blas, or vecdot routines.
numpy · charris · Nov 25, 2024 · Jan 22, 2024 · Nov 25, 2024 · 391e5def375db6f0d1098fedf410f69cb7b8f5f6
commit 391e5def375db6f0d1098fedf410f69cb7b8f5f6
diff --git a/benchmarks/benchmarks/bench_ufunc.py b/benchmarks/benchmarks/bench_ufunc.py
@@ -16,12 +16,12 @@
          'isinf', 'isnan', 'isnat', 'lcm', 'ldexp', 'left_shift', 'less',
          'less_equal', 'log', 'log10', 'log1p', 'log2', 'logaddexp',
          'logaddexp2', 'logical_and', 'logical_not', 'logical_or',
-          'logical_xor', 'matmul', 'maximum', 'minimum', 'mod', 'modf',
-          'multiply', 'negative', 'nextafter', 'not_equal', 'positive',
+          'logical_xor', 'matmul', 'matvec', 'maximum', 'minimum', 'mod',
+          'modf', 'multiply', 'negative', 'nextafter', 'not_equal', 'positive',
          'power', 'rad2deg', 'radians', 'reciprocal', 'remainder',
          'right_shift', 'rint', 'sign', 'signbit', 'sin',
          'sinh', 'spacing', 'sqrt', 'square', 'subtract', 'tan', 'tanh',
-          'true_divide', 'trunc', 'vecdot']
+          'true_divide', 'trunc', 'vecdot', 'vecmat']
 arrayfuncdisp = ['real', 'round']

 for name in ufuncs:

diff --git a/doc/release/upcoming_changes/25675.new_feature.rst b/doc/release/upcoming_changes/25675.new_feature.rst
@@ -0,0 +1,20 @@
+New functions for matrix-vector and vector-matrix products
+----------------------------------------------------------
+
+Two new generalized ufuncs were defined:
+
+* `numpy.matvec` - matrix-vector product, treating the arguments as
+  stacks of matrices and column vectors, respectively.
+
+* `numpy.vecmat` - vector-matrix product, treating the arguments as
+  stacks of column vectors and matrices, respectively. For complex
+  vectors, the conjugate is taken.
+
+These add to the existing `numpy.matmul` as well as to `numpy.vecdot`,
+which was added in numpy 2.0.
+
+Note that `numpy.matmul` never takes a complex conjugate, also not
+when its left input is a vector, while both `numpy.vecdot` and
+`numpy.vecmat` do take the conjugate for complex vectors on the
+left-hand side (which are taken to be the ones that are transposed,
+following the physics convention).
diff --git a/doc/source/reference/routines.linalg.rst b/doc/source/reference/routines.linalg.rst
@@ -62,6 +62,8 @@ Matrix and vector products
   outer
   matmul
   linalg.matmul (Array API compatible location)
+   matvec
+   vecmat
   tensordot
   linalg.tensordot (Array API compatible location)
   einsum

diff --git a/numpy/__init__.py b/numpy/__init__.py
@@ -151,10 +151,10 @@
        left_shift, less, less_equal, lexsort, linspace, little_endian, log,
        log10, log1p, log2, logaddexp, logaddexp2, logical_and, logical_not,
        logical_or, logical_xor, logspace, long, longdouble, longlong, matmul,
-        matrix_transpose, max, maximum, may_share_memory, mean, memmap, min,
-        min_scalar_type, minimum, mod, modf, moveaxis, multiply, nan, ndarray,
-        ndim, nditer, negative, nested_iters, newaxis, nextafter, nonzero,
-        not_equal, number, object_, ones, ones_like, outer, partition,
+        matvec, matrix_transpose, max, maximum, may_share_memory, mean, memmap,
+        min, min_scalar_type, minimum, mod, modf, moveaxis, multiply, nan,
+        ndarray, ndim, nditer, negative, nested_iters, newaxis, nextafter,
+        nonzero, not_equal, number, object_, ones, ones_like, outer, partition,
        permute_dims, pi, positive, pow, power, printoptions, prod,
        promote_types, ptp, put, putmask, rad2deg, radians, ravel, recarray,
        reciprocal, record, remainder, repeat, require, reshape, resize,
@@ -165,8 +165,8 @@
        str_, subtract, sum, swapaxes, take, tan, tanh, tensordot,
        timedelta64, trace, transpose, true_divide, trunc, typecodes, ubyte,
        ufunc, uint, uint16, uint32, uint64, uint8, uintc, uintp, ulong,
-        ulonglong, unsignedinteger, unstack, ushort, var, vdot, vecdot, void,
-        vstack, where, zeros, zeros_like
+        ulonglong, unsignedinteger, unstack, ushort, var, vdot, vecdot,
+        vecmat, void, vstack, where, zeros, zeros_like
    )

    # NOTE: It's still under discussion whether these aliases

diff --git a/numpy/__init__.pyi b/numpy/__init__.pyi
@@ -4490,6 +4490,7 @@ logical_not: _UFunc_Nin1_Nout1[L['logical_not'], L[20], None]
 logical_or: _UFunc_Nin2_Nout1[L['logical_or'], L[20], L[False]]
 logical_xor: _UFunc_Nin2_Nout1[L['logical_xor'], L[19], L[False]]
 matmul: _GUFunc_Nin2_Nout1[L['matmul'], L[19], None, L["(n?,k),(k,m?)->(n?,m?)"]]
+matvec: _GUFunc_Nin2_Nout1[L['matvec'], L[19], None, L["(m,n),(n)->(m)"]]
 maximum: _UFunc_Nin2_Nout1[L['maximum'], L[21], None]
 minimum: _UFunc_Nin2_Nout1[L['minimum'], L[21], None]
 mod: _UFunc_Nin2_Nout1[L['remainder'], L[16], None]
@@ -4519,6 +4520,7 @@ tanh: _UFunc_Nin1_Nout1[L['tanh'], L[8], None]
 true_divide: _UFunc_Nin2_Nout1[L['true_divide'], L[11], None]
 trunc: _UFunc_Nin1_Nout1[L['trunc'], L[7], None]
 vecdot: _GUFunc_Nin2_Nout1[L['vecdot'], L[19], None, L["(n),(n)->()"]]
+vecmat: _GUFunc_Nin2_Nout1[L['vecmat'], L[19], None, L["(n),(n,m)->(m)"]]

 abs = absolute
 acos = arccos

diff --git a/numpy/_core/code_generators/generate_umath.py b/numpy/_core/code_generators/generate_umath.py
@@ -1152,6 +1152,22 @@ def english_upper(s):
          TD(O),
          signature='(n),(n)->()',
          ),
+'matvec':
+    Ufunc(2, 1, None,
+          docstrings.get('numpy._core.umath.matvec'),
+          "PyUFunc_SimpleUniformOperationTypeResolver",
+          TD(notimes_or_obj),
+          TD(O),
+          signature='(m,n),(n)->(m)',
+          ),
+'vecmat':
+    Ufunc(2, 1, None,
+          docstrings.get('numpy._core.umath.vecmat'),
+          "PyUFunc_SimpleUniformOperationTypeResolver",
+          TD(notimes_or_obj),
+          TD(O),
+          signature='(n),(n,m)->(m)',
+          ),
 'str_len':
    Ufunc(1, 1, Zero,
          docstrings.get('numpy._core.umath.str_len'),

diff --git a/numpy/_core/code_generators/ufunc_docstrings.py b/numpy/_core/code_generators/ufunc_docstrings.py
@@ -44,7 +44,7 @@ def add_newdoc(place, name, doc):

    skip = (
        # gufuncs do not use the OUT_SCALAR replacement strings
-        'matmul', 'vecdot',
+        'matmul', 'vecdot', 'matvec', 'vecmat',
        # clip has 3 inputs, which is not handled by this
        'clip',
    )
@@ -2793,7 +2793,9 @@ def add_newdoc(place, name, doc):

    See Also
    --------
-    vdot : Complex-conjugating dot product.
+    vecdot : Complex-conjugating dot product for stacks of vectors.
+    matvec : Matrix-vector product for stacks of matrices and vectors.
+    vecmat : Vector-matrix product for stacks of vectors and matrices.
    tensordot : Sum products over arbitrary axes.
    einsum : Einstein summation convention.
    dot : alternative matrix product with different broadcasting rules.
@@ -2808,10 +2810,10 @@ def add_newdoc(place, name, doc):
      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.
+      the prepended 1 is removed. (For stacks of vectors, use ``vecmat``.)
    - 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.
+      the appended 1 is removed. (For stacks of vectors, use ``matvec``.)

    ``matmul`` differs from ``dot`` in two important ways:

@@ -2910,8 +2912,8 @@ def add_newdoc(place, name, doc):
        Input arrays, scalars not allowed.
    out : ndarray, optional
        A location into which the result is stored. If provided, it must have
-        a shape that the broadcasted shape of `x1` and `x2` with the last axis
-        removed. If not provided or None, a freshly-allocated array is used.
+        the broadcasted shape of `x1` and `x2` with the last axis removed.
+        If not provided or None, a freshly-allocated array is used.
    **kwargs
        For other keyword-only arguments, see the
        :ref:`ufunc docs <ufuncs.kwargs>`.
@@ -2933,6 +2935,9 @@ def add_newdoc(place, name, doc):
    See Also
    --------
    vdot : same but flattens arguments first
+    matmul : Matrix-matrix product.
+    vecmat : Vector-matrix product.
+    matvec : Matrix-vector product.
    einsum : Einstein summation convention.

    Examples
@@ -2949,6 +2954,135 @@ def add_newdoc(place, name, doc):
    .. versionadded:: 2.0.0
    """)

+add_newdoc('numpy._core.umath', 'matvec',
+    """
+    Matrix-vector dot product of two arrays.
+
+    Given a matrix (or stack of matrices) :math:`\\mathbf{A}` in ``x1`` and
+    a vector (or stack of vectors) :math:`\\mathbf{v}` in ``x2``, the
+    matrix-vector product is defined as:
+
+    .. math::
+       \\mathbf{A} \\cdot \\mathbf{b} = \\sum_{j=0}^{n-1} A_{ij} v_j
+
+    where the sum is over the last dimensions in ``x1`` and ``x2``
+    (unless ``axes`` is specified).  (For a matrix-vector product with the
+    vector conjugated, use ``np.vecmat(x2, x1.mT)``.)
+
+    Parameters
+    ----------
+    x1, x2 : array_like
+        Input arrays, scalars not allowed.
+    out : ndarray, optional
+        A location into which the result is stored. If provided, it must have
+        the broadcasted shape of ``x1`` and ``x2`` with the summation axis
+        removed. If not provided or None, a freshly-allocated array is used.
+    **kwargs
+        For other keyword-only arguments, see the
+        :ref:`ufunc docs <ufuncs.kwargs>`.
+
+    Returns
+    -------
+    y : ndarray
+        The matrix-vector product of the inputs.
+
+    Raises
+    ------
+    ValueError
+        If the last dimensions of ``x1`` and ``x2`` are not the same size.
+
+        If a scalar value is passed in.
+
+    See Also
+    --------
+    vecdot : Vector-vector product.
+    vecmat : Vector-matrix product.
+    matmul : Matrix-matrix product.
+    einsum : Einstein summation convention.
+
+    Examples
+    --------
+    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.]])
+    >>> np.matvec(a, v)
+    array([[ 0., -1.,  0.],
+           [ 1.,  0.,  0.],
+           [ 0.,  0.,  1.],
+           [ 6.,  0.,  8.]])
+
+    .. versionadded:: 2.1.0
+    """)
+
+add_newdoc('numpy._core.umath', 'vecmat',
+    """
+    Vector-matrix dot product of two arrays.
+
+    Given a vector (or stack of vector) :math:`\\mathbf{v}` in ``x1`` and
+    a matrix (or stack of matrices) :math:`\\mathbf{A}` in ``x2``, the
+    vector-matrix product is defined as:
+
+    .. math::
+       \\mathbf{b} \\cdot \\mathbf{A} = \\sum_{i=0}^{n-1} \\overline{v_i}A_{ij}
+
+    where the sum is over the last dimension of ``x1`` and the one-but-last
+    dimensions in ``x2`` (unless `axes` is specified) and where
+    :math:`\\overline{v_i}` denotes the complex conjugate if :math:`v`
+    is complex and the identity otherwise. (For a non-conjugated vector-matrix
+    product, use ``np.matvec(x2.mT, x1)``.)
+
+    Parameters
+    ----------
+    x1, x2 : array_like
+        Input arrays, scalars not allowed.
+    out : ndarray, optional
+        A location into which the result is stored. If provided, it must have
+        the broadcasted shape of ``x1`` and ``x2`` with the summation axis
+        removed. If not provided or None, a freshly-allocated array is used.
+    **kwargs
+        For other keyword-only arguments, see the
+        :ref:`ufunc docs <ufuncs.kwargs>`.
+
+    Returns
+    -------
+    y : ndarray
+        The vector-matrix product of the inputs.
+
+    Raises
+    ------
+    ValueError
+        If the last dimensions of ``x1`` and the one-but-last dimension of
+        ``x2`` are not the same size.
+
+        If a scalar value is passed in.
+
+    See Also
+    --------
+    vecdot : Vector-vector product.
+    matvec : Matrix-vector product.
+    matmul : Matrix-matrix product.
+    einsum : Einstein summation convention.
+
+    Examples
+    --------
+    Project a vector along X and Y.
+
+    >>> v = np.array([0., 4., 2.])
+    >>> a = np.array([[1., 0., 0.],
+    ...               [0., 1., 0.],
+    ...               [0., 0., 0.]])
+    >>> np.vecmat(v, a)
+    array([ 0.,  4., 0.])
+
+    .. versionadded:: 2.1.0
+    """)
+
 add_newdoc('numpy._core.umath', 'modf',
    """
    Return the fractional and integral parts of an array, element-wise.

diff --git a/numpy/_core/multiarray.py b/numpy/_core/multiarray.py
@@ -83,11 +83,11 @@ def _override___module__():
        'isfinite', 'isinf', 'isnan', 'isnat', 'lcm', 'ldexp', 'less',
        'less_equal', 'log', 'log10', 'log1p', 'log2', 'logaddexp',
        'logaddexp2', 'logical_and', 'logical_not', 'logical_or',
-        'logical_xor', 'matmul', 'maximum', 'minimum', 'remainder', 'modf',
-        'multiply', 'negative', 'nextafter', 'not_equal', 'positive', 'power',
-        'rad2deg', 'radians', 'reciprocal', 'rint', 'sign', 'signbit', 'sin',
-        'sinh', 'spacing', 'sqrt', 'square', 'subtract', 'tan', 'tanh',
-        'trunc', 'vecdot',
+        'logical_xor', 'matmul', 'matvec', 'maximum', 'minimum', 'remainder',
+        'modf', 'multiply', 'negative', 'nextafter', 'not_equal', 'positive',
+        'power', 'rad2deg', 'radians', 'reciprocal', 'rint', 'sign', 'signbit',
+        'sin', 'sinh', 'spacing', 'sqrt', 'square', 'subtract', 'tan', 'tanh',
+        'trunc', 'vecdot', 'vecmat',
    ]:
        ufunc = namespace_names[ufunc_name]
        ufunc.__module__ = "numpy"