{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing a sparse solution of a set of linear inequalities\n",
    "\n",
    "A derivative work by Judson Wilson, 5/11/2014.<br>\n",
    "Adapted from the CVX example of the same name, by Almir Mutapcic, 2/28/2006.\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* Section 6.2, Boyd & Vandenberghe \"Convex Optimization\" <br>\n",
    "* \"Just relax: Convex programming methods for subset selection and sparse approximation\" by J. A. Tropp\n",
    "\n",
    "\n",
    "## Introduction\n",
    "\n",
    "\n",
    "We consider a set of linear inequalities \n",
    "$Ax \\preceq b$ \n",
    "which are feasible. We apply two heuristics to find a sparse point $x$ that satisfies these inequalities.\n",
    "\n",
    "The (standard) $\\ell_1$-norm heuristic for finding a sparse solution is:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\|x\\|_1 \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b.\n",
    "    \\end{array}\n",
    "\n",
    "The log-based heuristic is an iterative method for finding\n",
    "a sparse solution, by finding a local optimal point for the problem:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\sum_i \\log \\left( \\delta + \\left|x_i\\right| \\right) \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b,\n",
    "    \\end{array}\n",
    "where $\\delta$ is a small threshold value (which determines if a value is close to zero).\n",
    "We cannot solve this problem since it is a minimization of a concave\n",
    "function and thus it is not a convex problem. However, we can apply\n",
    "a heuristic in which we linearize the objective, solve, and re-iterate.\n",
    "This becomes a weighted $\\ell_1$-norm heuristic:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   &  \\sum_i W_i \\left|x_i\\right| \\\\\n",
    "    \\mbox{subject to} & Ax \\preceq b,\n",
    "    \\end{array}\n",
    "which in each iteration re-adjusts the weights $W_i$ based on the rule:\n",
    "    $$W_i = 1/(\\delta + \\left|x_i\\right|),$$\n",
    "where $\\delta$ is a small threshold value.\n",
    "\n",
    "This algorithm is described in papers:\n",
    "\n",
    "* \"An affine scaling methodology for best basis selection\"<br>\n",
    "  by B. D. Rao and K. Kreutz-Delgado\n",
    "* \"Portfolio optimization with linear and fixed transaction costs\"<br>\n",
    "  by M. S. Lobo, M. Fazel, and S. Boyd\n",
    "\n",
    "## Generate problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "# Fix random number generator so we can repeat the experiment.\n",
    "np.random.seed(1)\n",
    "\n",
    "# The threshold value below which we consider an element to be zero.\n",
    "delta = 1e-8\n",
    "\n",
    "# Problem dimensions (m inequalities in n-dimensional space).\n",
    "m = 100\n",
    "n = 50\n",
    "\n",
    "# Construct a feasible set of inequalities.\n",
    "# (This system is feasible for the x0 point.)\n",
    "A  = np.random.randn(m, n)\n",
    "x0 = np.random.randn(n)\n",
    "b  = A.dot(x0) + np.random.random(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\ell_1$-norm heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "status: optimal\n",
      "Found a feasible x in R^50 that has 40 nonzeros.\n",
      "optimal objective value: 28.582394099513873\n"
     ]
    }
   ],
   "source": [
    "# Create variable.\n",
    "x_l1 = cp.Variable(shape=n)\n",
    "\n",
    "# Create constraint.\n",
    "constraints = [A*x_l1 <= b]\n",
    "\n",
    "# Form objective.\n",
    "obj = cp.Minimize(cp.norm(x_l1, 1))\n",
    "\n",
    "# Form and solve problem.\n",
    "prob = cp.Problem(obj, constraints)\n",
    "prob.solve()\n",
    "print(\"status: {}\".format(prob.status))\n",
    "\n",
    "# Number of nonzero elements in the solution (its cardinality or diversity).\n",
    "nnz_l1 = (np.absolute(x_l1.value) > delta).sum()\n",
    "print('Found a feasible x in R^{} that has {} nonzeros.'.format(n, nnz_l1))\n",
    "print(\"optimal objective value: {}\".format(obj.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative log heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1: Found a feasible x in R^50 with 48 nonzeros...\n",
      "Iteration 2: Found a feasible x in R^50 with 36 nonzeros...\n",
      "Iteration 3: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 4: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 5: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 6: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 7: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 8: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 9: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 10: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 11: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 12: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 13: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 14: Found a feasible x in R^50 with 33 nonzeros...\n",
      "Iteration 15: Found a feasible x in R^50 with 33 nonzeros...\n"
     ]
    }
   ],
   "source": [
    "# Do 15 iterations, allocate variable to hold number of non-zeros\n",
    "# (cardinality of x) for each run.\n",
    "NUM_RUNS = 15\n",
    "nnzs_log = np.array(())\n",
    "\n",
    "# Store W as a positive parameter for simple modification of the problem.\n",
    "W = cp.Parameter(shape=n, nonneg=True); \n",
    "x_log = cp.Variable(shape=n)\n",
    "\n",
    "# Initial weights.\n",
    "W.value = np.ones(n);\n",
    "\n",
    "# Setup the problem.\n",
    "obj = cp.Minimize( W.T*cp.abs(x_log) ) # sum of elementwise product\n",
    "constraints = [A*x_log <= b]\n",
    "prob = cp.Problem(obj, constraints)\n",
    "\n",
    "# Do the iterations of the problem, solving and updating W.\n",
    "for k in range(1, NUM_RUNS+1):\n",
    "    # Solve problem.\n",
    "    # The ECOS solver has known numerical issues with this problem\n",
    "    # so force a different solver.\n",
    "    prob.solve(solver=cp.CVXOPT)\n",
    "    \n",
    "    # Check for error.\n",
    "    if prob.status != cp.OPTIMAL:\n",
    "        raise Exception(\"Solver did not converge!\")\n",
    "\n",
    "    # Display new number of nonzeros in the solution vector.\n",
    "    nnz = (np.absolute(x_log.value) > delta).sum()\n",
    "    nnzs_log = np.append(nnzs_log, nnz);\n",
    "    print('Iteration {}: Found a feasible x in R^{}'\n",
    "          ' with {} nonzeros...'.format(k, n, nnz))\n",
    "\n",
    "    # Adjust the weights elementwise and re-iterate\n",
    "    W.value = np.ones(n)/(delta*np.ones(n) + np.absolute(x_log.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Result plots\n",
    "\n",
    "The following code plots the result of the $\\ell_1$-norm heuristic, as well as the result for each iteration of the log heuristic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13704 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 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",
    "plt.figure(figsize=(6,6))\n",
    "\n",
    "# Plot the two data series.\n",
    "plt.plot(range(1,1+NUM_RUNS), nnzs_log, label='log heuristic')\n",
    "plt.plot((1, NUM_RUNS), (nnz_l1, nnz_l1), linestyle='--', label='l1-norm heuristic')\n",
    "\n",
    "# Format and show plot.\n",
    "plt.xlabel('iteration', fontsize=16)\n",
    "plt.ylabel('number of non-zeros (cardinality)', fontsize=16)\n",
    "plt.ylim(0,n)\n",
    "plt.xlim(1,NUM_RUNS)\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "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
}
