Sparse covariance estimation for Gaussian variables =================================================== A derivative work by Judson Wilson, 5/22/2014. Adapted (with significant improvements and fixes) from the CVX example of the same name, by Joelle Skaf, 4/24/2008. Topic References: - Section 7.1.1, Boyd & Vandenberghe “Convex Optimization” Introduction ------------ Suppose :math:y \in \mathbf{\mbox{R}}^n is a Gaussian random variable with zero mean and covariance matrix :math:R = \mathbf{\mbox{E}}[yy^T], with sparse inverse :math:S = R^{-1} (:math:S_{ij} = 0 means that :math:y_i and :math:y_j are conditionally independent). We want to estimate the covariance matrix :math:R based on :math:N independent samples :math:y_1,\dots,y_N drawn from the distribution, and using prior knowledge that :math:S is sparse A good heuristic for estimating :math:R is to solve the problem .. math:: \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} where :math:Y is the sample covariance of :math:y_1,\dots,y_N, and :math:\alpha is a sparsity parameter to be chosen or tuned. Generate problem data --------------------- .. code:: python import cvxpy as cp import numpy as np import scipy as scipy # Fix random number generator so we can repeat the experiment. np.random.seed(0) # Dimension of matrix. n = 10 # Number of samples, y_i N = 1000 # Create sparse, symmetric PSD matrix S A = np.random.randn(n, n) # Unit normal gaussian distribution. A[scipy.sparse.rand(n, n, 0.85).todense().nonzero()] = 0 # Sparsen the matrix. Strue = A.dot(A.T) + 0.05 * np.eye(n) # Force strict pos. def. # Create the covariance matrix associated with S. R = np.linalg.inv(Strue) # Create samples y_i from the distribution with covariance R. y_sample = scipy.linalg.sqrtm(R).dot(np.random.randn(n, N)) # Calculate the sample covariance matrix. Y = np.cov(y_sample) Solve for several :math:\alpha values --------------------------------------- .. code:: python # The alpha values for each attempt at generating a sparse inverse cov. matrix. alphas = [10, 2, 1] # Empty list of result matrixes S Ss = [] # Solve the optimization problem for each value of alpha. for alpha in alphas: # Create a variable that is constrained to the positive semidefinite cone. S = cp.Variable(shape=(n,n), PSD=True) # Form the logdet(S) - tr(SY) objective. Note the use of a set # comprehension to form a set of the diagonal elements of S*Y, and the # native sum function, which is compatible with cvxpy, to compute the trace. # TODO: If a cvxpy trace operator becomes available, use it! obj = cp.Maximize(cp.log_det(S) - sum([(S*Y)[i, i] for i in range(n)])) # Set constraint. constraints = [cp.sum(cp.abs(S)) <= alpha] # Form and solve optimization problem prob = cp.Problem(obj, constraints) prob.solve(solver=cp.CVXOPT) if prob.status != cp.OPTIMAL: raise Exception('CVXPY Error') # If the covariance matrix R is desired, here is how it to create it. R_hat = np.linalg.inv(S.value) # Threshold S element values to enforce exact zeros: S = S.value S[abs(S) <= 1e-4] = 0 # Store this S in the list of results for later plotting. Ss += [S] print('Completed optimization parameterized by alpha = {}, obj value = {}'.format(alpha, obj.value)) .. parsed-literal:: 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 Result plots ------------ .. code:: python import matplotlib.pyplot as plt # Show plot inline in ipython. %matplotlib inline # Plot properties. plt.rc('text', usetex=True) plt.rc('font', family='serif') # Create figure. plt.figure() plt.figure(figsize=(12, 12)) # Plot sparsity pattern for the true covariance matrix. plt.subplot(2, 2, 1) plt.spy(Strue) plt.title('Inverse of true covariance matrix', fontsize=16) # Plot sparsity pattern for each result, corresponding to a specific alpha. for i in range(len(alphas)): plt.subplot(2, 2, 2+i) plt.spy(Ss[i]) plt.title('Estimated inv. cov matrix, $\\alpha$={}'.format(alphas[i]), fontsize=16) .. parsed-literal::