RFC: SciPy array types & libraries support

[rgommers]

This is a long read and a proposal with quite a large scope and a decent amount of backwards-compatibility impact. I hope it'll make SciPy's behavior much more consistent and predictable, as well as yield significant performance gains. I'll post this on the mailing list too. I'd suggest to discuss the big picture first; will ask for process suggestions on the mailing list too in case we'd like to split up the discussion, open a doc/proposal PR for detailed review, or anything like that.

The basic design principle I propose aiming for in all of SciPy for array libraries is: container type in == container type out.

Python sequences (lists, tuples, etc.) will continue to be converted to NumPy arrays, as will other unknown types that are coercible with np.asarray.

The scope of this proposal is: how to treat array inputs. This includes different kinds of NumPy arrays and ndarrays with particular dtypes, as well as array/tensor objects coming from other libraries.

Out of scope for this proposal are (a) dataframe libraries, and (b) implementation of a dispatch mechanism for when non-numpy array types hit C/C++/Cython/Fortran code inside SciPy. Both of these topics are touched upon in the Appendix section.

I'll dive straight into the design here; for context on why we'd want/need this design or what has been done and discussed before, see the Context and Problems & opportunities sections further down.

array/tensor types support

Array types and how to treat them:

Input type	Return type	When	Behavior	Notes
numpy.ndarray	numpy.ndarray	always
torch.Tensor (CPU)	torch.Tensor	always
torch.Tensor (GPU)	torch.Tensor	when possible	raise TypeError otherwise	for pure Python implementation for now. See text below for more details
cupy.ndarray	cupy.ndarray	when possible	raise TypeError otherwise
has __array_namespace__	same as input type	when possible	raise TypeError otherwise
array with object dtype	n/a		raise TypeError
numpy.matrix	n/a	for NumPy >= 2.0 (or always?)	raise TypeError
numpy.ma.MaskedArray	numpy.ma.MaskedArray	only for stats.mstats	raise TypeError otherwise	future nan_policy/mstats plans are unaffected by this
numpy.ndarray subclasses	same subclass type	always		so always use np.asanyarray checks
torch.Tensor subclasses	same subclass type	same as for torch.Tensor on CPU/GPU
PyData/Sparse arrays	same sparse array type	once it has array API support	TypeError otherwise	TypeError is already the current state; no auto-densification
dask.Array	dask.Array	once it has array API support	TypeError otherwise (?)	bc-breaking, but mixing with current conversion to numpy array may not be healthy
jax.numpy.Array	jax.numpy.Array	once added to array-api-compat	same as PyTorch for CPU and GPU

When a non-NumPy array type sees compiled code in SciPy (which tends to use the NumPy C API), we have a couple of options:

1. dispatch back to the other library (PyTorch, CuPy, etc.).
2. convert to a NumPy array when possible (i.e., on CPU via the buffer protocol, DLPack, or __array__), use the compiled code in question, then convert back.
3. don't support this at all

I'll note that (1) is the long-term goal; how to implement this is out of scope for this proposal - for more on that, see the Appendix section. For now we choose to do (2) when possible, and (3) otherwise. Switching from that approach to (1) in the future will be backwards-compatible.

