Skip to content

DOC: NEP: array in -> array out #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
231 changes: 231 additions & 0 deletions doc/neps/nep-0057-array-in-array-out.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,231 @@
==============================
NEP 57 — Array In -> Array Out
==============================

:Author: Matt Haberland <[email protected]>, Add Your Name Here <[email protected]>
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LMK if you want your name added.

:Status: Draft
:Type: Standards Track
:Created: 2025-05-14
:Resolution:

Abstract
--------

:ref:`NEP56` proposed adding nearly full support for the array API standard,
but many operations involving higher rank arrays still return scalars instead
of zero-rank arrays. This NEP would redefine the result of these operations
Comment on lines +15 to +16

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This suggests that numpy scalars are not compatible with what the array-api defines as an "array object". But, is that really the case? Because on most situations, you can get away with treating a scalar as if it's a full-fledged array, even if you index it as one:

>>> duck = np.array(1)
>>> ducknt = np.int_(-1)
>>> duck[()]
np.int64(1)
>>> ducknt[()]
np.int64(-1)

I noticed that it also has an __array_namespace__. But maybe I'm missing something?

Copy link
Owner Author

@mdhaber mdhaber May 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think most of this discussion has already played out in scipy/scipy#22947 (comment). So yes, we can play devil's advocate and say that ndarray and generic work together to implement a weird array protocol where arrays might be immutable when they are 0d. But the standard itself calls out NumPy explicitly in several places for producing scalars instead of arrays, so with the standard as written, it's not a very fun game.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So array-api objects are mutable, and numpy scalars are no. Got it.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hehe, even trickier. Array API is either mutable or unmutable. But the abstract class of numpy array|scalar would have undefined mutability.

to be zero-rank arrays.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The term "rank" (in the "no. of dimensions" sense) isn't used very often in the numpy docs. Maybe because in linear algebra (and in the UK) "rank" has a different meaning?

Either way, I suppose it would be more numpy-esque to e.g. refer to "zero-rank" as "0-D", "zero-dimensional", or casual "0d".

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

NEP 27 proposed the addition of "Zero-rank arrays" to NumPy. That seemed like a fairly appropriate source to borrow the term from.

So it's fine if the majority of authors would prefer a different term, but I don't immediately see a reason to change it.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would've said the same thing if I would've been around when reviewing NEP 27 🤷🏻.


Motivation and scope
--------------------

The 2024.12 version of the array API standard [1]_ states:

Apart from array object attributes, such as ``ndim``, ``device``, and
``dtype``, all operations in this standard return arrays (or tuples of
arrays)...

Beginning with :ref:`NEP56` and NumPy 2.0.0, NumPy added nearly full support
for the standard, but explicitly deferred compliance with this aspect.

We note that one NumPy-specific behavior that remains is returning array
scalars rather than 0-D arrays in most cases where the standard, and other
array libraries, return 0-D arrays (e.g., indexing and reductions)...
There have been multiple discussions over the past year about the
feasibility of removing array scalars from NumPy, or at least no longer
returning them by default. However, this would be a large effort with some
uncertainty about technical risks and impact of the change, and no one has
taken it on.

