[Numpy-discussion] Two questions about PEP 465 dot product

[Alexander Belopolsky]

1. Is there a simple expression using existing numpy functions that
implements PEP 465 semantics for @?

2. Suppose I have a function that takes two vectors x and y, and a matrix M
and returns x.dot(M.dot(y)). I would like to "vectorize" this function so
that it works with x and y of any ndim >= 1 and M of any ndim >= 2 treating
multi-dimensional x and y as arrays of vectors and M as an array of
matrices (broadcasting as necessary). The result should be an array of xMy
products. How would I achieve that using PEP 465's @?

[Nathaniel Smith]

Post by Alexander Belopolsky
1. Is there a simple expression using existing numpy functions that

Not yet.

Post by Alexander Belopolsky
2. Suppose I have a function that takes two vectors x and y, and a matrix M
and returns x.dot(M.dot(y)). I would like to "vectorize" this function so
that it works with x and y of any ndim >= 1 and M of any ndim >= 2 treating
multi-dimensional x and y as arrays of vectors and M as an array of matrices
(broadcasting as necessary). The result should be an array of xMy products.

(x[..., np.newaxis, :] @ M @ y[..., :, np.newaxis])[..., 0, 0]

Alternatively, you might prefer something like this (though it won't
yet take advantage of BLAS):

np.einsum("...i,...ij,...j", x, M, y)

Alternatively, there's been some discussion of the possibility of
adding specialized gufuncs for broadcasted vector-vector,
vector-matrix, matrix-vector multiplication, which wouldn't do the
magic vector promotion that dot and @ do.

-n

[Alexander Belopolsky]

.. there's been some discussion of the possibility of
adding specialized gufuncs for broadcasted vector-vector,
vector-matrix, matrix-vector multiplication, which wouldn't do the

This would be nice. What I would like to see is some consistency between
multi-matrix
support in linalg methods and dot.

For example, when A is a matrix and b is a vector and

a = linalg.solve(A, b)

then

dot(A, a) returns b, but if either or both A and b are stacks, this
invariant does not hold. I would like
to see a function (say xdot) that I can use instead of dot and have xdot(A,
a) return b whenever a = linalg.solve(A, b).

Similarly, if w,v = linalg.eig(A), then dot(A,v) returns w * v, but only
if A is 2d.

[Nathaniel Smith]

Post by Alexander Belopolsky

.. there's been some discussion of the possibility of
adding specialized gufuncs for broadcasted vector-vector,
vector-matrix, matrix-vector multiplication, which wouldn't do the

This would be nice. What I would like to see is some consistency between

multi-matrix

Post by Alexander Belopolsky
support in linalg methods and dot.
For example, when A is a matrix and b is a vector and
a = linalg.solve(A, b)
then
dot(A, a) returns b, but if either or both A and b are stacks, this

invariant does not hold. I would like

Post by Alexander Belopolsky
to see a function (say xdot) that I can use instead of dot and have

xdot(A, a) return b whenever a = linalg.solve(A, b).

I believe this equivalence holds if xdot(x, y) = x @ y, because solve()
does follow the pep 465 semantics for shape handling. Or at least, it's
intended to. Of course we will also expose pep 465 matmul semantics under
some name that doesn't require the new syntax (probably not "xdot" though
;-)).

Post by Alexander Belopolsky
Similarly, if w,v = linalg.eig(A), then dot(A,v) returns w * v, but only

if A is 2d.

Again A @ v I believe does the right thing, though I'm not positive -- you
might need a swapaxes or matvec or something. Let us know if you work it
out :-).

Note that it still won't be equivalent to w * v because w * v doesn't
broadcast the way you want :-). You need w[..., np.newaxis, :] * v, I think.

-n

Post by Alexander Belopolsky
_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Benjamin Root]

At some point, someone is going to make a single documentation page
describing all of this, right? Tables, mathtex, and such? I get woozy
whenever I see this discussion go on.

Ben Root

Post by Alexander Belopolsky

.. there's been some discussion of the possibility of
adding specialized gufuncs for broadcasted vector-vector,
vector-matrix, matrix-vector multiplication, which wouldn't do the

This would be nice. What I would like to see is some consistency

between multi-matrix

Post by Alexander Belopolsky
support in linalg methods and dot.
For example, when A is a matrix and b is a vector and
a = linalg.solve(A, b)
then
dot(A, a) returns b, but if either or both A and b are stacks, this

invariant does not hold. I would like

Post by Alexander Belopolsky
to see a function (say xdot) that I can use instead of dot and have

xdot(A, a) return b whenever a = linalg.solve(A, b).
does follow the pep 465 semantics for shape handling. Or at least, it's
intended to. Of course we will also expose pep 465 matmul semantics under
some name that doesn't require the new syntax (probably not "xdot" though
;-)).

Post by Alexander Belopolsky
Similarly, if w,v = linalg.eig(A), then dot(A,v) returns w * v, but

only if A is 2d.
might need a swapaxes or matvec or something. Let us know if you work it
out :-).
Note that it still won't be equivalent to w * v because w * v doesn't
broadcast the way you want :-). You need w[..., np.newaxis, :] * v, I think.
-n

