CuPy Accelerates NumPy on the GPU? Hold My Cider, Here's Clojure | Hacker News Reader

CuPy Accelerates NumPy on the GPU? Hold My Cider, Here's Clojure

dragan.rocks

120 pointssigna1146 comments

The benchmark is flawed. Cupy follows Numpy in doing a type-promotion of the data to float64 before computing the covariance.

https://github.com/cupy/cupy/blob/b19577199d3f75c76e23c75194...

The following function in PyTorch completes the float32 computation in 27 ms (1080 Ti). I'm fairly sure that a similar function in Cupy would achieve very similar results.

  def corr(dat):
      cov = (dat @ dat.T).div_(dat.shape[-1]-1)
      sigma = torch.diag(cov).sqrt_()    
      return cov.div_(sigma[None]).div_(sigma[None].T)

The heart of the computation SGEMM is in anycase provided by CUBLAS, and the only areas in which PyT and Cupy differ are in the tricks used to unroll the loop for the broadcasting & L1 op. The performance difference b/w PyT, Cupy and others should be negligible IMO.

Note: np.corrcoef takes about 1.12 s and the above float32 version takes 325ms on the CPU.

So, you first state that the benchmark is flawed, then you confirm the results of the benchmark, and repeat everything that the OP describes as the solution, just with a different tool (PyTorch instead of Clojure)?

BTW did you check that your corr function gives the correct result? It looks to me that you forgot a few steps.

The PyT version above and the both the CPU/GPU versions on Neantherdal work completely in float32. I confirm that the performance of Neantherdal GPU is similar to that in the above PyT version. Yes.

The part that is flawed is that you're comparing this to Numpy (CPU)/ Cupy (GPU) both of which coerce the input array to float64 (for precision reasons) before computing the covariance and correlation. You only need to check the output type of the result to verify this (if the pointer to the code is not sufficient).

What is flawed there? The point of the article is to show that CuPy very often does not accelerate NumPy, especially on consumer-grade hardware. This is something that most users of NumPy/CuPy do not know, and they are led by the docs to think it does.

The reason for that is that CuPy is poorly implemented. And CuPy is poorly implemented because it is constrained by what NumPy does, which, in turn, does stuff that is OK on the CPU, and often translates poorly to the GPU.

I didn't get that from the article at all. There was no "that's because..." or "you need to do this to make it fast..." in the article. Instead, it's "clojure is faster without additional work". While that's super neat and thank you for showing this, bashing python because the library doesn't automatically do that and then hiding behind this silly argument isn't that enlightening.

You can check in the article that I made many checks to make sure that NumPy/CuPy get the data in float32. What do you suggest to do to instruct NumPy/CuPy to "automatically do that". Is there a way to say to python "I want you to use float32 precision for this computation" other than, well, providing everything as float32?