This NEP represents an effort to "take it on". It is a worthwile undertaking:
scalars "basically duck type 0-D arrays", but they do not *fully* duck type
zero-rank arrays, with the most fundamental difference being that scalars are
immutable and zero-rank arrays are not.
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe not surprising, I don't like the Array API focus, I would much prefer to only mention it as an additional argument.
First, all of the discussion is much older than it and all of the arguments work just as well without it. (In fact, to some degree better, because scalars aren't as bad at array-api as they at numpy.)

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are "all of the arguments"?

It sounds like we have two completely different perspectives about why 0d arrays should be returned instead of scalars. I think we could include both in the NEP so that we have a convincing arguments both with and without Array API, but I am less familiar with why scalars are "bad... at numpy".

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do what you prefer. I would expect all of your problems while doing Array API work to not actually be related to array api adoption itself?
Having to sprinkling in asarray() was always annoying and problematic maybe I am missing the argument why it is much worse now than what it was?

So, I would prefer to start with the issues below, rather than starting with "violating Array API".

The 0-D array->scalar conversion creates difficulties when writing functions that work with N-D inputs that may have 0-D intermediates. I.e. we have to add that np.asarray() call, but even if we do things may misbehave for object arrays.
The object array problem isn't huge for downstream typically, but has no good work-around.

Copy link
Owner Author

@mdhaber mdhaber May 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I'll need to research this "object array problem" more. I haven't run into it, personally. I assume this is "It changes the dtype of 0d object arrays containing array-likes" from numpy#13105.

maybe I am missing the argument why it is much worse now than what it was?

Probably the fundamental reasons are the same. I've just seen more frequent reminders of the inconsistency.

  • Our private array API test assertions (xp_assert_close) have spurred new arguments among SciPy maintainers about whether to return 0d arrays or scalars.
    • stats chose years ago to follow the lead of NumPy and special, and so we find lots of [()] shortly before functions return. But with the array API, all these are getting converted from something like return x[()] to the more verbose return x[()] if x.ndim == 0 else x. Also, the array API prompts us to update our documentation about return types. Having to write "array or scalar" all over the place makes me think how much simpler it would be if we could follow the array in -> array out rule.
  • The array API prompts us to rewrite old functions in a vectorized way. if/else branching is often converted to logic based on boolean arrays, and sometimes we are tempted to use boolean index assignment. Scalar immutability causes issues.
  • With the array API, we are now more careful about dtypes. We still support NumPy 1.26, though, so I frequently run into dtype issues due to the old scalar dtype promotion rules.

These are just a few examples off the top of my head. These particular things may not be relevant to everyone in the future, so I wouldn't mention them in the NEP. The point is that that the standard prompts code rewrites and reminds us to think about details that might have been ignored before. This inevitably exacerbates the old scalar/0d array problems.

One might say that by the time this NEP would take effect, all array API rewrites will be done, so it won't matter any more. I think the same argument can be made about any bug or undesirable API, though - once you've worked around the problem, it stops bothering you... for a while, at least, until the next time you encounted the same problem. I say we take the cue and fix it once and for all.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Object 0-D is the same as your mutability issue in the end. It leads to some code in NumPy wrapping np.asarray() around it, I have already fixed at least two bugs with out=... that would have been very hard to fix without and this addresses similar issues.
I.e. res = arr0d + arr0d, means that res is the scalar. And in the worst case np.asarray(res) doesn't even round-trip (e.g. because res is itself an array) if:

arr0d = np.array(None)
arr0d[()] = np.arange(10)

But of course typically, we just want to avoid the np.asarray() because it is easy to forget and may need to be sprinkled in multiple times...

For me the issues around this are the main problem I want to solve.

With the array API, we are now more careful about dtypes. We still support NumPy 1.26, though, so I frequently run into dtype issues due to the old scalar dtype promotion rules.

Yeah, although it seems unrelated? Scalars and 0-D always behaved the same in promotion, only Python scalars changed behavior.

stats chose years ago to follow the lead of NumPy and special, and so we find lots of [()] shortly before functions return. But with the array API, all these are getting converted from something like return x[()] to the more verbose return x[()] if x.ndim == 0 else x. Also, the array API prompts us to update our documentation about return types. Having to write "array or scalar" all over the place makes me think how much simpler it would be if we could follow the array in -> array out rule.

This is something where I am honestly not sure what you actually want ☹. If you don't like the arr[()] return, it seems to me that you can choose to stop using it. Is the argument that NumPy leading the way for ufuncs makes that easier?

But if NumPy preserves scalars in -> scalars out, that seems not really easier for on you, unless you are considering something where np.add(np.float64(1), 1) does return an array, but maybe np.float64(1) + 1 doesn't?

So in that case the question is what additional places you want an array return that I bracketed out in the minimal version.
But I suppose something that is short of the maximal change with the env var?

My view has always been that the feasible change to me is the one where we don't implicitly convert 0-D arrays to scalars, but arr0d[()] and arr.sum() are explicit enough to me.

One might say that by the time this NEP would take effect, all array API rewrites will be done, so it won't matter any more.

These type of things never go away, and they do have real, if rare, user impact. But of course the question is how much library author pain is OK if the alternative is for some user code to silently start doing the wrong thing.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the argument that NumPy leading the way for ufuncs makes that easier?

Yes. We should not have to fight NumPy's behavior to make this happen. We should be able to run the code and not have to either [()] or xp.asarray on the output of every function.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, but these are two distinct annoyances:

  1. The need for you to have to call np.asarray() when working with N-D arrays, since suddenly facing a scalar is very annoying and bug-prone.
  2. You want to remove the if arr.ndim == 0: return arr[()] that matches NumPy logic currently. To me, that one is more your choice than mine. Right now you choose to align with NumPy ufuncs so you ensure a scalar result with arr[()].
    In the future you would probably just remove it and always return arrays. That makes sense, but will still diverge from NumPy, which will return scalars for scalar only inputs.

So this NEP solves the first issue fully. The second one is still the downstream libraries problem at least in the sense that you must change your code to make it happen.
And in principle you could have made that choice earlier to say "our functions only work with and return arrays".

Now, transitioning SciPy depending on the NumPy version having implemented this NEP makes sense to me as it makes the transition clearer and cleaner.
But I suppose the fact that libraries are unlikely to fully align with NumPy makes me want this point to be clear, at least between us. Right now, I am not quite sure we are on the same page that this NEP only makes problem 2 go away in the very long term and only by diverging from NumPy by only returning arrays.

(And of course the reason why NumPy doesn't change, is that scalar + scalar -> array is very awkward unless scalars effectively don't exist, and I think the tries showed that this is a too impactful change probably.)

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right now, I am not quite sure we are on the same page that this NEP only makes problem 2 go away in the very long term and only by diverging from NumPy by only returning arrays.

I think we're on the same page.

I understand that with this change, I will generally no longer need to use asarray within code written for N-d arrays because (unless axis=None or I fully index a single element) it will no long suddenly turn into a scalar. That is 1.

I understand that removing [()] has always been a choice I could make independently of NumPy.
"transitioning SciPy depending on the NumPy version having implemented this NEP makes sense to me as it makes the transition clearer and cleaner" - yes!
"only by diverging from NumPy by only returning arrays." - I think see why you would write this, but I'm not sure it's necessarily true. IIUC, special uses the ufunc machinery from NumPy, so it would follow NumPy ufunc machinery behavior exactly, whatever that is. These will stay consistent with NumPy.
The functions that have [()] at the end now are mostly reducing functions, so they don't accept scalars or 0d arrays, and therefore there's never a reason for them to return scalars. These will stay consistent with NumPy reducing functions (except maybe that they will not accept 0d objects).
For other functions which do accept scalars and don't use ufunc machinery... I don't know right now what we'd do. In the long run, if SciPy is all-in on the array API, we might choose to only accept arrays (and therefore we'd never return scalars). If we allow scalars, we might choose to explain that this will behave identically to using np.asarray before passing in input. I understand that this would be inconsistent with NumPy, so instead, we might choose to accept scalars and treat them differently from 0d arrays. If internally we use xp.asarray on the input, we'd need to make sure that we keep track of whether the original input was a 0d array or scalar so we could return the appropriate thing at the end.

So if all that makes sense, I think we're on the same page w.r.t. 2, although I'm not certain that means we will diverge from NumPy by only returning arrays. There are a few possibilities.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are "all of the arguments"?

It would also help a lot for static typing. Because, currently, the simplest numpy-esque function is currently pretty difficult to properly type:

@overload
def noop[...](a: ndarray[Shape0D, dtype[ScalarT]]) -> ScalarT: ...
@overload
def noop[ShapeT: AtLeast1D, ...](a: ndarray[ShapeT, DTypeT]) -> ndarray[ShapeT, DTypeT]: ...
@overload
def noop(a: ArrayLike) -> TheShapeIsNotKnownSoEitherArrayOrScalar: ...

If it always returned an array, this could be simplified a lot. It would also help the users, because unknown "array-likes" will retun "array of unknown dtype" instead of "array of unknown dtype or some unknown scalar". The latter often requires users to manually cast it to the "actual" type. That's why currently the numpy-stubs will often return Any in such situations, even though Any is horribly type-unsafe.

tldr; This NEP would help a lot for static typing.


It may be argued that if instances of array base class ``np.ndarray`` and
scalar base class ``np.generic`` were fully interoperable, together, they
would implement a protocol compatible with the array API standard. Even if this
were the case, this design is complex and leads to confusion and errors due to
self-inconsistency (zero-rank array-like scalars are immutable, but arrays of
other rank are mutable) and inconsistency with all other array API compatible
libraries. In particular, it leads to difficulties in working with vectorized
reducing functions, which begin with arrays of rank :math:`N` and return
objects of rank :math:`M < N`: when :math:`M = 0`, the rules change. This
prompts an unfortunate pattern of calling ``asarray`` on the results of
intermediate array operations to ensure that operations like boolean mask
assignment still work. The inconsistency also presents downstream library
authors with an unfortunate choice: should they maintain consistency with
NumPy and prefer to return scalars when possible (e.g. ``scipy.stats``, which
explicitly uses empty-tuple indexing on all results to *ensure* consistency,
and ``scipy.special``, which relies on NumPy ufunc machinery), or should they
follow the lead of the array API standard and prefer zero-rank arrays (e.g.
``scipy.interpolate``).

