Skip to content

Arithmetic between real arrays and Python complex scalars (revisited) #841

Closed
#871
Closed
Arithmetic between real arrays and Python complex scalars (revisited)#841
#871
Labels
API changeChanges to existing functions or objects in the API.Needs DiscussionNeeds further discussion.topic: Complex Data TypesComplex number data types.topic: Type PromotionType promotion.
Milestone
v2024

Description

@mdhaber
mdhaber
opened on Sep 11, 2024

The array API standard is now clear (in Type Promotion Rules):

image

This means that following does not have defined behavior:

x = xp.asarray([1., 2., 3.])
x + 1j

However, the behavior is defined consistently in NumPy, CuPy, PyTorch, jax.numpy1, Dask array, and Tensorflow2 in the intuitive way: the real and imaginary component dtypes of the complex output array match the dtype of the input array.

This issue was discussed explicitly in gh-478 beginning with #478 (comment). I see several comments that seem supportive of allowing this operation, e.g.

I would hope it at least works for Python scalars. x + 0j seems like the simplest and most obvious way to convert a real array to a complex one.

The comment I see against it (#478 (comment)) is

I know it's very common to do this (also * 1.0 to convert integer to real). However, it doesn't seem like great practice irrespective of the decision here - there's an explicit astype function for converting to other dtypes.

It looked like the summary comments suggested that this would be allowed, but it was not specifically addressed in the PR that closed the issue. (See postscript for more information.)

I thought it might help to add some perspective on how this impacts a developer translating code to the standard: when tests including an operation like

x + 1j

begin to fail with array_api_strict only, the developer is faced with the choice of skipping the array_api_strict tests or figuring out what is wrong and changing the code to something like:

xp.astype(x, xp.result_type(x.dtype, xp.complex64)) + 1j

If they also want it to be able to preserve the lower precision types supported by some libraries (e.g. for NumPy), there would be additional hurdles. I tested for a few libraries, and array API standard compatible version of the code also increases the execution time notably for small arrays. Individually, the inconvenience, complexity, and overhead are small, but they add up. I am excited about adding array API support to SciPy, but not everyone is supportive, and performance regressions and complicated-looking diffs can make it difficult to garner support.

I would suggest that the operation be defined, even if it means adding an exception to the simply stated rules for array/Python arithmetic.

The more specific language about this not being defined was added in gh-513, which provided the justification:

Given that current guidance requires converting a scalar to the array data type, to accommodate complex scalars and real-valued arrays, we'd need to specify casting rules for complex to real,

I think that is referring to the language:
image

That guidance was added in gh-74, but I can't trace the origin further. Perhaps there can be exceptions to that rule?

Footnotes

  1. when higher precision support is enabled, see changelog ↩

  2. With NumPy experimental behavior (https://github.com/data-apis/array-api/issues/478#issuecomment-1272631595). Vanilla Tensorflow doesn't seem to support the other array-scalar operations that are defined, but already this was not deemed to be sufficient reason to prevent their inclusion in the standard (https://github.com/data-apis/array-api/issues/478#issuecomment-1270409660). ↩

Activity

  1. kgryte
    added
    API changeChanges to existing functions or objects in the API.
    topic: Complex Data TypesComplex number data types.
    topic: Type PromotionType promotion.
    Needs DiscussionNeeds further discussion.
    on Sep 11, 2024
  2. kgryte
    added this to the v2024 milestone on Sep 11, 2024

  3. asmeurer commented on Sep 11, 2024

    @asmeurer
    asmeurer
    on Sep 11, 2024
    Member

    I still agree that this should be supported. We should be careful to not confuse downcasting a complex array (or scalar) to real float, which is indeed not supported and for very good reasons, with upcasting a real float to complex.

    In fact, real and complex floating-point arrays are already supported to promote together:

    image 
    >>> import array_api_strict as xp
>>> xp.asarray([0., 1.]) + xp.asarray([0j, 1j])
Array([0.+0.j, 1.+1.j], dtype=array_api_strict.complex128)

    So it seems odd that this doesn't work for Python scalars.

    The piece of text that @mdhaber quoted at the top of the issue seems to imply that there is ambiguity whether the resulting array should be complex64 or complex128. But I think this is already pretty well spelled out for the Python float case https://data-apis.org/array-api/latest/API_specification/type_promotion.html#mixing-arrays-with-python-scalars

    1. Convert the scalar to zero-dimensional array with the same data type as that of the array used in the expression.

    2. Execute the operation for array 0-D array (or 0-D array array if scalar was the left-hand argument).

    i.e., float32_array + python_float gives float32 and float64_array + python_float gives float64. In other words, the corresponding rule should be that a float32_array + python_complex gives a complex64 array and float64_array + python_complex gives a complex128 array.

    I think the real question is, what is the idea behind the rule "convert the scalar to zero-dimensional array with the same data type as that of the array"? If the idea is that "when adding a Python scalar to an array, the dtype doesn't change", then indeed real_array + complex_scalar shouldn't be supported. But if it's rather that "the precision doesn't change", then it should be supported.

  4. rgommers commented on Oct 3, 2024

    @rgommers
    rgommers
    on Oct 3, 2024
    Member

    xp.astype(x, xp.result_type(x.dtype, xp.complex64)) + 1j

    This is quite clumsy indeed. Independent of the decision on allowing x + 1j, we should also make this nicer by supporting dtype kinds in astype I think, so the above becomes:

    xp.astype(x, 'complex')  # or maybe `'complex floating', if we'd use the `isdtype` kind specifiers

    Doing that would also address the problem for unsigned -> signed integers, which cannot be done with a Python scalar equivalent like + 0j.

  5. ev-br commented on Nov 25, 2024

    @ev-br
    ev-br
    on Nov 25, 2024
    Member

    One rather annoying example I tripped over in scipy is

    >>> n = 5
>>> x = xp.arange(n, dtype=xp.float64)
>>> xp.exp(1j * xp.pi * x / n)
...
TypeError: Python complex scalars can only be promoted with complex floating-point arrays.

    An "obvious" workaround,

    >>> I = xp.asarray(1j)
>>> xp.exp(I * xp.pi * x / n)

    is actually a footgun because if x.dtype == xp.float32, the exponential ends up being complex128 (at least with numpy). So the correct incantation is something like

    >>> I = xp.asarray(1j, dtype=xp.complex64 if x.dtype == xp.float32 else xp.complex128)

    which looks like something 1j * x should do under the hood.

  6. kgryte
    closed this as completedin #871on Feb 6, 2025
  7. kgryte
    added a commit that references this issue on Feb 6, 2025

    feat: define complex scalar <op> fp array

    c7476ba
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    API changeChanges to existing functions or objects in the API.Needs DiscussionNeeds further discussion.topic: Complex Data TypesComplex number data types.topic: Type PromotionType promotion.

    Type

    No type

    Projects

    No projects

      Milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions