As far as I can tell this article doesn't explain why this happens on the GPU (for example, why Tensorflow's reduce_sum is non-deterministic). My hypothesis is that this is entirely due to concurrency: if the same code can be run in two or more different interleavings, they can produce different results. This is corroborated by the first answer here [0].
If so, this exact same issue happens in CPU code as well: have two or more threads, run the program many times, observe different interleavings that expose race conditions which (depending on the algorithm) may or may not produce different results. This can happen even if you don't use floating point, and has nothing to do with floating point non-determinism itself. For example, have a thread print "Hello" and another thread print "World"; even without tearing, you may see either Hello World or World Hello on the screen.
Now, proper floating point non-determinism happens in two cases. One is that when you run the same code in two different architectures you could have different answers (because of rounding modes, or because some architecture doesn't support subnormal numbers or signaling nans, because transcedental functions like sine are implemented with different accuracy, etc). In this case it's deterministic when run the same in the same machine, but may run differently in another machine with a different architecture.
The other case is that some "optimizations" actually break your code if applied carelessly (you enable those broken optimizations with -ffast-math in C for example). Among other things, this may break numerical stability of algorithms like Kahan summation. And, if you let the compiler decide which exact optimizations will be applied and in what order, you get non-determinism between different compilers. So in this case it's deterministic when compiled with the same compiler, but may run differently with another compiler.
[0] https://stackoverflow.com/questions/50744565/how-to-handle-n...