Consider adding the complex multivariate normal distribution

[NeilGirdhar]

Please consider adding the complex multivariate normal distribution.

A helpful reference: Picinbono, B. (1996). Second-order complex random vectors and normal distributions. IEEE Transactions on Signal Processing, 44(10), 2637–2640. http://doi.org/10.1109/78.539051

The Wikipedia page is enough to fill in most of what is needed for a distribution.

I have added code that implements rvs and pdf below along with some tests.

[NeilGirdhar]

import numpy as np
from scipy.stats import multivariate_normal

__all__ = ['ScipyComplexNormal']


class ScipyComplexNormal:

    def __init__(self, mean, covariance, pseudo_covariance):
        self.mean = mean
        self.covariance = covariance
        self.pseudo_covariance = pseudo_covariance
        #assert mean.shape == covariance.shape == pseudo_covariance.shape == ()

    @classmethod
    def init_using_angle(cls, mean, covariance, angle, polarization):
        r = polarization *  np.exp(1j * 2 * np.pi * angle * 2)
        return cls(mean, covariance, r * covariance)

    def scale_and_natural_parameters(self):
        r = self.pseudo_covariance.conjugate() / self.covariance
        p = self.covariance.conjugate() - r * self.pseudo_covariance
        p_inv_c = (1.0 / p).conjugate()
        scale = 1.0 / (np.pi
                       * np.sqrt(self.covariance.real * p.real))
        return (scale,
                2.0 * (p_inv_c * self.mean.conjugate()
                       - r * p_inv_c * self.mean),
                -p_inv_c,
                r * p_inv_c)

    def pdf(self, z, out=None):
        it = np.nditer([z, out], [],
                       [['readonly'], ['writeonly', 'allocate']])
        scale, eta, H, J = self.scale_and_natural_parameters()
        for (z, out_i) in it:
            out_i[...] = self._unbroadcasted_density(z, scale, eta, H, J)
        return it.operands[1].real

    def _unbroadcasted_density(self, z, scale, eta, H, J):
        return scale * np.exp(
            (eta * z).real
            + (z.conjugate() * H * z).real
            + (z * J * z).real)

    def rvs(self, size=1, random_state=None):
        cov_sum = self.covariance + self.pseudo_covariance
        cov_diff = self.covariance - self.pseudo_covariance
        xx = 0.5 * cov_sum.real
        xy = 0.5 * -cov_diff.imag
        yx = 0.5 * cov_sum.imag
        yy = 0.5 * cov_diff.real
        cov = np.array([[xx, xy], [yx, yy]])

        xy_rvs = multivariate_normal.rvs(size=size, random_state=random_state,
                                         cov=cov)
        return xy_rvs[..., 0] + 1j * xy_rvs[..., 1] + self.mean

[NeilGirdhar]

from math import log, pi

import numpy as np
from numpy.linalg import det, inv

__all__ = ['ScipyComplexMultivariateNormal']


