# Covariance and gradient support

**URL:** <https://discuss.pytorch.org/t/covariance-and-gradient-support/16217>\
**Category:** autograd\
**Created:** [April 10, 2018, 3:27pm UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217 "2018-04-10T15:27:06Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![dem123456789](https://discuss.pytorch.org/letter_avatar_proxy/v4/letter/d/ce7236/32.png) [@dem123456789](https://discuss.pytorch.org/u/dem123456789)\
**Post date:** [April 10, 2018, 3:27pm UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/1 "2018-04-10T15:27:06Z")

</div>

Is there a proposal for torch.cov() and its gradient support?

---

<div class="post-metadata">

**Author:** ![ModarTensai](https://discuss.pytorch.org/user_avatar/discuss.pytorch.org/modartensai/32/10307_2.png) [@ModarTensai](https://discuss.pytorch.org/u/ModarTensai)\
**Post date:** [July 30, 2018, 2:17am UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/2 "2018-07-30T02:17:27Z")

</div>

The covariance of `X` of shape `[*, n]` can be easily [estimated as](https://en.wikipedia.org/wiki/Covariance_matrix#Estimation):  
`cov(X) = MM^T/(n-1)` where `M = X - columns_means(X)`.  
If we looked into the implementation of [`np.cov`](https://github.com/numpy/numpy/blob/master/numpy/lib/function_base.py#L2122-L2328), we can optimize it as follows:

```python
import warnings
import numpy as np

def np_cov(m, rowvar=False):
    # Handles complex arrays too
    m = np.asarray(m)
    if m.ndim > 2:
        raise ValueError('m has more than 2 dimensions')
    dtype = np.result_type(m, np.float64)
    m = np.array(m, ndmin=2, dtype=dtype)
    if not rowvar and m.shape[0] != 1:
        m = m.T
    if m.shape[0] == 0:
        return np.array([], dtype=dtype).reshape(0, 0)

    # Determine the normalization
    fact = m.shape[1] - 1
    if fact <= 0:
        warnings.warn("Degrees of freedom <= 0 for slice",
                      RuntimeWarning, stacklevel=2)
        fact = 0.0

    m -= np.mean(m, axis=1, keepdims=1)
    c = np.dot(m, m.T.conj())
    c *= np.true_divide(1, fact)
    return c.squeeze()

```

This can be easily tested as follows:

```python
r = False
x = np.random.rand(3, 1000)
np_c = np.cov(x, rowvar=r)
our_c = np_cov(x, rowvar=r)
print(np.allclose(np_c, our_c))

```

To port it to pytorch, I did the following:

```python
import torch

def cov(m, rowvar=False):
    '''Estimate a covariance matrix given data.

    Covariance indicates the level to which two variables vary together.
    If we examine N-dimensional samples, `X = [x_1, x_2, ... x_N]^T`,
    then the covariance matrix element `C_{ij}` is the covariance of
    `x_i` and `x_j`. The element `C_{ii}` is the variance of `x_i`.

    Args:
        m: A 1-D or 2-D array containing multiple variables and observations.
            Each row of `m` represents a variable, and each column a single
            observation of all those variables.
        rowvar: If `rowvar` is True, then each row represents a
            variable, with observations in the columns. Otherwise, the
            relationship is transposed: each column represents a variable,
            while the rows contain observations.

    Returns:
        The covariance matrix of the variables.
    '''
    if m.dim() > 2:
        raise ValueError('m has more than 2 dimensions')
    if m.dim() < 2:
        m = m.view(1, -1)
    if not rowvar and m.size(0) != 1:
        m = m.t()
    # m = m.type(torch.double) # uncomment this line if desired
    fact = 1.0 / (m.size(1) - 1)
    m -= torch.mean(m, dim=1, keepdim=True)
    mt = m.t() # if complex: mt = m.t().conj()
    return fact * m.matmul(mt).squeeze()

```

It is maybe important here to note that currently pytorch [doesn’t support complex tensors](https://github.com/pytorch/pytorch/issues/755), so the conjugate of a real matrix is itself. Another thing is that I like the pytorch convention where each function that consumes tensors and produces others outputs tensors of the same dtype. If you don’t desire this behaviour just uncomment the type casting line. Otherwise, cast your tensor to double before passing it to the function. Finally, I prefer `rowvar=False` by default which is not the default behaviour of `np.cov`.

Here is how to test this implementation:

```python
r = False
x = np.random.rand(3, 1000)
xt = torch.from_numpy(x).type(torch.double)
np_c = np.cov(x, rowvar=r)
our_c = cov(xt, rowvar=r).numpy()
print(np.allclose(np_c, our_c))

```

Because no special functions were used, we can acquire the gradient of this function with auto-differentiation.

---

<div class="post-metadata">

**Author:** ![ModarTensai](https://discuss.pytorch.org/user_avatar/discuss.pytorch.org/modartensai/32/10307_2.png) [@ModarTensai](https://discuss.pytorch.org/u/ModarTensai)\
**Post date:** [July 31, 2018, 3:55am UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/3 "2018-07-31T03:55:11Z")

</div>

I came across [this on StackOverflow](https://stackoverflow.com/questions/51416825/calculate-covariance-matrix-for-complex-data-in-two-channels-no-complex-data-ty/51422002) which asks for computing the covariance matrix for complex tensors in torch.

---

<div class="post-metadata">

**Author:** ![themightyoarfish](https://discuss.pytorch.org/letter_avatar_proxy/v4/letter/t/f07891/32.png) [@themightyoarfish](https://discuss.pytorch.org/u/themightyoarfish)\
**Post date:** [January 18, 2019, 10:30am UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/4 "2019-01-18T10:30:56Z")

</div>

Yeah the problem with that is that you can’t always work with the resulting matrices, i.e. I would have had to do fourier transform on the result, which of course doesn’t work since PyTorch doesn’t have complex yet.

---

<div class="post-metadata">

**Author:** ![ModarTensai](https://discuss.pytorch.org/user_avatar/discuss.pytorch.org/modartensai/32/10307_2.png) [@ModarTensai](https://discuss.pytorch.org/u/ModarTensai)\
**Post date:** [July 31, 2019, 11:04pm UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/5 "2019-07-31T23:04:42Z")

</div>

This is the same code with `inplace` flag.

```auto
import torch

def cov(m, rowvar=True, inplace=False):
    '''Estimate a covariance matrix given data.

    Covariance indicates the level to which two variables vary together.
    If we examine N-dimensional samples, `X = [x_1, x_2, ... x_N]^T`,
    then the covariance matrix element `C_{ij}` is the covariance of
    `x_i` and `x_j`. The element `C_{ii}` is the variance of `x_i`.

    Args:
        m: A 1-D or 2-D array containing multiple variables and observations.
            Each row of `m` represents a variable, and each column a single
            observation of all those variables.
        rowvar: If `rowvar` is True, then each row represents a
            variable, with observations in the columns. Otherwise, the
            relationship is transposed: each column represents a variable,
            while the rows contain observations.

    Returns:
        The covariance matrix of the variables.
    '''
    if m.dim() > 2:
        raise ValueError('m has more than 2 dimensions')
    if m.dim() < 2:
        m = m.view(1, -1)
    if not rowvar and m.size(0) != 1:
        m = m.t()
    # m = m.type(torch.double) # uncomment this line if desired
    fact = 1.0 / (m.size(1) - 1)
    if inplace:
        m -= torch.mean(m, dim=1, keepdim=True)
    else:
        m = m - torch.mean(m, dim=1, keepdim=True)
    mt = m.t() # if complex: mt = m.t().conj()
    return fact * m.matmul(mt).squeeze()

```

---

<div class="post-metadata">

**Author:** ![neurosurg](https://discuss.pytorch.org/user_avatar/discuss.pytorch.org/neurosurg/32/21232_2.png) [@neurosurg](https://discuss.pytorch.org/u/neurosurg)\
**Post date:** [May 12, 2020, 8:52am UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/6 "2020-05-12T08:52:02Z")

</div>

cor() and cov() for a 2D data matrix is easy to write.

However, TBH, it is not that straight-forward e.g. for a OHE-encoded tensor built from that 2D data matrix.

So, what about nice-written cor() and cov() for a tensor?

---

<div class="post-metadata">

**Author:** ![ModarTensai](https://discuss.pytorch.org/user_avatar/discuss.pytorch.org/modartensai/32/10307_2.png) [@ModarTensai](https://discuss.pytorch.org/u/ModarTensai)\
**Post date:** [April 6, 2021, 11:31pm UTC](https://discuss.pytorch.org/t/covariance-and-gradient-support/16217/7 "2021-04-06T23:31:58Z")

</div>

This topic has gone stale. Please, check the latest updates on this [issue](https://github.com/pytorch/pytorch/issues/19037).

As a sneak peek, check this implementation that supports any tensor type:

```python
def cov(tensor, rowvar=True, bias=False):
    """Estimate a covariance matrix (np.cov)
    https://gist.github.com/ModarTensai/5ab449acba9df1a26c12060240773110
    """
    tensor = tensor if rowvar else tensor.transpose(-1, -2)
    tensor = tensor - tensor.mean(dim=-1, keepdim=True)
    factor = 1 / (tensor.shape[-1] - int(not bool(bias)))
    return factor * tensor @ tensor.transpose(-1, -2).conj()

```