计算科学 - 白化变换不返回单位协方差矩阵 - 吾爱随笔录

白化变换不返回单位协方差矩阵

计算科学线性代数 Python svd

2021-12-01 05:17:50

认为 $X$ 是具有非奇异协方差矩阵的随机（列）向量 $\Sigma$ 并且均值为 0。那么变换 $Y=WX$ 带有美白矩阵 $W$ 满足条件 $W^TW=\Sigma^{-1}$ 产生白化的随机向量 $Y$ 具有单位对角协方差。

根据定义，我期望协方差矩阵 $Y$ 为单位矩阵。然而，这远非事实！

这是复制品：

import numpy as np
# random matrix
dim1 = 512 # dimentionality_of_features
dim2 = 100 # no_of_samples

X = np.random.rand(dim1, dim2)
# centering to have mean 0
X = X - np.mean(X, axis=1, keepdims=True)

# covariance of X
Xcov = np.dot(X, X.T) / (X.shape[1] - 1)

# SVD decomposition
# Eigenvecors and eigenvalues
Ec, wc, _ = np.linalg.svd(Xcov)
# get only the first positive ones (for numerical stability)
k_c = (wc > 1e-5).sum()
# Diagonal Matrix of eigenvalues
Dc = np.diag((wc[:k_c]+1e-6)**-0.5)
# E D ET should be the whitening matrix
W = Ec[:,:k_c].dot(Dc).dot(Ec[:,:k_c].T)

# SVD decomposition End

Y = W.dot(X)
# Now apply the same to the whitened X
Ycov = np.dot(Y, Y.T) / (Y.shape[1] - 1)
print(Ycov)

>> [[ 0.19935189 -0.00740203 -0.00152036 ...  0.00133161 -0.03035149
      0.02638468]  ...

似乎它不会给我一个单位对角矩阵，除非，dim2 >> dim1。

如果我接受dim2=1，那么我会得到一个向量（尽管在示例中由于除以 0 会导致错误），并且根据 Wiki 的定义，这是不正确的吗？

1个回答

正如评论所指出的，您可能对协方差和样本协方差有些混淆。但是，这不是导致您的错误的原因。

首先，忘记获得协方差 $Y$ 是单位矩阵（“单位对角线”）。出于排名原因，只有当 $XX^T$ 是满秩的，但事实并非如此，因为您处于一个不寻常的设置中，其中样本数小于维数： $XX^T$ 最多dim2是 , 的等级也是 $YY^T=(WX)(WX)^T$ .

我不知道你在哪里找到了白化矩阵的算法，但我觉得它很混乱。有一些步骤可以解决排名不足（整个k_c业务），但我不确定它们是否真的有用；在我看来，只要k_c不等于dim1你就会得到错误的结果：如果我没记错的话，公式简化为 $YY^T/(n-1) = EDE^T$ ，在哪里 $E$ 是正交的，并且 $D$ 是一个对角线，有k_c一个，其余的都是零；什么时候 $D$ 包含零，这个矩阵不是对角线（除非 $E$ 很特别）。

然而，即使 $XX^T$ 没有满级，你仍然可以得到 $Y$ 对角线协方差等于 $D$ , 对角线上有 1 和 0。为此，您只需使用类似的公式，但没有第一个因素 $E$ ：

Dc = np.zeros_like(Ec)
Dc[:k_c,:k_c] = np.diag((wc[:k_c])**-0.5)
W = Dc.dot(Ec.T)

我还删除了+1e-6，这似乎是多余的（甚至是有害的），因为无论如何您已经截断了小特征值。

最后一点，作为从事数值线性代数工作的人，我不能避免指出svd（ $XX^T$ ) 不是计算奇异值和向量的最数值稳定的方法 $X$ ，但这又不是这里的主要问题。改为使用Ec, wc, _ = np.linalg.svd(X); wc = wc / math.sqrt(X.shape[1]-1)。

其它你可能感兴趣的问题

上一篇多重网格方法如何处理吉布斯现象？下一篇显式中点法的稳定区域