首页 > 解决方案 > 样本协方差的特征值远离协方差的特征值

问题描述

我编写了以下代码来创建具有指定特征值的随机矩阵 Sigma,然后从具有零均值和协方差 Sigma 的多元正态分布中采样向量。

def generate_data(N, d):
    eigenvalues = [0] * (d + 1)
    for k in range(1, d + 2):
        eigenvalues[k - 1] = k

    random_matrix = numpy.random.randn(d + 1, d + 1)
    random_orthogonal = numpy.linalg.qr(random_matrix)[0]
    sqrt_cov = random_orthogonal @ numpy.diag(numpy.sqrt(eigenvalues))

    X = numpy.zeros((N, d + 1))
    for i in range(N):
        vec = numpy.random.standard_normal(d + 1)
        X[i] = sqrt_cov @ vec

在此代码之后,X 应该是从所需分布中采样的 N x d+1 矩阵。

现在想知道X的样本协方差矩阵的特征值是多少。如果我没记错的话,应该和Sigma差不多

def get_sample_covariance(data):
    data = data - data.mean(axis=0)
    sample_cov = data.T @ data / (data.shape[0] - 1)
    return sample_cov

然后我绘制了 sample_cov 的特征值

我期望一个大致线性的函数,从 d(500)到 1。我得到了这个

在此处输入图像描述

是什么赋予了?哪里错了?

标签: pythonnumpycovarianceeigenvalue

解决方案


样本的生成似乎是正确的。但是仅当 N >> d 时,来自 N 个随机样本的协方差估计才是正确的。特别是最高特征值倾向于系统地关闭。但是,如果你从所有特征值中取出平均特征值,它是非常准确的。

def get_max_mean_eig_run(N, d):
    """Return highest and mean eigenvalue for random samples"""

    X = generate_data(N, d)
    cov = get_sample_covariance(X)
    eig = np.linalg.eigvals(cov)
    return eig.max(), eig.mean()

plt.close('all')
np.random.seed(1)
ds = np.arange(1, 501, 25)

# shape (N, 2) for max and mean eigenvalues
eigs = np.array([get_max_mean_eig_run(1000, d) for d in ds])

plt.close('all')
plt.xlabel('d')
plt.ylabel('Eigenvalue')
plt.plot(ds, eigs[:, 0], 'r-', label='Max')
plt.plot(ds, eigs[:, 1], 'r--', label='Mean')
plt.scatter(ds, ds+1, color='k', marker='o', label='max expected')
plt.scatter(ds, (ds+2)/2, color='k', marker='*', label='mean expected')
plt.legend()

推断的特征值

对于这个特定的特征值输入谱,您需要N > 100*d使最大特征值接近 (+5%) 到最大输入特征值,但对于更实际的情况,这可能会有所不同。

这是特征值的直方图(对于 N=1000,d=500)

plt.close('all')
X = generate_data(1000, 500)
eigs = np.linalg.eigvals(get_sample_covariance(X))
hist, bin_edges = np.histogram(eigs, bins=20)
bin_centers = (bin_edges[1:] + bin_edges[:-1])/2
plt.step(bin_centers, hist, where='mid')
plt.xlabel('Eigenvalue')
plt.ylabel('Count (per bin)')

特征值直方图


推荐阅读