稀疏协方差估计的高斯变量

Judson Wilson 於 5/22/2014 创作的衍生作品。经过大量改进和修复，改编自 Joelle Skaf 於 4/24/2008 所做的 CVX 示例。

主题参考：

  • 书籍《Convex Optimization》第 7.1.1 节，作者 Boyd & Vandenberghe

介绍

假设 \(y \in \mathbf{\mbox{R}}^n\) 是一个均值为零、协方差矩阵为 \(R = \mathbf{\mbox{E}}[yy^T]\) 的高斯随机变量，且具有稀疏逆 \(S = R^{-1}`（如果 :math:`S_{ij} = 0\)，则表示 \(y_i\) 和 \(y_j\) 在条件独立下）。我们想要根据从分布中抽取的 \(N\) 个独立样本 \(y_1,\dots,y_N\) 来估计协方差矩阵 \(R\)，同时利用先验知识辅助估计 \(S\) 的稀疏性。

估计 \(R\) 的一个良好启发式方法是求解以下问题：

\[\begin{split}\begin{array}{ll} \mbox{maximize} & \log \det(S) - \mbox{tr}(SY) \\ \mbox{subject to} & \sum_{i=1}^n \sum_{j=1}^n |S_{ij}| \le \alpha \\ & S \succeq 0, \end{array}\end{split}\]

其中 \(Y\) 是 \(y_1,\dots,y_N\) 的样本协方差，\(\alpha\) 是要选择或调整的稀疏参数。

生成问题数据

import cvxpy as cp
import numpy as np
import scipy as scipy

# 修复随机数生成器以便重复实验。
np.random.seed(0)

# 矩阵的维度。
n = 10

# 样本数，y_i
N = 1000

# 创建稀疏、对称PSD矩阵S
A = np.random.randn(n, n)  # 单位正态分布。
A[scipy.sparse.rand(n, n, 0.85).todense().nonzero()] = 0  # 稀疏化矩阵。
Strue = A.dot(A.T) + 0.05 * np.eye(n)  # 强制严格正定。

# 创建与S相关的协方差矩阵R。
R = np.linalg.inv(Strue)

# 从具有协方差R的分布中创建样本y_i。
y_sample = scipy.linalg.sqrtm(R).dot(np.random.randn(n, N))

# 计算样本的协方差矩阵。
Y = np.cov(y_sample)

对于几个 \(\alpha\) 值求解

# 每个尝试生成稀疏逆协方差矩阵的alpha值。
alphas = [10, 2, 1]

# 空的结果矩阵S列表
Ss = []

# 对每个alpha值解决优化问题。
for alpha in alphas:
    # 创建一个被限制在正半定锥中的变量。
    S = cp.Variable(shape=(n,n), PSD=True)

    # 形成logdet(S) - tr(SY)目标函数。注意使用集合推导来形成S*Y的对角线元素的集合，
    # 并且使用原生的sum函数来计算迹。TODO：如果cvxpy中提供了迹运算符，请使用它！
    obj = cp.Maximize(cp.log_det(S) - sum([(S*Y)[i, i] for i in range(n)]))

    # 设置约束条件。
    constraints = [cp.sum(cp.abs(S)) <= alpha]

    # 形成和解决优化问题
    prob = cp.Problem(obj, constraints)
    prob.solve(solver=cp.CVXOPT)
    if prob.status != cp.OPTIMAL:
        raise Exception('CVXPY错误')

    # 如果需要协方差矩阵R，可以这样创建。
    R_hat = np.linalg.inv(S.value)

    # 将S元素阈值化以确保精确零值：
    S = S.value
    S[abs(S) <= 1e-4] = 0

    # 将此S存储在结果列表中以备后续绘图使用。
    Ss += [S]

    print('以alpha = {} 为参数的优化完成，目标值为 {}'.format(alpha, obj.value))

Completed optimization parameterized by alpha = 10, obj value = -16.167608186713004
Completed optimization parameterized by alpha = 2, obj value = -22.545759632606043
Completed optimization parameterized by alpha = 1, obj value = -26.989407069609157

结果图

import matplotlib.pyplot as plt

# 在 IPython 中内联显示图形。
%matplotlib inline

# 图形的属性。
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

# 创建图形。
plt.figure()
plt.figure(figsize=(12, 12))

# 绘制真实协方差矩阵的稀疏模式。
plt.subplot(2, 2, 1)
plt.spy(Strue)
plt.title('真实协方差矩阵的逆', fontsize=16)

# 对于每个特定的 alpha 值，绘制相应的稀疏模式。
for i in range(len(alphas)):
    plt.subplot(2, 2, 2+i)
    plt.spy(Ss[i])
    plt.title('估计的逆协方差矩阵，$\\alpha$={}'.format(alphas[i]), fontsize=16)

<Figure size 432x288 with 0 Axes>
../../_images/sparse_covariance_est_5_1.png