I just did this:
def cov(x: Tensor, y: Optional[Tensor] = None) -> Tensor:
    """
    Compute covariance of input

    :param x: [M, D]
    :param y: [M, D]
    :return:
    """
    if y is not None:
        y = x
    else:
        assert x.size(0) == y.size(0)
    # - center first
    x = center(x, dim=0)
    y = center(y, dim=0)
    # - conv = E[XY] is outer product of X^T Y or X Y^T depending on shapes
    sigma_xy: Tensor = x.T @ y
    return sigma_xy

should be the same no?