Student seminar talk on matrix 2-norm as a dual linear program

cvxopt_spectral_norm.ipynb

rfp_low_rank.ipynb

sdp_matrix_dual.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Find the dual program for the Nuclear norm\n",
    "\n",
    "We seek the check a [claim][1] involving the dual program corresponding to a primal SDP for computing the nuclear norm\n",
    "$$\\|X\\|_{\\ast} = \\sigma_1 + \\sigma_2 + \\cdots + \\sigma_r.$$\n",
    "\n",
    "[1]: http://www.eecs.berkeley.edu/~brecht/papers/07.rfp.lowrank.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SDP Duality\n",
    "\n",
    "From El Ghaoui [notes][1], we write the primal in standard form\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{maximize}}_{X} & \\operatorname{Tr}\\left(C^T X\\right) \\\\\n",
    "\\mbox{subject to} & X = X^T \\succeq 0 \\\\\n",
    "    & C^T = C, \\\\\n",
    "    & \\operatorname{Tr}\\left(A_i^T X\\right) = b_i, \\quad i = 1, \\ldots, m \\\\\n",
    "    & A_i^T = A_i, \\quad i = 1, \\ldots, m\n",
    "\\end{array}$$\n",
    "\n",
    "[1]: http://www.eecs.berkeley.edu/~wainwrig/ee227a/SDP_Duality.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SDP Duality: Developing the dual\n",
    "\n",
    "Glossing over some details, the dual is\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{\\nu} & \\nu^T b \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\sum_{i = 1}^m \\nu_i A_i \\succeq C\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying Duality\n",
    "\n",
    "In [showing][1] that the operator norm $\\|X\\|_{\\ast}$ is dual to the spectral norm $\\|X\\|_2$, we represent the dual norm as\n",
    "$$\\|X\\|_{2,\\text{dual}} = \\max_Y \\left\\{\\left\\langle X, Y\\right\\rangle_F \\mid  \\|Y\\|_2 \\leq 1\\right\\}$$\n",
    "where (as we've used above)\n",
    "$$\\left\\langle X, Y\\right\\rangle_F = \\operatorname{Tr}\\left(X^T Y\\right).$$\n",
    "\n",
    "This can be recast as the SDP\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{maximize}}_{Y} & \\operatorname{Tr}\\left(X^T Y\\right) \\\\\n",
    "\\mbox{subject to} \n",
    "    & \\left[ \\begin{array}{c c} I_m & Y \\\\ Y^{T} & I_n \\end{array}\\right] \\succeq 0\n",
    "\\end{array}$$\n",
    "As we've shown [elsewhere][2], we can realize a Schur complement via the [matrix congruence][3]\n",
    "$$\\left[ \\begin{array}{c c}\n",
    "I_m & Y \\\\\n",
    "Y^T & I_n \\end{array}\\right] = \\left[ \\begin{array}{c c}\n",
    "I_m & 0 \\\\\n",
    "Y^T & I_n \\end{array} \\right]^T \\left[ \\begin{array}{c c}\n",
    "I_m - Y Y^T & 0 \\\\\n",
    "0 & I_n \\end{array} \\right] \\left[ \\begin{array}{c c}\n",
    "I_m & 0 \\\\\n",
    "Y^T & I_n \\end{array}\\right]$$\n",
    "so the matrix in question is positive semidefinite provided the block diagonal is which occurs precisely when $\\|Y\\|_2 \\leq 1$.\n",
    "\n",
    "[1]: http://www.eecs.berkeley.edu/~brecht/papers/07.rfp.lowrank.pdf\n",
    "[2]: rfp_low_rank.ipynb\n",
    "<!---\n",
    "The link in [2] is a relative path, i.e. the file is expected to be in the same directory as this one.\n",
    "-->\n",
    "[3]: http://en.wikipedia.org/wiki/Matrix_congruence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying Duality: A Generalized Dual\n",
    "\n",
    "The dual of the previous SDP is\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{W_1, W_2} & \\frac{1}{2}\\left(\\operatorname{Tr}\\left(W_1\\right) +\n",
    "\\operatorname{Tr}\\left(W_2\\right)\\right) \\\\\n",
    "\\mbox{subject to} \n",
    "    & \\left[ \\begin{array}{c c} W_1 & X \\\\ X^{T} & W_2 \\end{array}\\right] \\succeq 0 \\\\\n",
    "    & W_1 = W_1^T, \\quad W_2 = W_2^T\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How does the dual above fit the typical framework?\n",
    "\n",
    "We recast the primal as\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{maximize}}_{\\overline{Y}} & \\frac{1}{2} \\left\\langle \\overline{X}, \\overline{Y} \\right\\rangle_F \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\overline{Y} = \\overline{Y}^T \\succeq 0 \\\\\n",
    "    & \\overline{Y}_{1:m, 1:m} = I_m \\\\\n",
    "    & \\overline{Y}_{(m + 1):(m + n), (m + 1):(m + n)} = I_n \n",
    "\\end{array}$$\n",
    "We don't need the condition that the upper-right and lower-left blocks are $Y$ and $Y^T$ since $\\overline{Y}$ symmetric gives us that for free.\n",
    "\n",
    "This requires that\n",
    "$$\\overline{X} = \\left[ \\begin{array}{c c} 0 & X \\\\ X^{T} & 0 \\end{array}\\right]$$\n",
    "which leads to\n",
    "$$\\left\\langle \\overline{X}, \\overline{Y} \\right\\rangle_F = \n",
    "\\operatorname{Tr}\\left(\\left[ \\begin{array}{c c} \n",
    "X Y^T & X \\\\ X^{T} & X^T Y \\end{array}\\right]\\right) =\n",
    "2 \\left\\langle X, Y \\right\\rangle_F.$$\n",
    "\n",
    "Most importantly for this \"translation\", the on the blocks of $\\overline{Y}$ are all linear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aside\n",
    "\n",
    "Letting $E_{ij}$ be the matrix with a $1$ in row $i$, column $j$, we see that\n",
    "$$\\left\\langle E_{ij}, M \\right\\rangle_F = \\operatorname{Tr}\\left(E_{ji}M\\right) =\n",
    "\\operatorname{Tr}\\left(\n",
    "\\left[ \\begin{array}{c c c c c}\n",
    "\\mathbf{0} & \\cdots & e_j & \\cdots & \\mathbf{0}\n",
    "\\end{array}\\right]\n",
    "\\left[ \\begin{array}{c}\n",
    "M_{1}^T \\\\ \\vdots \\\\ M_{m}^T\n",
    "\\end{array}\\right]\n",
    "\\right) = \\operatorname{Tr}\\left(e_j M_i^T\\right) = M_{ij}$$\n",
    "hence $\\left\\langle \\cdot, \\cdot \\right\\rangle_F$ acts as a standard dot product when viewing matrices as vectors in $\\mathbf{R}^{mn}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Back to Primal\n",
    "\n",
    "The conditions on the blocks of $\\overline{Y}$ can be seen as\n",
    "$$\\begin{array}{l}\n",
    "\\left\\langle E_{ii}, \\overline{Y} \\right\\rangle_F = 1 \\quad \\text{for} \\quad 1 \\leq i \\leq m + n \\\\\n",
    "\\left\\langle E_{ij} + E_{ji}, \\overline{Y} \\right\\rangle_F = 0 \\quad \\text{for} \\quad 1 \\leq i, j \\leq m, \\; i \\neq j \\\\\n",
    "\\left\\langle E_{ij} + E_{ji}, \\overline{Y} \\right\\rangle_F = 0 \\quad \\text{for} \\quad m < i, j \\leq m + n, \\; i \\neq j\n",
    "\\end{array}$$\n",
    "Since $\\overline{Y}$ is symmetric, a condition like\n",
    "$$0 = \\left\\langle E_{ij} + E_{ji}, \\overline{Y} \\right\\rangle_F = \\overline{Y}_{ij} + \\overline{Y}_{ji}$$\n",
    "actually forces\n",
    "$$\\overline{Y}_{ij} = \\overline{Y}_{ji} = 0.$$\n",
    "Also, due to symmetry, WLOG we need only consider $i < j$ for the off-diagonal conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now Find the Dual\n",
    "\n",
    "Having recast the primal as\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{maximize}}_{\\overline{Y}} & \\frac{1}{2} \\left\\langle \\overline{X}, \\overline{Y} \\right\\rangle_F \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\overline{Y} \\succeq 0 \\\\\n",
    "    & \\left\\langle E_{ii}, \\overline{Y} \\right\\rangle_F = 1 \\quad \\text{for} \\quad i = 1, \\ldots, m + n \\\\\n",
    "    & \\left\\langle E_{ij} + E_{ji}, \\overline{Y} \\right\\rangle_F = 0 \\quad \\text{for} \\quad 1 \\leq i < j \\leq m \\\\\n",
    "    & \\left\\langle E_{ij} + E_{ji}, \\overline{Y} \\right\\rangle_F = 0 \\quad \\text{for} \\quad m < i < j \\leq m + n\n",
    "\\end{array}$$\n",
    "we find the dual\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{\\overline{W}} & \\frac{1}{2}\\left(\\overline{W}_{11} + \\cdots + \\overline{W}_{m + n, m + n}\\right) \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\sum_{i = 1}^{m + n} \\overline{W}_{ii} E_{ii} + \\sum_{1 \\leq i < j \\leq m} \\overline{W}_{ij} \n",
    "    \\left(E_{ij} + E_{ji}\\right) + \\sum_{m < i < j \\leq m + n} \\overline{W}_{ij} \n",
    "    \\left(E_{ij} + E_{ji}\\right) \\succeq \\overline{X}\n",
    "\\end{array}$$\n",
    "since the only non-zero conditions (i.e. non-zero entries of the generic set of conditions $b$) correspond to the diagonal conditions $E_{ii}$.\n",
    "\n",
    "We have happily swapped out the vector $\\nu$ for a multi-indexed $\\overline{W}$ and with this, we can re-write the above as\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{W_1, W_2} & \\frac{1}{2}\\left(\\operatorname{Tr}\\left(W_1\\right) + \\operatorname{Tr}\\left(W_2\\right)\\right) \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\left[ \\begin{array}{c c} W_1 & 0 \\\\ 0 & W_2 \\end{array}\\right] \\succeq \n",
    "    \\left[ \\begin{array}{c c} 0 & X \\\\ X^{T} & 0 \\end{array}\\right] \\\\\n",
    "    & W_1 = W_1^T, \\quad W_2 = W_2^T\n",
    "\\end{array}$$\n",
    "The symmetrized $E_{ij} + E_{ji}$ relationships among the $\\overline{W}_{ij}$ is hiding in the fact that $W_1$ and $W_2$ are symmetric."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Re-writing the dual\n",
    "\n",
    "To make the final step from\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{W_1, W_2} & \\frac{1}{2}\\left(\\operatorname{Tr}\\left(W_1\\right) + \\operatorname{Tr}\\left(W_2\\right)\\right) \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\left[ \\begin{array}{c c} W_1 & -X \\\\ -X^T & W_2 \\end{array}\\right] \\succeq 0\n",
    "\\end{array}$$\n",
    "to the stated dual from the paper\n",
    "$$\\begin{array}{ll}\n",
    "\\mathop{\\text{minimize}}_{W_1, W_2} & \\frac{1}{2}\\left(\\operatorname{Tr}\\left(W_1\\right) + \\operatorname{Tr}\\left(W_2\\right)\\right) \\\\\n",
    "\\mbox{subject to}\n",
    "    & \\left[ \\begin{array}{c c} W_1 & X \\\\ X^T & W_2 \\end{array}\\right] \\succeq 0\n",
    "\\end{array}$$\n",
    "is not so difficult.\n",
    "\n",
    "Swapping $X$ out with $-X$, notice that\n",
    "$$\\begin{align*}\n",
    "\\left[ \\begin{array}{c c} -v_1^T & v_2^T \\end{array}\\right]\n",
    "\\left[ \\begin{array}{c c} W_1 & -X \\\\ -X^T & W_2 \\end{array}\\right]\n",
    "\\left[ \\begin{array}{c} -v_1 \\\\ v_2 \\end{array}\\right] &= (-v_1)^T W_1 (-v_1) + \n",
    "v_2^T W_2 v_2 - (-v_1)^T X v_2 - v_2^T X^T (-v_1) \\\\\n",
    "&= v_1^T W_1 v_1 +  v_2^T W_2 v_2 + v_1^T X v_2 + v_2^T X^T v_1.\n",
    "\\end{align*}$$\n",
    "This is equivalent to the matrix congruence\n",
    "$$\\begin{align*}\n",
    "\\left[ \\begin{array}{c c} W_1 & -X \\\\ -X^T & W_2 \\end{array}\\right] &\\sim \n",
    "\\left[ \\begin{array}{c c} -I_m & 0 \\\\ 0 & I_n \\end{array}\\right]^T\n",
    "\\left[ \\begin{array}{c c} W_1 & -X \\\\ -X^T & W_2 \\end{array}\\right]\n",
    "\\left[ \\begin{array}{c c} -I_m & 0 \\\\ 0 & I_n \\end{array}\\right] \\\\\n",
    "&= \\left[ \\begin{array}{c c} -W_1 & X \\\\ -X^T & W_2 \\end{array}\\right]\n",
    "\\left[ \\begin{array}{c c} -I_m & 0 \\\\ 0 & I_n \\end{array}\\right] \\\\\n",
    "&= \\left[ \\begin{array}{c c} W_1 & X \\\\ X^T & W_2 \\end{array}\\right].\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def get_random_X(M, N):\n",
    "    return np.random.random((M, N))\n",
    "\n",
    "\n",
    "def get_X_bar(X, M, N):\n",
    "    X_bar = np.zeros((M + N, M + N))\n",
    "    X_bar[:M, M:] = X\n",
    "    X_bar[M:, :M] = X.T\n",
    "    return X_bar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use `cvxopt` (see `cvxopt_spectral_norm.ipynb` for more details) to implement the dual program (actually the one using $-X$) and compare it to the out-of-the-box computation of the nuclear norm via summing singular values.\n",
    "\n",
    "It is doubtful we'll do better since $W_1$ and $W_2$ contain $m^2 + n^2$ values, and after accounting for symmetry\n",
    "$$\\frac{m^2 + n^2 + m + n}{2}$$\n",
    "variables in total.\n",
    "\n",
    "We'll order the variables first by row and then by column within each row:\n",
    "$$\\overline{W}_{11}, \\overline{W}_{12}, \\ldots, \\overline{W}_{1m}, \\overline{W}_{22}, \\ldots, \\overline{W}_{2m}, \\ldots, \\overline{W}_{mm}, \\overline{W}_{m + 1, m + 1}, \\ldots, \\overline{W}_{m + 1, m + n}, \\ldots$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making `cvxopt` happy\n",
    "\n",
    "To use the SDP solver from `cvxopt` we need to write\n",
    "$$\\sum \\nu_i A_i \\preceq C$$\n",
    "so we use\n",
    "$$\\left[ \\begin{array}{c c} W_1 & X \\\\ X^T & W_2 \\end{array}\\right] \\succeq 0 \\Longleftrightarrow\n",
    "\\left[ \\begin{array}{c c} W_1 & 0 \\\\ 0 & W_2 \\end{array}\\right] \\succeq \n",
    "\\left[ \\begin{array}{c c} 0 & -X \\\\ -X^T & 0 \\end{array}\\right] \\Longleftrightarrow\n",
    "\\left[ \\begin{array}{c c} -W_1 & 0 \\\\ 0 & -W_2 \\end{array}\\right] \\preceq \n",
    "\\left[ \\begin{array}{c c} 0 & X \\\\ X^T & 0 \\end{array}\\right].$$\n",
    "We represent this using the variables defined above\n",
    "$$\\sum_{i} \\overline{W}_{ii} \\left(-E_{ii}\\right) + \\sum_{i < j} \\overline{W}_{ij} \n",
    "\\left(-E_{ij} - E_{ji}\\right) \\preceq \\overline{X}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_G_row(zeros_mat, i, j):\n",
    "    # NOTE: Attempted speedup via index computations and using an already\n",
    "    #       flattened zeros matrix. Speedup was negligible.\n",
    "    # NOTE: We should use cvxopt.spmatrix.\n",
    "    G = zeros_mat.copy()\n",
    "    if i == j:\n",
    "        G[i, i] = -1.0\n",
    "    else:\n",
    "        G[i, j] = G[j, i] = -1.0\n",
    "    return G.flatten()\n",
    "\n",
    "\n",
    "# Would be useful to memoize results since the output is expensive. \n",
    "def get_coefficients(M, N):\n",
    "    all_zeros = np.zeros((M + N, M + N))\n",
    "    G_rows = []  # List of rows in single G.\n",
    "    optimized_coefficients = []  # List of c's.\n",
    "    \n",
    "    # W1, top left.\n",
    "    for i in xrange(M):\n",
    "        G_rows.append(make_G_row(all_zeros, i, i))\n",
    "        optimized_coefficients.append(0.5)\n",
    "        for j in xrange(i + 1, M):\n",
    "            G_rows.append(make_G_row(all_zeros, i, j))\n",
    "            optimized_coefficients.append(0.0)\n",
    "\n",
    "    # W2, bottom right.\n",
    "    for i in xrange(M, M + N):\n",
    "        G_rows.append(make_G_row(all_zeros, i, i))\n",
    "        optimized_coefficients.append(0.5)\n",
    "        for j in xrange(i + 1, M + N):\n",
    "            G_rows.append(make_G_row(all_zeros, i, j))\n",
    "            optimized_coefficients.append(0.0)\n",
    "\n",
    "    optimized_coefficients = cvxopt.matrix(optimized_coefficients)\n",
    "    G = cvxopt.matrix(np.vstack(G_rows).T)\n",
    "    constraint_coefficients = [G]\n",
    "    \n",
    "    return constraint_coefficients, optimized_coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this in place, we can write a solver to find the nuclear norm $\\|X\\|_{\\ast}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxopt\n",
    "import cvxopt.solvers\n",
    "# Suppress cvxopt logging.\n",
    "cvxopt.solvers.options['show_progress'] = False\n",
    "\n",
    "\n",
    "def get_solution(X, get_coefficients_method=get_coefficients):\n",
    "    M, N = X.shape\n",
    "    \n",
    "    Gs, optimized_coefficients = get_coefficients_method(M, N)\n",
    "    \n",
    "    X_bar = get_X_bar(X, M, N)\n",
    "    # List of h's, i.e. scalar RHS in the SDP.\n",
    "    hs = [cvxopt.matrix(X_bar)]\n",
    "\n",
    "    solution = cvxopt.solvers.sdp(optimized_coefficients, Gs=Gs, hs=hs)\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we put it to work and compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Primal objective: 14.2920671924\n",
      "        Duration: 0.10688996315\n",
      "----------------------------------------------------------------------\n",
      "     NumPy ||X||: 14.2920672228\n",
      "        Duration: 0.00019383430481\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg\n",
    "import time\n",
    "\n",
    "\n",
    "def compare(X=None, M=2, N=2, get_coefficients_method=get_coefficients):\n",
    "    if X is None:\n",
    "        X = get_random_X(M, N)\n",
    "\n",
    "    start = time.time()\n",
    "    solution = get_solution(X, get_coefficients_method=get_coefficients_method)\n",
    "    sdp_optimum = solution['primal objective']\n",
    "    duration = time.time() - start\n",
    "    print 'Primal objective:', sdp_optimum\n",
    "    print '        Duration:', duration\n",
    "\n",
    "    print '-' * 70\n",
    "\n",
    "    # Use SciPy to compute nuclear norm.\n",
    "    start = time.time()\n",
    "    nuclear_norm = np.sum(scipy.linalg.svdvals(X))\n",
    "    duration = time.time() - start\n",
    "    print '     NumPy ||X||:', nuclear_norm\n",
    "    print '        Duration:', duration\n",
    "    return sdp_optimum, nuclear_norm, X\n",
    "\n",
    "\n",
    "M, N = 10, 15\n",
    "_, _, X = compare(M=M, N=N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BIIIIG Slowdown\n",
    "\n",
    "This is quite slow, so we attempt to speed it up using sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_coefficients_sparse(M, N):\n",
    "    G_sparse_I = []  # List of rows in single G.\n",
    "    # I: Rows is equal to (M + N)^2.\n",
    "    G_sparse_J = []  # List of columns in single G.\n",
    "    # J: Columns is equal to number of coefficients.\n",
    "\n",
    "    # This is always 1.0.\n",
    "    coefficients_sparse_I = []  # List of rows for single c.\n",
    "    # There is only 1 column, so this will just be\n",
    "    #     [0] * len(coefficients_sparse_I).\n",
    "\n",
    "    curr_coeff_index = 0\n",
    "\n",
    "    # W1, top left.\n",
    "    total_rows = M + N\n",
    "    for i in xrange(M):\n",
    "        # Handle case i = j.\n",
    "        coefficients_sparse_I.append(curr_coeff_index)\n",
    "\n",
    "        G_sparse_I.append(total_rows * i + i)\n",
    "        G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "        # Update the coefficient index.\n",
    "        curr_coeff_index += 1\n",
    "\n",
    "        for j in xrange(i + 1, M):\n",
    "            G_sparse_I.append(total_rows * i + j)\n",
    "            G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "            G_sparse_I.append(total_rows * j + i)\n",
    "            G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "            # Update the coefficient index.\n",
    "            curr_coeff_index += 1\n",
    "\n",
    "    # W2, bottom right.\n",
    "    for i in xrange(M, M + N):\n",
    "        # Handle case i = j.\n",
    "        coefficients_sparse_I.append(curr_coeff_index)\n",
    "\n",
    "        G_sparse_I.append(total_rows * i + i)\n",
    "        G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "        # Update the coefficient index.\n",
    "        curr_coeff_index += 1\n",
    "\n",
    "        for j in xrange(i + 1, M + N):\n",
    "            G_sparse_I.append(total_rows * i + j)\n",
    "            G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "            G_sparse_I.append(total_rows * j + i)\n",
    "            G_sparse_J.append(curr_coeff_index)\n",
    "\n",
    "            # Update the coefficient index.\n",
    "            curr_coeff_index += 1\n",
    "\n",
    "    # Use the I, J to create our sparse matrices.\n",
    "    G_sparse = cvxopt.spmatrix(-1.0, G_sparse_I, G_sparse_J)\n",
    "    G = [G_sparse]\n",
    "\n",
    "    coefficients_sparse_J = [0] * len(coefficients_sparse_I)\n",
    "    coefficients_sparse = cvxopt.spmatrix(\n",
    "        0.5, coefficients_sparse_I, coefficients_sparse_J)\n",
    "    \n",
    "    coefficients_dense = cvxopt.matrix(coefficients_sparse)\n",
    "\n",
    "    return G, coefficients_dense"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Primal objective: 14.2920671924\n",
      "        Duration: 0.123693943024\n",
      "----------------------------------------------------------------------\n",
      "     NumPy ||X||: 14.2920672228\n",
      "        Duration: 0.000292062759399\n"
     ]
    }
   ],
   "source": [
    "# Use previous X.\n",
    "_, _, _ = compare(X, get_coefficients_method=get_coefficients_sparse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aw Man\n",
    "\n",
    "Sadly, we don't get any sort of speed up."
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}