Difference between `*` and `*=` operators

[ds-steventondeur]

Hi,

Is the following expected behavior

import quaternionic, numpy as np

q1 = quaternionic.array(np.random.rand(4))
q2 = quaternionic.array(np.random.rand(4))
print(q1 * q2)
q3 = np.copy(q1)
q3 *= q2
print(q3)

prints

[-0.04188505  0.95070465  0.78841118  0.75945989]
[0.01233402 0.00490972 0.00292491 0.00742983]

I would've expected both results to be equal.

Kind regards, steven

[moble]

That code does not work for me; I get a ValueError on the *= line, which is actually what I expect. Which versions of quaternionic and numpy are you running? Mine are the current versions on PyPI:

>>> np.__version__, quaternionic.__version__
('1.20.2', '0.3.4')

More generally, I don't expect in-place operators to work with quaternionic arrays. As far as I can tell, this is essentially a limitation of numpy itself that I can't get around. They are working on making numpy friendlier to custom types like this, but currently I don't see how it's possible. So my advice is to just use the standard operators and store the result as a new array. (If you overwrite the old variable name, it should get garbage collected soon, which frees up the memory as needed.)

If you have such huge arrays that you can't spare the extra storage, or you need to squeeze out a little extra speed, you could start writing your own numba functions using things like quaternionic.algebra_ufuncs.multiply to do the quaternion math — but that's obviously pretty advanced usage.

[ds-steventondeur]

Interesting, it seems that from numpy version >= 1.20 the code indeed throws a ValueError. numpy versions 1.19 and 1.18 don't throw an error but produce these unexpected results. I could reproduce this in a google colab.

I'm working with tensorflow, which still uses numpy 1.18 (tf 2.3) or numpy 1.19 (tf 2.4)

I can certainly use the regular * operator (I'm working on small arrays, memory and/or performance is not a constraint). This is just one of these unexpected gotcha's I guess. A bit of a dangerous one; I'd rather have a crash than a successful run + wrong result, but I understand there's no working around that. I hope tensorflow will soon update their numpy requirements.

Thanks again for the quick reply!

[moble]

Aha! Now I understand what's going on. I should actually be able to make it work (there's a PR going through the pipeline now), but you'll need to avoid np.copy.

There's an unfortunate choice of default value in np.copy. I actually ran into this in a different repo using sub-classed arrays, and complained about it here. It was "fixed" in 1.19, in the sense that you can figure it out if you pay more attention to the documentation and optional arguments. Of course, none of us do, so it's still pretty subtle. To see what's happening, change your code to

q3 = np.copy(q1)
print(q3 * q2)
q3 *= q2
print(q3)

You would expect those two print statements to output the same thing, but they're definitely not! When you just use q3 = np.copy(q1), you get back a plain numpy array of 4 floats that doesn't know anything about quaternions, and q3 * q2 does what it should do — even though this isn't what you meant for it to do — which is to just scalar-multiply the quaternion q2 by each of those floats in turn. Thus, you get a 4x4 quaternionic array. Again, this isn't actually wrong, it's just a surprise the numpy devs have given us.

Unfortunately, numpy<1.20 also seems to have an actual bug where q3 *= q2 tries to "reduce" the result in some way, which is why it gives you back a set of four numbers. (I can't understand what it's actually doing in that reduction; it's not addition or multiplication of any kind...)

To fix this, you should use what the numpy documentation calls the "preferred method for creating an array copy":

q3 = q1.copy()

(For numpy>=1.19, you could also use q3 = np.copy(q1, subok=True), though I don't see why that would be better.) This way, q3 will also be a quaternionic.array, and regular multiplication will work as you expect. In-place multiplication still won't work currently, but I think I've figured out how to enable that.

I've added a note to the README about this np.copy issue, which is the best I can do. Do this, and update to the new version in a little while, and your code should work.

[moble]

Closed by #23