And even if your argument stands, that does not change the fact that CuPy does not accelerate NumPy (in this particular case, but I'd say often).

I didn't read it as bashing Python, only as bashing NumPy/CuPy.

If you look at the article, it is not even using Python anywhere, it uses NumPy/CuPy directly from Clojure. So it really is comparing Neanderthal vs NumPy/CuPy, and uses Clojure in both cases.

> What is /flawed/ there?

You're comparing float32 vs float64 computation. I don't need to tell you how much slower DGEMM is vs SGEMM esp. on the GPU (you mention this in the post yourself!).

Numpy does this for precision reasons, and CuPy simply follows its behavior. This is precisely why I noted that the float32 version runs 3x faster on the CPU.

> The reason for that is that CuPy is poorly implemented.

It's a cheap shot to call something 'poorly implemented' when you don't understand what you're benchmarking.

Can you please provide a benchmark for corrcoef where CuPy is noticeably faster than NumPy on your (and mine) GPU, Nvidia GTX 1080Ti?

I still think I'd have to disagree with you. What was benchmarked was NumPy/CuPy, and the numbers in the article are not flawed. It isn't that they are using NumPy/CuPy wrongly, that's what you'd do, and even if you try really hard to specify everything as float32 it still will have the same performance timing as in the article.

It would be interesting to compare it against a float64 version in Neanderthal as well I agree with that.

That said, a flawed benchmark would mean to me that it isn't indicative of the performance one can expect when actually using the library on real world use case, but for now this benchmark for NumPy/CuPy does seem to be indicative of what you'd expect.

Now, the next question is, for model accuracy vs scale, is going with float64 coercion always the ideal trade off? What if you still needed to squeeze more performance? Is it really a bad idea to do so by going down to float32? Especially considering how much faster GPU can accelerate that?

The flaw, in my opinion, is not necessarily in the benchmark but in the claims and observations around the benchmark. For instance the blog flat out says:

> Nope, we work with float32.

It seems that is not true. The blog should make clear that it isn't straightforward or possible to get cupy to do this.

This is a reply to bearzoo, but we reached the depth level of the comment thread.

>> Nope, we work with float32.

>It seems that is not true.

Well, it is true. We work with float32. I explicitly checked that and NumPy/CuPy answered that the array is indeed float32. The fact, which I didn't know then, is that NumPy/CuPy internally decides on its own to use float64 without warning or possibility for us to order it not to do that. But, it is not I (or "us") that work with float64.

I would say that your comment would stand if I mistakenly ordered NumPy/CuPy to use float64 while claiming that I work with float32 (which could happen hypothetically if there was a non-obvious, but documented option to use float32 that I missed).

Yeah, you passed in float32. The function always promotes to float64 instead. You can argue that NumPy/CuPy should provide a purely float32 option, instead you went something like “look at how my expf beats their exp in benchmarks! What are you talking about ‘not the same function’, I made sure to pass float not double!” (a libc example to illustrate the point). And it’s not clear to me you realized the difference before akssri pointed it out, which renders the benchmarks pretty meaningless.

Where did you get that? Even the main title is refferring to the main point being that CuPy does not accelerate NumPy even in the case where it should be absolutely expected to. Then I used my implemetation to demonstrate that indeed GPU implementation for such a huge matrix should be many times faster.

I never claimed that my library aims for being a replacement for CuPy, or to have any compatibility with NumPy.

It would be more valuable to CuPy developers if I debugged CuPy to discover why that problem exists, but why should I be obligated to? I was writing this for a perspective of a user of these libraries.

Would it be an idea to benchmark Neanderthal with float64 as well, just to gather some data on it?

I agree with both of you, you’re both looking at this thing from a different perspective. It’s perhaps just better to gather timing measurements on a few variants with the trade-offs that each library has made, and how that affects implementation / speed.

It's a lazy benchmark that compares two functions without clarifying that they do different things. If this article was pointing out that the numpy function had poor UX, that would be valid, but presenting a table that suggests Clojure does the same work 200 times faster is blatantly misleading. If you're going to write a custom implementation for one side of the benchmark you either need to (a) make sure it does the same thing as the other side of the benchmark or (b) write a custom implementation for both sides.

> I would love to improve any part of this article, if possible!

You should state that the reason the neanderthal version runs faster is that it is doing the computation in lower precision than numpy / cupy.

The article walks through an investigation of why cupy's result is underwhelming (is the input data accidentally fp64? is the input data on the cpu? is the computation happening on the cpu?), so you should finish it by explaining that numpy and cupy do the computation in fp64.

Thanks for sharing! Why does CuPy switch to float64 in this case? Is it a bug?

Numpy/Cupy do explicit coercing to float64. There is no documentation for why this is done, but since GEMM is used to compute the covariates (summation over data points), it makes sense to increase the precision.

You could get away using a double to only accumulate the sum, but it's a pain to write such mixed-precision (slow) function in C and then wrap it in Python, esp. for something like this.

(See, https://github.com/numpy/numpy/blob/v1.17.0/numpy/lib/functi...)

CuPy is a drop in replacement for NumPy. Ask Numpy why they do that.

I'm not sure it is flawed. It just seems that you found the root cause of why NumPy/CuPy are slower, but it seems there is no way as a user of NumPy/CuPy to perform the float32 computation, even though you were explicit and said `:dtype "float32"`.

All these libraries under the hood delegate to the same things for the actual computation, so the only advantage of one over the other would be this type of stuff.

That NumPy/CuPy coerce to a float64 under the hood, and gives us no way to force it not to do so, is a disadvantage that they have in my opinion.

In any case, like the article says, the performance of NumPy/CuPy is still good enough, and maybe NumPy thinks it knows better then it's users that they should be using higher precision always.

Edit: I know I'll face the brigade of NumPy users. But I invite people to engage in the discussion unbiased. With a switch from CPU to GPU you'd expect more than 29% performance boost. If CuPy went with a float32 proc, it would give you a 96% performance increase! So ya, I think we can start to argue CuPy isn't doing the right thing here.

> user of NumPy/CuPy to perform the float32 computation

This is just getting tiring.

Numpy and Cupy are perfectly capable of doing float32 computation - their only "fault" is that they coerce data to float64 in this one fairly unimportant function (which you can implement to your liking in 3-4 LOC). Hell, PFNs entire deep learning system, Chainer (which mind you, both predated and inspired PyTorch, and is still quite competitive with it), is built entirely on top of Cupy!

Benchmarks make sense only when the outputs are the same - in this case, they certainly are not. It's the responsibility of the author to make sure that the outputs are the same (a real benchmark), or to argue that Numpy/Cupy are wasting time by using float64 in xp.cov (an issue of implementation). He does neither.

> Numpy and Cupy are perfectly capable of doing float32 computation

You're arguing that they would perfectly be capable of an implementation which would use float32, ya sure, if you change their source code you can have them be faster, that's an obvious to me? The question is, what kind of performance can you expect as a user using them as a library, and this benchmark is revealing of that in this case. What is wrong with it?

It's like saying that you can change the implementation of Python and it will make Python code run faster, like sure you can, but you don't benchmark hypothetical future versions, the current version of CuPy uses float64 for this function, and that makes it that it barely runs faster then the CPU version. Point in case, I don't know what else you're trying to disagree about here.

> their only fault is that they coerce data to float64 in this one fairly unimportant function

I don't know that, all I know is when picking one function and benchmarking it, we found a flaw which result in performance gains from CPU to GPU to be extremely underwhelming. How many more if we started to benchmark the rest?

> if you change their source code you can have them be faster, that's an obvious to me? The question is, what kind of performance can you expect as a user using them as a library, and this benchmark is revealing of that in this case. What is wrong with it?

You don't need to change the numpy or cupy source code, just write your own version of the function in a few lines of code using other primitives that numpy/cupy provides.

Note that is exactly what is done for the neanderthal version.

How am I supposed to know that corrcoef is slow and I should rewrite it? I didn't know before, now thanks to this article I do.

> Power users even know that Nvidia provides a drop-in replacement, CuPy, that brings the power of GPU acceleration.

> but Nvidia provides CuPy

> It's Nvidia's stuff

Is CuPy actually developed by Nvidia? It looks like the original authors were at Preferred Networks, and none of the top few contributors on github work at Nvidia.

If anyone's looking for NumPy that works on distributed machines with GPUs, see also Legate, which is produced directly by NVIDIA:

https://research.nvidia.com/publication/2019-11_Legate-NumPy...

https://developer.nvidia.com/legate-early-access

Disclaimer: I work on Legion, the runtime system that Legate sits on top of, but otherwise don't have anything to do with the Legate project.

CuPy is not developed by NVIDIA, but by Preferred Networks, the company behind Chainer.

Uses free Clojure library Neanderthal: https://github.com/uncomplicate

Also uses the book Numerical Linear Algebra for Programmers: https://aiprobook.com/numerical-linear-algebra-for-programme...

One of the things I never understood about the hype around these Clojure GPU libraries is that there is a lot of marketing jumbo about how high-level it all is. Even the book for these libraries has the words "no C++!" half a dozen times on its ordering page [0].

However, these are just high-level libraries and you cannot write the GPU shaders in Clojure. You still must use C or C++ for that (depending on if you are going the CL or Cuda route), and your C/C++ code is embedded or called from the Clojure side. This is no different than using another high-level language like Python -- I mean it's focusing on a rather misleading and irrelevant detail. You can use Clojure or Python (or many other languages) to call OpenCL or Cuda shaders, or you can just use the libraries that call them for you. Clojure is not special in this regard.

[0] https://aiprobook.com/deep-learning-for-programmers/

The book you've linked contains exactly 0 lines of C or C++ code. 100% of the code in the book is Clojure.

I think you missed the point of my comment. You don’t have to write your own C++ code when using high-level python GPU libraries either.

Good for these python GPU libraries. You don't have to do that in Clojure too.

Hum.. I'm not sure there's that much hype. It's just saying that Clojure can also be used as a high level layer for fast linear algebra. I do see some claims that as a high level layer, Clojure works even better, it gives you a full REPL, and the power of macros and all that.

It also goes to show you that there's barely any overhead added by the Clojure layer, which is something you want from a higher level layer.

Maybe you need some additional context, but in the Clojure world, there are many offerings for that layering, core.matrix, neanderthal, NumPy, CuPy and others can all be used. So I guess a lot of what draganj is trying to show is he believes neanderthal is the layer with the least overhead, which everyone can agree to disagree on.

You're never going to beat an assembly level implementation, or 50 year's worth of optimization done to a Fortran implementation, or straight openCL/cuda implementation in a higher level language, especially not in Clojure. So in the market of high level languages like Julia, Clojure, Python, R, Swift, etc. and their corresponding libraries, they all are just competing on being more user friendly, having better more useful/flexible abstractions, a better development flow, and the least amount of overhead. So that's what you'll be hearing them all fight over.

Common Lispers have had the ability to write CUDA kernels (with shaders) on-the-fly with a DSL in cl-cuda for many years now. See,

https://github.com/takagi/cl-cuda

Cupy (a project in which takagi is/was actively involved) also has the ability to compile kernels in the Python REPL, but the kernel code needs to be written in C and passed in as a string (pyOpencl can do the same for OpenCL). In fact, not too many years ago, the kernels for gradient descent, max pooling etc. were all entirely in Python in Chainer. Projects like JAX take this to another level by having the ability to transform the bytecode of a restricted class of Python functions straight into CUDA kernels.

I can't speak much about Clojure/Neanderthal, but I'd advise the people to stop dissing projects they are not familiar with, least of all using silly benchmarks like the above.

Good advice. And, yet, you've taken it pretty seriously to diss Clojure/Neanderthal and my blog post, mostly by talking about unrelated stuff and projects that the post didn't even mention. And while introducing these themes left and right you didn't even bother to show some code related to this topic, just a suggestion of great projects by cool people.

Yes, you showed the PyTorch code related to the blog post that confirms what the blog post says, but when I pointed out that the code has incorrect functionality (by missing some calculations) you didn't even bother to correct it, or to confirm that the code is good and that I'm wrong.

So, it seems that your standard is that it is enough for one side to throw bits and pieces around and call it a day, and for the other to run around and prove that their stuff is better than everything that could possibly be done in every technology.

I choose to stick to the theme. The theme is CuPy, NumPy, Clojure & Neanderthal. The related theme could be code in another technology. Great - write about it. But, even if every other technology were a million times better than what I describe in the article, it does not change the fact that CuPy and NumPy have the issue I've described.

> And, yet, you've taken it pretty seriously to diss Clojure/Neanderthal and my blog post

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.

I still don't get how the fact that someone could implement the same thing that I did in Clojure in PyTorch has anything to do with NumPy and CuPy, or my benchmark?

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...

- Function in Cupy takes 29.4ms, Numpy takes 427 ms. Happy ?

- 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.

Thank you so much, that's phenomenal news for me! (Since I can make neanderthal code go at 23ms (GPU) and 3XX ms (CPU) when I implement it as NumPy/CuPy/PyTorch does (sans float64 conversion, of course) You saved me from having to fiddle with Python (which I don't particularly enjoy). Thanks again!

Can you please post your implementation of this function, here, so I can try it on my machine and compare it to Neanderthal?

Great article. I learned a bunch of things I didn’t know about Numpy and Clojure. As far as I’m concerned the Dragan does rock.

Or just use Julia where the speed is out of the box.

Oh? I don't suppose that julia wraps lapack functions for no good reason. What makes you say this

A lot of the wrappers are largely legacy at this point. BLAS1 and BLAS2 operations are handled in pure Julia through the broadcast system, there's no performance degradation there and the fusion of operations actually makes for a lot of performance advantages over simple BLAS calls. For BLAS3, tools like LoopVectorization let you write a pure BLAS3 in Julia and get good performance, Gaius.jl is a good example of this by matching OpenBLAS with pure Julia code. Of course you still have to write out how to do blocking, but that's about it. At this point it's mostly needing more coverage before deprecating BLAS bindings in a Julia 2.0 down the line.

https://github.com/MasonProtter/Gaius.jl

Thanks, that is pretty impressive

I'm the dev of Gaius.jl. I'm not sure it will ever end up Julia's Base, but I think of it as a prototype implementation that might inspire (or possibly evolve into) a real implementation that could end up being bundled in the language eventually.

Julia really does have the tools to solve this problem natively, and they're maturing quite a bit more quickly than I ever would have guessed a year ago.

Last year if you had asked me I would have predicted the development of something like LoopVectorization.jl as appearing maybe by 2022 or something. Instead it's here now, albeit with a lot of work ahead to make it more robust.

Reply on news.ycombinator.com