[Numpy-discussion] Support of '@='?

[Hans Larsen]

I have Python 3-5-1 and NumPy 1-11! windows 64bits!
When will by side 'M=***@P' be supported with 'M@=P'???:-(

--
Hans Larsen Galgebakken Sønder 4-11A 2620 Albertslund Danmark/Danio

[Sebastian Berg]

On Mi, 2016-06-22 at 02:38 +0200, Hans Larsen wrote:
> I have Python 3-5-1 and NumPy 1-11! windows 64bits!
> When will by side 'M=***@P' be supported with 'M@=P'???:-(
>

When someone gets around to making it a well defined operation? ;) Just
to be clear, `M = M @ P` is probably not what `M @= P` is, because the
result of that should probably be `temp = M @ P; M[...] = temp`.
Now this operation needs copy back to the original array from a
temporary array (you can't do it in-place, because you still need the
values in M after overwriting them if you would).

Just if you are curious why it is an error at the moment. We can't have
it be filled in by python to be not in-place (`M = M @ P` meaning), but
copying over the result is a bit annoying and nobody was quite sure
about it, so it was delayed.

- Sebastian

[Marten van Kerkwijk]

>
> Just if you are curious why it is an error at the moment. We can't have
> it be filled in by python to be not in-place (`M = M @ P` meaning), but
> copying over the result is a bit annoying and nobody was quite sure
> about it, so it was delayed.


The problem with using out in-place is clear from trying `np.matmul(a, a,
out=a)`:
```
In [487]: a
array([[ 1. , 0. , 0. ],
[ 0. , 0.8660254, 0.5 ],
[ 0. , -0.5 , 0.8660254]])

In [488]: np.matmul(a, a)
Out[488]:
array([[ 1. , 0. , 0. ],
[ 0. , 0.5 , 0.8660254],
[ 0. , -0.8660254, 0.5 ]])

In [489]: np.matmul(a, a, out=a)
Out[489]:
array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]])
```

It would seem hard to avoid doing the copying (though obviously one should
iterate over higher dimensiones, ie., temp.shape = M.shape[-2:]). Not
dissimilar from cumsum etc which are also not true ufuncs (but where things
can be made to work by ensuring operations are doing in the right order).

-- Marten

[Nathaniel Smith]

To repeat and (hopefully) clarify/summarize the other answers:

It's been left out on purpose so far.

Why was it left out? A few reasons:

- Usually in-place operations like "a += b" are preferred over the
out-of-place equivalents like "a[...] = a + b" because they avoid some
copies and potentially large temporary arrays. But for @= this is
impossible -- you have to make a temporary copy of the whole matrix,
because otherwise you find yourself writing output elements on top of input
elements that you're still using. So it's probably better style to write
this as "a[...] = a @ b": this makes it more clear to the reader that a
potentially large temporary array is being allocated.

- The one place where this doesn't apply, and where "a @= b" really could
be a performance win, is when working with higher dimensional stacks of
matrices. In this case we still have to make a temporary copy of each
matrix, but only of one matrix at a time, not the whole stack together.

- But, not that many people are using matrix stacks yet, and in any case "a
@= b" is limited to cases where both matrices are square. And making it
efficient in the stacked case may require some non-trivial surgery on the
internals. So there hasn't been much urgency to fix this.

My guess is that eventually it will be supported because the stacked matrix
use case is somewhat compelling, but it will take a bit until someone
(maybe you!) decides they care enough and have the time/energy to fix it.

-n
On Jun 21, 2016 17:39, "Hans Larsen" <***@mail.dk> wrote:

> I have Python 3-5-1 and NumPy 1-11! windows 64bits!
> When will by side 'M=***@P' be supported with 'M@=P'???:-(
>
> --
> Hans Larsen Galgebakken SÞnder 4-11A 2620 Albertslund Danmark/Danio
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-***@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>

[Hans Larsen]

I thank you for the answer!!!

Bu this with "temp" is also in core-Python bihind the scene!;-)



Den 22-06-2016 kl. 21:39 skrev Nathaniel Smith:
>
> To repeat and (hopefully) clarify/summarize the other answers:
>
> It's been left out on purpose so far.
>
> Why was it left out? A few reasons:
>
> - Usually in-place operations like "a += b" are preferred over the
> out-of-place equivalents like "a[...] = a + b" because they avoid some
> copies and potentially large temporary arrays. But for @= this is
> impossible -- you have to make a temporary copy of the whole matrix,
> because otherwise you find yourself writing output elements on top of
> input elements that you're still using. So it's probably better style
> to write this as "a[...] = a @ b": this makes it more clear to the
> reader that a potentially large temporary array is being allocated.
>
> - The one place where this doesn't apply, and where "a @= b" really
> could be a performance win, is when working with higher dimensional
> stacks of matrices. In this case we still have to make a temporary
> copy of each matrix, but only of one matrix at a time, not the whole
> stack together.
>
> - But, not that many people are using matrix stacks yet, and in any
> case "a @= b" is limited to cases where both matrices are square. And
> making it efficient in the stacked case may require some non-trivial
> surgery on the internals. So there hasn't been much urgency to fix this.
>
> My guess is that eventually it will be supported because the stacked
> matrix use case is somewhat compelling, but it will take a bit until
> someone (maybe you!) decides they care enough and have the time/energy
> to fix it.
>
> -n
>
> On Jun 21, 2016 17:39, "Hans Larsen" <***@mail.dk
> <mailto:***@mail.dk>> wrote:
>
> I have Python 3-5-1 and NumPy 1-11! windows 64bits!
> When will by side 'M=***@P' be supported with 'M@=P'???:-(
>
> --
> Hans Larsen Galgebakken Sønder 4-11A 2620 Albertslund Danmark/Danio
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-***@scipy.org <mailto:NumPy-***@scipy.org>
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-***@scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion

--
Hans Larsen Galgebakken Sønder 4-11A 2620 Albertslund Danmark/Danio