`jax.numpy.cov` and `numpy.cov` have different results for empty arrays

[adriagarp]

Description

While working on data-apis/array-api-extra#451, we discovered different behavior between jax.numpy.cov and numpy.cov for empty arrays:

import jax.numpy as jnp
import numpy as np
print(repr(np.cov(np.array([]))))
print(repr(jnp.cov(jnp.array([]))))

Output:

array(nan)
Array(-0., dtype=float32)

Is this intended? Is this something that you would consider changing?

Thanks!

System info (python version, jaxlib version, accelerator, etc.)

jax:    0.7.1
jaxlib: 0.7.1
numpy:  2.3.3
python: 3.13.7 | packaged by conda-forge | (main, Sep  3 2025, 14:24:46) [Clang 19.1.7 ]
device info: cpu-1, 1 local devices"
process_count: 1
platform: uname_result(system='Darwin', node='xxxxxxx', release='25.0.0', version='Darwin Kernel Version 25.0.0: Mon Aug 25 21:17:56 PDT 2025; root:xnu-12377.1.9~3/RELEASE_ARM64_T6030', machine='arm64')

[mohsinm-dev]

I reproduced this issue, and if I am not wrong, JAX only handles one case of empty arrays

if X.shape[0] == 0:  # Only handled 0 variables
      return array([]).reshape(0, 0)

But there is no check for the case where X.shape[1] == 0 (0 observations), which can happen when

• 1D empty array [] → atleast_2d() → (1, 0)
• 2D arrays like (2, 0) with variables but no observations

[jakevdp]

Thanks for the report! I don't think this is intended, but JAX's output does follow directly from the implementation of the covariance formula; e.g.

>>> x = jnp.array([])
>>> jnp.sum((x - x.mean()) ** 2) / (len(x) - 1)
Array(-0., dtype=float32)

Here x.mean() is NaN, but since there are zero terms in the sum, that NaN never actually enters the computation. It looks like NumPy special-cases this, and JAX should probably do the same for consistency.

[mohsinm-dev]

Hi @jakevdp, thanks for the clarification!

Just to check if I am understanding this correctly: the -0.0 result happens only because the current formula runs all the way through with n=0 (so it ends up doing 0.0 / -1), while NumPy adds a special case and returns NaN.

I was thinking the cleanest way is to add an early check right after shaping X, but before the main covariance computation, and return NaN when there are variables but no observations. Something like this:

# Handle case where we have variables but no observations
if X.shape[1] == 0:
    result_dtype = dtypes.result_type(X.dtype, dtypes.float_)
    if X.shape[0] == 1:
        # Single variable, no observations -> scalar NaN
        return array(np.nan, dtype=result_dtype)
    else:
        # Multiple variables, no observations -> matrix of NaNs
        return array_creation.full((X.shape[0], X.shape[0]), np.nan, dtype=result_dtype)

Does this look like the right approach?
If yes, I’d be happy to open a PR and add the appropriate tests to match NumPy’s behavior.

Thanks!

[jakevdp]

Hey @mohsinm-dev sorry I missed your comment here. I ended up sending a fix that looks pretty similar to what you proposed.