{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Douglas Rachford Proximal Splitting\n", "===================================\n", "\n*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_matlab/) for details about how to install the toolboxes.\n", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n", "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$\n", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n", "$\\newcommand{\\CC}{\\mathbb{C}}$\n", "$\\newcommand{\\RR}{\\mathbb{R}}$\n", "$\\newcommand{\\EE}{\\mathbb{E}}$\n", "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n", "$\\newcommand{\\Ww}{\\mathcal{W}}$\n", "$\\newcommand{\\Vv}{\\mathcal{V}}$\n", "$\\newcommand{\\Nn}{\\mathcal{N}}$\n", "$\\newcommand{\\NN}{\\mathcal{N}}$\n", "$\\newcommand{\\Hh}{\\mathcal{H}}$\n", "$\\newcommand{\\Bb}{\\mathcal{B}}$\n", "$\\newcommand{\\Ee}{\\mathcal{E}}$\n", "$\\newcommand{\\Cc}{\\mathcal{C}}$\n", "$\\newcommand{\\Gg}{\\mathcal{G}}$\n", "$\\newcommand{\\Ss}{\\mathcal{S}}$\n", "$\\newcommand{\\Pp}{\\mathcal{P}}$\n", "$\\newcommand{\\Ff}{\\mathcal{F}}$\n", "$\\newcommand{\\Xx}{\\mathcal{X}}$\n", "$\\newcommand{\\Mm}{\\mathcal{M}}$\n", "$\\newcommand{\\Ii}{\\mathcal{I}}$\n", "$\\newcommand{\\Dd}{\\mathcal{D}}$\n", "$\\newcommand{\\Ll}{\\mathcal{L}}$\n", "$\\newcommand{\\Tt}{\\mathcal{T}}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\al}{\\alpha}$\n", "$\\newcommand{\\la}{\\lambda}$\n", "$\\newcommand{\\ga}{\\gamma}$\n", "$\\newcommand{\\Ga}{\\Gamma}$\n", "$\\newcommand{\\La}{\\Lambda}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\Si}{\\Sigma}$\n", "$\\newcommand{\\be}{\\beta}$\n", "$\\newcommand{\\de}{\\delta}$\n", "$\\newcommand{\\De}{\\Delta}$\n", "$\\newcommand{\\phi}{\\varphi}$\n", "$\\newcommand{\\th}{\\theta}$\n", "$\\newcommand{\\om}{\\omega}$\n", "$\\newcommand{\\Om}{\\Omega}$\n"], "metadata": {}, "cell_type": "markdown"}, {"source": ["This numerical tour presents the Douglas-Rachford (DR) algorithm to\n", "minimize the sum of two simple functions. It shows an\n", "application to\n", "reconstruction of exactly sparse signal from noiseless measurement using\n", "$\\ell^1$ minimization."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 2, "cell_type": "code", "language": "python", "metadata": {}, "input": ["addpath('toolbox_signal')\n", "addpath('toolbox_general')\n", "addpath('solutions/optim_4b_dr')"]}, {"source": ["Douglas-Rachford Algorithm\n", "--------------------------\n", "The Douglas-Rachford (DR) algorithm is an iterative scheme to minimize\n", "functionals of the form\n", "$$ \\umin{x} F(x) + G(x) $$\n", "where $F$ and $G$ are convex functions for which one is able to\n", "comptue the proximal mappings $ \\text{prox}_{\\gamma F} $ and\n", "$ \\text{prox}_{\\gamma G} $ which are defined as\n", "$$ \\text{prox}_{\\gamma F}(x) = \\text{argmin}_{y} \\frac{1}{2}\\norm{x-y}^2 + \\ga F(y) $$\n", "(the same definition applies also for $G$).\n", "\n", "\n", "The important point is that $F$ and $G$ do not need to be smooth.\n", "One only needs them to be \"proximable\".\n", "\n", "\n", "This algorithm was introduced in:\n", "\n", "\n", "P. L. Lions and B. Mercier\n", "_Splitting Algorithms for the Sum of Two Nonlinear Operators_\n", "SIAM Journal on Numerical Analysis\n", "Vol. 16, No. 6 (Dec., 1979), pp. 964-979\n", "\n", "\n", "as a generalization of an algorithm introduced by Douglas and Rachford in\n", "the case of quadratic minimization (which corresponds to the solution of\n", "a linear system).\n", "\n", "\n", "To learn more about this algorithm, you can read:\n", "\n", "\n", "Patrick L. Combettes and Jean-Christophe Pesquet,\n", "_Proximal Splitting Methods in Signal Processing_,\n", "in: Fixed-Point Algorithms for Inverse\n", "Problems in Science and Engineering, New York: Springer-Verlag, 2010.\n", "\n", "\n", "\n", "A DR iteration reads\n", "$$ \\tilde x_{k+1} = \\pa{1-\\frac{\\mu}{2}} \\tilde x_k +\n", "  \\frac{\\mu}{2} \\text{rprox}_{\\gamma G}( \\text{rprox}_{\\gamma F}(\\tilde x_k)  )\n", "  \\qandq x_{k+1} = \\text{prox}_{\\gamma F}(\\tilde x_{k+1},) $$\n", "\n", "\n", "\n", "We have used the following shortcuts:\n", "$$   \\text{rprox}_{\\gamma F}(x) = 2\\text{prox}_{\\gamma F}(x)-x $$\n", "\n", "\n", "It is of course possible to inter-change the roles of $F$ and $G$,\n", "which defines another set of iterations.\n", "\n", "\n", "One can show that for any value of $\\gamma>0$, any $ 0 < \\mu < 2 $,\n", "and any $\\tilde x_0$, $x_k \\rightarrow x^\\star$\n", "which is a minimizer of the minimization of $F+G$.\n", "\n", "Compressed Sensing Acquisition\n", "------------------------------\n", "Compressed sensing acquisition corresponds to a random projection\n", "$y=Ax_0$ of a signal $x_0$ on a\n", "few linear vectors (the lines of $A$). For the recovery of $x_0$ to be possible, it is assumed\n", "to be sparsely represented in an orthogonal basis. Up to a change of\n", "basis, we suppose that the vector $x$ itself is sparse.\n", "\n", "\n", "Initialize the random number generator to have always the same signals."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["set_rand_seeds(123456,123456);"]}, {"source": ["Dimension of the problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["n = 400;"]}, {"source": ["Number of measurements."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["p = round(n/4);"]}, {"source": ["Create a random Gaussian measurement matrix $A$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["A = randn(p,n) / sqrt(p);"]}, {"source": ["Sparsity of the signal."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["s = 17;"]}, {"source": ["We begin by generating a $s$-sparse signal $x_0$ with $s$ randomized values.\n", "Since the measurement matrix is random, one does not care about the sign\n", "of the Diracs, so we use +1 values."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sel = randperm(n);\n", "x0 = zeros(n,1); x0(sel(1:s))=1;"]}, {"source": ["We perform random measurements $y=Ax_0$ without noise."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["y = A*x0;"]}, {"source": ["Compressed Sensing Recovery using DR\n", "------------------------------------\n", "Compressed sensing recovery corresponds\n", "to solving the inverse problem $y=A x_0$, which is ill posed because\n", "$x_0$ is\n", "higher dimensional than $y$.\n", "\n", "\n", "The reconstruction can be perform with $\\ell^1$ minimization,\n", "which regularizes the problems by using the sparsity of the solution.\n", "$$ x^\\star \\in \\uargmin{ A x = y } \\norm{x}_1 $$\n", "where the $\\ell^1$ norm is defined as\n", "$$ \\norm{x}_1 = \\sum_i \\abs{x_i}. $$\n", "\n", "\n", "This is the minimization of a non-smooth function under affine\n", "constraints. This can be shown to be equivalent to a linear programming\n", "problem, for wich various algorithms can be used (simplex, interior\n", "points). We propose here to use the DR iterative algorithm.\n", "\n", "\n", "It is possible to recast this problem as the minimization of $F+G$\n", "where $G(x) = \\norm{x}_1$ and $F(x)=\\iota_{H}$ where $H =\n", "\\enscond{x}{Ax=y}$ is an affine space, and $\\iota_H$ is the indicator\n", "function\n", "$$ \\iota_H(x) = \\choice{ 0 \\qifq x \\in H, \\\\ +\\infty \\qifq x \\notin H. } $$\n", "\n", "\n", "The proximal operator of the $\\ell^1$ norm is the soft thresholding:\n", "$$ \\text{prox}_{\\gamma \\norm{\\cdot}_1}(x)_i = \\max\\pa{ 0, 1-\\frac{\\ga}{\\abs{x_i}} } x_i. $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["proxG = @(x,gamma)max(0,1-gamma./max(1e-15,abs(x))).*x;"]}, {"source": ["Display the 1-D curve of the thresholding."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAANl0lEQVR4nO3d3Xbi\nuBKAUemsvP8r+1y4wxBwCD+WVCXtfdGLNkxTY4w/bJzuum1bAYBs/jd6AAB4h4ABkJKAAZCSgAGQ\nkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQkoABkJKAAZCS\ngAGQkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQkoABkJKAAZCSgAGQ0tfoAdqq\ntY4eASClbdtGj/CHyQNWYrwGtVZjxJnBGNFmMEa0GUqST/9OIQKQkoABkJKAAZBSiJOt7QQ5mwyQ\nS4qdpyMwAFISMABSEjAAUhIwAFISMABSEjAAUhIwAFISMABSEjAAUhIwAFISMABSEjAAUhIwAFIS\nMABSEjAAUhIwAH6odfQEzxEwAP5Tawn/L1n+I2AA/JOoXkXAANjlqlcRMABKwnoVAQMgY72KgAEs\nLmm9ioABrCxvvYqAASwrdb2KgAGsKXu9ioABLGiCehUBA1jNHPUqAgawlGnqVQQMYB0z1asIGMAi\nJqtXETCAFcxXryJgANObsl5FwADmNmu9Silfowd4Wf3+t0K3u9fkwV0AC5q4XiVdwGqtlzhd3y7f\n9dqX3NwFsKC561XSBexP1xkDWNb09SqTBezw4OxyXvHmMQCzeq9eN3vL+KYK2CHFApby9rHX/Zcy\nwbkKEWAeK5w5vEh2BLZt2/2lhvsJw8O7ANaxVL1KuoCVozhdlugWsKzV6lWcQgSYwIL1KgIGkN2a\n9SoCBpDasvUqAgaQ18r1KgIGkNTi9SoCBpCRehUBA0hHvXYCBpCJel0IGEAa6nVNwAByUK8bAgaQ\ngHrdEzCA6NTrkIABhKZevxEwgLjU6wEBAwhKvR4TMICI1OtPAgYQjno9Q8AAYlGvJwkYQCDq9TwB\nA4hCvV4iYAAhqNerBAxgPPV6g4ABDKZe7xEwgJHU620CBjCMen1CwADGUK8PCRjAAOr1OQED6E29\nTiFgAF2p11kEDKAf9TqRgAF0ol7nEjCAHtTrdAIG0Jx6tSBgAG2pVyMCBtCQerUjYACtqFdTAgbQ\nhHq1li9g9duDB/ScB+CeenXwNXqA19Rat++N4vr29QO6DwXwg3r1ke8I7IHDpAH0pF7dJDsCe8PN\nMZnCAe2krle6M1jzBGxf9ZdfL6FSLKCP1PUqP/eWKWI2T8AefzcG0FT2emWULGDbtl0+FygWEIR6\nDZEsYOXolODNEjEDelKvUaa6ChGgM/UaSMAA3qReYwkYwDvUazgBA3iZekUgYACvUa8gBAzgBeoV\nh4ABPEu9QhEwgKeoVzQCBvA39QpIwAD+oF4xCRjAI+oVloAB/Eq9IhMwgGPqFZyAARxQr/gEDOCW\neqUgYAA/qFcWAgbwH/VKRMAA/lGvXAQMoBT1SkjAANQrJQEDVqdeSQkYsDT1ykvAgHWpV2oCBixK\nvbITMGBF6jUBAQOWo15zEDBgLeo1DQED1qJe0xAwYBW1jp6AUwkYsARnDucjYMD81GtKAgZMTr1m\nJWDAzNRrYgIGTEu95iZgwJzUa3oCBkxIvVYgYMBs1GsRX6MHeFn9/lnE7W4LfXAXsAj1WkeygNVa\nL3G6vn2xLzm8C5ieei0lWcAeO4xW/fm3xwgbzEq9PlSz/V1bUwVsd3P4pViwAvX63PXeMkXMpgrY\nvsYVC1ajXmua7SpE9YLVqNeykh2Bbdt2f6nhfs5wX+5CRFiKeq0sWcDKUZn2JYoFq1Gvxc12ChFY\nhHohYEA+6kURMCAd9WInYEAm6sWFgAFpqBfXBAzIQb24IWBAAurFPQEDolMvDgkYEJp68RsBA+JS\nLx4QMCAo9eIxAQMiUi/+JGBAOOrFMwQMiEW9eJKAAYGoF88TMCAK9eIlAgaEoF68SsCA8dSLNwgY\nMJh68R4BA0ZSL94mYMAw6sUnBAwYQ734kIABA6gXnxMwoDf14hQCBnSlXpxFwIB+1IsTCRjQiXpx\nLgEDelAvTidgQHPqRQsCBrSlXjQiYEBD6kU7Aga0ol40JWBAE+pFawIGnE+96EDAgJOpF30IGHAm\n9aKbr9EDnKnWut/YvIFgBPWip3kCVmu9dOv6NtDB/unR246e5gkYMIR0Mcr8AbucV9w5MpvDz1eV\nkbylplGzva/mD5hizccXLdDC9d4yRczmDxgzcbYKuJgnYNu2uQpxbg68gGvzBKzo1tTUC7jhB5lJ\nQL2AewJGdOoFHBIwQlMv4DcCRlzqBTwgYASlXsBjAkZE6gX8ScAIR72AZwgYsagX8CQBIxD1Ap4n\nYEShXsBLBIwQ1At4lYAxnnoBbxAwBlMv4D0CxkjqBbxNwBhGvYBPCBhjqBfwIQFjAPUCPidg9KZe\nwCkEjK7UCziLgNGPegEnEjA6US/gXAJGD+oFnE7AaE69gBYEjLbUC2hEwGhIvYB2BIxW1AtoSsBo\nQr2A1gSM86kX0IGAcTL1AvoQMM6kXkA3AsZp1AvoScA4h3oBnQkYJ1AvoD8B41PqBQwhYHxEvYBR\nBIz3qRcwkIDxJvUCxvoaPcDLaq37je1u9/ngLs6lXsBwyQJWa73E6fr2xb7k8C7Ool5ABFOdQhSt\nDtQLCCLZEdgzbg6/LucVdyL3CfWCid3sLeMLHbBX27M//uZhinUW9YK5PfjoH1PogL3RHrlqRL2A\naEIH7N62bfeXGu7nDPflLkRsQb2AgJIFrByVaV+iWI2oFxDTVFchcjr1AsISMH6lXkBkAsYx9QKC\nEzAOqBcQn4BxS72AFASMH9QLyELA+I96AYkIGP+oF5CLgFGKegEJCRjqBaQkYKtTLyApAVuaegF5\nCdi61AtITcAWpV5AdgK2IvUCJiBgy1EvYA4Cthb1AqYhYAtRL2AmArYK9QImI2BLUC9gPgI2P/UC\npiRgk1MvYFYCNjP1AiYmYNNSL2BuAjYn9QKmJ2ATUi9gBQI2G/UCFiFgU1EvYB0CNg/1ApYiYJNQ\nL2A1AjYD9QIWJGDpqRewJgHLTb2AZQlYYuoFrEzAslIvYHEClpJ6AQhYPuoFUDIGrH578ICe83Sm\nXgC7r9EDvKbWun3vv69vXz+g+1D9qBfARbKAPbYn7aZhN7+9b14W6gU0le4AYKqAHcpbrGvqBbR2\nvbdMEbPQAXvp4Gl/8OXXObq1Uy+Ae6ED9lKEHn83lpd6ARwKHbB7119xzVqsa+oF8JtkAStHh2U3\nS6aJmXoBPJDv58AWoV4AjwlYROoF8CcBC0e9AJ4hYLGoF8CTBCwQ9QJ4noBFoV4ALxGwENQL4FUC\nNp56AbxBwAZTL4D3CNhI6gXwNgEbRr0APiFgY6gXwIcEbAD1AvicgPWmXgCnELCu1AvgLALWj3oB\nnEjAOlEvgHMJWA/qBXA6AWtOvQBaELC21AugEQFrSL0A2hGwVtQLoCkBa0K9AFoTsPOpF0AHAnYy\n9QLoQ8DOpF4A3QjYadQLoCcBO4d6AXQmYCdQL4D+BOxT6gUwhIB9RL0ARhGw96kXwEAC9ib1AhhL\nwN6hXgDDCdjL1AsgAgF7jXoBBPE1eoCX1Vr3G9tdSR7cddJTqxdAFMkCVmu9xOn6dvmu177k5q6T\nnlq9AAJJFrA/XWfs1D9WvQBimSpghwdnl/OKN495nnoBK7jZW8YXOmCft+ft/+pqBvUClnD/pUxw\noQPW6FqM56kXQFihA3Zv27b7Sw33E4aHd31CvQAiSxawchSny5ITj9jUCyA4P8h8QL0A4hOwW+oF\nkIKA/aBeAFkI2H/UCyARAftHvQByEbBS1AsgIQFTL4CUVg+YegEktXTA1Asgr3UDpl4AqS0aMPUC\nyG7FgKkXwASWC5h6AcxhrYCpF8A0FgqYegHMZJWAqRfAZJYImHoBzGf+gKkXwJTmD5h6AUxp/oAB\nMCUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAeuh\n1jp6hFJijBFhhmKMYDMUYwSbIQsBAyAlAQMgJQEDIKW6bdvoGRpyNhngPfHrMHnAAJiVU4gApCRg\nAKQkYACkJGAApPQ1eoDmah1wocrl6sf7p97v6jDSnzMMH6N0fHV+G6PnqojwvA+euvNIETaM+O+R\nbruLQ0N2ni+ZOWCjrqG/ftVvtoDLb1tvGY9nKN9viYFjlI4v0OMx+qyK5+cZ+NTdVkWEDSP+e6Tb\n7uJwqp5P97aZTyFu2xb848NAtdaxH6+CfLiLMEMQQVZFkA2jBHiPDJRl5znzEVhMY88J7EZ9/A/L\neriwKnZB3iMRdheRTRKwmwPezudknnzqdieOBv7vvzfG/sjLr+cO/NLasIO4iLAqmm4YuQw8z5zI\nJAEb+OpG2LAizFBeGaPpO/PVPzDI2otg+Kqwy+Yl828lEa5CvP4y9np5zxmCjPHb19RDxohw4B5k\n4+y/KiJsGMHfIwM3kstgwQMRfT4AODTzVYgATEzAAEhJwABIScAASEnAAEhJwABIScAASEnAAEhJ\nwABIScAASEnAAEhJwABIScAASEnAAEhJwABIScAASEnAAEhJwABIScAASEnAAEhJwABIScAASEnA\nAEhJwABIScAASEnAAEhJwABIScAASEnAAEhJwABIScAASEnAAEhJwABIScAASOn/aC+Gvj33lHwA\nAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["t = linspace(-1,1);\n", "plot(t, proxG(t,.3));\n", "axis('equal');"]}, {"source": ["The proximal operator of the indicator function of $H$ is the\n", "projector, and does not depends on $\\gamma$.\n", "$$ \\text{prox}_{\\gamma \\iota_H}(x)_i = \\text{prox}_F(x) = x + A^* (A A^*)^{-1} (y-Ax). $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["pA = A'*(A*A')^(-1);\n", "proxF = @(x,y)x + pA*(y-A*x);"]}, {"source": ["There are two kinds of Douglas-Rachford iterations, depending on wether\n", "you first apply the projection or the thresholding.\n", "\n", "\n", "The first algorithm, (DR1), reads:\n", "$$ \\tilde x_{k+1} = \\pa{1-\\frac{\\mu}{2}} \\tilde x_k + \\frac{\\mu}{2} \\text{rprox}_F( \\text{rprox}_\\gamma(\\tilde x_k)  )\n", "\\qandq x_k = \\text{prox}_\\gamma(\\tilde x_k) $$\n", "\n", "\n", "The first algorithm, (DR2), reads:\n", "$$ \\tilde x_{k+1} = \\pa{1-\\frac{\\mu}{2}} \\tilde x_k + \\frac{\\mu}{2} \\text{rprox}_{\\gamma G}( \\text{rprox}_F(\\tilde x_k)  )\n", "\\qandq x_k = \\text{Prox}_F(\\tilde x_k) $$\n", "\n", "\n", "The advantage of (DR2) is the the iterates $x_k$ always satisfy\n", "$Ax_k=y$, so that one can only keep track of the evolution of the\n", "$\\ell^1$ norm during the iterations. We will use only (DR2) in the\n", "following.\n", "\n", "\n", "Set the value of $\\mu$ and $\\gamma$.\n", "You might consider using your own value to speed up the convergence."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["mu = 1;\n", "gamma = 1;"]}, {"source": ["Define the rprox operators."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["rproxG = @(x,tau)2*proxG(x,tau)-x;\n", "rproxF = @(x,y)2*proxF(x,y)-x;"]}, {"source": ["Number of iterations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["niter = 500;"]}, {"source": ["__Exercise 1__\n", "\n", "Implement the DR iterative algorithm on |niter| iterations.\n", "Keep track of the evolution of the $\\ell^1$ norm."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAMAUlEQVR4nO3d0VLj\nuBqFUekU7//K/7lw4zZ2gDRDLG9rrQsKZ5iK2oR8kaIkvaoaAKT53+gBAMBPCBgAkQQMgEgCBkAk\nAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIj0NnoAr9V7a62P\nHgVAnqoaPYRv3Dxg7QK/g9778DEYxgWHcYUxGMYFh3GFMSzDGD2E71lCBCCSgAEQScAAiHSJxdbX\n6b3d+t8H8BIXeSrua2ZgAES6fcAittIA8M9uHzAA7ik4YP3d7sJR4wHgTMEBa61VVVWt0VIvgHkE\nB2y3QyZizwwAvyX7raSenHJtf0zkAB6KW8S6w6zleNLXf1TvXbMA/lXEmlbqEuIuWvWu6RXAHFKX\nELd7NxQLYEKpAWufdEvMACaRuoQIwOQEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACK9jR7Az/Xel2+q6ngIwL1lz8Cqqqq26doeAnBj\nwQHbzbRMvACmEryE2DbLhttLjiXb/pjOATwUt3yVHbDts1/L14d9Ei2Ab23vKiNilrqE+PDkChXA\nPFJnYLu9G9tJWFMygAmkBqx9rJRiAcwmdQkRgMkJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQS\nMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACI\nJGAARBIwACIJGACRBAyASAIGQKT7B6yq9T56EAD8tvsHDIBbEjAAIgkYAJHeRg/g5/r7U1tVdTwE\n4N6CA9Y+pmvtVu9dwwBuL3gJUaUAZpY9A+tPbJDvvbdWFhgBvvbMPeqlZAds9+zXZz/Tu24BfGN7\nPxkRs9QlxIiTC8DrpM7AqvarghYJAaaSGrB2CJVuAUwldQkRgMkJGACRBAyASAIGQCQBAyCSgAEQ\nScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAER6Gz2An+u9L99U1faS\n9RCAG8uegVVVVS3d6r1vDwG4t+CAmWkBzCx4CXGxTLzW79shbL331uq43gjAVtzyVXDAtrnaZWxb\nqarqXbcAvrG9n4yIWfASYnt6OlXVEn4XAPyD1BnYunFjOdzu3TDZAphBasCOldItgKlkLyECMC0B\nAyCSgAEQScAAiDRLwOykB7iZWQIGwM0IGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiTRQw\nb8YBcCcTBQyAOxEwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEhzBcz7\n+QLcxlwBA+A2BAyASAIGQCQBAyDS2+gB/Fx/349RVcdDAO4tOGDtk3T13jUM4PaCA3as1FIy9QKY\nQXDAFst8azvr2s3A+v6VX/IG8MDh3vLqggP25Hxr9wNpvyCAk3z50P+KsnchWi0EmFbqDGx5dLDd\nvmEXIsBUUgN2rJRuAUwlewkRgGkJGACRBAyASAIGQCQBAyDSdAHzocwA9zBdwAC4BwEDIJKAARBJ\nwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMg\nkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIFJwwPq73YWjxgPAmYID1lqrqqpao6VeAPMI\nDlhVbQ9777tLALixt9ED+K++7dZuWiZyAA/FLWIFB2w510uQlu/Xr9tKKRbAM7b3lhExCw5Y25zu\n9ZtnFhKrWu9N1wCipQZsO+VqplkA80kN2GfFUjKASQTvQgRgZgIGQCQBAyCSgAEQScAAiCRgLeHl\negDszR4w9QIINXvAAAglYABEEjAAIgnYn/f2BSCLgAEQScAAiCRgAEQSMAAiCRgAkQQMgEhTB6z3\n5gOcAUJNHTAAcgkYAJEEDIBIAgZAJAEDINKkAfMGvgDpJg0YAOkEDIBIAgZAJAFrzVNiAIEEDIBI\nAgZAJAEDINLUAfNW9AC53kYP4Of6+76LqjoeAnBv2TOwqqqqbbq2h9/9v68cGQAvFhyw3UzLxAtg\nKsFLiIve+zZdu8O2WVpc6BzAQ08uX11HcMCWc70GaXe4UiyAZ+wmAwNH8qTgJcRmFRFgYqkzsOXR\nwe4xgo2IAPNIDZhEAUwuewkRgGkJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgC9kdVS/j0AAD+EDAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgC9ldV6330IAB4ztvoAfxcf69NVR0PAbi34IC1\nj+lau9V71zCA2wteQlQpgJllz8Da+3yrf/7k1e4/yR7AQ1/ckV5TcMB2K4efUSyAZ2zvLSNiFryE\n2MQJYGKpM7Dl0cF256FdiABTSQ3YsVK6BTCV7CVEAKYlYABEEjAAIgkYAJEEbC/hxQ8ACNhHdjIC\npBAwACIJGACRBAyASAIGQCQBAyCSgAEQScD2qrwUDCCAgAEQScAAiCRgAEQSsAc8DQZwfQIGQCQB\nAyCSgD1mFRHg4gTsUxoGcGUC9hUNA7gsAfuGhgFck4ABEEnAvmcSBnBBAvYsDQO4FAF7StXf73sX\nM4Dx3kYPIMnSrSVmvX+oGgAnE7BnyRXApVhCBCCSgP2QrYkAYwkYAJHiA9Y386DeezctAphD8CaO\nXat671W1/QaAGwuegVWVUAFMK3gGdrTMyXZV203UfrF5yz4ODQXuIe4pmJsEbLtsuFtCNEsDeMb2\n3jIiZsFLiADM7CYzsKpaHy+YcgHMID5ga650C2AqlhD/E+/HATCKgAEQScAAiCRgv8AqIsD5BOy/\nOn5Ys54BnCB+F+J1bN+Vwzt0ALyaGdgvOL6nlN2JAK8mYL/j6/mWdUWAX2cJ8VW2k7Alb9YVAX6R\ngL3QLlfevR7gF1lCPJXnxgB+i4C93OEDycY07CIfjmAYlxpDM4yPrjCMK4whhSXEAR42zNIiwD8R\nsDGOufrsUZewATwkYFfxWaieXE7QOWA2/d4fozXNavKdf4nAENePw80DBsBd2YUIQCQBAyCSgAEQ\nScAAiGQb/Utsdz8u22TWS07bNdP73x06u2s/czDrMEadk2//7ScM44sr/WJgZw7jCr+Rk4exDmbg\nDWM3hjb6hvHZlZ5/9/UkAXuV7W96dwM94e9hd7j78zhnMMfXMIw6J9s/wt2VnjaMdQzrVQw5G9th\njLphtMPZGHXDWG+ixys9bRgX+TP54krbiTeMf2UJ8VV676NehVZVV7iRHYcx5Jxc5FQcLzz/bFzh\nVLRHwxhyw7jCffHDMYw6G4mvmjUDe5XjI24GnpP1YfXJ13scw3o45Gxc5E7q8A7X/lj+GnI2Qn8F\nAvYSWTeCc4w6J7s1kIuMYdR4jk87DR/GkFOxXPX26/AxHFdTT5N7f2UJ8fcNv3e4oLHn5Ap/n7vn\nGM4fwEVulscnaIcMo961oQ8mdmOY+YbxM2ETxhQDdzRtx3DNXYhDtv+tV3T+MD5bLjt5GM9c6Tz7\nQreDGT6ML/5UzxnGRW4YPyBgAESyhAhAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyDS/wHryy+35m8+8gAAAABJRU5E\nrkJggg==\n", "output_type": "display_data"}], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["We display the convergence speed of the $\\ell^1$ norm on the first half iterations, in log\n", "scales."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAOe0lEQVR4nO3da4+b\nyBaGUXOU//+XOR/oIQRf2sZc6t21lqJRpqMeOYzbj3dR4GEcxxsApPnf1Q8AALYQMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJ\nGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEOnP1Q/gWMNw\nu92Gqx8FQJ5xHK9+CL8oHrDbbUj4v3C4YRgchZmjseRozByKpWEIeOtvCRGASPUDNo63hHcSAHym\nfsAAKEnAumBlf8nRWHI0Zg5FnC4CZhURoJ4uAgZAPQIGQKSCAXt4+YJVRIBiSl3I/OuVd8Nwc5oW\noIZSAZs2ET3L2JSu+Q+VDCBaqYA9tOrZvFPWNAawFHH7qKX6AXt2bcd8VkzGAG7/vlpGxKzgJo73\njaPNHQCpug4YALkKBuzT+8EYwgASFQzYBhoGEEfAAIgkYABEErAfVhEBstS/DuxL7twB0CYT2F/3\nQ9h0tw6XiwE0SMDepWEATRGwfywr5WaJAC0TsMce1ssQBtAOAVubumX2AmicgH3GEAbQCAEDIJKA\nARBJwD5mFRGgBQIGQCQBAyCSgG1hFRHgcgIGQCQBAyCSgG1kFRHgWgIGQCQBAyCSgG1nFRHgQgIG\nQCQB+4ohDOAqAgZAJAH7liEM4BICtgMNAzifgO1DwwBOJmAARBKw3RjCAM705+oHUMqyYeN46UMB\nqE7AdjZ3S8kADiVgR1EygEMJ2OGWJdMwgL3YxHGe6QyZjR4AuzCBnWqawKaGPdy1aEQDeJOAXWDO\n2H2u5rYB8JolxMs8rJR0AbxJwJrjgmiAdwhYizQM4FelzoEN/73qj1biAKqrE7BhGOZuLX8PQEmW\nEBtlFRHgtToT2DPDvx0wmQE8NKS9a64fsNxiTUNY7MMHwixfLSNiZgkRgEh1JrBxHO1CBOhHnYDd\nKnbLKiLAM5YQAYgkYABEErDWuSAM4CEBAyCSgAEQScACWEUEuCdgAEQqdR1Yh5aTmcvFgK4IWIb7\nK5qndM1fscYI9MYSYqQpZsueOU8G9MYEFmOZqIerhe47BXRFwJL8GicNA/phCRGASAJWjZNhQCcE\nrKBnDRuGn18ABTgH1oWHe+6dKgOimcBqmoewaeS633MPkE7Aypp3JL7Ycw+QS8AqM2kBhQlYvwxh\nQDQBAyCSgHXNEAbkEjAAIgkYAJEErHdWEYFQAgZAJAEDIJKAYRURiCRgAEQSMAAiCRi3m1VEIJCA\nARDJB1ryY/74ldlqJnNve6ApJjD+MUdr/tTmZ58oBnAtExh/TaGa07X6o9V8BnAtExhrPsQZiCBg\nfEDDgHYIGACRBIzPGMKARtjEwcdsuAdaIGBsMTfs4ZZF+xWBExRcQhyscJ1ibth9qywzAicoNYFJ\n18lejFmuGwOOVmoCG8dx9JIJ0IdSE9hDq7FM4U5jCIMscYtYqQF7P0uKBfCO5atlRMxSAyZLEQxh\nwHFKnQMDoB8FA2Y4a4ot9cBBCgaM1mgYcAQB4wymYmB3AsZJDGHAvlJ3IRJnuZBoIAO+J2CcZ+6W\nvfXA9ywhcgHbOoDvCRjX0DDgSwLGZTQM+IaAcSUNAzYTMC6mYcA2AgZAJNvoud6Lm9YvhzM774El\nExhNeLiQOFVt+gWwImC0YmrYlLHpN8tuOVUGrFhCpCFTsZ4tJ/p4TGDJBEZzJAp4h4CRxEIiMBMw\nACIJGGEMYcBEwMijYcBNwAAIJWBEenbhM9APASPVsmHzhc/zpdBAeS5kJtjcsPnSsdeXQgOVmMDI\n9vBOiXZ5QA8EDIBIAkZNhjAoT8AAiCRgAEQSMMqyigi1CRgAkQQMgEgCRmW/riJOd+6w0giJ3ImD\nfi1v2OHmHRDHBEanVsWy4wPiCBjFKRNUJWD06OGCodRBFgEDIJKAARBJwKhvtTb4YsOhVUQIImAA\nRCp1Hdjw35vn0RU9POF6Lyij2gQ2juM4joNlIP6119qgZxa0o1TADF68MDXs1+fIi9RN3+7WU9CI\nUkuIk2EYliVbTWMi17Nv/ufP8bPRg6ri1q5SA/YwS9MXV4lSLHb35jAHWV689W9TasCeZUmu+N59\nn+5zpWFwudSA3ZveL9iICNCJOgFTLA5i0oI2ldqFCGeymwOuJWDwwBwn4xc0S8AAiCRg8Nh81RfQ\nJgGD7ZwGgwvV2YUIDVrmzTAH+xIwONbcLftBYF+WEOEoigWHEjD4yvunwZwwg30JGBzC+AVHEzAA\nIgkYnMcqIuxIwACIJGCwPyfA4AQCBt/6aGHQKiLsRcAAiCRg0JZhMKLBWwQMdvbNCbDpey0zwjsE\nDFqxLJ+Gwa8EDM4mTrALAYMm3C886hy8JmCwA7GB8wkYtEsX4QUBgz25BwecRsDgerIHGwgY7MMN\npeBkAgZAJAGDpr0zq013nzLS0Zs/Vz8AqGPbeawvT4DN3y5g9MYEBsHcfYqeCRjUoWF0RcDgGmID\nXxIwaN2z1Ll6jM4JGFxp9wgZ7OiHgAEQScAgkvVDEDAAIgkYXGYcD5minAajEwIGQCQBgwD7DlXu\nnUgNAgY9sgGEAkrdzHf47y3l6KcTnpi3L05TnZ8VclWbwMZxHMdxsDhCaZvDs/pG2z2IVipgBi+Y\nKBM9KLWEeFusIj77isgBPBS3dpUasGdZmn6z/FPFgsnDhUdnwpgtXy0jYpYasPssDcOgVQD9SA3Y\nveXeDSUDKK9OwG66BZtYRSRUqV2IUNiXGwslinoEDMJckiK3nqJBAgY17Xgp2JRM15bRGgED3qVh\nNEXAgFdlcvcpmiVgAEQSMOApexdpmYABn7GKSCMEDIBIAgZJtq3pWQmkJAEDbrdHC4OyR+MEDPiY\n02C0QMAAiCRgUJY5idoEDGJknZGa7v+roByn1OeBAY1YbgCxGYSDmMAgybY99Bv8Wp33b58IBzGB\nQWUthMQnPnMQExgAkQQM+LHLrkXDFqcRMOBwNvRzBAEDIJKAAbuxfsiZBAzY6KOFQauI7E7AAIjk\nOjCgCfN8ZhGSNwkYcD23nmIDS4jAxVbFcraMNwkY8Nc38fjm9omwgYABaxbxiCBgQHPMarxDwIAr\nmfbYTMCA7YxKXEjAgPO8Hzxp5FcCBlzG+iHfEDAAIgkY8I9tI5FZivMJGNCod06DDYNTZf0SMGAt\nZZaaxj7bPbolYECk5aKlhvVJwIBTzbHZ97SZhnWoYMAGz2Kozp4RbvUCpl5QibmKF0p9oOUwDOM4\nrhq2+tfR2zYoaqqdH/HN4gaAUgF7SLHgUNuy8c1LpUodZPlqGRGz1IDdz1XTV+Z/6hY067ifTkNY\nV1IDdt+n+SvqBdCDaps4AOhEwYAZvwB6UDBgQCWrnfS/nuJ6vfPevRMrETDgWymrHu6dWIyAAV1w\n78R6BAzokYYVIGBANfdxcnFYSQIGxNAhlgQMaN1By31WEdMJGFCcua0qAQP6ZQiLJmBAQcrUAwED\nKrN+WFjq3eiBrpz/qWMP/zta2BQBA2p6Mza/pnH5p+a5plhCBHhKsVomYECGFkJib0hTBAyASAIG\n9O7ZXPVw/dAQ1g4BAyCSgAE88GL7hiGsEQIGQCQBAyCSgAF8zCpiCwQMYB0k1y9HEDAAIgkYAJEE\nDGCLF6fBhuHnF4dyN3qAPbl7/WlMYAC322Ki+qY6q++1WfFQAgZAJAED2OidzfeGsOMIGACRBAzg\nxzQtbTsB5ua/5xMwACIJGMBfn45fpqsLuQ4M4CvfLDwuLUPo6rF3mMAAvvVrb34d1KYETr94k4AB\nXMzlz9sIGEBzNOwdAgZwBk3anYABXMkNfzcrtQtxWLy9GT0jgGR7bW4srFTAbroF0I1qS4jDMAyW\nmYEQZqxv1JzAhmGYR7FVz4xoAA/FvftPDdjDLD2Mk2IBjfj0tNbJp8GWr5YRMUsN2H2WllMXAOWl\nBuzeOI7zWwYlAyivTsBuugVEeWd50Gb6F6rtQgSgE6UmMIDG7ThR+fgVAQPIs6pgn8uMlhABLtBn\ncvYlYABh7uPX563uBQygaX3G6R0CBnCqg4LUYecEDOACm0+AOXk2swsR4GwXFqjS5nsTGECM1+PX\nr6uI07dPvwoQMIDWHXF+q8A5MwED6EK9k2cCBhDgiLv6pg9hAgZQR3qTPiJgABm+Gb/qrR/eBAyg\nZ9ETm4ABEEnAAIoruX54EzAAQgkYQCmfntbKPQ0mYABEEjAAIgkYAJEEDKCyd7Yghp4GEzAAIgkY\nQDWhE9WnBAyASAIGQCQBAyCSgAGUVfUuiBMBAyCSgAEU1MMdEQUMgEgCBsC9gFNnAgZQVuEdHDcB\nAyCUgAHUtHn8Gobb7Rawo0PAAIgkYF0Y4rbHHsnRWHI0Zg5F3E56AQMgkoABEEnAAPgxrSKmbL4f\nxpRHuolFbYAPjdMWxPbrUDxgAFRlCRGASAIGQCQBAyCSgAEQ6c/VD4CjLHdgTlt15q/0tnNnGP5u\nVlodhA6PyXw0en6G3P9N3/lKVS/+7tMXmz0UAlbZ8tm2ehFv7Yl4kNV1FPcl6+qY3F9V0vMzZG72\n9AK9+rv3fDSWX7m1/cSwhFjZMAydXwk3jmNTP2/Xuj8a3T5DPCuWHh6NiOeGCayy+3dVsNT5M6Tb\nv/hDq6Nxv87cIAEry08mr/X8DFktIHfu/mikHBlLiDU1/r6Jy3mGpLxGn2N1NvTCR/IRE3RZPe+q\nWrELcel+F2Jvz5DVC7SjsfzXh0ej2UMhYABEsoQIQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJ\nGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMg0v8B+elfg7n1C3YA\nAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["plot(log10(lun(1:end/2)-lun(end)));\n", "axis('tight');"]}, {"source": ["Display the original and the recovered signals.\n", "Since the original signal is highly sparse, it is perfectly recovered."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAATA0lEQVR4nO3d0ZKq\nPLcFUDm13/+VORfWZ/GDYKRB1gxjXHW3CouITBIiPYzj+ACANP93dQEAsIcAAyCSAAMgkgADIJIA\nAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAOOOhomPz2xZ2l8q2XjobZF/\nWd2OMqCsf1cXAD/1PFJP/w3e8i9TLf8w7/B/qve2yNev/ocfPOmBcTuzAKiZB8uqdJJgRoBxI9N+\nzNQ4jq94eP7wGribxsZ0QG/29+Vrl+N+jYOWb43j+Kq8vaQzKoE6DCHC3FrOPf7rGG1fuJomzfT5\ny79v1/D4ZgBzWdJRlUBZAgzm3h7TZ1eh1jJsrYe3o4BZJ+ntoOJGSYdUApUJMPiRb4fsZmOGB/aW\nDB7SBwEGp3sFxscRyLf2veqMSqAUkzi4kbWhvx9cCprOwvjo1Fz5qhKoTIBxO/u+F7ycqXjI2huf\n9jZi/1iS7hfpDCFyL8uhs6+6I69ZfF8d/afPf/683ed7PaelyK9K+rYSqMy+C3sUPO4XLAlOZQgR\nmsy+JlwhKgqWBL9kCBGafDts+AMFS4JfctYGQCRDiABEEmAAROoqwFwPALiPrgKMp2OD/LzTgu17\nuu9+7YEv+eMLGxd+7YnXJWs/aqXTfyXz4w1Z/lObX679ED1sQh+TOGZ3eLut5S64u0EOXNRXC/94\nS/UdhbXcpv3YF+5Y+OHLr7n2o9a7cbQ9e0O2D/T1D0Efg6r+Jrx00gNze7fH+afSB541tz+077ZP\nXy3z8BfuW/ixy9/hB2v/wS504Fr2LTyxHzMTtAmdBBjTPujaf+/dvbRDKny75GmdX21C4zbubpZj\n23PDJedea1t3yUrPOFyedAje+FyknECvbUJoH6CTIcSn5c0Igk4lAEqpnw7934mj/nvwOOI+QLM7\ntE4f2rHk5RKOusq4Uefb1c1Wunbj2uWJy+ymvd82y7HtubH8mUsuRH18fw+5T9WBTdp42+KvlvnV\nerfvgdJyUfaq41LjOX3KTV4MIXaifRLEDgfuyht1/mUex47VnfTCCGsDFadu3S+b9KQNOXXY8zf6\n2IFfBBhvZO3lWdU+va05cUMuUbmhKtfWKGgTBFg/jr2wvFzUUbv1Rp0fLyzv2MbdF6vPvsp97SyA\n897f7ZUetdVrLzx7Q5Z7RUtVpWzXHLEJLwKsN9MwqLOorxbePr7XXtjhB8pjXXLUOPX9/cF6LzzU\nvs2woEP/rObETXjqKsDiWv8lpXJ1Hkud96Q9j9JVgAFwHwIMgEgCDIBIAgyASAIMgEgCDIBIAgyA\nSAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBI\nAgyASAIMgEgCDIBIAgyASP+uLuCDYRheP4/j2PK07WcC0IfqAfb4L41mETU1DIPEArib0kOI02Qa\nx3EjwwC4m4AeWIvGkUYAutFJgD0mI42zEUWXxwBaxI1y9RBg00xajjRKLIAWG6f+NZW+BgYAa3oI\nsIgzBQCOVTrApuOBjXPlTakHuImAa2BrHaxXVk1zTnoB3ET1ANsIpNncjZ+UA0AVpYcQAWCNAAMg\nkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCS\nAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADINK/qws4xjAMr5/HcbywEgB+o58e2DiOogvg\nPnoIsGEYXtE1juO0NwZArzoZQmy0MdI4i72NR9sf2v3CQ5a58ardazxjA6fnH7NHNx56q2WNH5e5\n74W7d6f2TWhf41fFPNbf37OL2d2kG2MtBy7z1N3pjGK2H/q42Ja/VzNv/USzfWj6a8rbAFBQ8YDo\nYQiRUBufjTM+NruXWfwznOjwJv39JXB7RQX998A62ECAH4s4eOqBARCphwCbzjyMOGsA4O/6mYVo\nvgbArXQSYHpdAHfTwxAiADckwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACI\nJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIj0\n7+oCPhiG4fXzOI4tT9t+JgB9qB5gj//SaBZRU8MwSCyAuyk9hDhNpnEcNzIMgLsJ6IG1aBxpBKAb\nnQTYYzLSOBtRdHkMoEXcKFcPATbNpOVIo8QCaLFx6l/T9QH2tpmkDgDbrg+wv2eVWYgAN1R6FuJ0\nPLAxpYQZwE1c3wP7aG0o9pVV05yTXgA3UT3ANgJpNnfjJ+UAUEXpIUQAWCPAAIgkwACIJMAAiCTA\nAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAA\niCTAAIgkwACIJMAAiCTAAIjUVYANw3B1CQD8SFcBlisletV5LHUeS5138+/qAo5hhwC4m056YOM4\njuN4dRUA/E4nAQbA3Qw9dVyGYb45hhYB9qmfDp1cA1tT/w0AYJ+YAHvbl5JPALcVE2CyCoApkzgA\niCTAAIjU1SxEAO4j5hoYv7T8QkJN9eucTj6qXGpcnZWLfOTUma7nAKv8mZxNqnyVl1LztLZqNddp\n2+can+sahmEWt3Xas73OR6X23PjSZ6n2bGnMR73PUYSeA+wx2YeuLuR/bPcbrq15e71rtf2+5rV1\nVWvbVzHjOLa3W506i7fnrLay7bkWum9fVe14VVm3kzim+8rbY0dBFWpeu6vkWm1X1fzt3S8rtG1L\nPdXqXFOtTu15T533wMoyXHCeOm2b8s5u11mzPSsf4jfqrNOYfei2B1bfqw9R+aMYqlrbLi8s1bRW\nZ6n2fBb5KB8Aa3WWasx0emAXmO7QhguOVbBtI462j5U6C7bndBZJ5VZd1lmwMdPpgcGJotOrvpQM\nSKkzjgC7gL35PKXaNiUVNuos2J71bUyR/XEl3es2wNZmIlXTMnuqgpaZXdVqrlDn2oqqtWfLiirU\nOfPVzMMKdb7MiilSZ5zOr4HVPOWZ7qxv9+yfV9Sq/qlltbbd+OLq8tGPfz/PWp2l2nN23ai9niJ1\nbjfmo9LnKIWoByBSt0OIAPRNgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQ\nSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgHGK4Z2ri9pprfK/\nbNHH1/69uVqW8O1aDnkTc/cEqvl3dQF0axzH6a/PDJv9kfO0NPW3b8fut2/61tsHOIoAg1p0UKCR\nAOMyryP1sq/W8vfnr9PnzHp4y+U8nzB74dvVrdXwcSs+VvXRdDkbK31tyPMlaxu79vzZc1oWOHv+\n27JnD02f/GqZj2/TRsEwJcD4nbXD+tufZ0OOy5+nB7jH4pD6cfkfH2oPkvaqDjRb9bTgt4f7jS3d\nscCNNnwsAnIZWi1VfSwYTOLgLLMZHBudktdBf+1KydrzZzaiaLs/0VjGxhpbqjrWbLEf1/LtExrL\nnm3yM65aXvj49LZKLD7SA+MsywHAM6YM1Dw3r1nV4ZbjpdO/w9kEGNc45Bg37TPtW/4Zh9q3Pbn+\nvE2vt1cE4SQCjGsc3rvaNwLWccD8wNvW06T8jGtgXGBtpsP077PLIbPnv53G9nH5fyxj6auqetKy\ndR+fs/22wkd6YPxOy/zptb+vTcTYnuH2dvmNT5tO+1577Y6qluudvmrj+XWsnWQ8Fg0ynb/zdsjR\nXHn+wikPnEJ/As5mCBGO1/f4IRQhwOBga5PLgWMZ5QAgkh4YAJEEGACRugowV84B7sP3wHpz7Ldq\nzpuPsFHnx03YsY27m+XsbyldO+Pjku9gHbjSq75L9/aOWVnzCTrYhEc3kzh8F/LpwM/zqYeGtYW3\n3M3928J23yH+7FvLX/tF5ssP/X9c7/Zwy3nb0nJ7kZNWfZQONuGlkyHEr/6JQ6/OHkE9avl/WU43\no8RvN+TarfvB2ivsQmerXFujoE3oJMCY9kGncb5vX1xb2t9t1Nm+Ce3buLtZjm3PDZece/1s665a\n6UkbMtuE6UMpJ9AdbMJUJ0OIT8ub9wSdSgCUUj8d+p/EUf89eBxx37y3d1J//WX30pbLPK/Ot6ub\nrXTtBrvLE5fZRbVvm+XY9txY/swlF6I+vr+H3NfxwCZtOTE9oyVb9t6WVV94n8zGc/rtDazDEGIn\n1gZkTj3a7rBRZ/smtNezu1lObc/LrQ1UnLp1v2zSkzbkN2Otp+pjB34RYLyRtZdnVfv0tubEDblE\n5YaqXFujoE3ofwjxPpaDZocv7ZATz406P27Cjm3c3SzHtufb5S//n9bPnL11Z690bYzr7A3ZHpqO\nOPR3sAkvemC9SRmQ2Vhm+/jeIas76YUF11Jk7Uet7sJGW551Pa5+E/d51py7CV0FWFzrv6RUrs5j\nqfOetOdRugowAO5DgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmA\nARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARDp39UFfDAM\nw+vncRxbnrb9TAD6UD3AHv+l0SyipoZhkFgAd1N6CHGaTOM4bmQYAHcT0ANr0TjSCEA3Ogmwx2Sk\ncTai6PIYQIu4Ua4eAmyaScuRRokF0GLj1L+m0tfAAGBNDwEWcaYAwLFKB9h0PLBxrrwp9QA3EXAN\nbK2D9cqqac5JL4CbqB5gG4E0m7vxk3IAqKL0ECIArBFgAEQSYABEEmAARBJgAEQSYABEEmAARBJg\nAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAA\nRBJgAEQSYABE+nd1AccYhuH18ziOF1YCwG/00wMbx1F0AdxHDwE2DMMrusZxnPbGAOjV0EGvZRpg\ns19fP6+l2kbgbbzwjGWWKubjMksVk97apYqxgfteuCGltZcvrJ8O/QfYRUUBxCseED0MIRJq47Nx\nxsdm9zKLf4YTHd6kv78Ebq+ooP8eWAcbCPBjEQdPPTAAIvUQYNMLjxFnDQD8XSdfZH6YrwFwM50E\nmF4XwN30MIQIwA0JMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkw\nACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACL9u7qAD4Zh\neP08jmPL07afCUAfqgfY4780mkXU1DAMEgvgbkoPIU6TaRzHjQwD4G4CemAtGkcaAehGJwH2mIw0\nzkYUXR4DaBE3ytVDgE0zaTnSKLEAWmyc+td0fYC9bSapA8C26wPs71llFiLADZWehTgdD2xMKWEG\ncBPX98A+WhuKfWXVNOekF8BNVA+wjUCazd34STkAVFF6CBEA1ggwACIJMAAiCTAAIgkwACIJMAAi\nCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJ\nMAAiCTAAIgkwACJ1FWDDMFxdAgA/0lWA5UqJXnUeS53HUufd/Lu6gGPYIQDuppMe2DiO4zheXQUA\nv9NJgAFwN0NPHZdhmG+OoUWAfeqnQyfXwNbUfwMA2CcmwN72peQTwG3FBJisAmDKJA4AIgkwACJ1\nNQsRgPuIuQbGLy2/kFBT/Tqnk48qlxpXZ+UiHzl1pus5wCp/JmeTKl/lpdQ8ra1azXXa9rnG57qG\nYZjFbZ32bK/zUak9N770Wao9WxrzUe9zFKHnAHtM9qGrC/kf2/2Ga2veXu9abb+veW1d1dr2Vcw4\nju3tVqfO4u05q61se66F7ttXVTteVdbtJI7pvvL22FFQhZrX7iq5VttVNX9798sKbdtST7U611Sr\nU3veU+c9sLIMF5ynTtumvLPbddZsz8qH+I066zRmH7rtgdX36kNU/iiGqta2ywtLNa3VWao9n0U+\nygfAWp2lGjOdHtgFpju04YJjFWzbiKPtY6XOgu05nUVSuVWXdRZszHR6YHCi6PSqLyUDUuqMI8Au\nYG8+T6m2TUmFjToLtmd9G1Nkf1xJ97oNsLWZSNW0zJ6qoGVmV7WaK9S5tqJq7dmyogp1znw187BC\nnS+zYorUGafza2A1T3mmO+vbPfvnFbWqf2pZrW03vri6fPTj38+zVmep9pxdN2qvp0id2435qPQ5\nSiHqAYjU7RAiAH0TYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJg\nAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAA\nRBJgAEQSYABE+n9456pZj0OJRAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "subplot(2,1,1);\n", "plot_sparse_diracs(x0);\n", "set_graphic_sizes([], 15);\n", "title('Original Signal');\n", "subplot(2,1,2);\n", "plot_sparse_diracs(x);\n", "set_graphic_sizes([], 15);\n", "title('Recovered by L1 minimization');"]}, {"source": ["__Exercise 2__\n", "\n", "Test the recovery of a less sparse signal.\n", "What do you observe ?"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAdUUlEQVR4nO3d2bba\nuhIFULgj///L3AdGfBQ3Qu6r5DmfEjaYQm6WJQvz/nw+LwDI5n93FwAAWwgwAFISYACkJMAASEmA\nAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQHGE70LP5/ZsrQ9lVT+NFvk\nnrfbUAaE9efuAuBS3yN1+TN400dKLT+Yd/iP6s0WOfzXb/jBlx4YjzMKgJh5MK1KJwlGBBgPUvZj\nSp/PZ4iH7z+GgbsyNsoBvdHj09dOx/0aBy1nfT6fofL2ks6oBOIwhAhjSzn3+tsxql+4KpOmfP70\n8XoNrzUDmNOSjqoEwhJgMDZ7TB9dhVrKsKUe3oYCRp2k2UHFSkmHVAKRCTC4yNohu9GY4YG9JYOH\n9EGAwemGwPg5Ajlr26vOqARCMYmDB1ka+rvgUlA5C+OnU3NlVSUQmQDjcbZ9L3g6U/GQd2982mzE\n7ixJ94vsDCHyLNOhs1XdkWEW36qjf/n877/rfb7hOS1FrippbSUQmW0Xtgh43A9YEpzKECI0GX1N\nOEJUBCwJrmQIEZqsHTa8QMCS4ErO2gBIyRAiACkJMABS6irAXA8AeI6uAmzWqp+NmP5QxTlFtSp/\nFGPVS86rp/zxjtG3aGd/uaNxmbOPLz1/bdktNbQ8Z1Wd7Us40O1b7DXe/1r1qqOeNjx51UuWNqGd\nK+7KnSWUTiZxjO7wNnpw8PPboEuub6Wleto/wrE179wTlu7vPvu02cdXrc1GjctcVeeeN9rjgreI\noLIdrtq7V63oVZWctAk1lnHUYrPopAc2vb3bsSf1F5/IVN6u/SMcWPP+RbWXd29f87I6jy2++xPt\nr/rH3Lmmjjpi7Hm79pcvPflpXbH+vwfW8guES8/M8j2baQd0GNm48cyrsZiy2SsfZPbxPdYus7HO\n+g2iDql8JObaD2i0Bjf/ZulajVv12nc8abG5dLWJv//9Cdp7iwFILX469N8DG5n9WdvhT9tG2Pcb\nnQ+uvSBX+QiHlL3tbKBezLt6C9qf62L2qme58J8ffFpb5Upqe52VtTN6yfcfB17JO2ntv8J05lq2\nw6X2/9k4ldW0tpLGJU9fVXnHSg2NO8vOq2sxdXINbI8OOm1LH+Heg05LMY1D9ktHhD0fcNUyt9XZ\n8pL9Yq79aDYP0O3Xcnl+wzs2HrjyHtZadBtglcstRy3tFrdUctmbxmnnulV1HrspPlm90XaulKNW\nU/tLTtoGHrVphRgZOMp0oGM0I+Pny2cfv6CJZodoNsyFPaRr8nPhP42uHjfORV4aRJ0906x/wPYh\nr8ZlrqqzfQnHDs2t2trXLjnOgWJp/O1VrXPnim6spH3a/Yat+ueSt22ZlcXGWelLuu2BjbSsieE5\nn8+n/PeJZbXV017JqdWWC59tq+mf2pfZspAzVkrjMlfVueeN9oiw0V6gsfO09JxtK7qlkrM3ocYy\njlpsCl0FWN5VlbfymLK0pzqPpc6n6SrAAHgOAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIA\nAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQk\nwABISYABkJIAAyClP3cX8MP7/R7+/fl8Wp5WfyYAfYgeYK+/aTSKqNL7/ZZYAE8TegixTKbP51PJ\nMACeJkEPrEXjSCMA3egkwF7FSOPsiOKQcOINYFa6Ua4eAqzMpKWRRrkFUFceJ1OEWehrYACwpIcA\nS3GmAMCxQgdYOR7YOFfelHqAh0hwDWypgzVkVZlz0gvgIaIHWCWQRnM3LikHgChCDyECwBIBBkBK\nAgyAlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQ\nkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQ0p+7CzjG+/0e/v35fG6sBIBr9NMD+3w+\nogvgOXoIsPf7PUTX5/Mpe2MA9KqTIcRGlZHG8k9lIq564bQLuO2F9QxueWH9I5xUzNKffhbTvszy\nCUe19ugJ7aW2/+nizSlUMRdsh+3LrLzwmtZ+tW1p24qZ7miNL6z/KbKZD5zOaLWV/82yGgACCh4Q\nPQwhHij42mpx9kdYtfwzismyTE5llZUe2xr998A6+IAAF0tx8NQDAyClHgKsnHmY4qwBgP36mYVo\nvgbAo3QSYHpdAE/TwxAiAA8kwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYAB\nkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJg\nAKT05+4Cfni/38O/P59Py9PqzwSgD9ED7PU3jUYRVXq/3xIL4GlCDyGWyfT5fCoZBsDTJOiBtWgc\naQSgG50E2KsYaRyNKLo8BtAi3ShXDwFWZtJ0pFFiAbSonPrHdH+AzTaT1AGg7v4A259VZiECPFDo\nWYjleGBjSgkzgIe4vwf209JQ7JBVZc5JL4CHiB5glUAazd24pBwAogg9hAgASwQYACkJMABSEmAA\npCTAAEhJgAGQkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlAQY\nACkJMABSEmAApCTAAEhJgAGQkgADIKWuAuz9ft9dAgAX6SrA8soSveo8ljqPpc6n+XN3AcewQQA8\nTSc9sM/n8/l87q4CgOt0EmAAPM27p47L+z3+OIYWAbaJnw6dXANbEn8FALBNmgCb7UvJJ4DHShNg\nsgqAkkkcAKQkwABIqatZiAA8R5prYFxp+oWEmOLXWU4+ilxqujojF/nKU2d2PQdY5H1yNKlyKC9L\nzWVt0WqO07bfd/y+1/v9HsVtnPZsr/MVqT0rX/oM1Z4tjfmKtx+l0HOAvYpt6O5C/lHvN9xbc/19\nl2q7vual94rWtkMxn8+nvd3i1Bm8PUe1hW3PpdCdfVW041Vk3U7iKLeV2WNHQBFqXrqr5FJtd9W8\n9u6XEdq2pZ5odS6JVqf2fKbOe2BhGS44T5y2zbJm63XGbM/Ih/hKnXEasw/d9sDiG/oQkXfFpKK1\n7fTCUkxLdYZqz2+Rr/ABsFRnqMbMTg/sBuUGbbjgWAHbNsXR9rVQZ8D2LGeRRG7VaZ0BGzM7PTA4\nUer0ii9LBmSpMx0BdgNb83lCtW2WVKjUGbA946tMkb24ku51G2BLM5GiaZk9FUHLzK5oNUeoc+mN\norVnyxtFqHNk1czDCHUORsUEqTOdzq+BxTzlKTfW2S378opaxT+1jNa2lS+uTv/68/HzLNUZqj1H\n143a6wlSZ70xX5H2oyxEPQApdTuECEDfBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICU\nBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAOMU\n7zl3F7XRUuV7PtHP1+5vrpYlrH2XQ1Zi3i2BaP7cXQDd+nw+5X+/GTZ6kPO0NPXa1bF59ZWr3jbA\nUQQYxKKDAo0EGLcZjtTTvlrL49//ls8Z9fCmy/k+YfTC2bdbquHnp/hZ1U/lcipvOnyQ70uWPuzS\n80fPaVng6PmzZY/+VD55aJmfq6lSMJQEGNdZOqzP/ns05Dj9d3mAe00OqT+X//NP7UHSXtWBRm9d\nFjx7uK980g0LrLThaxKQ09BqqepnwWASB2cZzeCodEqGg/7SlZKl549Uoqjen2gso/KOLVUda7TY\nn++y9gmNZY8+8jeuWl74+rVaJRY/6YFxlukA4BlTBmKem8es6nDT8dLycTibAOMehxzjyj7TtuWf\ncaid7cn1Zza9Zq8IwkkEGPc4vHe1bQSs44C5wGzraVIu4xoYN1ia6VA+ProcMnr+7DS2n8vfWcbU\nqqp60vLpfj6nvlrhJz0wrtMyf3rp8aWJGPUZbrPLb3xaOe176bUbqpq+b/mqyvPjWDrJeE0apJy/\nMzvkaK48ezjlgVPoT8DZDCHC8foeP4QgBBgcbGlyOXAsoxwApKQHBkBKAgyAlLoKMFfOAZ6jqwAj\no+9Xhe6ugpRsPA/XySQO34XMKO83ebmdjYdXNz2wVT/iQASzJ87Opmlh4+GrkwAjLycfbGbjebhO\nhhC/pjfvcVIGsE38dOj/Zr7x18Erz33zDqxz6Va5hyz/ge15qmh11m+vfGUl20Rrz1kpzv4NIXKn\nFDsJEJMA4x5+C5HNbDx8CTBus+03lOFl4+H1egkw7jUcdxyAWMvGQ1cBlnc7zlK5Oo+lzmOp82m6\nCjAAnkOAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFLq/2a+qfnVPjjKcP/ouwvhMAIsrqVf7bMHwirl\nrlSPsfKZdrT4BNi8UNvxtwA3bocN2nec0TOdLMbnGtiM6XZ8VyUvN3yDg5Q/3zzdqZcGPE4vix30\nwMaGTbbs99x7LmYvYq0zfiY0telPp+5/JrfTA5sXp99jX2Ktxu4Fs7RVIgJs0fv9DjJ4+PNBGJRD\nCLaWweZ9WZ5FZgjxP9OpSq+7N9/RaEbk41GWOh+iHEJwCJ6yifZBD+z1auhs3bi5xxnMrAg17YWv\n24cQgpj2RI1tdEOA/WbLrjN9K6AgQwhxbDgRtOPHZwhxPO3wZTRsk8oEZa5kzPAomjE+PbB/GHXZ\nQ+sFMT3rch7WwtBiOgKMY4iuUMoRM4fgdi1Xy4hDgBn7AmZIr/gEGHsZeAFuYRLH6zW5W66b565V\ntpjoAq6hB/YfR16ARAQYAClFH0Js/EqWe28/gbUMIw/fKaIH2KvhipTfnXuC2TusW+/slPpSt50i\n9BBiuTLcX+BR3n8N//3+w7eaOFDqQ8rS6FTqD7VWgh5Yi9l1OfwW5fRPm5fv6Hm2+n2Bh4mOEX5o\nlG6knnj8Lf6Zp/idBNir2AS/B7WlO8yuPd4tHUB3VArcrzwrzX7oP6r+dO3QQ4CVcTLdFlOfW5Wi\n9QLPyPL6jZUrT4OnKcchRo/vWebw7xTHzNDXwA4x9Mn2LOTz17DMI0rb6PYNa3R16qR63BqYUwXZ\nndmjhx5YSz7t30bv6gocNRZ6lMZ6Dp/du9S3huw2H0bsFKF7YOXqaTxkX3ZkH/UPLuguRJuAV6ln\ndnZv4zJnnz9cpv751hG8/3V3OSxq/LHm80wPI2uXUH6EyDvFSUIH2NfSUWB0DXY0VHjsHWYjHIYO\nGQs90FI9O2f3xvmA22wOb+6yMwA2n6zUD2u0iD6EWNmkRnM3Zp+wdDA9yvXD6Im2782ze8uX5Mqz\n0faWaGWxzXSW8oYt1tW4zaIH2FE2HwfrN6ofzcQL1UO6Rn2X27BDJtqHK5cfRuH9wA3jCY66Pj06\njNAuwRBiBI3b1qmdvLPfqNHPeurXseoSjb/VSw149cs1uZOMpigH8ZALsc8KsKNWan3G3dlu3FVm\n2+2y2I6pUmqQgdDpTIGOj2i3OOT69IErZf+J4DBycFRJJ3lQgJ19dn/NQerG2fOzzdXYQVzV/br3\n/KDFqNREcTsI2KotAnYsRhcU2p00uLLzwn+u050HBdjgkG+bj2avZjyEtdvQSjtn90beZ0Z+7uFB\nto2YI12rJNoqGo3WyLFrZ3R21V/rvZ4TYB2cMkdwzREw/p7WeFCIEBiJzqbbRWjY19FdqMM/0YbV\nvXScDLvZPCXASt3sxtc7temWduAIh6ptImxmeVtvSahPFH8kZs+F2AgbcN1TAiz+qQSvrRfPblGJ\n25g1j7b8mEU2cg5akXrNrvWUAFuSfWX3tyeHDYCpaamfvzeFuaukWVnas1G05g1o84ngnu/A3OIp\nX2R+dXfjy9G1jdfyx9n/qadNx6ywrfQpvoB/zWYf5CsEoVzZJsM+u/a9ch0nHxRgryN2451TVI8y\ne6Cc/SbK7JcH9meYSGPJaMM4MDWjbXXtxZzXJme4/nRns2cF2LEi7EvlpjZ1bOREO3wMIpxVzDbO\n0hc2gh8UIpvd4Pd/iXitmDvCAz39GtgGcY4+7VdcskyK3WB6LKt8ulO/Blu5ejd7Ap5XpfXKUbKy\nQQ78yEvDD0ctf7Mb2+SxBNg6Sb9GFnBmwSFWfajZ0dRDy2l938ve+nDlhvRzozp7q7t9NxwF0k+9\n7oY3MoS42iGb4CGXSddWEvm6+s4B91VntZ/1P5R6iPp4b3ztV16vcVdLTudPcRc9sHXOSK8Ni228\nBUDAoKrbfIra/sINfej3xIYKh/HedCtlqtLtSDcPe4MN8/rK/0Y+j0xHD+xOe6ZFTKcFNr5FqN3m\n2OG1xhf+nBJ2xml12edbepc9q+bKtTy02MP7H5H3rIfQA9tu57jf/okVw2vrlZRPC7uPbaut8SWj\nKzdLz6k/IfLl97su780arcqfa3Zt17Zx+OFUGzphP9vEFbINBNg6jZnRyCY7OHB4bWnEprGGpWXu\nGTys//XAXLzyOF6vtv3con2ZO9/rWJtPB+vb56rNrBL8DzmwCLDbPGQLW2VbN7TSxWwZZf252P2W\nijzqNKhc2nlH88O7Pk8+8n7t+bBL3e5RpPXdngLsaj+vfnOGysF9dvLLqEd43lHg2F742X36UVsd\ntdF+/jpkaUmtaoGlM7OnnRMIsBt0MF1wp2Ho48rrGe2DLS3Xw44qaf9CLr4y97P7uOqCVvmSA4rL\nbFu36eEnxALsHoePI2Xxc3zj2NY44+B+8RhdLnumk8iw0Whw40tmd6KHNKYAI5CYR/Azzjam40V7\nlhywT99SwO1F3u68wYYzFhtQJwG2diYuO21r5/IEc8M1j6UBwG3T4SqPV2aFHOvAXIzQp1/qQKza\nWqSa7G/XzxeZP8nv05PUtrnv75U3jpqdbD0dJ1labLltVEZXDj8obBgR6kbjnjhaNefWtGDpkuc1\nq6/cINdOlB0VvHQM7Hjz6yHAymPod6VGXmH9RezZDb7UYse25C0foVdHfd7ZA3T7S9Y+v/KcmIeU\n2bOxUYzFrPwoQVfMKqPNa9Xs581XO+sn8tMdr6WS17+9hPbaVhUzfcf2U7nKYodHthVzoNFJTGMx\noxPhn43Wssy1NVfe8TW3pjaPOlReuPYDlut91Qtn/3TGFrK50W5s7aUjQL21yy1/w2paemHTh7zJ\n0wMMgCXBA6KTSRxHqayt4Cvyeqdeal6aPVF5fE8xVu7hzp7/MrzLSYu6eJPY83Z7doqfS965hLP1\n3wM74wPunIMw8nPsZWV1rbZd7F0afDjvOvz+8brZqkYF76l/+trZhe9cfuMSZj/IeVvXgeu9ZSjy\nq75C60Pfs6+afeFobPNn/Qda2n5Klc8y/euG1RT2yl9JD4z/bL5+c3glx77pXSPJ779uefep+Mej\nRpUm3XzYnS4zzoqbFby8a/QQYOWJUuSzhlCFbd76P4VjS/r5vmUBEUpaK0J58Vvpq15k/a+rhham\nf2p55EoPmU+4TQ/T6L+Sno8sbZRDKl8Tye3DcTfa0A53fa7ROFiQQ0/8VVxau+5anj+7N0VLrIqW\n8dVH6STAwm5wI5+5ibnT51xYUfSmq/St693uCLt68LYNXt5rIZN+dqcqc9/jf+Sy+NFVPbk11UmA\n9eTi72Hs3yvO6yOOaqt8QaISb88xe+z7/nfpyB7ftm736++HzbUZrO1xnldJFgLsXNNjShA7x9am\n+bG7oh9vMTy4VPn3T2vP2bvR+M3cCyoJK+nHX/pG8131hCLAVmufMBIqsaZaxjOz6/JDtYs8p+ls\nZWc9VyM89iRsgx5mIV7p7J7HxRPqtr3dlSeDo0Gwxj+98ky3+6llG0t0jed6Gdukm633bHpgK1RG\nripPTjHBL6yle7st/amzpq586tJn8pXbztrhUTasu8eubgF2hV6PKdecJNa/tbr03/PuDHIZBzJe\nzfdeeSZDiKs1DruVQ/Chbsew3zUfrfIdg/ZYahx/K69rNheYQKL87nUVrDI9qmxegw9pRgG22oF7\n2kM2ss3K/Xm0by/9ae0On+IGQi1XRA489l1v2pMOuBYus+3KdKLVfSBDiFs07l29jhy+5r4fet7H\nrOyZZV/wwHdMN29tkLRsDtHxAWeJHtgKlUGt+qs29w+Cu3LC5M9Kyn+vLSnCR5gaDdXeW8x5lia1\ndvyRz5Dxi9v76YGts/NIV5lTN3ra9CW02DzFP+bZ67Rra2NgSfdf65zSA7vaz7PpFFdlAtqzx0Y+\ne43TzT3b9DLnjcUk8uSG0gO7yOxVolWz6Z68mTYqe1F7ZnNo6it93Kx2t3KLfVQz6oFdp+Ww6JYK\n+63qr0yfqeUjsBZooQd2qVUXqGNelenP/rHH4d/nHXZHX97o7Phe7hedfTROpQcWyyjhpFdkl12t\nfM5lUenFKgLsaqP7PjTeJ8mOncXh0dJrVsF+Aux+P+dJS6+A7rpaGeprUu7/xL1cA7vUNKtm93yJ\nlcs1VytDbRWNt8nndp3d4XpEDywEJ7DpXH+1cnTr5ICHIZtxNN2vEQF2D9/ZvFH7Zcj2BQ7/PnxV\nBtw2IocoP/WUagLsHqf+EAkVs22+bUVcdrUy8mXRh9wlJKMnnGe4BnapU+/azvUuOzSEOgaV984I\nVRilzwN+p1uAXW36QyS3lvNQw77tlGInX/aIr+MtXIDdw65+r8afBWCJ4M+rp4NP9ABrPL/zkxOs\n5fi7k70svulZWmdrLXqAvRpGbw3E027nzwJAOh1v3qFnIZbJZMiCo4z25453b+hbgh5Yi477yPFl\n/Kp//AqBnzoJsFcx0jgaUXR57FRLX6vSzpBOulGuHgJsdEsLiQWwQbpfdr4/wNzNtgOuUALXuz/A\n9meVAasIfK0KuNj9AVbxPa8fjowtKSXMbiG3gOuFDrCvpYPjkFXTW35xGV+rAu4SPcAqh0I/RxLE\nKMOsC+Aa0QOMFIQWcL3Qd+IAgCUCDICUBBgAKQkwAFISYCEMs/h8oQqgkQC73yi0ZBhACwF2s6W7\nuV9fCZDd08ZyBBhADx44liPAQih/ePreSoCMnjmW404cIZT3LL67FoAc9MCi+P6W9N1VAIk9bSxH\ngN1sdjt7yMYHHOs7lvOcH5g1hHi/cvBQdAF7PCS6vvTAovieN91dBZDSM8dy9MAAevDAsRwBBtCP\nh0TXlyFEAFISYACkJMAASEmAAZCSAAMgJQEGQEpdBdijvoIO8HBdBVheWaJXncdS57HU+TSdfJHZ\nBgHwNJ30wNxIEOBpOgkwAJ7m3VPHZfhd4/KRu4oBSC1+OnRyDWxJ/BUAwDZpAmy2LyWfAB4rTYDJ\nKgBKJnEAkJIAAyClrmYhAvAcaa6BcaXpFxJiil9nOfkocqnp6oxc5CtPndn1HGCR98nRpMqhvCw1\nl7VFqzlO237f8fte7/d7FLdx2rO9zlek9qx86TNUe7Y05ivefpRCzwH2Krahuwv5R73fcG/N9fdd\nqu36mpfeK1rbDsV8Pp/2dotTZ/D2HNUWtj2XQnf2VdGOV5F1O4mj3FZmjx0BRah56a6SS7XdVfPa\nu19GaNuWeqLVuSRandrzmTrvgYVluOA8cdo2y5qt1xmzPSMf4it1xmnMPnTbA4tv6ENE3hWTita2\n0wtLMS3VGao9v0W+wgfAUp2hGjM7PbAblBu04YJjBWzbFEfb10KdAduznEUSuVWndQZszOz0wOBE\nqdMrviwZkKXOdATYDWzN5wnVtllSoVJnwPaMrzJF9uJKutdtgC3NRIqmZfZUBC0zu6LVHKHOpTeK\n1p4tbxShzpFVMw8j1DkYFROkznQ6vwYW85Sn3Fhnt+zLK2oV/9QyWttWvrg6/evPx8+zVGeo9hxd\nN2qvJ0id9cZ8RdqPshD1AKTU7RAiAH0TYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFIS\nYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICU\nBBgAKQkwAFISYACkJMAASEmAAZCSAAMgpf8DeAOJw+3UssgAAAAASUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Evaluation of the CS Recovery Probability\n", "-----------------------------------------\n", "In order to bench in a randomized manner the efficiency of compressed\n", "sensing, we compute the probability $p_s$ for a $s$-sparse signal\n", "with random non-zero coefficient locations to be recovered by $\\ell^1$\n", "minimization.\n", "\n", "\n", "Put formally, if we call $ x^\\star(y) $ the resolution of the $\\ell^1$\n", "problem using measurements $y$, then we want to compute with\n", "Monte-Carlo sampling\n", "$$ p_s = \\mathbb{P}( x^\\star(Ax)=x \\:\\backslash\\: \\norm{x}_0=s ) $$\n", "\n", "\n", "An important feature of the DR algorithm is that it can be run on many\n", "signal at once.\n", "\n", "\n", "Number of signals."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["q = 1000;"]}, {"source": ["List of benched sparsity."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["slist = 14:2:42;"]}, {"source": ["List of sparsity of each signal"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Slist = slist(mod(0:q-1,length(slist))+1);"]}, {"source": ["per\n", "Genetate signals so that |x0(:,j)| has sparsity |Slist(i)|."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["U = rand(n,q);\n", "v = sort(U);\n", "v = v( (0:q-1)*n + Slist );\n", "x0 = U<=repmat( v, [n 1] );"]}, {"source": ["Measurements."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "cell_type": "code", "language": "python", "metadata": {}, "input": ["y = A*x0;"]}, {"source": ["__Exercise 3__\n", "\n", "Perform DR on the set of signals |x0|. Note that the proximal mappings\n", "operate in parallel on all the signals in |x0|.\n", "Each |i|, count the average number |proba(i)|\n", "of recovered vectors of sparsity |slist(i)| (up to a given, small, precision)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo3()"]}, {"collapsed": false, "outputs": [], "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}], "metadata": {}}], "nbformat": 3, "metadata": {"kernelspec": {"name": "matlab_kernel", "language": "matlab", "display_name": "Matlab"}, "language_info": {"mimetype": "text/x-matlab", "name": "matlab", "file_extension": ".m", "help_links": [{"url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md", "text": "MetaKernel Magics"}]}}}