Specify the expected behaviors for handling complex numbers

[leofang]

Follow-up of #102. Currently we state that complex numbers will be supported in the next version (v2) of the standard. However, there is nothing forbidding the complaint libraries from implementing them, if they choose so (ex: CuPy has full support already).

In today's call, people seem to be happy about the idea that we specify the expected behaviors of APIs that could accept (or generate) complex numbers, so that we can provide a guideline for those libraries choosing to support it, and also earn a long enough time to potential changes or iterations before v2 is finalized.

Therefore, off top of my head a few linalg PRs will need some changes to set up the expectations, including

• Eigensolvers (#113)
• SVD (#114)
• QR (#126)

In addition, this would also unblock FFT (#159).

[leofang]

We can also have separate PRs to address this. Maybe easier to review.

[asmeurer]

I would expect most elementwise functions that accept floating-point inputs would also accept complex inputs, and likely many special cases will need to be added.

Another question: for floats there is IEEE 754. Is there a standard for complex numbers we can reference, so we don't need to specify basic semantics?

[leofang]

I don't think any IEEE standards specify it. Our best shot is the C++ (std::complex) or C99 (complex.h) ISO standards.

[rgommers]

Thinking about the outline here:

Maybe a Future API header below the API sections?

I think we need a separate section somewhere for v2 content that's clearly separate from v1. That way we can merge PRs and can reference them without people being confused about what's in/out of scope for v1.

Our best shot is the C++ (std::complex) or C99 (complex.h) ISO standards.

Do you know off the top of your head if these two standards are 100% compatible regarding numerical behavior?

[leofang]

Maybe a Future API header below the API sections?

I think we need a separate section somewhere for v2 content that's clearly separate from v1. That way we can merge PRs and can reference them without people being confused about what's in/out of scope for v1.

Yes, we should definitely have a Future API section. But, for the case of complex numbers, I am afraid not much clean separation we can do. For example, do we want to duplicate the docstrings for eigen, svd, qr, ..., just to add the mention of complex numbers? I doubt we want to do this.

As I said earlier it's fine to have separate PRs. The first is to focus on real floats (as is being done for many linalg functions), and a follow-up PR to add docstring modifications in place with a big v2 warning marker (and perhaps add to the list in Future API).

Our best shot is the C++ (std::complex) or C99 (complex.h) ISO standards.

Do you know off the top of your head if these two standards are 100% compatible regarding numerical behavior?

In terms of memory layout, yes. In terms of math function behavior, I am unaware of a strong guarantee for the compatibility and we need to evaluate case by case.

[mhvk]

(Hoping this is the right place) For complex numbers, it would be important to think through exact definitions for some of the elementwise functions. While most are clear, some are less so: in particular, for numpy, an arguably poor choice was made for sign - best would be to use z/|z| (e.g., as done in julia); for discussion, see https://github.com/numpy/numpy/issues/3621#issuecomment-22701337; more general issue about useful complex elementary functions, https://github.com/numpy/numpy/issues/13179).

[rgommers]

Thanks for pointing that out @mhvk! (and yes, this is a great place for this topic)

In addition to thinking carefully about all the known challenges with functions like on the numpy tracking issue you opened, I think it would be useful to have complex numbers in the test suite, so we can run it against multiple libraries easily and compare behavior. Detecting diverging behavior will help find pain points.

Also, PyTorch just completed its complex dtype support, and carefully compared every function to its numpy equivalent. The PyTorch issue tracker with the module: numpy plus module: complex labels will allow easily finding all the deviations (link).

The final thing that comes to mind is this warning that numpy has a habit of emitting:

>>> x = np.array([2+2.j])
>>> np.asarray(x, dtype=float)
...
ComplexWarning: Casting complex values to real discards the imaginary part

That should raise an error, it's been the cause of too many issues and confused users.

[asmeurer]

The behavior of nans with complex numbers is very odd, even with Python complex. I think some of this behavior might also happen in C++, based on limited testing

>>> complex(float('inf'), float('inf'))
(inf+infj)
>>> float('inf') + float('inf')*1j
(nan+infj)
>>> complex(0, float('nan'))
nanj
>>> float('nan')*1j
(nan+nanj)

Are these complex numbers that can apparently only be constructed directly via the complex constructor valid?

[asmeurer]

There are also some complex specific functions that should be added. Off the top of my head, real, imag, and angle (we should check the API comparison to see if there are any others or if they are spelled differently).

[mhvk]

@asmeurer - yes. Complementary to the argument (angle) is of course the modulus, which is captured by abs. The latter is not complex-specific, but it would also be really good to have the equivalent of abs(z)**2 (modulus squared, or power), which is most useful for complex numbers (in numpy, square does a direct square, but one wants z * z.conj()).

[asmeurer]

conjugate is another one.

@mhvk is there an existing function that does z*conjugate(z) in one of the existing Python libraries, or in the C or C++ standards?

[mhvk]

@asmeurer - Julia has abs2. C++ has (the rather badly named) std::norm.

[rgommers]

This is all completed now (see gh-373). The nan/inf special cases, z/|z| definition, and error for complex->real were all specified in line with the discussion here.

Thanks everyone!