Usage and impact
----------------

Currently, most operations in NumPy involving zero-rank arrays return scalars,
reducing operations that would naturally result in a zero-rank array actually
produce a scalar, and indexing operations that would naturally result in a
zero-rank array actually produce a scalar. The proposal is for these operations
to return zero-rank arrays instead of scalars.
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is confusing, the below examples are almost right, but speaking too generally about reductions isn't?
I have been very clear that in my branch both:

  • arr.sum()
  • arr1d[0]`
    return scalars.

You could try more for either points (I don't really believe in it from a BC stand-point, because those are the two cases where it makes a lot of sense to have code that continues relying on scalars).

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have been very clear that in my branch both:

Must have missed that. I saw:

The simple rule is, if any input is already an array, the output should also be an array. The one extension is that axis=None for reductions should return scalars.

So I suppose there is a second exception, and the easiest way to summarize it is that indexing behavior is unchanged.

Copy link

@seberg seberg May 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤷 fair, I wanted to refer to ufuncs/functions there. It does continue with:

Since that doesn't work at all, it adds an arr.to_scalar() method, that does the same as arr[()] did previously.

In the context of doing the bigger change with NUMPY_DISLIKE_SCALARS=1


If you feel strong about changing indexing, we could still try it out. As you can tell, I am focused on the minimal useful change.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. And I was hoping to go a little further and follow the array API more closely. We're coming from different perspectives. I think it's fine now that we understand that.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As for whether we should try it out - yeah, I think there's value in seeing how much additional trouble it causes. I think the best decision would come from understanding the tradeoff between switching behaviors and downstream pain. Right now, we only know that NUMPY_DISLIKE_SCALARS=1 - which is extremely aggressive, and more so than I have any motivation to propose - looks quite a bit more painful than NUMPY_DISLIKE_SCALARS=0. You made an educated guess that leaving axis=None as a special case would be much less painful - but I don't think we have any data. Maybe the transition from the painful NUMPY_DISLIKE_SCALARS=1 is very sharp and happens mostly because scalar-scalar math suddenly results in 0d arrays or something. It would be good to know!

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, TBH, I just don't expect arr.sum() to be a huge gain either, so I want to err on the safer side.
The next question is how much you want to try arr[()], which is an orthogonal knob (and there may be more?).


.. code:: python

import numpy as np
x = np.asarray(1.)
np.isscalar(x + x) # True (main), False (NEP)
np.isscalar(np.exp(x)) # True (main), False (NEP)
y = np.ones(10)
np.isscalar(np.sum(y), axis=-1)) # True (main), False (NEP)
np.isscalar(y[0]) # True (main), False (NEP)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The last one doesn't change. In fact, indexing doesn't change at all.


For exceptions to these rules, ask Sebastian.
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One exception is axis=None -> scalar.
If I understood correctly, indexing an array like:

x = np.asarray([1, 2, 3])
x[0]

produces a scalar, too?

axis=None is one thing, but this would make selling the proposal even harder, since it's a half-measure.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

since it's a half-measure.

It makes it much easier because it avoids a swaths of BC issues since arr1d[0] will in practice be one of the most common paths to get scalars in user code.
(While it will actually be not that common in major libraries.)

You can call it a half-measure if you like, but I don't see it as such at all. The point is x[0, ...] or x[..., 0] is a simple way to get what you want already.

To home in on what IMO, is the core of all of this: Make scalar returns predictable and avoidable when working with N-D arrays.
Yes, scalars have more annoyances, but I think focusing on the above makes it easier to see that this is the minimal way to achieve it and avoiding the temptation to change things that really don't need changing.

Copy link
Owner Author

@mdhaber mdhaber May 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make scalar returns predictable and avoidable

I don't see how to sell that, though. Individual NumPy operations are pretty predictable right now - operations other than array creation functions will always return a scalar instead of a 0d array. There may be exceptions to that rule, but doesn't that sum it up, or is there a ton of inconsistencies that this PR resolves? (Can you provide examples of surprises we get right now that are not summed up by the rule above?)

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, the problem isn't really them being unpredictable, the problem is that there is no good way to avoid the 0-D array -> scalar conversion.
The issue is preserving the 0-D arrays to be able to work with N-D arrays without hassle.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there is no good way to avoid the 0-D array -> scalar conversion.

OK, so do I understand that somewhere in the code, the results of ufuncs and reducing operations is a 0d array, and then it gets converted (unavoidably) to a scalar. You are reversing that choice except for reducing operations when axis is None.

When we fully index an array (e.g. one integer for a 1d array, 2 integers for a 2d array), is the result naturally a 0-d array or a scalar? Is it like ufuncs in that the result is naturally 0d and it gets converted to a scalar, or is it the other way arround - the result is naturally a scalar and if we wanted it to return a 0d array, it would not be so different from just calling asarray on the result?


Empty-tuple indexing may still be used to cast any resulting zero-rank array
to the corresponding scalar.

.. code:: python

import numpy as np
x = np.asarray(1.)
np.isscalar(x + x) # True (main), False (NEP)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please replace all isscalar with isinstance(..., np.ndarray). isscalar is very awkward API and much less precise anyway.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, I can't replace isscalar(x) with literally isinstance(x, np.ndarray) because these would basically be opposites. Do you mean change it to not isinstance(..., np.ndarray)?
I've often found isscalar useful and less awkward than trying to acheive the same thing as isinstance(... np.ndarray) because in the past SciPy has treated lists and DataFrames and such just like ndarrays, but they are not instances of ndarray. Why is it so bad?

Copy link

@seberg seberg May 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because it has a very narrow meaning of "scalar" that is just not quite right in this context (honestly it is never quite right, but does a decent enough job for numeric types).
For example isscalar(None) is False even though clearly it's not an array and np.asarray(None)[()] round-trips, so it clearly is a scalar (unless you define it very narrowly).

So yes, of course I mean not isinstance, because I think in this context is clear and has no caveats that isscalar doesn't have a precise/correct definition.


The main impact to users is more predictable results due to improved
consistency within NumPy, between NumPy and the array API standard, and
between NumPy and other array libraries. Working with the results of reducing
functions, in particular, will be easier because return values of any rank
will support boolean indexing assignment.

There is a secondary impact on performance. On typical hardware, execution
time of conversion from zero-rank arrays to scalars and elementary arithmetic
operations involving only scalars is on the order of tens of nanoseconds,
whereas operations involving only zero-rank arrays is on the order of hundreds
of nanoseconds. Consequently, some elementary arithmetic calculations will be
slower. On the other hand, conversion from scalars to zero-rank arrays takes a
few hundred nanoseconds, and many operations, such such as ufuncs and
operations involving both scalars and rank-zero arrays require conversion from
scalars to zero-rank arrays. These operations will be faster. We will not
speculate as to whether this will have a net positive or net negative impact on
user applications, but the net impact is expected to be small since impact on
downstream library test suites has been minimal in testing.

Backward compatibility
----------------------

The motivation of this proposal is to eliminate the surprises associated with
a scalar result when a zero-rank array would be expected. However, existing
code may rely on the current behavior, and this presents backward compatibility
concerns.

The main concern for user code is that the mutable zero-rank arrays that
replace immutable scalars are no longer hashable. For instance, they cannot
by used directly as keys of dictionaries, the argument of an ``lru_cache``
-decorated function, etc. In all circumstances, tbe patch is simple: convert
the zero-rank array to a scalar with empty-tuple indexing.
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is OK for the hashability but needs to be expanded on the mutability. While less common, it does need to be pointed out explicitly because it can lead to incorrect results.

(Unlike hashability, which should mostly lead to a relatively clear error.)

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought about this, but I don't really understand how existing code could run into this problem.

If you try to mutate a scalar now, it raises an error. Presumably this is not what was intended, so this will not be part of a typical code path.

Do you mean that people might be relying on mutating a scalar to raise an error, but with the change it would works as expected, and that would be undesirable?

I was assuming this would hold true:

Going from error to non-error (nan or whatever) should always be fine, the reverse is typically not fine.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you try to mutate a scalar now, it raises an error.

It doesn't raise for in-place operators. I actually have problems of thinking of an example (which is good!). But a silly one:

def my_metric(x):
    x += 1
    return x

arr = np.random.random()

s = arr.sum(axis=0)
metric = my_metric(s)  # modifies s

metric_s_values.append((s, metric))

The point is, having a 0-D array where code is clearly written for scalars (using in-place operators) is bad news. I could imagine it also for a simplistic integration_step function?

Is that common? Hopefully not, and to me if arr.sum() and arr[0] keep returning scalars it seems like it shouldn't be common to suddenly pass 0-D arrays to a helper function that uses in-place arithmetic.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah right. I realized that at the conference and recently forgot.


Running the test suites of dependent libraries against a branch of NumPy that
implements these changes has revealed a few other issues. <Library maintainers,
please summaries these issues here.>

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One worry I still have -- not sure how big it is -- is the fact that it isn't necessarily trivial for downstream libraries to follow NumPy behavior.
I.e. numpy does scalar in -> scalar out, but for downstream that is hard to impossible to copy.
Unless downstream libraries want to be -> array out always on newer NumPy versions, they will API exposed and use that API likely.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but for downstream that is hard to impossible to copy.

I don't see how it is impossible. Even a decorator could be written to detect input type (scalar or array) and ensure that the output types are consistent. It might require an explicit conversion, but aside from the tens to hundreds of nanoseconds, a user would not notice the difference.

Unless downstream libraries want to be -> array out always on newer NumPy versions

They might. This is something SciPy will certainly be discussing for SciPy 2.0.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see how it is impossible. Even a decorator could be written to detect input type (scalar or array) and ensure that the output types are consistent

Yes, but we will need public API from NumPy to do this correctly, and that machinery will again not fit well into your Array API code.

Detailed description
--------------------

The new functionality will be used in much the same was as old functionality,
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

s/same was/same way/

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The new functionality will be used in much the same was as old functionality,
The new functionality will be used in much the same way as old functionality,

except that sometimes authors will need to convert zero-rank arrays to scalars
rather than converting scalars to zero-rank arrays. For instance, consider the
following:

.. code:: python

import numpy as np
rng = np.random.default_rng(85878653462722874072976519960992129768)
x = rng.standard_normal(size=10)
y = np.sum(x, axis=-1)
z = {y: 'a duck'} # use scalar result as dictionary key
y = np.asarray(y) # convert to array to allow mutation
y[y < 0] = np.nan

The use of ``x`` as a dictionary key would need to become ``x[()]``, but ``z``

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The x and y here should both (all 3 cases) be y?

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oops; I was messing with the code and forgot to update the text. Thanks!

no longer needs to be explicitly converted to an array.

.. code:: python

z = {y[()]: 'a duck'} # extract scalar to use as dictionary key
y[y < 0] = np.nan # no conversion to array required

Realistic examples can be found throughout the codebases of dependent
libraries.

Related work
------------

All known libraries that attempt to implement the array API standard
(e.g. ``cupy``, ``torch``, ``jax.numpy``, ``dask.array``) return
zero-rank arrays as specified by the standard.

Implementation
--------------

To implement the NEP, the branch prepared by Sebastian needs to be merged,
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should I link https://github.com/seberg/numpy/tree/try-0d-preservation-rebase? (Or will there be a different branch, do you think?)

and dependent libraries will need to prepare releases that adapt to (and take
advantage of) the changes. Branches of libraries including SciPy and Matplotlib
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@j-bowhay @pfackeldey @ksunden can you provide links to your branches?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the link to Awkward Array with the NumPy branch provided by @seberg: https://github.com/scikit-hep/awkward/tree/pfackeldey/test_NEP57. There are only 2 failures and they are related to pandas. For Awkward Array the migration seems to be simple.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just catching up on notification now, happy to provide a branch for SciPy and read over this but will probably need a week or two to do so

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the status here? I am still happy to work on a branch but don't want to invest the time if things are blocked elsewhere

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have worked on it here: numpy#29067 but there are still some issues around polynomials mostly I think, where the NumPy python code needs to properly deal with the changes (which is a bit tedious).

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's ok to wait, @j-bowhay. We'll ping you when things are closer.

have already been prepared without much difficulty. Initially, users will have
the option of opting into this behavior using an environment variable, so these
releases need to be compatible with both old and new behaviors. To make the new
behavior the default and only behavior, NumPy will need to advise users of the
pending change and give the appropriate notice. The initial draft of this
document does not specify the appropriate timeline and procedures; this line
will be updated to reflect the consensus of the maintainer team.
Comment on lines +177 to +179
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@seberg should I propose something? What would you suggest?

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Propose something. I would say:

  1. Create a branch and use it to adapt and test downstream.
  2. The final NEP acceptance will happen after this with potentially more discussion.

Or in other words, I would just not say much about the precise path of achieving step 2. That makes sense IMO, because without 1 it is even harder to judge how big of a change it is (not that 1 will be a silver bullet).

How exactly you want to write it, I don't care. You could even add a .. note:: early on that a final decision is deferred until after further testing.


Alternatives
------------

There are two main alternatives to this proposal.

The alternative suggested by :ref:`NEP56` is to maintain the current behavior
and make NumPy scalars more fully duck-type zero-rank arrays, such as adding
missing behaviors.

While this work would still be valuable, we propose the behavior change to more
fully comply with the standard and to resolve the problems mentioned in the
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this somewhat leaves out a huge part of it which is the non-array API part? The whole "just add methods" works very well for array-api numeric types, but badly for general scalars like object/string dtypes.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this somewhat leaves out a huge part of it which is the non-array API part

Could you be more specific about what the "non-array API part" is? From context, I think you are referring to an important argument in favor of the change, but I don't know what argument that is.

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you reply to this part @seberg?

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think our discussion above covered this, but I'll try to see if there is more (it may make sense to look at NumPy bugs around this, but SciPy fixes seem actually more interesting for the NEP).
I.e. I mainly mean the fact that sometimes you have to sprinkle in np.asarray() because of the magic 0-D -> scalar conversion.

Motivation.

More extreme alternatives are also available, such as eliminating NumPy
scalars entirely. We do not take this approach for two reasons:

1. Scalars still have some advantages as hashable, dtyped objects that support
very fast elementary arithmetic.
2. The backward compatibility concerns with eliminating scalars entirely are
much more severe.

A variant of this proposal is to eliminate the exceptional behavior associated
with reducing operations and ``axis=None``. I cannot present any justification
for why we shouldn't do this; I would certainly prefer it because it would be
even more consistent and fully compliant with the standard. Ask Sebastian.
Comment on lines +202 to +205
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please replace the sentences after the first for me : )

Copy link

@seberg seberg May 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minimal change, minimal change, minimal change. And maybe also why not?

If you change arr.sum() to return 0-D arrays you will create a class of places where you return 0-D arrays for, honestly, unclear reasons and you will break user code. Not as much as arr1d[0], but still...
If you change arr1d.sum(-1) I think that is much less, because who would write that in a non N-D context?!

I keep repeating myself, but the point is making the return predictable for N-D arrays, changing arr.sum(), frankly, to me smells like "I don't like scalars", and not "I want to work with N-D arrays in a predictable way".

Now, "I don't like scalars" is a valid opinion, but I think we established well enough that we won't achieve that easily, so I don't see why we should try to make inroads towards it that we don't need to fix the actual core issue (don't get random scalars when working with N-D arrays).

EDIT: To be clear, I am not in love with this rule, also because it is slightly different from the rule for indexing (where adding ... is the magic). But it seems like a rule that should have a very good hit-to-miss ratio in existing code and I suppose I also can't think of a neat rule that would look the same as indexing (relying on out=... or adding ... to the axis seems not so nice).

And yes, this is an opinion that you can reasonably disagree on and the desire for minimal change is at the core of it.


Discussion
----------

This section may just be a bullet list including links to any discussions
regarding the NEP:

- https://github.com/numpy/numpy/issues/24897
- https://github.com/scientific-python/summit-2025/issues/38
- https://github.com/scipy/scipy/pull/22947#discussion_r2080108060


References and footnotes
------------------------

.. [1] `Python array API standard 2014.12 — Array Object`_
.. [2] Each NEP must either be explicitly labeled as placed in the public domain (see
this NEP as an example) or licensed under the `Open Publication License`_.

.. _Open Publication License: https://www.opencontent.org/openpub/
.. _Python array API standard 2014.12 — Array Object: https://data-apis.org/array-api/latest/API_specification/array_object.html

Copyright
---------

This document has been placed in the public domain. [2]_
Loading