class ScipyComplexMultivariateNormal:

    def __init__(self, mean, variance, pseudo_variance):
        self.size = mean.shape[0]
        self.mean = mean
        self.variance = variance
        self.pseudo_variance = pseudo_variance
        if mean.shape != (self.size,):
            raise ValueError(f"Mean has shape {mean.shape} "
                             f"instead of {(self.size,)}.")
        if variance.shape != (self.size, self.size):
            raise ValueError(
                f"The variance has shape "
                f"{variance.shape} "
                f"instead of {(self.size, self.size)}.")
        if pseudo_variance.shape != (self.size, self.size):
            raise ValueError(
                f"The pseudo-variance has shape "
                f"{pseudo_variance.shape} "
                f"instead of {(self.size, self.size)}.")
        if not np.all(np.linalg.eigvals(variance) >= 0):
            raise ValueError(
                "The variance is not positive semidefinite.")
        if not np.allclose(variance, variance.T.conj()):
            raise ValueError(
                "The variance is not Hermitian.")
        if not np.allclose(pseudo_variance, pseudo_variance.T):
            raise ValueError(
                "The pseudo-variance is not symmetric.")

    # New methods -------------------------------------------------------------
    def log_normalizer(self):
        mu = self.mean
        _, precision, pseudo_precision = self.natural_parameters()
        det_s = det(self.variance).real
        det_h = det(-precision).real
        return (-mu.conj().dot(precision.dot(mu)).real
                - mu.dot(pseudo_precision.dot(mu)).real
                + 0.5 * log(det_s)
                - 0.5 * log(det_h)
                + self.size * log(pi)).real

    def natural_parameters(self):
        r, p_c = self._r_and_p_c()
        p_inv_c = inv(p_c)
        precision = -p_inv_c
        pseudo_precision = r.T @ p_inv_c
        eta = -2.0 * (precision.T.dot(self.mean.conj())
                      + pseudo_precision.dot(self.mean))
        return (eta, precision, pseudo_precision)

    def natural_to_sample(self, eta, precision, pseudo_precision):
        inv_precision = inv(precision)
        r = -(pseudo_precision @ inv_precision).T
        s = inv(r.conj() @ r - np.eye(self.size)) @ inv_precision
        u = (r @ s).conj()
        k = inv_precision.T @ pseudo_precision
        l_eta = 0.5 * inv(k @ k.conj()
                          - np.eye(self.size)).dot(inv_precision.T.dot(eta))
        mu = (l_eta.conj()
              - (inv_precision.T @ pseudo_precision).conj().dot(l_eta))
        return mu, s, u

    def pdf(self, z):
        log_normalizer = self.log_normalizer()
        eta, precision, pseudo_precision = self.natural_parameters()
        return np.exp((eta.dot(z)
                      + z.conj().dot(precision).dot(z)
                      + z.dot(pseudo_precision).dot(z)).real
                      - log_normalizer)

    def rvs(self, size=(), random_state=None):
        if random_state is None:
            random_state = np.random
        xy_rvs = random_state.multivariate_normal(
            mean=self._multivariate_normal_mean(),
            cov=self._multivariate_normal_cov(),
            size=size)
        return xy_rvs[..., :self.size] + 1j * xy_rvs[..., self.size:]

    # Private methods ---------------------------------------------------------
    def _r_and_p_c(self):
        r = self.pseudo_variance.conj().T @ inv(self.variance)
        p_c = self.variance - (r @ self.pseudo_variance).conj()
        return r, p_c

    def _multivariate_normal_mean(self):
        return np.concatenate([self.mean.real, self.mean.imag])

    def _multivariate_normal_cov(self):
        cov_sum = self.variance + self.pseudo_variance
        cov_diff = self.variance - self.pseudo_variance
        xx = 0.5 * cov_sum.real
        xy = 0.5 * -cov_diff.imag
        yx = 0.5 * cov_sum.imag
        yy = 0.5 * cov_diff.real
        return np.block([[xx, xy], [yx, yy]])

[NeilGirdhar]

import pytest
import numpy as np
from numpy.testing import assert_allclose

from . import ScipyComplexMultivariateNormal, ScipyComplexNormal


class Test:

    def setup_method(self):
        self.rng = np.random.RandomState(123)

    def test_complex_normal_rvs(self):
        mean = 2+1j
        covariance = 1.9
        angle = 0.05
        polarization = 0.98
        dist = ScipyComplexNormal.init_using_angle(
            mean, covariance, angle, polarization)

        size = (200, 100)
        rvs = dist.rvs(random_state=self.rng, size=size)

        assert rvs.shape == size

        estimated_mean = np.average(rvs)

        centered_rvs = rvs - dist.mean

        estimated_covariance = np.average(
            centered_rvs.conjugate() * centered_rvs).real
        estimated_pseudo_covariance = np.average(np.square(centered_rvs))

        assert_allclose(estimated_mean, dist.mean, rtol=1e-2)
        assert_allclose(estimated_covariance, dist.covariance, rtol=1e-2)
        assert_allclose(estimated_pseudo_covariance, dist.pseudo_covariance,
                        rtol=1e-2)

    def test_multivariate_complex_normal_rvs(self):
        mean = np.array([2+1j, 1+3j])
        zs = np.array([
            [1+3j, 2-1j],
            [3-7j, 6.0],
            [9+3j, -3-2j],
                       ])
        weights = np.array([1, 4, 2])

        covariance = np.average([z.conjugate() * z[:, np.newaxis]
                                 for z in zs],
                                weights=weights,
                                axis=0)
        pseudo_covariance = np.average([z * z[:, np.newaxis]
                                        for z in zs],
                                       weights=weights,
                                       axis=0) * 0.98

        dist = ScipyComplexMultivariateNormal(
            mean, covariance, pseudo_covariance)

        size = (500, 700)
        axis = tuple(range(2))
        rvs = dist.rvs(random_state=self.rng, size=size)

        assert rvs.shape == size + mean.shape

        estimated_mean = np.average(rvs, axis=axis)

        centered_rvs = rvs - dist.mean

        estimated_covariance = np.average(
            centered_rvs[..., np.newaxis]
            * centered_rvs.conjugate()[..., np.newaxis, :],
            axis=axis)
        estimated_pseudo_covariance = np.average(
            centered_rvs[..., np.newaxis]
            * centered_rvs[..., np.newaxis, :],
            axis=axis)

        assert_allclose(estimated_mean, dist.mean, rtol=1e-2, atol=1e-2)
        assert_allclose(estimated_covariance, dist.covariance, rtol=1e-2)
        assert_allclose(estimated_pseudo_covariance, dist.pseudo_covariance,
                        rtol=1e-2)