Post by Alexander Belopolsky
_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Nathaniel Smith]

Post by Benjamin Root
At some point, someone is going to make a single documentation page

describing all of this, right? Tables, mathtex, and such? I get woozy
whenever I see this discussion go on.

That does seem like a good idea, doesn't it. Following the principle that
recently-confused users write the best docs, any interest in taking a shot
at writing such a thing?

-n

[Benjamin Root]

That assumes that the said recently-confused ever get to the point of
understanding it... and I personally don't do much matrix math work, so I
don't have the proper mental context. I just know that coworkers are going
to be coming to me asking questions because I am the de facto "python guy".
So, having a page I can point them to would be extremely valuable.

Ben Root

Post by Benjamin Root

Post by Benjamin Root
At some point, someone is going to make a single documentation page

describing all of this, right? Tables, mathtex, and such? I get woozy
whenever I see this discussion go on.
That does seem like a good idea, doesn't it. Following the principle that
recently-confused users write the best docs, any interest in taking a shot
at writing such a thing?
-n
_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Nathaniel Smith]

Post by Benjamin Root
That assumes that the said recently-confused ever get to the point of

understanding it...

Well, I don't think it's that complicated really. For whatever that's worth
:-). My best attempt is here, anyway:

https://www.python.org/dev/peps/pep-0465/#semantics

The short version is, for 1d and 2d inputs it acts just like dot(). For
higher dimension inputs like (i, j, n, m) it acts like any other gufunc
(e.g., everything in np.linalg) -- it treats this as an i-by-j stack of
n-by-m matrices and is vectorized over the i, j dimensions. And 0d inputs
are an error.

-n

[Benjamin Root]

Then add in broadcasting behavior...

Post by Benjamin Root

Post by Benjamin Root
That assumes that the said recently-confused ever get to the point of

understanding it...
Well, I don't think it's that complicated really. For whatever that's
https://www.python.org/dev/peps/pep-0465/#semantics
The short version is, for 1d and 2d inputs it acts just like dot(). For
higher dimension inputs like (i, j, n, m) it acts like any other gufunc
(e.g., everything in np.linalg) -- it treats this as an i-by-j stack of
n-by-m matrices and is vectorized over the i, j dimensions. And 0d inputs
are an error.
-n
_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Nathaniel Smith]

Post by Benjamin Root
Then add in broadcasting behavior...

Vectorized functions broadcast over the vectorized dimensions, there's
nothing special about @ in this regard.

-n

Post by Benjamin Root

Post by Nathaniel Smith

Post by Benjamin Root
That assumes that the said recently-confused ever get to the point of

understanding it...

Post by Benjamin Root

Post by Nathaniel Smith
Well, I don't think it's that complicated really. For whatever that's
https://www.python.org/dev/peps/pep-0465/#semantics
The short version is, for 1d and 2d inputs it acts just like dot(). For

higher dimension inputs like (i, j, n, m) it acts like any other gufunc
(e.g., everything in np.linalg) -- it treats this as an i-by-j stack of
n-by-m matrices and is vectorized over the i, j dimensions. And 0d inputs
are an error.

Post by Benjamin Root

Post by Nathaniel Smith
-n
_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
http://mail.scipy.org/mailman/listinfo/numpy-discussion

[Alexander Belopolsky]

For higher dimension inputs like (i, j, n, m) it acts like any other
gufunc (e.g., everything in np.linalg)

Unfortunately, not everything in linalg acts the same way. For example,
matrix_rank and lstsq don't.

[Charles R Harris]

Post by Alexander Belopolsky
1. Is there a simple expression using existing numpy functions that
2. Suppose I have a function that takes two vectors x and y, and a matrix
M and returns x.dot(M.dot(y)). I would like to "vectorize" this function
so that it works with x and y of any ndim >= 1 and M of any ndim >= 2
treating multi-dimensional x and y as arrays of vectors and M as an array
of matrices (broadcasting as necessary). The result should be an array of

If you are willing to run Python 3.5 (use 3.6.0a3, a4 crawls with the
bugs), you can use gh-5878 <https://github.com/numpy/numpy/pull/5878>. The
override mechanisms are still in process in Nathaniel's PR, so that may
change. I'd welcome any feedback.

Chuck

[Charles R Harris]

On Fri, May 22, 2015 at 12:02 AM, Charles R Harris <

Post by Charles R Harris

Post by Alexander Belopolsky
1. Is there a simple expression using existing numpy functions that
2. Suppose I have a function that takes two vectors x and y, and a matrix
M and returns x.dot(M.dot(y)). I would like to "vectorize" this function
so that it works with x and y of any ndim >= 1 and M of any ndim >= 2
treating multi-dimensional x and y as arrays of vectors and M as an array
of matrices (broadcasting as necessary). The result should be an array of

If you are willing to run Python 3.5 (use 3.6.0a3, a4 crawls with the
bugs), you can use gh-5878 <https://github.com/numpy/numpy/pull/5878>.
The override mechanisms are still in process in Nathaniel's PR, so that may
change. I'd welcome any feedback.

Oops, make the 3.5.0a3.

Chuck