[Tracker] Implement all scipy.* APIs in CuPy

[kmaehashi]

Implement GPU version of scipy.* functions in cupyx.scipy.* namespace.

This is a tracker issue that summarizes the implementation status of each SciPy public module in CuPy. See the comparison table for details.

Legends

• 🥚 Module not yet available in CuPy.
• 🐥 Module available in CuPy but coverage needs improvement.
• Module available in CuPy and major features are implemented (may still have a few unimplemented APIs).
• ❌ Out of scope.

List of Modules

• 🥚 Clustering package (scipy.cluster)
  • 🥚 K-means clustering and vector quantization (scipy.cluster.vq) (Equivalent of scipy.cluster.vq.kmeans #5947)
  • 🥚 Hierarchical clustering (scipy.cluster.hierarchy) (Equivalent of scipy.cluster.hierarchy.linkage #3434)
• 🥚 Constants (scipy.constants)
• Discrete Fourier transforms (scipy.fft)
• Legacy discrete Fourier transforms (scipy.fftpack) (note: deprecated, no further development planned)
• 🥚 Integration and ODEs (scipy.integrate) (Numerical Differential Equation Integration #7019)
• Interpolation (scipy.interpolate) (Tracker: cupyx.scipy.interpolate #7186)
• 🥚 Input and output (scipy.io)
  • 🥚 scipy.io.arff
  • 🥚 scipy.io.matlab
  • 🥚 scipy.io.wavfile
• 🐥 Linear algebra (scipy.linalg)
  • 🥚 Low-level BLAS functions (scipy.linalg.blas)
  • ❌ scipy.linalg.cython_blas
  • 🥚 Low-level LAPACK functions (scipy.linalg.lapack) c.f. cupyx.lapack
  • ❌ scipy.linalg.cython_lapack
  • 🥚 Interpolative matrix decomposition (scipy.linalg.interpolative)
• Multidimensional image processing (scipy.ndimage)
• 🥚 Orthogonal distance regression (scipy.odr)
• 🥚 Optimization and root finding (scipy.optimize) (Any plans for scipy.optimize? #6112)
• Signal processing (scipy.signal) (Tracker: cupyx.scipy.signal #7403)
  • scipy.signal.windows (Tracker: cupyx.scipy.signal.windows #7404)
• Sparse matrices (scipy.sparse)
  • Sparse linear algebra (scipy.sparse.linalg)
  • 🐥 Compressed sparse graph routines (scipy.sparse.csgraph) (Compressed Sparse Graph Routines (scipy.sparse.csgraph) #2431)
• Spatial algorithms and data structures (scipy.spatial) (cupyx.scipy.spatial and cupyx.scipy.spatial.distance packages #5946)
  • 🐥 Distance computations (scipy.spatial.distance) (cupyx.scipy.spatial and cupyx.scipy.spatial.distance packages #5946)
  • 🥚 Spatial transformations (scipy.spatial.transform)
• 🐥 Special functions (scipy.special) (note: planned to be backed by SciPy in future Add complex support for the digamma function #8163)
• 🐥 Statistical functions (scipy.stats)
  • 🥚 Contingency table functions (scipy.stats.contingency) (note: RAFT could be used)
  • 🥚 scipy.stats.distributions
  • 🥚 Statistical functions for masked arrays (scipy.stats.mstats) (note: masked arrays unsupported in CuPy)
  • 🥚 Quasi-Monte Carlo submodule (scipy.stats.qmc)
  • 🥚 Random Number Generators (scipy.stats.sampling) (note: RAFT could be used)

Note: Several modules in SciPy such as scipy.cluster.hierarchy accepts CuPy ndarrays as inputs through the Array API standard.

Starter Task

If you are new to CuPy and wondering where to start, here are some good first things to try. They are relatively simple and independent.

• APIs under Statistical Functions (scipy.stats)
• scipy.sparse.block_diag (cupyx.scipy.sparse.block_diag working with cupy sparse matrix on GPU #7058)
• scipy.sparse.load_npz
• scipy.sparse.save_npz
• scipy.constants.* (hint: most things except for methods can just be borrowed from SciPy)

Steps to Contribute

1. Fork and star ⭐ the CuPy repository 😉

2. Pick a function you want to work on from any of the modules listed above. You can find the function in the SciPy API Reference to understand what should be implemented. If you want to work on new modules (🥚), first discuss with us on this issue.

3. Implement a function in your branch. If you need help, join Gitter or just ask for help in this issue.

4. Don't forget to write a test code!

5. Build CuPy and run tests to confirm that the function runs fine:
   pip install -e . && pytest tests/cupyx_tests/scipy_tests/PATH_TO_YOUR_TEST
   See the Contribution Guide for details.

6. Submit a pull-request to the main branch.

Note that you will need a GPU environment to develop CuPy.

See also: NumPy API Tracker Issue #6078

[leofang]

cc: @cjnolet @quasiben

[rgommers]

This looks great, thanks @kmaehashi!

Miscellaneous routines (scipy.misc)

Note that we're quite likely deprecating this module in SciPy 1.9.0 (see scipy/scipy#15608), so I suggest not working on that.

[kmaehashi]

Thanks for the info, @rgommers! Updated the list.

[cjnolet]

@kmaehashi, RAFT contains hierarchical clustering code as well: #3434.

RAFT also contains functions for computing contingency table and sampling from various distributions. Wondering if those might be useful here as well.

[kmaehashi]

Thanks @cjnolet! Added a link to the issue.

[Pranavchiku]

This is brilliant and interesting! I am Pranav, developer @lfortran, we are trying to compile scipy using LFortran, we are very close, almost got scipy.special.specfun working. I wish to ask, if we (cupy contributors) are translating fortran implementation present in SciPy to make it use cupy apis, right?

[kmaehashi]

Hi @Pranavchiku! Yes, we are rewriting SciPy code into CUDA C/C++. You'll get the idea by exploring our code base under cupyx/scipy:
https://github.com/cupy/cupy/tree/main/cupyx/scipy/special

[Pranavchiku]

Note that you will need a GPU environment to develop CuPy.

Do we need explicit GPU environment, or integrated GPUs work? By integrated I mean Apple GPU that is available with mac m1, m2, etc. chips.

[kmaehashi]

Discrete GPUs are required, currently we support CUDA and HIP (for AMD).
https://docs.cupy.dev/en/latest/overview.html

[Zhyrek]

Hey, I need non-batched BLAS functions for what I'm working on, so I will be working on implementing those in the near future. Just wondering if there are any thoughts on a best-practices for implementing these.

Roughly how I imagine using these functions (collapsible code section):

@jit.rawkernel(device=True)
def square(a):
  matrix2 = cupyx.jit.local_memory(float, (2, 2)) #I'll also try implementing multi-dimensional local/shared arrays
  matrix2[0][0] = a
  matrix2[1][1] = a
  return cupyx.scipy.linalg.det(matrix2) #can be used in device functions

@jit.rawkernel()
def elementwise_2_x_squared(x, size): #computes 2*x^2, toy example which uses linalg functions
  tid = jit.grid(1)
  ntid = jit.gridsize(1)

  matrix = cupyx.jit.local_memory(float, (2, 2))

  for i in range(tid, size, ntid):
      stored_val = x[i]
      x[i] = square(stored_val)
      matrix[0][0] = stored_val
      matrix[1][1] = stored_val
      x[i] += cupyx.scipy.linalg.det(matrix) #can be used inside raw kernel

c = cp.array([1., 2.])
elementwise_2_x_squared((256,), (256,), (c, 2))
print(c) #should return [2., 8.]

The main use case I have is roughly: I have a large array of scalar values (could be hundreds of millions). For each value, I need to generate a (potentially not small) matrix, solve it, and then use the result to compute a scalar output value. For batched code, I would need to store these matrices altogether (potentially hundreds of GB memory - not feasible), while for non-batched, dynamically generated matrices, the memory requirements would just be the scalar values.

[rcremese]

Hello,
I'm quite new to the open-source community but I would like to contribute to this project.
Is anybody on the integration of scipy.constants module ?

[ev-br]

It seems possible to mark scipy.interpolate, scipy.signal and scipy.signal.windows as now?

[ev-br]

Also, I wonder if the comparison table (super-useful BTW!) , https://docs.cupy.dev/en/stable/reference/comparison.html#scipy-cupy-apis is updated automatically on a release, or is it manual?

Would be great to remove SciPy deprecated items (signal.cmplx_sort, windows from scipy.signal namespace, linalg.tri{l,u} etc)

[Pranavchiku]

Hey all, I was always willing to contribute here and I have my first PR #8305 opened that ports scipy.stats.kstat. I do not have access to GPUs hence have not tested locally, I followed boxcox_llf and it was fun and easy :)

[kmaehashi]

It seems possible to mark scipy.interpolate, scipy.signal and scipy.signal.windows as now?

Good catch! Updated the list

Also, I wonder if the comparison table (super-useful BTW!) , https://docs.cupy.dev/en/stable/reference/comparison.html#scipy-cupy-apis is updated automatically on a release, or is it manual?

latest page is updated automatically for each commit merged into main branch. stable page is updated automatically for each stable (v13.*) release.

Would be great to remove SciPy deprecated items (signal.cmplx_sort, windows from scipy.signal namespace, linalg.tri{l,u} etc)

The page is generated with this tool. Deprecated items are currently listed in footnotes section to describe why we don't implement them.

cupy/docs/source/_comparison_generator.py

Lines 206 to 207 in 1aa70f0

	'numpy', 'cupy', klass='ndarray', footnotes={
	'tostring': _deprecated,

[kmaehashi]

Hey all, I was always willing to contribute here and I have my first PR #8305 opened that ports scipy.stats.kstat. I do not have access to GPUs hence have not tested locally, I followed boxcox_llf and it was fun and easy :)

Hi @Pranavchiku, thank you for your interest in contributing! We, however, would like to kindly ask all contributors to build and test the branch before submitting a pull request. I understand that there are situations where access to GPUs is difficult but in that case please consider using offerings such as Google Colab.

[ev-br]

Thanks @kmaehashi A couple of comments, now that I looked at this tracker a bit more:

• With KDTree and Delaunay implemented, maybe scipy.spatial is fair to be marked as , too. Submodules (spatial.transform and spatial.distance might be a separate story, don't have an opinion)
• scipy.linalg.{blas, lapack} are unlikely to be implemented at all, correct?
• even more so for scipy.linalg.cython_{blas,lapack}. Maybe even add an "out of scope" type symbol, ❌ or some such?
• scipy.io --- is there an appetite for ever adding functionality for reading from matlab to CuPy? I suppose using scipy.io followed by cupy.asarray is the way, is it?
• Several SciPy subpackages support CuPy array inputs via Array API, e.g. scipy.cluster. It might be nice to show these too here?

[kmaehashi]

Thanks for the suggestions @ev-br!

• With KDTree and Delaunay implemented, maybe scipy.spatial is fair to be marked as , too. Submodules (spatial.transform and spatial.distance might be a separate story, don't have an opinion)

Absolutely, marked !

• scipy.linalg.{blas, lapack} are unlikely to be implemented at all, correct?

Can't really say they are unlikely at all, as we have non-public (undocumented) cupyx.lapack module. I think we can say they're considered low-priority, though.

• even more so for scipy.linalg.cython_{blas,lapack}. Maybe even add an "out of scope" type symbol, ❌ or some such?

Sounds great, marked ❌.

• scipy.io --- is there an appetite for ever adding functionality for reading from matlab to CuPy? I suppose using scipy.io followed by cupy.asarray is the way, is it?

I began to wonder if I could delete it from the list, as the audience is very limited compared to NumPy's I/O routines (e.g., cupy.load which we implemented by wrapping numpy.load).

• Several SciPy subpackages support CuPy array inputs via Array API, e.g. scipy.cluster. It might be nice to show these too here?

Nice, added a link to the Array API support docs!