[rgommers]

Thanks @NeilGirdhar , looks like an interesting distribution to add. Could you submit your code as a pull quest?

The code looks good to me from a quick browse (I didn't test it). Only a couple of comments:

• Class should be named something like multivariate_complex_normal_gen for consistency with multivariate_normal_gen.
• Class should inherit from multi_rv_generic
• Parameter covariance should be called cov and density should be called pdf, also for consistency.
• Don't use any plain asserts. For testing use functions from numpy.testing, and in the main code use if ... checks where needed and remove the rest.
• The circularly-symmetric normal distribution has some properties that should hold that would be good to test for.

[NeilGirdhar]

@rgommers Thanks for the feedback. I'm very busy with school and was initially thinking that I would leave this code in scipy's capable hands. However, I think I can find some time to get this in if you don't mind helping minimize the back and forth. (When I worked on PEP 448, it took a lot longer than I anticipated although not through anyone's fault.)

Would you mind commenting on my use of nditer? What is the correct way to write pdf? Also, I couldn't figure out what to do about dtypes. I want nditer to allocate out to be np.float64 when z is np.complex128 (or np.float32 when z is np.complex64, etc.)

How would you test circularly-symmetric normal? It sounds like you know, but just so we're on the same page, it is the special case of multivariate normal with zero pseudo-covariance and mean.

[rgommers]

A SciPy PR is nowhere near as bad as a PEP. Not much in Python land is as bad as a PEP to be fair:)

• Would you mind commenting on my use of nditer?

That's not a common patterns at all, but it looks better than some of the approaches to handling broadcasting in the stats distributions. It's very similar to np.broadcast_to, could you use that instead? If so, use scipy._lib._numpy_compat.broadcast_to for compatibility with numpy <1.10.

The out=None argument to density complicates things, and probably shouldn't be there in the first place because none of the other distributions have it.

• What is the correct way to write pdf?

It should have all parameters that __init__ has, so it looks like the signature should be pdf(self, z, mean, cov, pseudo_cov).

• I want nditer to allocate out to be np.float64 when z is np.complex128 (or np.float32 when z is np.complex64, etc.)

Maybe make your life easy and always return float64?

How would you test circularly-symmetric normal. It sounds like you know, but just so we're on the same page, it is the special case of multivariate normal with zero pseudo-covariance and mean.

Yes, I know. I meant that in testing you can create that distribution and then check for properties that it should have, e.g. "The density function depends only on the magnitude of z but not on its argument"
and "the standard complex normal distribution has density ... [easy to implement expression in test]".

[NeilGirdhar]

Thank you. Good idea regarding the circular test.

Just to confirm. I'm writing two distributions: multivariate_complex_normal_gen and complex_normal_gen? I don't see how to roll them into one without explicitly having a special case in most methods.

[rgommers]

I'm writing two distributions

Ah, I missed that so far, but yes indeed. Then I have to amend by comment about inheritance: the multivariate distribution should inherit from multi_rv_generic, but the univariate one I'm now unsure about. Normally you'd inherit from rv_continuous, but you'd have to check which methods make sense to add for a distribution with complex input. Some of the generic methods may break.

[NeilGirdhar]

Okay, I think I have a first draft for one of the distributions. I'm having trouble making it work with complex arguments.

