BTW, your PyTorch code is still incorrect (or so it seems to me although I don't use PyTorch so I can't try it on the computer). The formula for correlation requires division by sigma_x * sigma_y (which has dimension n x n), and you are dividing by (sigma_x)^2 (which has dimension n). So you still forgot at least one L2 operation that computes all combinations of sigma_x_y. A couple operations here, a couple operations there, an edge case here, and edge case there, it adds up. That's why people use NumPy/CuPy after all...
- Broadcasting semantics + division takes care of the outer-product normalization. This is 2 L1 ops in size of the matrix & the input (~ xSCAL).
Pedantry is still not an argument.
Can you please post your implementation of this function, here, so I can try it on my machine and compare it to Neanderthal?
I have not - all I've said so far is that your benchmark is flawed.
The fact that the code fragment above assumes zero mean data (thus using 2 fewer L1 ops) doesn't change a single thing in anything that has been written; to wit, the timings change to 28.6ms (GPU) and 333 ms (CPU). Pedantry is not an argument.