I would like to add a complex random normal generator to
mtrand/RandomState.
A scalar complex normal is a (double) bivariate normal. The main
motivation is to simplify the construction of complex normals where are
generally parameterized in terms of three values: location, covariance and
relation. location is the same as in a standard normal. The covariance
and the relation jointly determine the variance of the real and imaginary
parts as well as the covariance between the two.
#1 The initial implementation in the PR has followed the standard template
for scalar RV generators with three paths, scalar->scalar, scalar->array
and array->array. It is bulky since the existing array fillers that handle
the scalar->array and array->array for double rvs cannot be used. It
supports broadcasting and has a similar API to other scalar RV generators
(e.g. normal).
#2 The PR discussion has moved towards exploiting the relationship with the
multivariate normal. Currently the MV normal doesn't broadcast, and so
following this path would only allow the scalar->scalar and the
scalar->array paths. This could theoretically be extended to allow
broadcasting if multivariate_normal was extended to allow broadcasting.
#3 If broadcasting is off the table, then it might make more sense to skip
a scalar complex normal and just move directly to a
multivariate_complex_normal since this is also just a higher dimension
(double) multivariate normal. This function could just wrap
multivariate_normal and would be relatively straightforward forward. The
only down side of this path is that it would not easily support a
scalar->scalar path, although this could be added.
This probably isn't much of a performance hit for following #2 or #3. I
checked how normal and multivariate normal perform for large draws:
%timeit np.random.normal(2.0,4.0,size=1000000)
30.8 ms ± 125 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit np.random.multivariate_normal([2.0],[[4.0]],size=1000000)
32.2 ms ± 308 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
For smaller draws the performance difference is larger:
%timeit np.random.normal(2.0,4.0,size=10)
2.95 µs ± 16.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit np.random.multivariate_normal([2.0],[[4.0]],size=10)
49.4 µs ± 249 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
And for scalars the scalar path is about 30x faster than the multivariate
path. It is also worth noting that multivariate_normal will only return a
vector even if the inputs only generate a single scalar.
%timeit np.random.normal(2.0,4.0)
1.42 µs ± 3.05 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit np.random.multivariate_normal([2.0],[[4.0]])
47.9 µs ± 167 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
It would be helpful do determine which path,
#1 Clone standard scalar generator with 3 paths including broadcasting
#2 Scalar generator using multivariate normal, excluding broadcasting
#3 Multivariate generator using multivariate normal, excluding broadcasting
is the preferred one.
Kevin