python-计算两个通道中的复杂数据的协方差矩阵(无复杂数据类型)

def cov(m, y=None):

    if m.ndimension() > 2:
        raise ValueError("m has more than 2 dimensions")

    if y.ndimension() > 2:
        raise ValueError('y has more than 2 dimensions')

    X = m
    if X.shape[0] == 0:
        return torch.tensor([]).reshape(0, 0)
    if y is not None:
        X = torch.cat((X, y), dim=0)

    ddof = 1

    avg = torch.mean(X, dim=1)

    fact = X.shape[1] - ddof

    if fact <= 0:
        import warnings
        warnings.warn("Degrees of freedom <= 0 for slice",
                      RuntimeWarning, stacklevel=2)
        fact = 0.0

    X -= avg[:, None]
    X_T = X.t()
    c = dot(X, X_T)
    c *= 1. / fact
    return c.squeeze()

现在以numpy格式,它可以透明地处理复数,但是我不能简单地使用最后一个维度为(real,imag)的3-d数组,并希望它可以工作.

如何调整计算以获得具有实数通道和虚数通道的复协方差矩阵？

cov()操作的描述

令M为HxW矩阵,其中H行中的每行对应于W个复杂观测值的变量.

现在,使cov(M)为M的H变量的HxH协方差矩阵(定义为numpy.cov()).可以如下计算(忽略边缘情况)：

cov(M) = 1 / (W - 1) . M * M.T

使用*矩阵乘法运算符,以及M的M.T转置.

注意：为了澄清下一个方程,让cov_prod(X,Y)= 1 /(W-1). X * Y.T,具有X,Y HxW矩阵.因此,我们有cov(M)= cov_prod(M,M).

到目前为止,没有什么新鲜的东西,它与您编写的代码相对应(减去y权重和边缘情况的数据检查).让我们仔细检查一下此公式的Pytorch实现是否对应于Numpy(对于实值数据)：

import torch
import numpy as np

def cov(m, y=None):
    if y is not None:
        m = torch.cat((m, y), dim=0)
    m_exp = torch.mean(m, dim=1)
    x = m - m_exp[:, None]
    cov = 1 / (x.size(1) - 1) * x.mm(x.t())
    return cov

# Real-valued matrix:
M_np = np.random.rand(3, 2)  
# Same matrix as torch.Tensor:
M = torch.from_numpy(M_np)

cov_real_np = np.cov(M_np)
cov_real = cov(M)
eq = np.allclose(cov_real_np, cov_real.numpy())
print("Numpy & Torch real covariance results equal? > {}".format(eq))
# Numpy & PyTorch real covariance results equal? > True

扩展到复杂矩阵

现在,这对复杂矩阵如何起作用？
从这里开始,令M为复杂值,即由W个复杂观测值的H行变量组成.
此外,令A和B为实值矩阵,使得M = Ai.B.

我将不进行数学演示,通过@zimzam可以找到here,但是在那种情况下,cov(M)可以分解为：

cov(M) = [cov_prod(A, A) + cov_prod(B, B)] + i.[-cov_prod(A, B) + cov_prod(B, A)]

给定M(A和B)的实部和虚部,这使得直接计算cov(M)的实部和虚部非常简单.

实施与验证方式

在下面找到优化的实现：

import torch
import numpy as np

def cov_complex(m_comp):
    # (adding further parameters such as `y` is left for exercise)
    # Supposing real and img are stored separately in the last dim:
    real, img = m_comp[..., 0], m_comp[..., 1] 
    x_real = real - torch.mean(real, dim=1)[:, None]
    x_img = img - torch.mean(img, dim=1)[:, None]
    x_real_T = x_real.t()
    x_img_T = x_img.t()
    frac = 1 / (x_real.size(1) - 1)

    cov_real = frac * (x_real.mm(x_real_T) + x_img.mm(x_img_T))
    cov_img = frac * (-x_real.mm(x_img_T) + x_img.mm(x_real_T))
    return torch.stack((cov_real, cov_img), dim=-1)

# Matrix with real/img values stored separately in last dimension:
M_np = np.random.rand(3, 2, 2)  
# Same matrix converted to np.complex format:                  
M_comp_np = M_np.view(dtype=np.complex128)[...,0]
# Same matrix as torch.Tensor:
M = torch.from_numpy(M_np)

cov_com_np = np.cov(M_comp_np)
cov_com = cov_complex(M)
eq = np.allclose(cov_com_np, cov_com.numpy().view(dtype=np.complex128)[...,0])
print("Numpy & Torch complex covariance results equal? > {}".format(eq))
# Numpy & PyTorch complex covariance results equal? > True