RFC: add support for computing the variance of complex array

[ntessore]

The standard should support taking the variance of a complex-valued array, which is statistically well-defined and unambiguous (var(z) = var(z.real) + var(z.imag)). As far as I'm aware, all array backends already support this.

[rgommers]

The tracking issue for adding complex number support to all relevant APIs was gh-373, and for both std and var it says only:

• No complex number support.
• Change to only real number data types.

No reason for it given unfortunately. Probably was because of the correction keyword and/or the inf/nan special cases. correction does seem well-defined as total variance; there's also pseudo/relational variance and generalized variance, in case numbers aren't simply uncorrelated points in the complex plane (hard to find a single concise reference, https://en.wikipedia.org/wiki/Variance#For_complex_variables has the latter).

If all libraries implement total variance though, I don't see a reason yet not to do the same in the standard.

@kgryte do you remember why you excluded std/var?

[kgryte]

Complex support was originally added in the v2022 revision, and, in that revision, complex number support in mean, std, and var were considered lower priority. In the v2024 revision, we added support for complex numbers in mean (ref: #850).

So, I think the answer is that we have added complex support to statistical APIs on an as requested basis since adding initial support in 2022.

[rgommers]

Then I think we should check all the statistical functions for whether they have universal complex dtype support; std is well-defined as well, so should be treated the same as `var.

[ev-br]

numpy, torch, jax and cupy seem to be in agreement:

In [1]: import numpy as np

In [2]: x = np.asarray([1 + 2j, 3 + 4j, 5 + 6j])

In [3]: np.var(x)
Out[3]: np.float64(5.333333333333333)

In [4]: np.var(x.real) + np.var(x.imag)
Out[4]: np.float64(5.333333333333333)

In [6]: import torch

In [7]: xt = torch.as_tensor(x)

In [8]: torch.var(xt, correction=0)
Out[8]: tensor(5.3333, dtype=torch.float64)

In [10]: torch.var(xt.real, correction=0) + torch.var(xt.imag, correction=0)
Out[10]: tensor(5.3333, dtype=torch.float64)

In [12]: torch.std(xt, correction=0)**2
Out[12]: tensor(5.3333, dtype=torch.float64)

In [13]: import jax.numpy as jnp

In [14]: xj = jnp.asarray(x)
An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.

In [15]: jnp.var(xj)
Out[15]: Array(5.3333335, dtype=float32)

In [16]: jnp.var(xj.real), jnp.var(xj.imag)
Out[16]: (Array(2.6666667, dtype=float32), Array(2.6666667, dtype=float32))

In [17]: jnp.std(xj)**2
Out[17]: Array(5.333333, dtype=float32)

In [18]: import cupy

In [19]: xc = cupy.asarray(x)

In [20]: cupy.var(xc)
Out[20]: array(5.33333333)

In [21]: cupy.var(xc.real), cupy.var(xc.imag)
Out[21]: (array(2.66666667), array(2.66666667))

In [22]: cupy.std(xc)**2
Out[22]: array(5.33333333)

[rgommers]

Discussed in today's meeting, everyone was happy to move forward with this.