A note on numpy.matrix: the only place this is needed is scipy.sparse, it can be vendored there and for NumPy >= 2.0 instances of the vendored matrix code can be returned. That allows deprecating it in NumPy. We need to support scipy.sparse.*_matrix long-term for backwards compatibility (they're too widely used to deprecate), however for new code we have sparse.*_array instances and PyData/Sparse.

Regarding array API support: when it's present in an array library, SciPy will require it to be complete (v2022.12), including complex number support and the linalg and fft submodules - supporting partial implementations without linalg support in particular seems unnecessary.

For as-yet-unsupported GPU execution when hitting compiled code, we will raise exceptions. The alternative considered was to transfer to CPU, execute, and transfer back (e.g., for PyTorch). A pro of doing that would be that everything works, and there may still be performance gains. A con is that it silently does device transfers, usually not a good idea. On balance, something like this is only a good idea if there's a well-defined plan to make GPU execution work for most functionality on a reasonable time scale (~12-18 months max). Which means both addressing the dispatch (uarray & co) problem that I am trying hard to avoid diving into here, and having time commitments for doing the work. Since we don't have that, raising exceptions is the way to go.

Development and introduction strategy

The following strategy is meant to allow implementing support in a gradual way, and have control over when to default to a change in behavior - which necessarily is going to have some backwards compatibility impact.

• develop behind an experimental flag (an environment variable and a global setting starting with an underscore)
• copy the scikit-learn approach, using array-api-compat as an optional dependency for now
• switch the new behavior of "input type in == input type out" by default once a submodule has complete support
  • at that point, vendor array-api-compat
  • may be a major release if it can all be done in one release cycle. otherwise do submodule by submodule
• explicit support can and should be tested in CI for: PyTorch CPU tensors, numpy.array_api for compliance testing in APIs that start supporting array API compliant input, and perhaps pandas Series/DataFrame. GPU CI may materialize at some point, but let's not count on or plan for that.
• implement new type coercion machinery centrally in scipy._lib and depend on that from other submodules. That means we'll have a single replacement for np.asarray/np.asanyarray in one central location.

Context

• For NumPy 2.0 we should have full array API standard support in the main numpy namespace
• The numpy.array_api module is useful for testing that code is actually using only standard APIs, because it's a minimal implementation. It's not going to be used as a layer for production usage though, it's for portability testing only.
• The array-api-compat library provides the compatibility we need right now between NumPy, CuPy and PyTorch. Other libraries can be added there too. This bridges any gaps between the standard and real-world implementations. The package is light-weight and can be depended on - but is also designed to be vendored to avoid a new runtime dependency.
• scikit-learn has had experimental array API support for a little while, and has now settled on array-api-compat as the way to go. This looks quite a bit nicer than the previous approach with numpy.array_api, and the work to support PyTorch that way has been quite well-received.
• As for the myriad of open NEPs on interoperability topics:
  • NEP 30 (__duckarray__), NEP 31 (uarray) and NEP 37 (__array_module__) are all effectively dead - I'll propose rejecting these.
  • NEP 47 (array API support via numpy.array_api) has been implemented, however I'm going to propose marking it superceded via a new NEP for the main numpy namespace
  • __array_function__ and __array_ufunc__ are fully implemented and will continue to exist and be supported in NumPy. We won't support those mechanisms in SciPy though, since we are coercing unknown input types to ndarray and error out if that fails. The exception here is ufuncs in scipy.special, which happen to work already because we're reusing the numpy ufunc machinery there. We can probably leave that as is, since it's not problematic.
• NumPy array subclasses are a long-standing problem, with masked arrays and np.matrix instances not being Liskov-substitutable (i.e. it's a drop-in replacement, no changes in behavior but only extensions), making it difficult to accept them. We can explicitly start rejecting those with clear exceptions, that will make regular subclasses a lot more useful.
• Missing data support is a broad and tricky subject:
  • numpy.ma has tons of issues and isn't well-maintained. There's a full rewrite floating around that is ~90% complete with some luck will make it into NumPy 2.0. However, it hasn't seen movement for several years, and the work on that is not planned. numpy.ma in its current form should be considered legacy.
  • scipy.stats.mstats is the only module that specifically supports masked arrays.
  • scipy.stats has a nan_policy keyword in many functions that is well-maintained at this point, and has a documented specification. That is probably not applicable to other submodules though.
  • Support for missing data in other places like interpolate and spatial.KDTree (see gh-18230) may make sense, but ideally that'd use solid numpy support (e.g., via a null-aware dtype) and that does not exist.
  • Dataframe libraries have much better support, but it's a difficult story to use non-numpy-backed dataframes at all.
  • For the purposes of this proposal, I'd prefer leaving it at "reject numpy.ma.MaskedArray instances".

Problems & opportunities

The motivation for all the effort on interoperability is because the current state of SciPy's behavior causes issues and because there's an opportunity/desire to gain a lot of performance by using other libraries (e.g., PyTorch, CuPy, JAX).

Problems include:

• pandas input with non-numpy dtypes being converted to object arrays and then support being very hit-and-miss; see for example gh-18254
• general pandas dataframe/series support being patchy; see for example gh-9616 and gh-12164
• masked arrays giving incorrect results because the mask is silently dropped; see, e.g., gh-1682 and gh-12898
• sparse array support being patchy; agreement that auto-densification needs to be avoided, but that means current behavior is inconsistent (xref gh-4239)
• object arrays, including use of Decimal and other random things that people stuff into object arrays on purpose or (more likely) by accident, being handled very inconsistently.

Opportunities include:

• Large performance gains. See for example mdhaber/scipy#63 (10x-50x for stats bootstrapping functionality with CuPy) and scikit-learn#25956 which shows how PyTorch improves over NumPy (2x-7x on CPU, 60x-100x on GPU)
• interest in support for Dask which we haven't really been able to do much with, see for example gh-10204, gh-8810 and gh-10087. Xarray is in the same boat, it builds on top of Dask (see, e.g., gh-14824)

References

• Scikit-learn PR for PyTorch support: scikit-learn#25956
• Scikit-learn docs for array API support: https://scikit-learn.org/stable/modules/array_api.html
• Tracking issue for SciPy array API standard support: gh-15354
• Initial PR for SciPy array API support (not even blocked, we just didn't push it forward because numpy.array_api is a bit cumbersome): gh-15395
• NumPy 2.0 developer meeting presentation on array API adoption

Appendix

Note: these topics are included for reference/context because they are related, but they are out of scope for this proposal. Please avoid diving into these (I suggest to ping me directly first in case you do see a reason to discuss these topics).

dataframe support

For dataframe library support, the situation is a little trickier. We have to think about pandas Series and DataFrame instances with numpy and non-numpy dtypes, the presence of nullable integers, and other dataframe types which may be completely backed by Apache Arrow, another non-numpy library (e.g., cuDF & CuPy), or have implemented things completely within their own library and may or may not have any NumPy compatibility.

This would be one option, which is somewhat similar to what scikit-learn does (except, it converts nullable integers to float64 with nans):

Input type	Return type	When	Behavior	Notes
pandas.Series	pandas.Series	if dtype is a numpy dtype	raise TypeError otherwise	non-NumPy dtypes get coerced to object arrays otherwise, avoid that
pandas.DataFrame	pandas.DataFrame	if dtype is a numpy dtype	raise TypeError otherwise
has __dataframe__	same as input type			needs a small extension to dataframe interchange protocol to reconstruct input type

There's a lot to learn from the effort scikit-learn went through to support pandas dataframes better. See, e.g., the scikit-learn 1.2 release highlights showing the set_output feature to request pandas dataframe as the return type.

Note: I'd like to not work this out in lots of detail here, because it will require time and that should not block progress on array library support. I just want to put it on the radar, because we do need to deal with it at some point; current treatment of dataframes is quite patchy.

Dispatching mechanism

For compiled code, other array types (whether CPU, GPU or distributed) are likely not going to work at all; the SciPy code is written for the NumPy C API. It's not impossible that some Cython code will work with other array types if those support the buffer protocol and the Cython code uses memoryviews - but that's uncommon (won't work at all on GPU, and PyTorch doesn't support the buffer protocol on CPU either).

There has been a lot of discussion on how this should work. The leading candidate is Uarray, which we already use in scipy.fft (as do matching fft APIs in CuPy and PyFFTW) and has other PRs pending in both SciPy and CuPy. However, there is also resistance to that because it involves a lot of complexity - perhaps too much. So significant work is needed to simplify that. Or switch to another mechanism. This is important work that has to be done, but I'd prefer not to mix that with this proposal.

Whatever the mechanism, it should work transparently such that scipy.xxx.yyy(x) where x is a non-numpy array should dispatch to the library which implements the array type of x.

We have a uarray label in the issue tracker. See gh-14353 for the tracker with completed and pending implementations. For more context and real-world examples, see:

• https://labs.quansight.org/blog/pydata-extensibility-vision (and in particular the PyData dispatching system section which covers alternative dispatching options)
• https://labs.quansight.org/blog/making-gpus-accessible-to-pydata-ecosystem-via-array-api
• https://quansight-labs.github.io/array-api-demo/GW_Demo_Array_API.html
• https://data-apis.org/array-api/latest/use_cases.html
• https://github.com/numpy/archive/blob/main/other_meetings/2020-04-21-array-protocols-discussion-notes.md
• Note that there's work to be done to write SPEC-0002 on this dispatching topic, see also Discourse posts on that: here, here and here.
• Scikit-learn is working on a different kind of system for this, with explicit plugins to extend scikit-learn with support for other libraries/devices where it's using compiled code - see Path for pluggable low-level computational routines scikit-learn/scikit-learn#22438.
• PyTorch has a more robust backend system (https://pytorch.org/docs/stable/dynamo/custom-backends.html).
• I'm sure there more examples of how to enable this kind of thing, inside and outside of Python. It's just a lot of engineering effort.

[tylerjereddy]

I read it over, sounds exciting. I have to admit that evaluating API designs on this scale is not one of my strengths. I'll briefly mention that the US National Labs will probably quietly test out their own array API compatible technologies behind the scenes as well, if this moves forward, to see if performance gains are observed on bleeding-edge supercomputer hardware. One question I have related to this--for nodes with multiples devices (like 4 x A100 is pretty common these days), do we leave the use of multiple devices up to the GPU library under the hood? Do we explicitly not encourage using kwargs for controlling concurrency that is specific to say GPU libs (i.e., too much risk of optional kwarg API bloat)?

You did briefly mention testing a few different CPU-based scenarios by default. I'm curious if we'd use a pytest decorator/magic to auto-expand eligible tests and/or if it might be helpful to be able to have i.e., a testing env variable so that we can request particular "backends" for testing (for example, could be useful to be able to first iteratively develop/test CPU-only, then see what kind of mess arises when you switch to testing with GPU locally/on a cluster node, etc.). That's perhaps more of a detail.

One thing that may be more general re: testing is how we might handle bug reports that are effectively isolated to a subset of pass-throughs--only library X with hardware Y. I suppose the array API is supposed to prevent that from happening, but there is also reality... I guess maybe occasional shims and regression tests that only run with the conditions are right, assuming that we can't just consider these upstream issues most of the time. I will say that I don't see any way for this overall improvement to make our testing situation any simpler, unfortunately/obviously.

Maybe after the first few PRs go in, it will be easier for those of us who aren't familiar with every corner of those NEPs to start contributing since we might be able to understand the changes needed as basic engineering tasks.

I wonder if spatial will unfortunately be one of the least-beneffiting modules--some of our most-used functionality relies heavily on carefully-written compiled backends, like Qhull, the C++ stuff driving KDTree, and the C++ backing much of the distance calculation code. I believe there are a few "Python/NumPy land" corners at least.

[WarrenWeckesser]

What is the policy for handling mixed inputs when a function accepts more than one argument? This is probably covered in one of the array API documents or discussions, but I think it we need it stated explicitly here, too. I expect the most common use-case is mixing $SOME_OTHER_ARRAY_TYPE and standard numpy arrays (e.g. multiplying a sparse matrix or sparse array by a numpy array, convolving some other array type with a kernel that is a numpy array, etc.).

[asmeurer]

What is the policy for handling mixed inputs when a function accepts more than one argument? This is probably covered in one of the array API documents or discussions, but I think it we need it stated explicitly here, too. I expect the most common use-case is mixing $SOME_OTHER_ARRAY_TYPE and standard numpy arrays (e.g. multiplying a sparse matrix or sparse array by a numpy array, convolving some other array type with a kernel that is a numpy array, etc.).

The array API spec says:

Non-goals for the API standard include:
...

• Making it possible to mix multiple array libraries in function calls.

  Most array libraries do not know about other libraries, and the functions they implement may try to convert “foreign” input, or raise an exception. This behaviour is hard to specify; ensuring only a single array type is used is best left to the end user.

https://data-apis.org/array-api/latest/purpose_and_scope.html#out-of-scope

See also the interchange section about DLPack https://data-apis.org/array-api/latest/design_topics/data_interchange.html

[AnirudhDagar]

Thanks Ralf, for the detailed RFC; after reading it a couple of times, I have a few implementation-specific comments/questions/clarifications (some arising from the past work on the Gravitational waves demo).

It seems array-api-compat library will be key.

This bridges any gaps between the standard and real-world implementations.

With my understanding, I'd try to expand the above for everyone with a real-world example taken from scipy.signal and PyTorch (note: PyTorch is just an example, it may apply to any other array/tensor lib for some other case). Specifically, the lines below are used in scipy.signal.welch:

scipy/scipy/signal/_spectral_py.py

Lines 1748 to 1752 in c8206b8

	if not same_data:
	y = np.asarray(y)
	outdtype = np.result_type(x, y, np.complex64)
	else:
	outdtype = np.result_type(x, np.complex64)

Here np.result_type is used; consider the scenario where scipy.signal.welch is to be made Array API compatible with, say torch tensors (which we tried in this prototype demo). It would require the behaviour of torch.result_type to follow the Array API spec. But the reality is that torch.result_type differs from how result_type is defined in array API. Difference being that, as of today, PyTorch only allows two tensors, and in contrast, array API specifies arbitrary number of arguments being either tensor and/or dtypes.

Now I'd expect at some point, PyTorch will update their implementation to support array API for these functions, which are still left, and there is probably effort around it already (see this result_type PR pytorch/pytorch#61168 ).

But for now, array-api-compat is a brilliant intermediary to solve such problems. See how array-api-compat fixes this specific compatibility issue below by wrapping around the current torch.result_type signature and behaviour:

https://github.com/data-apis/array-api-compat/blob/9ef7f72fa91ddc4561148f906d386d08e65def7d/array_api_compat/torch/_aliases.py#L103

---

Questions:

1. Are we considering cases in SciPy which at first might not seem portable with Array API due to the usage of NumPy methods that are outside the scope of Array API, but with slight refactoring (using an alternate NumPy method and making sure no performance regression in default behaviour), we might bring that particular SciPy method/class closer to Array API Spec or just closer to other array/tensor APIs? One close example, but not exactly, is this (MAINT: Replace np.rollaxis with np.moveaxis #14491). So basically, would it be okay to slightly refactor the current SciPy implementation of a function using some NumPy method to some other NumPy method for the sake of achieving interoperability through this work? Of course, BC is still maintained.

2. array-api-compat adds an extension on top of the Array API Spec. I'd like to know what happens to cases (if any, we will only find more such cases during this work) that are still outside array-api-compat. For example, NumPy and other array/tensor library, for a particular function, might have the same signature and behaviour. If such np.<some_method> (which is still out of the scope of Array API) is used within a particular SciPy method/class; will that SciPy method still be marked as incompatible with this work and be completely ignored for now? Even though it would be technically possible to keep interoperability in such cases. One example is np.iscomplex or cupy.iscomplex, (see use in https://github.com/scipy/scipy/blob/main/scipy/signal/_signaltools.py#L4546). I guess there will be other use cases of such functions, which might block the compatibility of a whole SciPy module/method just because that one NumPy method is still outside of Array API even though, it had an intersection with other array/tensor libraries.

3. Is there already an identified list of targeted submodules or methods in SciPy that can be made to work with array-api-compat and that we will focus on initially? What if a submodule is largely not compatible with Array API (and needs extra infra like dispatching with uarray etc.), but maybe a few methods/classes inside the submodule are compatible? Do we make only those methods/classes compatible in the largely incompatible submodule or completely ignore the whole submodule in such a case?

Apologies if my questions are a bit too specific with implementation details and are special examples, and if this RFC is about discussing the high-level design, feel free to ignore it.

Overall, excited that all of this work is coming together, and the current state looks pretty solid to me to start baking in this functionality in the "array consumer libraries" (scipy, sklearnetc.). Thanks for setting up the foundation tools like array-api-compat!

[rgommers]

Thanks for the feedback so far! Let me try to answer the high-level questions first.

for nodes with multiples devices (like 4 x A100 is pretty common these days), do we leave the use of multiple devices up to the GPU library under the hood? Do we explicitly not encourage using kwargs for controlling concurrency that is specific to say GPU libs

Yes, we should not carry anything specific to the execution model within SciPy code. That's a key design principle of the array API standard. You could use a multithreaded CPU library, a GPU library, or something even fancier like what you have here. The code should be unchanged. There's quite a bit of variation between different distributed libraries in particular, and that's out of scope to deal with I'd say.

I'm curious if we'd use a pytest decorator/magic to auto-expand eligible tests

Yes, that's what I'd like to aim for. This has worked pretty well for the array API test suite; a single switch allows you to select any compatible libraries. I imagine we'd have that as a decorator for relevant tests, or some such thing.

I wonder if spatial will unfortunately be one of the least-beneffiting modules

Probably. I imagine modules like stats and signal to be prime candidates early on.

What is the policy for handling mixed inputs when a function accepts more than one argument? This is probably covered in one of the array API documents or discussions, but I think it we need it stated explicitly here, too.

Good point, I'll add it to the issue description. @asmeurer gave a good answer already. I'd add that numpy arrays + sequences will work unchanged; non-numpy arrays + sequences should error; and in general it's a bad idea to mix array objects coming from different libraries. Rather get a clear exception there and have the user convert one kind of array to the other, to get the behavior they intended. That's how all libraries except NumPy work, and what NumPy does is an endless source of bugs and user confusion.

I expect the most common use-case is mixing $SOME_OTHER_ARRAY_TYPE and standard numpy arrays (e.g. multiplying a sparse matrix or sparse array by a numpy array, convolving some other array type with a kernel that is a numpy array, etc.).

The sparse-dense mix is the only one that in principle makes sense. I think that could be allowed in the future, but I'm not exactly sure under what circumstances. Right now all sparse array types forbid auto-densification, so where you can mix sparse-dense is quite restricted. PyTorch has the best story here, because it natively provides both dense and sparse tensors - but even there it's still a work in progress.

[rgommers]

So basically, would it be okay to slightly refactor the current SciPy implementation of a function using some NumPy method to some other NumPy method for the sake of achieving interoperability through this work?

Yes indeed, that is expected in many cases I'd say.

2. If such np.<some_method> (which is still out of the scope of Array API) is used within a particular SciPy method/class; will that SciPy method still be marked as incompatible with this work and be completely ignored for now?

I'd say "it depends". If it's simple to reimplement np.xxx with other functions that are part of the standard, then that can be done (special-casing numpy arrays for performance even if that makes a big difference). If it's not simple to do that, then we have a gap. If it's a big gap (hopefully won't happen often) then we should propose it for inclusion in the standard, and once that proposal is accepted implement it in array-api-compat. If it's niche functionality, then we have to deal with that on a case-by-case basis.

3. Is there already an identified list of targeted submodules or methods in SciPy that can be made to work with array-api-compat and that we will focus on initially? What if a submodule is largely not compatible with Array API (and needs extra infra like dispatching with uarray etc.), but maybe a few methods/classes inside the submodule are compatible? Do we make only those methods/classes compatible in the largely incompatible submodule or completely ignore the whole submodule in such a case?

I'd prefer for coherent sets of functionality to be upgraded at once. To stay with Tyler's example of spatial: if 90% of functionality will only work with NumPy, then there's not much point using the array API standard. In those cases we can still use the central infrastructure I proposed for input validation though (e.g., reject object and masked arrays) to gain consistency - there's nothing blocking that.

[seberg]

That's how all libraries except NumPy work, and what NumPy does is an endless source of bugs and user confusion.

Thanks for the proposal! I doubt its urgent to figure out (see end), but I need to push back a little on promotion/mixing story here.

The big issues are really not about mixing/promotion. They are about NumPy thinking that:

• any arbitrary object that doesn't shout "no" should be promoted to a NumPy array. (Even nested lists/sequences that contain objects which do shout "no". That last part could be something to allow changing even 🤔)
• we/libraries often call np.asarray(), which coerces and doesn't give proper promotion a chance to begin with (but you would stop doing this here!). (Proper promotion would mean e.g. JAX decides, not NumPy. The __duckarray__ proposal would actually fix it, and maybe that could even still relevant in some form as a simple improvement for __array_function__. But that is not important here.)

Unlike NumPy most array implementers do not pretend that its OK to promote almost everything to their array type.

But promotion/mixing is a different thing. Many do promote: JAX (as far as I can tell), dask.array, sparse arrays, masked arrays, quantities all promote/mix great with NumPy arrays. Also subclasses with their superclass, but if they don't customize the super class namespace, they are actually supported in the current state.

sparse/dense

I admit sparse/dense may be trickier. But I can see that for many things it either makes perfect sense or can be trusted to raise a good error. Some ops densify, some ops don't, and those that where it's not clear raise.
The main way this could end up problematic is if subtle code changes might somehow end up changing the result type. I doubt that is the case, but I don't know. My feeling is that this is something that sparse arrays need to figure out.

Is supporting limited promotion important enough to be needed? I don't know. I tend to think the answer is yes and I somewhat suspect e.g. Dask users may well be annoyed if its missing (I am assuming Dask has a good reason to support it very well, but trying to gauge it better. Would post the result on the array-api issue if there is a clearer one).

Overall my angle/answer would be: Fortunately, it should not be a big mess to add simple promotion/mixing support here. (I am confident about that, but I think its important to clarify since it is not obvious.)
So probably it need not be an urgent thing to figure this out. But it is an interesting question that I am sure will come up again and I expect it may well prove important enough sooner or later.

[rgommers]

Thanks for your thoughts Sebastian!

(Even nested lists/sequences that contain objects which do shout "no". That last part could be something to allow changing even thinking)

yes please:) I think our NumPy 2.0 plate is pretty full already, so I'm hesitant to invest a lot of time now - but if there's easy wins there to avoid promote-to-object-array issues, let's try and take them.

Proper promotion would mean e.g. JAX decides, not NumPy

This is not really "promotion" in any proper sense. JAX and NumPy arrays are two different types of arrays that are unrelated by class hierarchy. There happens to be an ordering though, through __array_priority__. I have to admit I'm not certain how future-proof that is. Best practice is clearly to explicitly convert to the desired kind here, mixing them is not great. But it may work. What my proposal is trying to say is to reject things that aren't related, rather than blindly stuffing them into np.asarray. It should be more precise in what relations are allowed/supported - I had not thought of __array_priority__.

I've worked on quite a few PyTorch issues with the same problem. PyTorch tensors and NumPy arrays are unrelated, so things like x + t working but t + x erroring out (or vice versa) is quite confusing. It should just error.

But promotion/mixing is a different thing. Many do promote: JAX (as far as I can tell), dask.array, sparse arrays, masked arrays, quantities all promote/mix great with NumPy arrays.

Again, only when there is a clear relation. scipy.sparse and pydata/sparse mostly errors out when mixed with ndarray. Masked arrays do work because they're a subclass (modulo silently dropping the mask of course). Some quantities objects work well with numpy because they're a subclass (see, e.g. AstroPy's Quantity) or because they use things like __array_function__. Those cases will continue to work. In the absence of that well-defined relationship, things will not work. So you cannot and should not mix Quantity with JAX/PyTorch/CuPy/sparse.

Overall my angle/answer would be: Fortunately, it should not be a big mess to add simple promotion/mixing support here.

If this is well-defined, I agree. We just need to be very explicit what the mechanisms are, and they should be introspectable. I want to avoid "just try it and see what rolls out, YMMV".

[seberg]

Collapsing everying into details, because... It doesn't matter much for the RFC, but I can't leave some of these things unanswered.

(We could get hung up and say I should use "implicit conversion" like the C-standard making the argument about implicit/explicit conversion – but Julia uses "promotion" the way I do as far as I understand, so I will change my habit when someone asks me to or points out why its very wrong.)

but if there's easy wins there to avoid promote-to-object-array issues, let's try and take them.

This is not about promote-to-object-array issues. Its about JAX sayings its higher priority then NumPy, but NumPy ignoring that when you write: numpy_arr + [jax_array]. Maybe also about a way to write np.asarray(jax_array) that raises an error because JAX says its the higher priority.

This is not really "promotion" in any proper sense

Aha? JAX and NumPy arrays are two different types of arrays that are unrelated by class hierarchy. Just as: int and float are two different types of numbers that are unrelated by class hierarchy.

Unless you also say 3 + 3.0 = 6.0 is not promotion, in that case, I agree.

So yes, there is an (abstract) base class (array and number). There is a clear expected hierarchy everyone agrees on and expects (JAX > NumPy, Dask > NumPy and float > int). And neither of these conversions when mixing the two are perfectly "safe" in some sense, but implicit int -> float conversions are still universally useful and numpy -> dask.array/JAX/masked array implicit conversions seem also useful (otherwise Dask would not do it, no?).

Sure "array" is a bit of a more complicated then "number".

There happens to be an ordering though, through array_priority

__array_priority__ is more of a remnant. The real thing is __array_ufunc__ which can just be None to tell NumPy: "No". It isn't relevant here.

I've worked on quite a few PyTorch issues with the same problem. PyTorch tensors and NumPy arrays are unrelated, so things like x + t working but t + x erroring out (or vice versa) is quite confusing. It should just error.

Fine its confusing and probably has things that need improving. But also __array_ufunc__ = None should guarantee you an error. Yeah, it seems like it doesn't end up with a great error but you could probably improve on that in your __add__ or by implementing a dummy __array_ufunc__.

But I don't even think that matters. If a type is comfortable with implementing __array_ufunc__ such as Dask (or even just the operators like JAX), then that type can mix with NumPy arrays fine in practice (not with np.<func>, but we know we won't use those here).

And there is no problem with allowing such a type to signal that it wants to do just that. That still allows Torch to signal that it doesn't want to mix!

So is it annoyingly hard to tell NumPy "no"? Yes, and maybe that needs to be improved. But frankly, I think it is misleading to use those difficulties as an argument that mixing is generally pretty much doomed to fail.

We have many examples of successfull implementations that mix fine. Like dask, masked arrays, or quantities...

What my proposal is trying to say is to reject things that aren't related, rather than blindly stuffing them into np.asarray

But because you stop stuffing things blindly into np.asarray() you suddenly can start reasoning about their relation. The problem is that right now we break any relation by explicit conversion via np.asarray().

Take np.asarray() out of the equation, and Dask and NumPy already have a well defined relation and mix great!

Masked arrays do work because they're a subclass

No, I would be very surprised if Allan's masked array is a subclass and I am sure it has better behavior than np.ma and mixes great with NumPy arrays.

The important point is that it wants to work. That can just as well be by implementing __array_ufunc__ (or even only operators). Allowing interop doesn't mean that np.func(my_object) needs to work (although mylib.func(numpy_array) must).

If this is well-defined, I agree. We just need to be very explicit what the mechanisms are, and they should be introspectable. I want to avoid "just try it and see what rolls out, YMMV".

This is a reason, unlike the "it's a mess" argument, IMO. We do have quite a bit of prior art, since every other protocol defines a mechanism. It still needs care of course (there are small variations after all).
For me the main point is conservative guidance similar and additionally to NEP 18, though.

But at that point it's: We still need to figure that out because we don't want to get it wrong. And not "its a bad idea".

[rgommers]

But at that point it's: We still need to figure that out because we don't want to get it wrong. And not "its a bad idea".

Maybe we're not disagreeing all that much. I said briefly "it's a bad idea" because (a) it's apparently so complex that after all these years we still haven't figured it out indeed, and (b) the gains are in many cases a minor amount of characters saved to deal with explicit conversions in the user's code instead. So my personal assessment is that this isn't worth the effort. Others may disagree of course, and spend that effort.

I think the appropriate place for discussing this topic further is either on the numpy repo or on pydata/duck-array-discussion, which was/is an effort specifically to sort out this hierarchy.

It doesn't matter much for the RFC

I agree. We're not closing any doors in this RFC, it can in principle be enabled in a future extension.