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