class complex_normal_gen(rv_continuous):
    """A complex normal continuous random variable.

    The location (loc) keyword specifies the mean.
    The variance keyword specifies the variance: E(abs(x)^2).
    The pseudo_variance keyword specifies the pseudo-variance: E(x^2).

    %(before_notes)s

    Notes
    -----
    Further reading about the complex normal distribution can be found in:
    Picinbono, B. (1996). Second-order complex random vectors and normal
    distributions. IEEE Transactions on Signal Processing, 44(10), 2637–2640.

    %(after_notes)s

    %(example)s

    """
    shapes = 'covariance pseudo_covariance'

    def __init__(self, covariance, pseudo_covariance, *args, **kwargs):
        super(complex_normal_gen, self).__init__(*args, **kwargs)
        self.covariance = covariance
        self.pseudo_covariance = pseudo_covariance
        self._ctor_param.update(dict(
            covariance=covariance,
            pseudo_covariance=pseudo_covariance))

    def _argcheck(self, covariance, pseudo_covariance):
        self.covariance = covariance
        self.pseudo_covariance = pseudo_covariance
        return True

    def _rvs(self, covariance, pseudo_covariance):
        cov = self._multivariate_normal_covariance(covariance,
                                                   pseudo_covariance)
        xy_rvs = self._random_state.multivariate_normal(mean=np.zeros(2),
                                                        cov=cov,
                                                        size=self._size)
        return xy_rvs[..., 0] + 1j * xy_rvs[..., 1]

    def _pdf(self, x, covariance, pseudo_covariance):
        ln, nat = self._log_normalizer_and_natural_parameters()
        ss = self._sufficient_statistics(x)
        return np.exp(sum((ss[i] * nat[i]).real for i in range(1, 3)) - ln)

    # When the multivariate_normal implements cdf, the cdf and related
    # functions will be easy to implement using the helper function.
    def _entropy(self, covariance, pseudo_covariance):
        from ._multivariate import multivariate_normal
        cov =self._multivariate_normal_covariance(covariance,
                                                  pseudo_covariance)
        return multivariate_normal(cov=cov).entropy()

    @inherit_docstring_from(rv_continuous)
    def fit(self, data, **kwds):
        data = np.asarray(data)
        axis = tuple(range(len(data.shape)))
        exp_parms = np.average(self._sufficient_statistics(data),
                               axis=axis)
        x = exp_parms[0]
        return (x,
                exp_parms[1] - x * x.conjugate(),
                exp_parms[2] - x * x)

    # Helpers
    def _multivariate_normal_covariance(self, covariance, pseudo_covariance):
        cov_sum = covariance + pseudo_covariance
        cov_diff = covariance - pseudo_covariance
        xx = 0.5 * cov_sum.real
        xy = 0.5 * -cov_diff.imag
        yx = 0.5 * cov_sum.imag
        yy = 0.5 * cov_diff.real
        return np.array([[xx, xy], [yx, yy]])

    def _log_normalizer_and_natural_parameters(self):
        r = self.pseudo_covariance.conjugate() / self.covariance
        p = self.covariance.conjugate() - r * self.pseudo_covariance
        p_inv_c = (1.0 / p).conjugate()
        log_normalizer = log(np.pi
                             * np.sqrt(self.covariance.real * p.real))
        return (log_normalizer,
                [0.0,
                 -p_inv_c,
                 r * p_inv_c])

    @staticmethod
    def _sufficient_statistics(x):
        return x, x * x.conjugate(), x * x


complex_normal = complex_normal_gen(covariance=1.0, pseudo_covariance=0.0,
                                    name='complex_normal')

[rgommers]

You put it in scipy/stats/_continuous_distns.py. Then you add it to:

• the __all__ dict in that file,
• the dists list at the top of scipy/stats/tests/test_distributions.py
• the distcont list in scipy/stats/_distr_params.py
• to the listing in the docstring of scipy/stats/__init__.py

The you can run the tests by from scipy import stats; stats.test(). You can also submit as a PR at that point, and the tests will be run automatically on TravisCI.

[rgommers]

That's quite a few places to add the distribution to, we should write that down somewhere ....

[NeilGirdhar]

@rgommers Thanks. It looks like it's starting to work, but it has trouble accepting complex-valued parameters. Would you mind taking a look at the distribution infrastructure to see how we could accept complex-valued parameters?

[rgommers]

Can you either point me to the branch you're working in, or send a WIP (work-in-progress) PR? Then I'm happy to have a look.

[NeilGirdhar]

I think it's this but I might have to put it into a new fork? NeilGirdhar@f4d6f04

[NeilGirdhar]

@rgommers I've separated out the task of improving the distribution infrastructure into a separate issue #7311.

[rgommers]

@rgommers I've separated out the task of improving the distribution infrastructure into a separate issue #7311.

thanks

I think it's this but I might have to put it into a new fork? NeilGirdhar/scipy@f4d6f04

Sorry, I've been distracted so missed that reply. Yes, you've been committing on your master branch. It's better to create a new branch and send a WIP PR from there.

[rgommers]

Looking at complex_normal_gen I'm a bit confused as to why it has a _multivariate_normal_covariance method. I'd expect that to only show up in the multivariate complex normal distribution, which this class inherits from rv_continuous and is hence univariate.

[NeilGirdhar]

I'd expect that to only show up in the multivariate complex normal distribution…

That method converts the univariate complex normal variance and pseudo-variance into a multivariate normal distribution's variance matrix. That's avoids reimplementing methods like entropy.