0

Say a and b are disjoint 1D complex NumPy arrays, and I do numpy.multiply(a, b, b).
Is b guaranteed to contain the same values as I would get via b[:] = numpy.multiply(a, b)?

I haven't actually been able to produce incorrect results, but I don't know if I'm just being lucky with my particular compilation or platform or if I can actually rely on that, hence the question.

Notice that with float (i.e. real numbers) the answer would obviously yes, since a reasonable implementation can't make it fail, but with complex numbers it's easy for the cross-multiply operation to give incorrect results by writing the real part and then reading an imaginary part:

# say the real part is at [0] and the imaginary part is at [1] and c is the product of a & b
c[0] = a[0] * b[0] - a[1] * b[1]
c[1] = a[0] * b[1] + a[1] * b[0]  # if c[0] overlaps a[0] or b[0] then this is wrong
user541686
  • 205,094
  • 128
  • 528
  • 886
  • I don't see how a reasonable implementation would would fail for complex numbers either, since you have to read both the real and imaginary parts before you can write either. – Dietrich Epp Jun 09 '17 at 02:25
  • `b*=a` is a vectorization of complex scalar mutiplication: `b[i]*=a[i]`. The latter requires fetching `a[i]` and `b[i]`, multiplying them, and storing the result back into `b[i]`. Both real and imaginary parts of both multiplicands are needed for multiplication and must be fetched _before_ the product is calculated. How could this possibly go wrong? – DYZ Jun 09 '17 at 02:26
  • @DietrichEpp, DYZ: See update. – user541686 Jun 09 '17 at 02:31
  • Well, I guess you can look into Numpy source code and check yourself: https://github.com/numpy/numpy/tree/master/numpy – DYZ Jun 09 '17 at 02:33
  • @DYZ: I don't know if I can rely on that, is the problem. What if they use another vectorized implementation in some other place or time? Is there a guarantee provided anywhere? Do they ever address aliasing at all? I haven't seen any evidence of it myself but I don't know. – user541686 Jun 09 '17 at 02:34
  • What better guarantee can you have than the source code??? – DYZ Jun 09 '17 at 02:35
  • To do vectorization properly, you calculate the right hand side then write to the left hand side. It doesn't matter whether the type is complex or not, whether you're doing b *= a or b = b*a or b[:] = b*a, you calculate b*a, then you write the result to b. I think it's highly unlikely that numpy got that wrong, but if you find that to be the case you should submit a bug. – OldGeeksGuide Jun 09 '17 at 02:36
  • @DYZ: Documentation? Or someone voicing the concern on a mailing list at least? I don't know. Source code can easily change; this can easily be an implementation detail that might change. They might use a different (say) vectorized implementation at some other place or time. If the API guarantees it then I can rely on that. – user541686 Jun 09 '17 at 02:36
  • @OldGeeksGuide: OK but then what about behavior on an implementation that *doesn't* vectorize? And I mean whether hypothetical or currently existing (lots of builds of NumPy exist for various platforms, I don't even know about most of them...) – user541686 Jun 09 '17 at 02:37
  • This seems like a weird question to ask. You're basically asking if this particular NumPy operation is buggy. Sure, there are obvious ways they could have gotten it wrong, but that's true for tons of things. – user2357112 Jun 09 '17 at 02:40
  • @user2357112: This *can't* be called a "bug" if there is no such guarantee provided for the case of aliasing; it would just be an implementation detail, hence the entire point of my question... – user541686 Jun 09 '17 at 02:40
  • It doesn't matter whether the underlying implementation uses vector machine instructions or not. When you say b = b*a, the semantics are that you calculate all of b*a before modifying b. The implementation may be able to cheat a bit if it knows that writing part of b won't affect the calculations of other parts, but the semantics don't change. Read up on data dependence analysis and automatic vectorization for more details - here's a good book on the subject: http://www.springer.com/us/book/9780792398097 – OldGeeksGuide Jun 09 '17 at 02:44
  • 1
    @OldGeeksGuide: No, there is no such semantics. You easily get into a slippery slop if you try that with *shifted* arrays; you can't expect NumPy to handle all cases of aliasing, so why this one. Data dependence analysis and automatic vectorization are only about not changing the semantics of the implementation; they provide no guarantee that the implementation is correct. (And there is no guarantee there's *automatic* vectorization occurring here, but I disgress.) If the code itself is implemented like I wrote above, the *correct* result for the compiler would be the *incorrect* result. – user541686 Jun 09 '17 at 02:47
  • They wouldn't provide this operation if they expected it not to work. – user2357112 Jun 09 '17 at 02:50
  • @user2357112: The point is "work" might not be the same thing you expect. Again: they can't reasonably handle *all* cases of aliasing, so an argument could be made fro why this one. – user541686 Jun 09 '17 at 02:51
  • The semantics of a vectorized operation `c[i]=a[i]*b[i]` is that the righthand side iis calculated before the assignment for each `i` (but not necessarily in any particular order of `i`s), regardless of the type of `a` or `b`. – DYZ Jun 09 '17 at 02:51
  • @DYZ: That's just wrong. **There is no guarantee of vectorization semantics in the first place**; this might not even be vectorized, so bringing that up is already wrong. Second, try this: `c = numpy.asarray([1+2j,3+4j,5+6j,7+8j]);` `a = c[1:]; b = c[:3];` `print(a * b);` `a *= b;` `print(a)` Do you honestly call this a bug? I certainly don't. The case of 0 offset is special and it could be debated whether that is a bug or not. – user541686 Jun 09 '17 at 02:56
  • Your example has nothing to do with complex number. Try `c = numpy.asarray([1,3,5,7]);` and then the rest of your code. As I said, the RHS for _any_ `i` is calculated before the assignment for the same `i`. – DYZ Jun 09 '17 at 02:58
  • @DYZ: *"Your example has nothing to do with complex number."* ...Wha...? *Everything* in my post and comments has been about complex numbers. My example *is* using complex numbers. What are you even looking at? – user541686 Jun 09 '17 at 02:59
  • c = numpy.asarray([1,3,5,7]); a = c[1:];b = c[:3];a*b; -> array([ 3, 15, 35]); a*=b; a->array([ 3, 15, 105])`. – DYZ Jun 09 '17 at 03:00
  • @DYZ: Oh, you mean it also fails in the real number case. Yes it does. But that doesn't stop it being a counterexample to your argument right? It makes it even more so. – user541686 Jun 09 '17 at 03:01
  • The point of data dependence analysis is exactly to guarantee that the implementation matches the semantics of the language. To put it simply, there is an implied (and often implemented) temporary variable, i.e. the semantics of `b *= a` is the same as `t = b*a; b = t`. The implementation may try to optimize this by eliminating that intermediate temporary if it can prove (with data dependence analysis) that there it is safe to do so. I guess the question you're asking is if such semantics are well defined and documented in Numpy, or if it's undefined and up to the implementation. – OldGeeksGuide Jun 09 '17 at 03:04
  • In your example, both a and b are aliases for c. When you do `a*b`, the value of `c` does not change, and you get what you expect. When you do `a*=b`, you change the value of `c` as you go. _However_, for any `i`, the RHS is calculated _atomically_ before the assignment, whether it is complex or float. – DYZ Jun 09 '17 at 03:05
  • @OldGeeksGuide; Data dependency analysis **does not guarantee anything about the implementation**. I thought I said this. It only guarantees the compiler output matches the implementation. – user541686 Jun 09 '17 at 03:05
  • I concede your point about numpy, it appears that it does not guarantee anything in the case of aliasing as your example illustrates. This tells me that they don't do dependence analysis, but rely on the developer to avoid aliasing, and that semantics are implementation dependent where there is overlap. This is true whether it's real or complex. My experience with data dependence is with traditional compiled languages automatically vectorizing or otherwise parallelizing, where DD is used to guarantee that the implementation conforms to the language semantics. – OldGeeksGuide Jun 09 '17 at 03:12
  • @OldGeeksGuide: Last time I tried, automatic vectorization was pretty suboptimal; as far as I'm aware people write their own vectorized code when they really care about speed. Automatic vectorization is just for getting icing on the cake with zero effort at the moment. – user541686 Jun 09 '17 at 03:14
  • @DYZ: If you feel strongly that that's the only acceptable semantics then by all means post it as an answer (it's a valid one, whether right or wrong) but I can't find myself agreeing in general without any extra support, because in my experience APIs that don't specify anything about aliasing don't really guarantee anything about that case, except for the sole exception of very formal "lawyer"y specs like the C++ standard. – user541686 Jun 09 '17 at 03:16
  • @DYZ: I actually made a slight edit to further illustrate my point. Now what do you think the answer is? If `a` and `b` are independent, and we do `numpy.multiply(a, b, b[some_slice])` what guarantees are there? You're essentially claiming there is a guarantee for the `some_slice = 0` case but not the `some_slice != 0` case, which is hardly convincing given I haven't seen them mention anything in the documentation regarding the aliasing of `out` with either array. – user541686 Jun 09 '17 at 03:20
  • Your original question is not about aliasing. It is about atomicity. What I claim is that for any operation `c[i]=a[i]*b[i]`, complex or real, aliased or not, the RHS for any `i` is calculated before the assignment for the same `i`, period. Therefore, all results are calculated correctly, but may be overwritten by the assignments. – DYZ Jun 09 '17 at 03:22
  • @DYZ: If I take your claim literally, it has absolutely nothing to do with my question, because yours is about out-of-place operations and mine is specifically about in-place operations. Maybe clarify what you mean because I'm not sure how to make it relevant to this case. – user541686 Jun 09 '17 at 03:24
  • @Mehrdad, I'm sure you are correct when it comes to overlapping vectors, I was mistaken in my thinking on that. Clearly numpy does not make any guarantees for that case. Back to your original question, about complex vs real arrays, I would not expect any differences if there are no cross iteration dependences, i.e. I would expect for each scalar value that it would calculate any individual complex number completely before storing it, rather than, storing the real part before calculating the imag part. I don't know if there is any guarantee in numpy though. – OldGeeksGuide Jun 09 '17 at 03:24
  • In your question, you have `c[0]=...` and `c[1]=...` and wonder if `c[0]` is overwritten before it is (or may be) used to calculate `c[1]`. I claim that because `c[0]` and `c[1]` are two parts of the same number, the assignments will not take place before both parts are calculated. – DYZ Jun 09 '17 at 03:26
  • @OldGeeksGuide: Okay so to clarify, what you're saying is that unless there are cross-iteration dependencies, it is most definitely a bug (unless there's an explicit disclaimer)? – user541686 Jun 09 '17 at 03:26
  • @DYZ: Uhm, `c[0]` is just an **example** of a very dumb NumPy implementation illustrating how the scenario in my question could potentially fail. It is irrelevant to the code in question itself. The question is about `numpy.multiply(a, b, b)` now (formerly `b *= a`). So your argument about `a[i]` has nothing to do with my question because I don't have that anywhere in the code I'm asking about. – user541686 Jun 09 '17 at 03:28
  • Then, you have to change your question, because (1) You explicitly cite the concern of `c[0]` being overwritten, and (2) You claim that }with float (i.e. real numbers) the answer would obviously yes", which it is not, because of aliasing. – DYZ Jun 09 '17 at 03:30
  • Well, after my previous insistence, I'm hesitant to say it's definitely anything :-) I would think it's probably a bug, but I don't know for sure. It's probably best to take up the question with the developers/creators of numpy to see for sure. – OldGeeksGuide Jun 09 '17 at 03:30
  • @DYZ: I think you're not understanding my argument. My question is about the particular case it's about. I'm just giving counterarguments to **your** argument for why **your** rationale fails in cases where I would expect it to hold. I *myself* don't care about the case where `float` would fail; I'm not asking about such a case. But my point is, if your rationale fails for that case, it might very well fail here too; why should I believe otherwise. – user541686 Jun 09 '17 at 03:33
  • @OldGeeksGuide: Haha okay. – user541686 Jun 09 '17 at 03:34

1 Answers1

1

Yes. Complex values ought to be treated atomically. If it doesn't work like that, then it's a bug that we'll fix.

Robert Kern
  • 13,118
  • 3
  • 35
  • 32