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