2_Reconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noiseless Sparse Reconstruction\n",
    "In this notebook we will review a few \"easy-to-implent\" Compressed Sensing (CS) reconstruction methods. The goal is to familiarize the reader with a few different techniques to gain an intuition of how these reconstruction approaches operate and help alleviate the \"black-magic\" view of CS reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn.linear_model as lm\n",
    "import scipy.linalg as la\n",
    "import pywt as pw\n",
    "from DepsRecon import *\n",
    "from scipy.stats import lomax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can start to review any CS signal reconstruction, we first require a signal to reconstruct! We will generate a target $K$-Sparse signal synthetically. We'll then take some measurements of this signal and from those measurements attempt to recover the original. Since we have an original target to compare to, this is often called a \"planted\" problem.\n",
    "\n",
    "We'll start by creating a function to generate a random orthonormal basis, $F$. Orthonormal, here, means that the basis vectors of $F$ are independent from one another and, additionally, that each of these vectors is normalized such that they have an $\\ell_2$ norm of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RandomOrthoProjection(M,N):\n",
    "    \"\"\" Create a random orthonormal projection matrix.\n",
    "    \n",
    "        Inputs\n",
    "        ------\n",
    "        M: int\n",
    "           Number of rows in the matrix.\n",
    "        N: int\n",
    "           Number of columns in the matrix.\n",
    "           \n",
    "        Returns\n",
    "        -------\n",
    "        F: array_like, float\n",
    "           The random orthonormal matrix\n",
    "    \"\"\"\n",
    "    \n",
    "    F = np.random.randn(M,N)   # Generate the random samples\n",
    "    F = la.orth(F.T).T         # Perform orthogonalization (Basis for Span(F))\n",
    "    \n",
    "    return F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll create particular realizations of our signal and projection matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 1024      # Signal Dimensionality\n",
    "M = 128       # Number of measurements\n",
    "K = 20        # Number of non-zeros\n",
    "x,xs = RandomKSparseSignal(K,N)\n",
    "F = RandomOrthoProjection(M,N)\n",
    "\n",
    "ShowProjection(F)\n",
    "ShowSignal(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, given $F$ and the true $x$, lets calculate the value of the observations for this planted problem. We accomplish this through the simple matrix multiplication,\n",
    "$$ y = F x$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = np.dot(F,x)\n",
    "ShowMeasurements(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal then is to find some reconstruction $\\hat{x}$, which should be close to $x$, using only the observations $y$ and the projection matrix $F$. The first thing that one might attempt to do is to find the Least-Squares solution to this problem. As we saw in lecture, we cannot say that this solution will be anything at all like our original one, since there is an entire space of possible $\\hat{x}$ which all map to the same $y$ in this case of under-sampling, $M<N$.\n",
    "\n",
    "We can calculate the least-squares solution by using the Pseudo-Inverse of the matrix $F$,\n",
    "\n",
    "$$\\hat{x} = F^{\\dagger} y = (F^T F)^{-1} F^T y$$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "FtF = np.dot(F.T,F)           # Calculate $(F^T F)$\n",
    "Fty = np.dot(F.T,y)           # Calculate $F^T y$\n",
    "xlsq = la.solve(FtF,Fty)      # Invert FtF and multiply against Fty\n",
    "\n",
    "ShowSignal(xlsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that this signal is nowhere near to the original one. However, as we mentioned in lecture, this signal does indeed match our observations, since it is one of the many vectors which lie in the space of possible solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ShowMeasurements(np.dot(F,xlsq))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, as mentioned earlier, one the of the easiest approaches to regularization is *Ridge Regression* or *Tikhonov Regularization*. Here, we will regularize the least-squares problem by solving instead\n",
    "$$\\arg \\min_{x}\\quad||y - Fx||_2^2 + \\lambda ||x||_2.$$\n",
    "\n",
    "This problem can be solved directly via\n",
    "$$\\hat{x} = (F^T F + \\lambda I)^{-1} F^T y.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lam = 0.2\n",
    "FtFI = FtF + lam * np.eye(np.shape(F)[1])\n",
    "\n",
    "xtik = la.solve(FtFI,Fty)\n",
    "ShowSignal(xtik)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we see that the $\\ell_2$ regularized solution is more correlated with our original signal, however, it still suffers from a lot of noise on what should be the non-zero values. This is because the $\\ell_2$ norm does not promote trunction to 0 for these small values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative Hard Thresholding (IHT)\n",
    "---\n",
    "One of the easiest-to-implement algorithms for CS reconstruction is iterative hard thresholding. The IHT approach seeks to solve the following problem,\n",
    "$$\\arg\\min_x\\quad ||y-Fx||_2^2\\quad s.t.\\quad||x||_0 = K.$$\n",
    "\n",
    "The IHT tackles this non-convex problem by iteratively bouncing between the space of satisfying solutions and the space of $K$-sparse signals. This sounds a bit complicated at first, but as we will see, this is really quite a simple procedure. In essence, the entire algorithm can be written in the following equation,\n",
    "\n",
    "$$ x^{(t+1)} =  H_K\\left(x^{(t)} + F^T(y - Fx^{(t)}) \\right),$$\n",
    "\n",
    "where $x^{(t)}$ refers to the approximation of the original signal at iteration $t$ and the operator $H_K(\\cdot)$ refers to a thresholding of all but the $K$ largest magnitude coefficients.\n",
    "\n",
    "---\n",
    "- T. Blumensath & M. E. Davies, \"Iterative Thresholding for Sparse Approximations,\" _Journal of Fourier Analysis and Applications_, vol. 14, no. 5, pg. 629-654, 2008.\n",
    "- T. Blumensath & M. E. Davies, \"Iterative Hard Thresholding for Comrpessed Sensing,\" _Applied and Computational Harmonic Analysis_, vol. 27, no. 3, pg. 265-274, 2009.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def IterativeHardThresholding(F,y,K,iters,xtrue):\n",
    "    \"\"\" Use the IHT alorithm to reconstruct a sparse signal.\n",
    "    \n",
    "        Attempt to solve the inverse problem\n",
    "        \n",
    "                    y = F x\n",
    "        \n",
    "        for known y, F, and unknown K-sparse x.\n",
    "                \n",
    "    \n",
    "        Inputs\n",
    "        -----\n",
    "        F: array_like, float\n",
    "           An (M,N) projection matrix\n",
    "        y: array_like, float\n",
    "           A vector of observations\n",
    "        K: int\n",
    "           Number of non-zeros in the recovered signal.\n",
    "        iters: int\n",
    "               Number of iterations to run for.\n",
    "               \n",
    "        Returns\n",
    "        -------\n",
    "        x: array_like, float\n",
    "           Recovered K-sparse signal\n",
    "        mseval: array_like, float\n",
    "                The MSE against ``xtrue`` at each iteration\n",
    "        resval: array_like, float\n",
    "                The MSE between ``y`` and the estimate ``F x`` at each\n",
    "                iteration.\n",
    "    \"\"\"\n",
    "    m,n = np.shape(F)                        # Find the dimensions of the problem\n",
    "    x = np.zeros((n,1))                      # Initial solution for ``x``, the starting point\n",
    "    mseval = np.zeros((iters,1))             # Initialize array for MSE values\n",
    "    resval = np.zeros((iters,1))             # Initialize array for Residual values\n",
    "    \n",
    "    for t in range(0,iters):                 # Loop over the number of iterations\n",
    "        ## Update Approximation\n",
    "        r = y - np.dot(F,x)                  # Find the gradient in Real^M, (y - F x)\n",
    "        z = x + np.dot(F.T,r)                # Take a step in Real^N by projecting gradient\n",
    "        x = KSparseApproximation(z,K)[0]     # Threshold the new value of x  \n",
    "        ## Measure Progress\n",
    "        mseval[t]   = MeanSquareError(xtrue,x) # Calculate how far from the true solution we are\n",
    "        resval[t] = MeanSquareError(y,np.dot(F,x)) # Calculate how far from the measurements we are\n",
    "        \n",
    "    return x,mseval,resval        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xiht,ihtmse,ihtres = IterativeHardThresholding(F,y,K,500,x)\n",
    "\n",
    "ShowRecovery(x,xiht)\n",
    "ShowHistory(ihtmse,ihtres)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthogonal Matching Pursuit (OMP)\n",
    "---\n",
    "For orthogonal matching pursuit, one is also attempting to solve the reconstruction problem of the form\n",
    "$$\\arg\\min_x\\quad ||y-Fx||_2^2\\quad s.t.\\quad||x||_0 = K.$$\n",
    "\n",
    "However, instead of an iteration based on bouncing between spaces, lets instead take a more support-oriented approach. OMP is an iterative algorithm, meaning that the approximation successively improves over multiple passes, however, the manner in which we do this is different than the IHT. Instead of attempting to locate **all** $K$ significant coefficients at **each** iteration, we will instead have each iteration _refine_ our estimation of the support.\n",
    "\n",
    "E.g. on each iteration (or pass) we will simply attempt to identify **one** element of the support, including this in a running list of support values. After we have completed $K$ passes, we will have identified up to $K$ support locations. The refinement is accomplished by looking at the difference $r^{(t)} = y - F x^{(t)}$ at each pass and seeing which column of $F$ is most correlated with this difference. The *most* correlated column is therefore the support-location which will move us the furthest towards the solution space and should therefore be included in the support list, $S^{(t)}$.\n",
    "\n",
    "Subsequently, we refine $r^{(t)}$ by finding the LSQ solution of the on-support (or sub-) problem, $F_S x_S = y$, which is an overdetermined. Finally, we update $r^{(t+1)} = y - F_S x_S$ and begin again.\n",
    "\n",
    "---\n",
    "- J. Tropp and A. Gilbert, \"Signal Recovery from Random Measurements via Orthogonal Matching Pursuit,\" _IEEE Transactions on Information Theory_, vol. 53, no. 12, pg. 4655-4666, 2007."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def OrthogonalMatchingPursuit(F,y,K,xtrue):\n",
    "    \"\"\" Use the OMP alorithm to reconstruct a sparse signal.\n",
    "    \n",
    "        Attempt to solve the inverse problem\n",
    "        \n",
    "                    y = F x\n",
    "        \n",
    "        for known y, F, and unknown K-sparse x.\n",
    "                \n",
    "    \n",
    "        Inputs\n",
    "        -----\n",
    "        F: array_like, float\n",
    "           An (M,N) projection matrix\n",
    "        y: array_like, float\n",
    "           A vector of observations\n",
    "        K: int\n",
    "           Number of non-zeros in the recovered signal.\n",
    "        xtrue: array_like, float\n",
    "               Original signal, used to track MSE progress\n",
    "               \n",
    "        Returns\n",
    "        -------\n",
    "        x: array_like, float\n",
    "           Recovered K-sparse signal\n",
    "        Support: array_like, int\n",
    "                List of support locations\n",
    "        mseval: array_like, float\n",
    "                The MSE against ``xtrue`` at each iteration\n",
    "        resval: array_like, float\n",
    "                The MSE between ``y`` and the estimate ``F x`` at each\n",
    "                iteration.        \n",
    "    \"\"\"\n",
    "    m,n = np.shape(F)                    # Find problem dimensions\n",
    "    Support = np.empty((0,1),dtype=int)  # Initialize list of support\n",
    "    r = y                                # Intialize ``r`` to ``y``\n",
    "    Fs = []                              # Initial sub-matrix is empty\n",
    "    mseval = np.zeros((K,1))             # Initialize array for MSE values\n",
    "    resval = np.zeros((K,1))             # Initialize array for Residual values\n",
    "    \n",
    "    ## Main Loop\n",
    "    for t in range(1,K+1):\n",
    "        ## Find most correlated column\n",
    "        response = np.abs(np.dot(F.T,r))              # Calculate the correlations between\n",
    "                                                      #   the columns of F and r.\n",
    "        maxcorrcol = response.argmax()                # Index of the maximum correlation\n",
    "        Support = np.append(Support,maxcorrcol)       # Add this index to the support list\n",
    "        Support = np.unique(Support)                  # Remove duplicate entries from the support list\n",
    "        NS = np.size(Support)                         # Record the number of on-support coefficients\n",
    "        ## Construct on-support sub-matrix \n",
    "        Fs = F[:,Support.astype(int)]                 # Sub-sample the columns of F according to the support list\n",
    "        Fs = np.reshape(Fs,(m,NS))                    # Ensure that A is in the right data-type (safety)\n",
    "        ## Solve LSQ for on-support sub-problem\n",
    "        if NS>1:\n",
    "            xs = la.lstsq(Fs,y)[0]                    # Solve A x = y using Least Squares\n",
    "        else:\n",
    "            # In the case that ``A`` consists of a single column,\n",
    "            # ``scipy.linalg`` will not allow us to use the ``lstsq``\n",
    "            # function. So, we perform a linear regression to find\n",
    "            # the best-fitting coefficient.\n",
    "            xs = lm.LinearRegression().fit(Fs,y).coef_[0]\n",
    "        ## Update residual\n",
    "        r = np.reshape(r,(m,1))                         # Ensure that r is in the data-type (safety)\n",
    "        r  = y - np.reshape(np.dot(Fs,xs),(m,1))        # Update r, ensure that A\n",
    "        ## Measure Progress\n",
    "        x = np.zeros((n,1))                             # Make a full version of x\n",
    "        x[Support] = xs                                 # Assign on-support coefficients\n",
    "        mseval[t-1]   = MeanSquareError(xtrue,x)        # Calculate how far from the true solution we are\n",
    "        resval[t-1] = MeanSquareError(y,np.dot(F,x))    # Calculate how far from the measurements we are\n",
    "        \n",
    "    return x, Support, mseval, resval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myxomp,Support,ompmse,ompres = OrthogonalMatchingPursuit(F,y,K,x)\n",
    "\n",
    "ShowRecovery(x,myxomp)\n",
    "ShowHistory(ompmse,ompres)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compressive Sampling Matching Pursuit (CoSaMP)\n",
    "---\n",
    "\n",
    "The CoSaMP approach is very similar to the original OMP approach we have just discussed. It too seekes to solve the same $\\ell_0$ regularized optimization problem as the IHT and OMP algorithms. \n",
    "\n",
    "As this approach is at its core a _Matching Pursuit_, the CoSaMP algorithm also retains a list of supports which are updated at each iteration. This support list is refined at each iteration, again, by looking at the difference $r^{(t)} = y - F x^{(t)}$. \n",
    "\n",
    "However, rather than add a single support at each iteartion, CoSaMP mixes the ideas $K$-Thresholding with a Matching Pursuit approach, retaining $K$ non-zeros at each iteration.\n",
    "\n",
    "At the beginning of every CoSaMP iteration, $2K$ possible support values are allowed. These possible locations are merged with the list at this time step, and this merged support list is used to form the on-support subproblem, much like OMP.\n",
    "\n",
    "Finally, after solving this reduced LSQ problem, the result is then thresholded to the top $K$ coefficeints, these $K$ coefficients are taken as the final support list for the iteration, and the process is repeated.\n",
    "\n",
    "\n",
    "---\n",
    "- D. Needell and J. Tropp, \"CoSaMP: Iterative Signal Recovery from Incomplete and Inaccurate Samples,\" _Applied and Computational Harmonic Analysis_, vol. 26, no. 3, pg. 301-321, 2009."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def CompSampMatchingPursuit(F,y,K, iters,xtrue):\n",
    "    \"\"\" Use the CoSaMP alorithm to reconstruct a sparse signal.\n",
    "    \n",
    "        Attempt to solve the inverse problem\n",
    "        \n",
    "                    y = F x\n",
    "        \n",
    "        for known y, F, and unknown K-sparse x.\n",
    "                \n",
    "    \n",
    "        Inputs\n",
    "        -----\n",
    "        F: array_like, float\n",
    "           An (M,N) projection matrix\n",
    "        y: array_like, float\n",
    "           A vector of observations\n",
    "        iters: int\n",
    "               Number of iterations to perform\n",
    "        xtrue: array_like, float\n",
    "               Original signal, used to track MSE progress\n",
    "               \n",
    "        Returns\n",
    "        -------\n",
    "        x: array_like, float\n",
    "           Recovered K-sparse signal\n",
    "        Support: array_like, int\n",
    "                List of support locations\n",
    "        mseval: array_like, float\n",
    "                The MSE against ``xtrue`` at each iteration\n",
    "        resval: array_like, float\n",
    "                The MSE between ``y`` and the estimate ``F x`` at each\n",
    "                iteration.        \n",
    "    \"\"\"\n",
    "    m,n = np.shape(F)                        # Obtain problem dimensions\n",
    "    r = y                                    # Initialize ``r`` to ``y``\n",
    "    S = np.empty((0,1),dtype=int)            # Start with an empty support list\n",
    "    mseval = np.zeros((iters,1))             # Initialize array for MSE values\n",
    "    resval = np.zeros((iters,1))             # Initialize array for Residual values\n",
    "    \n",
    "    for t in range(0,iters):                 # Loop over the number of iterations\n",
    "        ## Form Signal Proxy\n",
    "        q = np.dot(F.T,r)                    # Find q = F^T r, the gradient in Real^N\n",
    "        ## Identify large components\n",
    "        O = KSparseApproximation(q,2*K)[1]   #  Identify the locations of the 2K largest\n",
    "                                             #    magnitudes in the gradient.\n",
    "        ## Merge Supports (Union)\n",
    "        T = np.append(O,S)                   # Merge these 2K locations into the support list\n",
    "        T = np.unique(T)                     # Remove duplicates\n",
    "        ## Solve reduced system\n",
    "        b = np.zeros((n,1))                  # Initialize LSQ solution to zeros\n",
    "        b[T] = la.lstsq(F[:,T],y)[0]         # Solve over-determined on-support LSQ problem\n",
    "        ## Prune approximation\n",
    "        x,S = KSparseApproximation(b,K)      # Prune out all but the top $K$ magnitude coefficients\n",
    "                                             #     from the previous LSQ solution.\n",
    "        ## Update residual \n",
    "        r = y - np.dot(F,x)                  # Update r = y - F x        \n",
    "        ## Measure Progress\n",
    "        mseval[t]   = MeanSquareError(xtrue,x)        # Calculate how far from the true solution we are\n",
    "        resval[t] = MeanSquareError(y,np.dot(F,x))    # Calculate how far from the measurements we are\n",
    "    \n",
    "    return x,S,mseval,resval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xcosamp,Support,cosampmse,cosampres = CompSampMatchingPursuit(F,y,20,100,x)\n",
    "ShowRecovery(x,xcosamp)\n",
    "ShowHistory(cosampmse,cosampres)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## L1 Minimization (Lasso)\n",
    "---\n",
    "\n",
    "One can also attempt, instead, to solve the $\\ell_1$ regularized minimization problem\n",
    "$$\\arg\\min_x\\quad ||y-F x||_2^2 + \\lambda ||x||_1,$$\n",
    "which is commonly known as the \"Lasso\". One can write a mapping between the Lasso and other $\\ell_1$ penalized problems such as basis-pursuit denoising. So, we'll focus on attempting to solve this unconstrained optimization problem.\n",
    "\n",
    "The details of implementations of solvers for these kinds of problems are a little bit beyond the scope of this notebook, so we will rely on external packages to solve the Lasso problem. For example, `sckit-learn` already includes functions to solve the Lasso by a number of different methods. Below, we show how to use the Least Angle Regression (LARS) method to solve the Lasso, and thus, the $\\ell_1$ minimization problem for estimating sparse $x$ given $F$ and $y$.\n",
    "\n",
    "**Note:** Since we do not have access to the inner-workings of this algorithm, we cannot display the same historical progress of the MSE over the reconstruction iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def L1Minimization(F,y,lam=1):\n",
    "    return lm.LassoLars(alpha=lam).fit(F,y).coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xl1 = L1Minimization(F,y,0.000001)\n",
    "xl1 = np.reshape(xl1,(N,1))\n",
    "\n",
    "ShowRecovery(x,xl1)\n",
    "print \"Lars-Lasso MSE: %0.4e\" % MeanSquareError(x,xl1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Comparison\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def RandomCompressibleSignal(alpha,N):\n",
    "    x=lomax.rvs(c=alpha,size=(N,1))\n",
    "\n",
    "    x = x / la.norm(x)\n",
    "\n",
    "    signs = np.random.rand(N,1) > 0.5\n",
    "    x[signs] = -1.0 * x[signs]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N = 1024     # Signal dimensionality\n",
    "K = 20      # Number of non-zeros. Used both to\n",
    "             #    generate the signal as well as \n",
    "             #    a parameter for the algorithms\n",
    "             #    (oracle parameter setting).\n",
    "\n",
    "## Generate a truly sparse signal\n",
    "# x = RandomKSparseSignal(K,N)[0]\n",
    "## Generate a Random Compressible signal\n",
    "DecayStrength = 1.0\n",
    "x = RandomCompressibleSignal(DecayStrength,N)\n",
    "\n",
    "MToTest = range(2,786)[::32]\n",
    "Nexp = np.size(MToTest)\n",
    "\n",
    "mse_iht = np.zeros((Nexp,1))\n",
    "mse_omp = np.zeros((Nexp,1))\n",
    "mse_cosamp = np.zeros((Nexp,1))\n",
    "mse_l1 = np.zeros((Nexp,1))\n",
    "\n",
    "i = 0\n",
    "for m in MToTest:\n",
    "    print \"Testing M = %d\" % m\n",
    "    F = RandomOrthoProjection(m,N)\n",
    "    y = np.dot(F,x)\n",
    "    \n",
    "    xiht = IterativeHardThresholding(F,y,K,500,x)[0]\n",
    "    xomp = OrthogonalMatchingPursuit(F,y,K,x)[0]\n",
    "    xcosamp = CompSampMatchingPursuit(F,y,K,500,x)[0]\n",
    "    xl1 = L1Minimization(F,y,0.000001)\n",
    "    \n",
    "    \n",
    "    xiht = np.reshape(xiht,(N,1))\n",
    "    xomp = np.reshape(xomp,(N,1))\n",
    "    xcosamp = np.reshape(xcosamp,(N,1))\n",
    "    xl1 = np.reshape(xl1,(N,1))\n",
    "    \n",
    "    mse_iht[i] = MeanSquareError(x,xiht)\n",
    "    mse_omp[i] = MeanSquareError(x,xomp)\n",
    "    mse_cosamp[i] = MeanSquareError(x,xcosamp)\n",
    "    mse_l1[i] = MeanSquareError(x,xl1)   \n",
    "    \n",
    "    i = i + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xk = KSparseApproximation(x,K)[0]\n",
    "xkmse = MeanSquareError(x,xk)\n",
    "\n",
    "plt.plot(MToTest,xkmse*np.ones(np.shape(MToTest)),'-k')\n",
    "plt.plot(MToTest,mse_iht,'-*b',linewidth=1,label=\"IHT\")\n",
    "plt.plot(MToTest,mse_omp,'-xr',linewidth=1,label=\"OMP\")\n",
    "plt.plot(MToTest,mse_cosamp,'-ok',linewidth=1,label=\"CoSaMP\")\n",
    "plt.plot(MToTest,mse_l1,'-^g',linewidth=1,label=\"Lasso/LARS\")\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"M\",fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.axis('tight')\n",
    "plt.ylabel(\"$\\\\frac{1}{N}||x - \\hat{x}||_2^2$\",fontsize=20)\n",
    "plt.legend(loc=3)\n",
    "plt.ylim((1e-10,1e2))\n",
    "plt.title(\"Noiselss Reconstruction over $M$\",fontsize=20);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}