generalize multivarate_normal generator for complex numbers?

[cvanwynsberghe]

What's the need?
random.multivariate_normal works great for but is limited to sampling real vectors. Today, if ones would like to generate from complex multivariate distribution from (mu, Omega), sampling would need these steps:

1. get Omega^1/2 by cholesky decomposition (or svd)
2. generate a serie z of random complex gaussian numbers whose elements are in Nc(0, 1)
3. get the sample mu + Omega^1/2 * z

Featuring it? How?
Would it be possible to extend today's function to complex numbers? I imagine two possibilites
for implementing that:

1. in the case where the given input mean array and covariance matrix are complex type, get the dtype of mu and Omega, and infer the right distribution to sample
2. add a kwarg or dtype in function to enforce if resulting sampling is real or complex

[njsmith]

Probably the simplest API would be to make it a new, separate method.

[bashtage]

The first step should be to take a decision about allows standard_normal to return different dtypes, esp float32, float64, and possibly complex128 and complex64. If this was implemented then the rest could be cleanly implemented in the existing multivariate normal. If this isn't desirable, possible, it would seem strange to have a relatively special purpose generator for complex MV normals.

[eric-wieser]

+1 for a dtype argument to all of the random functions. I'd go a stage further and then add if dtype is np._NoValue and result_type(*args) != float64: warnings.warn("In future, the dtype will be inferred from the inputs")

[bashtage]

Unfortunately the current structure makes this really tedious for many distributions. For example the t is directly implemented so that the C code generates a normal and a chi square and then returns a single t rv. It would be much easier to reimplement types if it used the higher level normal and gamma and then combined after generation. In this hypothetical, changing styles would be as easy as changing the two core rv dtypes.

[cvanwynsberghe]

From my experience in complex distributions (somewhat limited, to signal processing), I actually mainly see the normal, normal multivariate, chi2 and whishart distributions, as extendable to complex numbers. For other ones, I don't know if such extension has been used or proposed, do you?

This question is just to guess how much an implementation of all the complex cases would be for already existing generators in numpy.

[njsmith]

Supporting dtype arguments is cool, but I'm not sure it should actually change what distribution we're sampling from? Usually these just change the storage / rounding.

[bashtage]

Yes, one could just cast to alternative types. This would be simple for float32 but the complex types would need additional calls/support since they require more bits than the defaults provide.

[njsmith]

Not sure what you mean. If sampling from a real distribution and storing the output in a complex dtype, you don't need any extra bits: the real part of the complex dtypes is identical to the corresponding real dtype.

(And to be clear: of course it also makes sense to support sampling from complex distributions, I just feel like maybe that should be a different method rather than something triggered by a dtype= argument.)

[eric-wieser]

There's already precedent for the dtype argument causing a very different result:

>>> np.sqrt(-1)
nan
>>> np.sqrt(-1, dtype=complex)
1j

Given that's the case, I think it would be very reasonable for the complex overload to behave in the same way.

[eric-wieser]

Today, if ones would like to generate from complex multivariate distribution from (mu, Omega)

Note that the wikipedia page claims that a complex normal distribution has three parameters, not two:

[bashtage]

I have cobbled an early version of a univariate complex_normal here:

https://github.com/bashtage/ng-numpy-randomstate/blob/complex-normal/randomstate/randomstate.pyx#L1690

It it doesn't have any test coverage, but it does pass interactive sanity tests. In particular, when generating from the same state, a single standard normal has the same value as complex_normal(0,1,1)

import randomstate as rs
st = rs.get_state()

rs.standard_normal()
Out[14]: 1.1538228506652775

rs.set_state(st)
rs.complex_normal(0,1,1)

Out[15]: (1.1538228506652775+0j)

rs.set_state(st)
rs.complex_normal(0,[1],[1])

Out[16]: array([ 1.15382285+0.j])

rs.set_state(st)
rs.complex_normal(0,[1],[-1])

Out[16]: array([ 1.15382285+0.j])

# Here the complex is the second std normal from a common state
rs.set_state(st)
rs.complex_normal(0,1,-1)

Out[23]: 1.9065545940505344j

rs.set_state(st)
rs.standard_normal(2)

Out[24]: array([ 1.15382285,  1.90655459])

FWIW, if multivariate_normal would broadcast inputs then it would be possible to implement virtually all of this function there.