CuPy: NumPy and SciPy for GPU

It's interesting to see hardware/software/API co-development in practice again.

The last time I think this happen at market-scale was early 3d accelerator APIs? Glide/opengl/directx. Which has been a minute! (To a lesser extent CPU vectorization extensions)

Curious how much of Nvidia's successful strategy was driven by people who were there during that period.

Powerful first mover flywheel: build high performing hardware that allows you to define an API -> people write useful software that targets your API, because you have the highest performance -> GOTO 10 (because now more software is standardized on your API, so you can build even more performant hardware to optimize its operations)

An excellent example of Array API usage can be found in scikit-learn. Estimators written in NumPy are now operable on various backends courtesy of Array API compatible libraries such as CuPy and PyTorch.

https://scikit-learn.org/stable/modules/array_api.html

Disclosure: I'm a CuPy maintainer.

And of course the native Python solution is memoryview. If you need to inter-operate with libraries like numpy but you cannot import numpy, use memoryview. It is specifically for fast low-level access which is why it has more C documentation than Python documentation: https://docs.python.org/3/c-api/memoryview.html

Jax variables are immutable.

Code written for CuPy looks similar to numpy but very different from Jax.

For those interested in the NumPy/SciPy API coverage in CuPy, here is the comparison table:

https://docs.cupy.dev/en/latest/reference/comparison.html

Indeed, has anyone so far successfully drop-in replaced numpy in a project with this cupy and achieved massive improvements? Because, you know, when dealing with GPU it is very important to actually understand how data flows back and forth to it, not only the algorithmic nature of the code written.

As a sidenote, it is funny how this gets released in 2024, and not in say 2014...

The moving-data-around cost is conventional wisdom in GP-GPU circles.

Is it changing though? Not only do PCIe interfaces keep doubling in performance, but CPU-GPU memory coherence is a thing.

I guess it depends on your target: 8x H100s across a PCIe bridge is going to have quite different costs vs an APU (which have gotten to be quite powerful, not even mentioning MI300a)

Exactly my experience. You end up having to juggle a whole different set of requirements and design factors in addition to whatever it is that you’re already doing. Usually after a while the results are worth it, but I found the “drop-in” idea to be slightly misleading. Just because the API is the same does not make it a drop-in replacement.

When I packaged the ROCm libraries that shipped in the Ubuntu 24.04 universe repository, I built and tested them with almost every discrete AMD GPU architecture from Vega to CDNA 2 and RDNA 3 (plus a few APUs). None of that is officially supported by AMD, but it is supported by me on a volunteer basis (for whatever that is worth).

I think that every library required to build cupy is available in the universe repositories, though I've never tried building it myself.

Fingers crossed that all future AMD parts ship with full ROCm support.

Go for it! There's nothing to lose.

You could checkout some of EuroCC's courses. That should get you up to speed. https://www.eurocc-access.eu/services/training/

You'll do fine :) PyTorch has an API that is somewhat similar to numpy, although if you've never programmed a GPU you might want to get up to speed on that first.

Yes it does: https://numba.readthedocs.io/en/stable/cuda/kernels.html

It does (for NVIDIA at least)

Mentioned above:

https://data-apis.org/array-api-compat/

It's typically just importing torch and s/np/torch, not too different from NumPy -> CuPy. Try it in your own code and see!

Personally, I prefer CuPy over your library. For example, your vectorAdd.cu implementation at https://github.com/eyalroz/cuda-api-wrappers/blob/master/exa... is much longer than a similar CuPy implementation:

    import cupy as cp

    vector_add = cp.RawKernel("""
    extern "C" __global__
    void vector_add(const float *A, const float *B, float *C, int num_elements) {
        int i = blockDim.x * blockIdx.x + threadIdx.x;

        if (i < num_elements) { C[i] = A[i] + B[i]; }
    }
    """, "vector_add")

    num_elements = 50_000

    block_size = 256
    # round up to next multiple of block_size
    grid_size = (num_elements + block_size - 1) // block_size

    a = cp.random.rand(num_elements, dtype=cp.float32)
    b = cp.random.rand(num_elements, dtype=cp.float32)
    c = cp.zeros(num_elements, dtype=cp.float32)

    args = (a, b, c, num_elements)

    print(f"[Vector addition of {num_elements} elements]")
    print(f"CUDA kernel launch with {grid_size} blocks of {block_size} threads each")

    vector_add((grid_size,), (block_size,), args)

    incorrect = cp.abs(a + b - c) > 1e-5

    if cp.any(incorrect):
        print("Result verification failed at element", cp.argmax(incorrect))

    print("Test PASSED")
    print("SUCCESS")

Although I have to concede that the automatic grid size computation in cuda-api-wrappers is nice.

A few marketing tips for your README:

* Put a code example directly at the top. You want to present the selling points of your library to the reader as fast as possible. For reference, look at the CuPy README https://github.com/cupy/cupy?tab=readme-ov-file#cupy--numpy-... which immediately shows reader what it is good for. Your README starts with lots of text, but nobody reads text anymore these days. A link to examples is almost at the end, and then the examples are deeply nested.

* The first links in the README should link to your own library, for example to documentation or examples. You do not want to lead the reader away from your GitHub page.

* Add syntax highlighting with "cpp" after triple backticks:

    ```cpp
    <code here>
    ```

Jax is, as far as an outsider can tell, Google research/DeepMind’s primary computation library now, not Tensorflow. So it’s safer than most Google projects, unless you think they’re secretly developing their third tensor computation library in 10 years (I suppose if anyone was going to do it, it would be Google).

this is not what I wanted to hear. NOT AT ALL. Please whisper sweet lies into my ears.

Actually yes it does....except I seem to remember that it doesn't go back that far in cuda versions. I can't seem to find it again right now.

The condos 200-employee threshold licence change is problematic for some.

I'm okay with containers generally, I think. Is this a situation where you put your code into the container and run it, or does the code make calls to the container's gpu?

The largest speedup I have seen was for a quantum mechanics simulation where I needed to repeatedly calculate all eigenvalues of Hermitian matrices (but not necessarily their eigenvectors).

This was basically the code needed:

    import scipy.linalg as la

    if cuda:
        import cupy as cp
        import cupy.linalg as cla

        ε = cp.asnumpy(cla.eigvalsh(cp.asarray(H)))
    else:
        ε = la.eigvalsh(H)

Not OP, but think about stuff like FFTs or Matmuls. It's not even a competition, GPUs win when the algorithm is somewhat suitable and you're dealing with FP32 or lower precision.

What about CPU-only loads? If one wants to write code that'll eventually run in both CPU and GPU but in the short-to-mid term will only be used in CPU? Since JAX natively support CPU (with numpy backend), but CuPy doesn't, this seems like a potential problem for some.

> CuPy does not support automatic gradient computation, so if you do deep learning, use JAX instead

DL is major use case; is CuPy planning on adding auto gradient comp?

Yeah, Jax provides a one-to-one reimplementation of the Numpy interface, and a decent chunk of the scipy interface. Random number handling is a bit different, but Numpy random number handling seeeeems to be trending in the Jax direction (explicitly passed RNG objects).

Jax also provides back-propagation wherever possible, so you can optimize.

yes

My systems are sparse, but might not fit on the GPU when factorized. Actually, usually I do CPU stuff with lots of ram, and Pardiso, so it isn’t an issue.

But I was hoping to try out something like ILU+bicgstab on the GPU and the python-verse seems like it has the lowest barrier-to-entry for just playing around.

Oh dang, that’s pretty awesome, thanks.

“array or tree of arrays” sounds very general, probably even better than an old fashioned 2D array.