{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Performance of Sparse Recovery Using L1 Minimization\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 tour explores theoritical garantees for the performance of recovery\n", "using $\\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/sparsity_6_l1_recovery')"]}, {"source": ["Sparse $\\ell^1$ Recovery\n", "--------------------------\n", "We consider the inverse problem of estimating an unknown signal $x_0 \\in\n", "\\RR^N$ from noisy measurements $y=\\Phi x_0 + w \\in \\RR^P$ where $\\Phi \\in \\RR^{P \\times N}$\n", "is a measurement matrix with $P \\leq N$, and $w$ is some noise.\n", "\n", "\n", "This tour is focused on recovery using $\\ell^1$ minimization\n", "$$ x^\\star \\in \\uargmin{x \\in \\RR^N} \\frac{1}{2}\\norm{y-\\Phi x}^2 + \\la \\normu{x}. $$\n", "\n", "\n", "Where there is no noise, we consider the problem $ \\Pp(y) $\n", "$$ x^\\star \\in \\uargmin{\\Phi x = y} \\normu{x}. $$\n", "\n", "\n", "We are not concerned here about the actual way to solve this convex\n", "problem (see the other numerical tours on sparse regularization) but\n", "rather on the theoritical analysis of wether $x^\\star$ is close to\n", "$x_0$.\n", "\n", "\n", "More precisely, we consider the following three key properties\n", "\n", "\n", "* *Noiseless identifiability*: $x_0$ is the unique solution of $\n", "\\Pp(y) $ for $y=\\Phi x_0$.\n", "* *Robustess to small noise:*: one has $\\norm{x^\\star - x_0} =\n", "O(\\norm{w})$ for $y=\\Phi x_0+w$ if $\\norm{w}$ is smaller than\n", "an arbitrary small constant that depends on $x_0$ if $\\la$ is well chosen according to $\\norm{w}$.\n", "* *Robustess to bounded noise:* same as above, but $\\norm{w}$ can be\n", "arbitrary.\n", "\n", "\n", "\n", "\n", "Note that noise robustness implies identifiability, but the converse\n", "is not true in general.\n", "\n", "\n", "Coherence Criteria\n", "------------------\n", "The simplest criteria for identifiality are based on the coherence of the\n", "matrix $\\Phi$ and depends only on the sparsity $\\norm{x_0}_0$ of the\n", "original signal. This criteria is thus not very precise and gives very pessimistic\n", "bounds.\n", "\n", "\n", "The coherence of the matrix $\\Phi = ( \\phi_i )_{i=1}^N \\in \\RR^{P \\times\n", "N}$ with unit norm colum $\\norm{\\phi_i}=1$ is\n", "$$ \\mu(\\Phi) = \\umax{i \\neq j} \\abs{\\dotp{\\phi_i}{\\phi_j}}. $$\n", "\n", "\n", "\n", "Compute the correlation matrix (remove the diagonal of 1's)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["remove_diag = @(C)C-diag(diag(C));\n", "Correlation = @(Phi)remove_diag(abs(Phi'*Phi));"]}, {"source": ["Compute the coherence $\\mu(\\Phi)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["maxall = @(C)max(C(:));\n", "mu = @(Phi)maxall(Correlation(Phi));"]}, {"source": ["The condition\n", "$$ \\normz{x_0} < \\frac{1}{2}\\pa{1 + \\frac{1}{\\mu(\\Phi)}} $$\n", "implies that $x_0$ is identifiable, and also implies to robustess to small and bounded noise.\n", "\n", "\n", "Equivalently, this condition can be written as $\\text{Coh}(\\normz{x_0})<1$\n", "where\n", "$$ \\text{Coh}(k) = \\frac{k \\mu(\\Phi)}{ 1 - (k-1)\\mu(\\Phi) } $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Coh = @(Phi,k)(k * mu(Phi)) / ( 1 - (k-1) * mu(Phi) );"]}, {"source": ["Generate a matrix with random unit columns in $\\RR^P$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["normalize = @(Phi) Phi ./ repmat(sqrt(sum(Phi.^2)), [size(Phi,1) 1]);\n", "PhiRand = @(P,N)normalize(randn(P,N));\n", "Phi = PhiRand(250,1000);"]}, {"source": ["Compute the coherence and the maximum possible sparsity to ensure\n", "recovery using the coherence bound."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["Coherence:    0.30\n", "Sparsity max: 2\n"], "output_type": "display_data"}], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["c = mu(Phi);\n", "fprintf('Coherence:    %.2f\\n', c);\n", "fprintf('Sparsity max: %d\\n', floor(1/2*(1+1/c)) );"]}, {"source": ["__Exercise 1__\n", "\n", "Display how the average coherence of a random matrix\n", "decays with the redundancy $\\eta = P/N$ of\n", "the matrix $\\Phi$. Can you derive an empirical law between\n", "$P$ and the maximal sparsity?\n", "= plot(eta_list, k_mean ); set(h, 'LineWidth', 2);\n", "abel('\\eta'); ylabel('1/2(1+1/E(\\mu))');"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAT1ElEQVR4nO3d65Kb\nuBaAUTw17//KnB8kHIabMRdpb7FWTaUSd6ejsbv9WUKGT9/3HQBk80/tAQDAGQIGQEoCBkBKAgZA\nSgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBK\nAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASv/WHsCz\nPp///7biMADS6fu+9hC+aDxgo1CPxOfzCTWeQcBRBRxSZ1SHBRxSF3JUAYfUdd3nk+BFf/NLiH8e\ngwyPBQA/aD5g3fjKRsMAWtJ+wKY0DKAZEddebzQuLk/T1fT/McANYh6Zm3nLDCz8AwHAb94SsM7B\nMIC2vChgnYYBNORdAZvSMIDUXhew6cEwDQPI63UB62zoAGjCGwPWORgGkN9LA9ZpGEByCU7mO55T\ncvauuum5Jvu+n/3xx3/CuiJAMtEDNn03+PKd4ft//Krv/z/90jCAXHIvIX4+n+nEa/bHI0QLIKno\nM7B9f89z+GdmNvvjYJa05SxtnIeZhAFvluIaYFOJA3Zw/fDIuqKGAey89I8p6xLi7M698b7O8KgB\nEH4GNt1eOF0hXL19+sff/yEbOgAyiR6wbnup8OIWxLV/yPQLII2sS4gP8e5mgCwEbE7DAFIQsD0a\nBhCWgK1wyRWA+ARsnV2IAMEJ2CYHwwAiE7A9GgYQloAdpWEAoQjYFzZ0AMQkYN/Z0AEQkIAd4mAY\nQDQJzoU4s3PS3uk56a+fGnHGJVcAQkkWsOnFKmcXrhzc3q2NYWgYQGWtLSF+Pp/nrsNmQwdAHMlm\nYF9Nrxk23DLr2cUpmkuuAK1KcRXmqaYCthonB8MAjpg+W6aIWTtLiCXvbpsSAapLNgPr+365C3FY\nMFz9UAHmYQBVJAtYtxan8ZaS3ZoeDNMwgPLaWUIsT7QAKhKwSxwMA6hFwK7SMIAqBOxOGgZQjIDd\nwBk6AMoTsHvY0AFQmIDdxsEwgJIE7E4aBlCMgD1FwwAeJWA3s6EDoAwBu58NHQAFJDgX4v4pesdL\nf03PRl/ypIirXHIF4GnRAza7NOWsTPderPJeGgbwqMRLiKs9i3kRtpCDAsgt+gzsJ+Na4s51RV1y\nBWBVzAnAjqwBG+7o8dfhgparn1l3XXHaMIDIdl76x5R1CbH/q/vvZZoD8u5mgCdEn4FN47S6Qrj1\naaHY0AFwu+gB6xZB2vpjwG6t0jCAW2RdQkzHGToA7iVg5WgYwI0ErKhZw2QM4DQBK212AEzGAM4R\nsAr6XsYArhKwamQM4AoBq0zGAM4RsBBkDOBXAhaIjAEcl+BMHG8zNGzarfH3TuEBMDIDC2o5G+tM\nyAAmQgTsp3PJf/665asFN2TMuiLAUrIlxOmp6FdPS99Svaa21hUtKgKvVS1g09LcdSL5IWmtNqyT\nMYCJagGbXfrzuYuhzHqW5aorO2QMeEK6V//JlhB3DHf9+Gu664T9SsaAe83mFRVHclCIJcTxj1di\ns39srFUyBrxWiCXEn/7WWL53FmuVjAEvlG8Jcdmq2S2vjZl3QAOvEuJ9YNzIO6CBl8g3A+OIsWHW\nFYFWmYE1zok8gFYJ2CvIGNAeAXsRGQNaImCvI2NAG2zieCl77oHszMBezZ57IC8zMFZmY50JGRCe\ngPHH6lvHOu8eA6KyhMicdUUghQQzsOXZe1dvf+IKmW+2cy6PzoQMCCB6wKZnml+edX5M1/Ab3XrC\nzhEy9zdQUfSA7Vjm6vpFxdhiQgZEkzhg3eKqmLMJ2c7ncJoJGbQqxVWYp3IHbHr0a6tMivUEEzJo\nz85L/5iy7kKc3bkp7usm2bII1BJ9Btb3/XK3Yd/3q7dP/0hJJmRAefN9fY1ZblykgNXpl8cBEknx\n5Bl9BkZGJmRAAVmPgZGCI2TAc8zAeJyTBQNPEDAKcbJg4F6WECnNuiJwCzMw6rDRA7jIDIzKTMiA\nc8zACGF/QtaZkwELAkYsq1sWl7foGSBgRLS1ZXGkZ4CAEdqsTHoGjPIFbOekvc7n27wTPfO9AK1K\nFrDp+SVXzzW5ek1LWjV9kE3O4G2SBWyfaL2ZxUZ4m6YC1q1d2XJ2i8i9hJ7Br9JdGbi1gM0ubtkp\nFl3XOXgGB0yfLVPErJ2AOe7FcUd6ZnIGwSULWN/3y62GQ7pWPwRH2AwCGSULWLcWp/EW3eK6cwfP\nln8ReFq+gEFJB3u2+iFJg0cJGPxg9cT5WyQNHiVgcImkQS0CBje7mLTVrwAsCRg87qekrX5U0mBJ\nwKACSYPrBAxCkDT4lYBBUBeTpmc0T8AgjSvbQ/SM9ggYJHY8aXpGewQMmuI6MrxHgoBtnaJ3drtL\nqMCSntGw6AGbXiRlecGUMV3Db3QL9ukZLYkesB3LXA2TMBmDg/SM1BIHbDCdls0mZOMnTD9f3mCL\nnr1ciqswTyUO2Gy+tVUmxYJzXBrtbXZe+seUOGDd4qLMdQcDbXNpNKKJHrC+75e7Dccbpx/a2qwI\nPOF4z3Y+6oeVKxqfuJiZQRUX15/81FaX4skz+gwMyGjrqe9g2MzYOELAgHKEjRsJGFCfsHGCgAFx\nPRS2na9MIgIG5HMxbMc/U+ciEzCgHdfDdvovSl15Aga070hdLm79N6UrT8AqiPkGi4CjCjikzqgO\nCzikbndUxwd7JXWmdDcSMIDfHEyLKd3TBAzgEYc79+m68xV685ROwAAqM6U7J+Ii9Y1SXBEA4G43\nPLHHj0PjM7C28wxwRfZX+I0HDIAt2V/h/1N7AABwhoABkJKAAZBSa8fAxm2Hq9s3ap0aYGtU+6Ot\nMqT9Dz0t4MN35I5a/eijAn5Hbf3rs53AQe6o6UfjfFPVffiWYp5LZaqpgE3v7tldX3E//c6our/f\nqYW/UfaHFHNUVR7Bg3dUYVujmj4dl3/q2RpVxR/Dr99Rob7Pv36zlZTlDUhvWULs+z7gS4mAQ+pC\njqr6z/OWz+cT7Ud9GFLAuyvmqFgV8wlzqakZWFIBf7CjPSmHVWu6syXOS/gUai0h7vPTd5yA1RTz\n56ebPC/XHkjX/R3G+GucuyvOSOIL9cB1wdbrBn3fj08IQX704nvLEmJYEX5ypgL+5PR/dZHuroB3\nFKkNHY3zHZ5CiJceN5pt41lu5YiwO2sYRrTdWdNj2lWGtPynIzx8Ke6obvFNFef7vKs6xUnx8EUY\n0qogc9Md0ccHAKssIQKQkoABkJKAAZCSgAGQkoABkJI3MkNpdd9BAc0wA4PSxjdlSxdcIWAApCRg\nAKQkYACkJGBQgaNfcJ2AAZCSk/kCkJIZGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkY\nACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgA\nKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAD/8fl0n0/tQRzwb+0BPOuT4kEACKGvPYDf\nNB6wruv6vv5D8vl8DCPOGAwj2hgMo+4Y8r7Obz9gACxtdWuo5+fziT8hEzCAF9nvVi4JAjYex1qd\nVo/T7enhruqLAAChtNStUf0F333T5eDl0vAQrTFgy/+XCCvaALWc7laKJ88EM7Atw/07nXhNewbw\nWjv7Mlp6gkwcsKXVqdhsJ728Aa262K107zvKGrDhjh5/7ft+q0yKBbTtrvnWzkv/mLIGbHZgLMVy\nLcBdXrJIuC96wKZHuXY2ayw/DaA9ujUVPWDdIkhbf9QtoFW6tSpBwABeSLS+EjCAQHTrOAEDCKHJ\nk2U8SsAAalrtlmgdIWAAdUjXRQIGUJp03ULAAArRrXsJGMDjpOsJAgbwFN16lIAB3E+6ChAwgNvo\nVkkCBnAD6SpPwAAuka5aBAzgpGW6dKukfAH7et0vF7cEniZdESQL2DROq6FKcRlsIC/piiNZwPYN\nSdMw4AnSFU1TAVs165nVReBXL0lXulf/7QRsuOvHX8dQKRZw2kvSNZg+W6aIWTsB2z82BvCTV6Ur\nqWQBmx7iUizgCdKVRbKAdWtLgrNbxAw4R7pyyRcwgNtJV0YCBryadOUlYMBLSVd2RQPmLVlAELN6\neTbKqFDAhnTNirV6I8CjpKsZhQK2WinpAoqxYNgeS4hAy7ZOKOHppwGWEIEG6dYbWEIE2qFbr1Jt\nG73zPwG32DnrrOeYtnkfGJCSblFtE4fpF3CCbjEqGrDZxWY0DDhIt1hKsIS4vH7K6u2md9Ae3WJH\nzfeBHdlGP52oLSdtY7qG3+gWtMFmQo6o+T6wi5ZfzRvLIDXd4idF3we2tRh42uqpPWYTNaf/gOB0\nK4jPzoptSHU2cdw1VZp2ceurKRYE5OBWQDsv/WOqs4njelSW0yyhgvh0ixsVPQY2cyQ5fd8vdxv2\nfb96+/EvC5SkWzyh3MSlyiTJzAxqEa3UUjx5JngfGJCIblFMuYDFjzlwmm5RnhkYcJ5uUVGFgKVY\nWgW2iBZBmIEBh+gW0dQ5F6JzPkEWukVYFc7EYQkR4tMt4quwhKheEJNokYtjYPB2ukVSAgbv5TTw\npBb9XIjA7XSLNpS7Hpi9G1CXbtEYS4jQvtV06RbZORciNMuUi7YVPQY2a9i5tzPvXPfLJcFgYMrF\nG5Q7BtYttnKcyMz0QNrqQTXvlebNdItXqXAmjrxfH8KSLl6oXMBmC4YPTZKW+/WvT/sgLN3iRqvv\nd4qs3DGwMot7y7VKxaJJ0sXtps+WKWJWehv9c28Ic9yLN9AtGBVdQpyuH577IkP/xt9Pv/Lqh6AN\nNsTDUtFdiNPfX2nY1i26RXtMuWBLtfeBiQ3sMOWCr+q8D0y9YIspFxxU531g587BAQ3TLfhVnZP5\nSheMpAvOcTZ6qEO34CIXtITSpAtu4YKWUIhuwb0sIcKzbIiHh7igJTzFlAsele+ClhDZzhlmfKfD\nvZJd0BIC2j8tmm9zeEhTF7SEYkQLqqu5icNVu0hHtyCOagGb7aq3yZ6wRAtiSrCNfusUwLPbzee4\nkWhBfNUCNr3+ZLednNllMGefNqZr+I1ucYVoQS41Z2AXe7P86/blc4JuQVIJlhD3TadlswnZ+AnT\nz5c3OtGCNasnrY2s5iaO2S2/pmU239r664rFQLRg385L/5gqHwO7axXRJkZWiRY0LPoxsOlej+kK\n4XDj9ENbmxV5Id2CNygasHPLhqs7D5d/UbdeTrTgbUrPwGaLfpb+uEK04M3S70LkbUQLGAgYCXzd\nD6Vb8EKlz0Y/22ph/ZAdJlvAjmrHwGCVaAEHlQ6Yze4siRZwQult9N56zMBhLeAimzgoymQLuIuA\nUcJOt0QLOKfyLkQaJlrAowoFbHkSKcfAWqVbQBmFAqZVzdMtoDDHwDhPtICKBIyf6RYQQb6A7WwD\nsUPkUboFhJIsYF/fCj296GXpwbVItICwkgVsn2jdRbeA+JoK2GA2/Zrt4Be5HboFb7Z8v1NwTQVs\nuPdniVKsfaIFDHZe+sfUVMA6uTpMt4DskgVs9WRUw5rhcLuNiPt0C2hGsoB1a2Vyced9ogU0KV/A\nOEi3gLYJWGt0C3gJAWuBaAEvJGCJ6RbwZgKWj24BdAKWhWgBzAhYaLoFsEXAItItgK8ELArRAviJ\ngFWmWwDnCFgdugVwkYCVI1oANxKwx+kWwBME7Cm6BfCoBAHbv8TXcDGw7r/XD611aRXRAigmesDG\nPs1+361d8Vq3AN4jesB2LCdew++LZUy3ACpKHLClMWk7E7WLeduJVqdbQGbLZa3g2gnYVplumZCZ\nbAHN2z9GE1AjAZvNuu77spsf0i2AuqIHrO/72S7E1VYtP+00i4QAKUQPWLcI0tYfn+uWaAEE9E/t\nAVT2+fz5b6nv//x3x78SYjU5wjAijKEzjGBj6Awj2BiySDADu51FQoAGvChgFgkBWtJ+wHQLoEmP\n7D6PY6NelpgBvohfh8YD1v1tWOv/lwCv0/4SonQBNOnt2+gBSErAAEhJwABIqf1jYA+d5/frPzr8\nZvZP33XCxitjmH601jDuvcbN6WHs3F5yDN1/75CnR3LwEudPy/KIRPhRjfMMFk3LAat1Rpatq0jv\nXF262Bhmt9QaRuGrNuw8CuNg6j4iXalniq/fGAXG8HUYoR6RisOI9gwWUMtLiH3fR77r6WL8eHw+\nnzjDqDuA6ndCF+Ylf4RheAb7quUZWFgRTtZZbAkxuDivNMu83k8hwg/IwMMRXCMBK39M5cQ/PV60\nbHr1ssJj6J5cH/j1UXjsMqTVvhlOj+GhQR4fxvCZ46/3jufcvVHxZ6R77BVehG/OljQSsIrfB8f/\n6edezQX5MTCMc2OI8I3x6Ew06c/IEyOJ8M3ZkvYnyBH28Ew3TUxvLzmGIMMotn/k+DCmt5ccQ5Bh\nLDcOVNmFGOGuGIdReJ4U5BHZGVJY7QcMgCa1vAsRgIYJGNzmln0Hq18kzsY8iEPAIIHbd65CAxrZ\nhQhxrG4NGD965E0Oq5sIhtsdtIaRgMGdlm+zW5686qvp3zLxgi2WEAFIScAASMkSItxpegRreuxq\n9mlbJ8UHjhMwuM1YrNXbu/+ecnC8ZfX8C8svYgcHzFhChMd9/toq3Mh2eTjODAwe99PMafXqjr9+\nEXgDMzAAUjIDg5qmV4mrPRZIxmFhAFKyhAhASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoC\nBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEoCBkBKAgZASgIGQEr/A9Vvsknlctn6\nAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Support and Sign-based Criteria\n", "-------------------------------\n", "In the following we will consider the support\n", "$$ \\text{supp}(x_0) = \\enscond{i}{x_0(i) \\neq 0} $$\n", "of the vector $x_0$. The co-support is its complementary $I^c$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["supp   = @(s)find(abs(s)>1e-5);\n", "cosupp = @(s)find(abs(s)<1e-5);"]}, {"source": ["Given some support $ I \\subset \\{0,\\ldots,N-1\\} $, we will denote as\n", "$ \\Phi = (\\phi_m)_{m \\in I} \\in \\RR^{N \\times \\abs{I}}$ the\n", "sub-matrix extracted from $\\Phi$ using the columns indexed by $I$.\n", "\n", "\n", "J.J. Fuchs introduces a criteria $F$ for identifiability that depends on the\n", "sign of $x_0$.\n", "\n", "\n", "J.J. Fuchs. _Recovery of exact sparse representations in the presence of\n", "bounded noise._ IEEE Trans. Inform. Theory, 51(10), p. 3601-3608, 2005\n", "\n", "\n", "Under the condition that $\\Phi_I$ has full rank, the $F$ measure\n", "of a sign vector $s \\in \\{+1,0,-1\\}^N$ with $\\text{supp}(s)=I$ reads\n", "$$ \\text{F}(s) = \\norm{ \\Psi_I s_I }_\\infty\n", "      \\qwhereq \\Psi_I = \\Phi_{I^c}^* \\Phi_I^{+,*} $$\n", "where $ A^+ = (A^* A)^{-1} A^* $ is the pseudo inverse of a\n", "matrix $A$.\n", "\n", "\n", "The condition\n", "$$ \\text{F}(\\text{sign}(x_0))<1 $$\n", "implies that $x_0$ is identifiable, and also implies to robustess to\n", "small noise. It does not however imply robustess to a bounded noise.\n", "\n", "\n", "Compute $\\Psi_I$ matrix."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["PsiI = @(Phi,I)Phi(:, setdiff(1:size(Phi,2),I) )' * pinv(Phi(:,I))';"]}, {"source": ["Compute $\\text{F}(s)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["F = @(Phi,s)norm(PsiI(Phi,supp(s))*s(supp(s)), 'inf');"]}, {"source": ["The Exact Recovery Criterion (ERC) of a support $I$,\n", "introduced by Tropp in\n", "\n", "\n", "J. A. Tropp. _Just relax: Convex programming methods for identifying\n", "sparse signals._ IEEE Trans. Inform. Theory, vol. 52, num. 3, pp. 1030-1051, Mar. 2006.\n", "\n", "\n", "Under the condition that $\\Phi_I$ has full rank, this condition reads\n", "$$ \\text{ERC}(I) = \\norm{\\Psi_{I}}_{\\infty,\\infty}$\n", "      =  \\umax{j \\in I^c} \\norm{ \\Phi_I^+ \\phi_j }_1. $$\n", "where $\\norm{A}_{\\infty,\\infty}$ is the $\\ell^\\infty-\\ell^\\infty$\n", "operator norm of a matrix $A$,\n", "computed with the Matlab command |norm(A,'inf')|."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["erc = @(Phi,I)norm(PsiI(Phi,I), 'inf');"]}, {"source": ["The condition\n", "$$ \\text{ERC}(\\text{supp}(x_0))<1 $$\n", "implies that $x_0$ is identifiable, and also implies to robustess to\n", "small and bounded noise.\n", "\n", "\n", "One can prove that the ERC is the maximum of the F criterion for all signs of the given\n", "support\n", "$$ \\text{ERC}(I) = \\umax{ s, \\text{supp}(s) \\subset I } \\text{F}(s). $$\n", "\n", "\n", "The weak-ERC is an approximation of the ERC using only the correlation\n", "matrix\n", "$$ \\text{w-ERC}(I) = \\frac{\n", "      \\umax{j \\in I^c} \\sum_{i \\in I} \\abs{\\dotp{\\phi_i}{\\phi_j}}$\n", " }{\n", "      1-\\umax{j \\in I} \\sum_{i \\neq j \\in I} \\abs{\\dotp{\\phi_i}{\\phi_j}}$\n", " }$$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["g = @(C,I)sum(C(:,I),2);\n", "werc_g = @(g,I,J)max(g(J)) / (1-max(g(I)));\n", "werc = @(Phi,I)werc_g( g(Correlation(Phi),I), I, setdiff(1:size(Phi,2),I) );"]}, {"source": ["One has, if $\\text{w-ERC}(I)>0$,  for $ I = \\text{supp}(s) $,\n", "$$ \\text{F}(s) \\leq \\text{ERC}(I) \\leq \\text{w-ERC}(I) \\leq\n", "      \\text{Coh}(\\abs{I}). $$\n", "\n", "\n", "This shows in particular that the condition\n", "$$ \\text{w-ERC}(\\text{supp}(x_0))<1 $$\n", "implies identifiability and robustess to small and bounded noise."], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 2__\n", "\n", "Show that this inequality holds on a given matrix.\n", "What can you conclude about the sharpness of these criteria ?"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "text": ["N=2000, P=1990, |I|=6\n", "F(s)    =0.19\n", "ERC(I)  =0.25\n", "w-ERC(s)=0.27\n", "Coh(|s|)=1.78\n"], "output_type": "display_data"}], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["__Exercise 3__\n", "\n", "For a given matrix $\\Phi$ generated using |PhiRand|, draw as a function of the sparsity $k$\n", "the probability that a random sign vector $s$ of sparsity\n", "$\\norm{s}_0=k$ satisfies the conditions $\\text{F}(x_0)<1$,\n", "$\\text{ERC}(x_0)<1$ and $\\text{w-ERC}(x_0)<1$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAASyUlEQVR4nO3d65Kb\nuBYGUHSq3/vMPDnzg5jQgDHmIrSltSo15didjoY2fGhrW0l933cAEM3/nh4AABwhwAAISYABEJIA\nAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIT0\n8/QAoGgppfFfLZ8+Xn7Z8GD6xdNnll/w7jtsfM2eb/vxL4JqCDD4YCO3xi/oJonS9/0s9rrfAfPu\nu+35mu53dM3+yM64hTooIcIHQzB8/LKUUobMkEkwMgOD78zCbEiUS+Y978qD0yeXz0OzBBh8Np2E\n3TcHmn3n7UqjGAMBBt95NzHK8BdZ1oIp5wNsua8LcfbddoaTLkQYCTB4hukUnOQUAiAkbfQAhCTA\nAAipwgDTXgzQgqra6EUXQDuqCrD1D3hekWrp/+e/R8H+fbKRRxcRcExVAbaUUrrk8tj/+/qGdSfZ\nE143GGbPUJbye9QrD7CuO32HP9sN4d/pb3Z97LTYN0H6Z/zE60MD+H1olyMp+eh1xQ+vK36EhnfG\n3cMLsSJTYRPHxfr+z6+llP7+Cqj/5+GTc3lcwx5L4AENzMCuMl5rVy+x45MF37KVaThg04OakqMI\nfFbhDOz2Wf84J6tlWjbWEh80O5yhjh/wDDOwc2YX3ZmhhcSEYrfFbKwv+eCVvEAyKHyEhQxvY7Gn\n8HWgy4dXyE9kPwF2ne0aI7v1/byi2CnNcqdwF+47FJ7WqyosIT5vo8BYmMf7ON7R3wF8JMAolxgD\nNggwuq6MPo53VmMMwBrYnWaLOZww6++wMEZ0szWnjXW4wj9S/SABRiSr/R0Z/lK4w8dYithYkZMS\nYuuK7eN4J0h/DIFNP8x57a/dAxi3eevNvTaYgRFSnpN6uIyU/Fk0QpsG1fQZobWTGRh/lNzH8RQf\n7WvTdLOda38t/qJ+nGOllIa1Lum1nwDLwvUvLBlGHkN0DTH29FjCEGCwlwsLdxtj7OmBxCDAiNfH\nkdn2hpdwzEapcPaSouI7Agw+k2FQIAHGX/o4NsgwKI0Ag71kGBRFgN1M8bouMgzK4YPM8B07XHKJ\n1b0Q322QuPzI88a3bafpwwyMrtOI+CUfDuMS/cTyyfGZ8QPO2x32DX6GzAwMjhjnYTaaqs993UwH\n7hRnM6rl7Gr8gtU5XN3MwPhFI+J+5mGclCaWT+784+0UDJfMwHJxo141P96aZKuor2bPct1rZv96\nWPXMwOA4TYncbbbuNVsPa6pguCTA+EMfxzEyjMNWS4hLY1Yt18Ma3zhRCRHOmjbWqyWy03b9cOPx\n9vdpqq5oBsacPo4DzMMgPwF2v5ZuiFomwyCzeCXEjQ4czTk8Sy0Rcgo2A9v4RPrOD6uzQR/HeUIL\nsok3A4PC2aSDC9kycUOwGdgey5nZzl5VRvo4TrJJBxfKtmViuEtlPTOw8Wc2++EVdJfhhhzoum4y\nAVo+OPkNB8e2TJzF4eHBZFNPgNUxI6YaComB3XftXnsrfLx2rd6R7w+Y8c5+7yDjCFZCXP1E+nTi\nZXfLk/RxQAbj9Wo6/Xq33rGzhDgzvR7WekmMNwNb/iTGZ2r9IQFZPXQlOXMFG/Jvelu/vMU/P8LS\nBJuBkY0+jvO0cnAJWya+E28GBlCBnaUjWyZuMAPLoor3CkBRBBjcyK0L3EeAMacR8Q41LkDAw6yB\nAa2rssGhBQKMt9I/yWzsvOkW9RSojnaGNikh5uVK1jA/fLiWAAMgJAHGCpXDa6lRwR0EGOSjiggX\nEmBssaEUUCwBBjmoIsLlBFguLmB0XaeKCNcRYKzTxwEUToBBJibhcC0Bxgf6OC6nigiXEGAAhCTA\nsnP73TBVRLiQAOMtfRz3cRsD5wkwAEISYHymj+NCqohwFQEGz1BFhJMEWEbuvQGuI8AgN3cycAkB\nxhaNiLdSRYQzBBgAIQkwdtGIeC1VRDgvXoCll41XMw/pa+WPkFy8F+Cwn6cH8J2UUv+6d50+nj2z\nfAmAysSbgX0r/fb0cOLRx3ETt1iUJtylMtgM7KPhuE+nX6ZiFC4lYUYRplfLEBlWzwxsKBsOQhz6\ncPRxAEWpJ8AgHBMvOCNYCXE6u5p1c6y+VJy+13bGkioiHBBvBjbWCafPvHuJS+jjAAoUL8CgJm63\n4DABxhf0cdxHaRm+JcAACEmAPcT9Ni+qiHCMAGMXfRwZuKuBrwgwAEISYHxHH8cdVBHhAAEGBVFF\nhP0EGAAhCbDsVItY430B3xJg7KURMQ9VRNhJgAEQkgDjaxoRb6KKCF8RYFAcVUTYQ4A9x1UK4AQB\nxhf0cdxNFRH2E2BQIvNz+EiAcYQ+DuBxAgzKoooIOwmwJ7hEsYMqImwTYHxHHwdQCAEGxTFFhz0E\nGAfp48hAFRE2CDAAQhJgj3KDzRtjFTElbxNYJ8D4mj6O/MQYLAkwKFTfz7s5xBhMxQuw9PLu+dVX\nuYM+jgxWYwzouu7n6QF8J6XUv87m6eOu66aPBRiVGd7d4/t6eKDbnsbFm4F9NAs2qMZsNqaiSOOC\nzcAOmM3GSsm2vnft4ZjZe8dsjKuEq13VNgNbTr/6354aWGU0Ij5Lfwd3CHeprC3AoB1ijMYFC7C+\n78c+w2k3x/ggyo1DNTQiPk6bIs2Ktwa2jKjxGelFs7Qp0qBgM7AKuVvmOqZiNEWAcZA+jjItK4pQ\nKwEGFZruBQy1EmAAhCTAAAhJgHGWTvoyWQmjegLsOS4wACcIMKicPg5qJcAACEmAARCSAAMgJAFW\nAGsU3EOfEHUTYBxnNyngQQIM6meST5UEGAAhCTAAQhJgj7LIzs28xaiYAOMCtkME8hNg0AR9HNRH\ngAEQkgADICQBVgb1HYAvCTConEZEaiXAOMVuUsBTBBi0QqGayggwAEISYACE9PP0AL6WXnWQfrE2\nvfFSufpeZYe7eZdRpWAzsJRS/5J+n5HDb1dfIgO7SQGZxZuBbRtj7OmBQIlS0lVPPaoKsDG3hona\n+Hj1awCYCle7qirAVkksgD2mV8sQYRZsDaxmEd4uAOUINgObNmjMCoarLwEDjYjUJ94MbOxCnD7z\n7iUysJsU8Ih4AQacYR5GNQQYACEJMABCEmAFsGhHFt5oVEaAcRm7SQE5CTBojj4O6iDAAAhJgAEQ\nkgAricoOwG4CDBqiEZGaCDAAQhJgXMB2iOEoV1MBAQZASAKsDJYmAL4kwKAtbpaohgADICQBxpVs\nhxiIPg6iE2AAhCTACuOuGGAfAQZASAIMmqMRkToIMABCEmBcw25SEVlyJTQBBkBIAqwY1iUAviHA\noEXul6iAAAMgpHgBll42Xlp9lTzsJhWLc4W4fp4ewHdSSv2r9jF9PFo+A0CV4s3AttUw/Yo+foAs\ngs3APhpmYLOJ2vILAJgJd/dfVYCthpPEglV9b7bPL9OrZYgwq6eEGOJwA3CVYAHW9/3YZzgrEq6+\nRE52kwrKvR9BxSshbnQeyi2AdgSbgVVOAAPsJsCgXW6ZCE2AARCSAON6dpMKRx8HEQkwAEISYEVy\nPwzwiQADICQBBk3TiEhcAgyAkAQY0HUWXglIgHEl2yEC2QiwwliRANhHgEHr3DURlAADICQBBvyh\nj4NYBFipgl9LbIcI3E2AARCSAAMgJAEGaEQkJAEGQEgCDPgrePMQbRFg5QlezbGbFJCHAAMgJAEG\ndF34mT8tEmAAhCTAAAhJgBUseEOY3aSCCv6+oyECDICQ4gVYetn4gpzjAeARP08P4Dsppf7VLDV9\nPP2C7IOCSvS9+iGRBAuwbUOkzTJs9ttl5gHQBZwAVBVgqyQWfCslHwtr0fRqGSLM4q2BvTMc7ul/\nAwt+8bCbFJBBPTOw7bUxACoTLMCmS1wSCy6nj4NAggVYt7amNXtGmAG0oJ41MACaIsDKFryaYzep\nuIK/9WiCAAMgJAEGQEgCDPhFFxRRCDAAQhJgwDp9HBROgJVKHQdgkwDjFrZDBO4mwAAISYABcwrY\nhCDAimclned491EyAQZsSUmMUSgBxr1shxjUrIo4xJgkoyjx/jkVII8hw2ahNf7WOhmPE2DAljGo\nJBmlEWDALttJJsbIzxpYwVwSKFLfr7w3LZKRnxkYcITSIo8zA+MudpNqxDAhezcng/uYgQHX0LVI\nZgIMuJLSItkIsAhSctITjiTjbtbAgHttLJLBGQKM29lNisEyxmQYZwgwIKvV2RgcIMCAB8gwzovX\nxJFeRYd+cQZsvASUSYsShwULsJTSGE7Tx6PhmdWXQup7qwQAq6oqIVYSWgDsUFWAdV2XUkq/pyzp\nt6cG1ia7SbHBDWdpwl0qg5UQPxpLiLNngGJZBivE9GoZIsPqmYGFONwAXCXYDKzv+2Wr4dCysfpS\nPdyjAvwWLMC6tXAan6kwt6Bq2mw5o54SIiWzmxTbxBgHCDAAQhJgwJMU/jlMgAFFUEXkWwKseG5Q\nAdYIMABCEmDAw1QZOEaAcS/bIbKfZTC+IsDicHIDTAgw4HmqiBwgwICCKDSwnwADICQBRia2QwSu\nJcCAIlgG41sCLAJnNi2xDMZOAgyAkAQYACEJMKAUiuV8RYBxO7tJ8S3LYOwhwEJxWgO8CDCgIKqI\n7CfAgBIpN/CRAAMgJAFGPnaTAi4kwICyWAZjJwEWhHOa9lgGY5sAAyCkn6cH8LX0uivrF5OSjZeA\nQPre9IvPggVYSmkMp+nj0fDM6ktAOCkpn/NWVSVEoVUsu0kBlws2A9tjNv1KvysR4UPOHSlwjxSt\nbltVgA1HfxZR4RMLmmQZLL+NW/8yVVVC7MQVVCfChZRnBJuB9X2/bDUcaobD8xoRARoRbwbWv0yf\nmT4/e5Wi2E2K/ZzHbIsXYADQCbBI3I7SKstgrBJgAIQkwIByqTuwQYABAagisiTAAnIqAwgwsrEd\nInAtAQYUzTIY7wgwIAa1c2YEGAAhCTCgdKqIrBJgoVRxHtsOkcNUEZkSYACEJMAACEmAAQFUUT7n\nYgIMiMQyGCMBFpOTGGieACMfu0lxhioiMwIMgJAEGBCMCjoDAQZASAIMCMMyGFMCLJoqzmC7SXGS\nKiKdAAMgKAEGQEgCDIikiiI61xBgQEiWwYgXYOll4wtyjucxjfxvArwRLMBSSv3LMqi2g40S2E2K\n81QRGQQLsG1DsD09CiAT96uN+3l6ALebzckkHNQkJROyy4SrYNUfYBIL6tP3f6dfMuwq06tliDCr\nqoQItGMaWhEutlwvWIANvRuD8WYhxJ3Claq427SbFOfJsMbFKyEuS4KzZ9QMoR1qiS0LNgMDmJnN\nw0zF2iHAgPD6XjmxRQIsMqcpTMiw1ggwoB4yrCkCDKiKDGuHACM32yFyN20djRBgQIW0dbRAgAHV\nkmF1E2Ax+bgm7CPDKibAgMrJsFoJMB5jO0Sy0dZRJQEGNEFbR33ibeYLcNhs899LviFPMQMLzm0k\nfOnayHEKPsgMDGjOJRk2m8mZiuVnBgZwxCyx9IbkJ8B4gN2kqMOsMaQTY3kJMIBTxNhTBBjABVZj\njFsJsLAsGUN5lp82E2P3EWAAFxNjeQgwnmQ3KSpmYexuAgzgLvo7biXAAO4lxm5iJ474xvNAWwcU\nbDhBp7k1fez0PUCAVUSSQfGWMTZYPuM8/kgJ8V7p1jLBsjDx+lv//tp07/BOM7yTCh9hy8Mbzt3t\niJqex2uBV/TRy0OAxTeeCtthVhi7SUH3+/Q9k2dtqqqEON6S9M3Ovcf/8e0KRbPHB8q2PDXfZ1U/\ne6nB07qeAEspjbk1fdyo7X96dkz6XMMBjtm/PVWDeVbPhX41wJSJZyr5YQMZFJ8O9czAVlUTz9cT\n7UBwlQcYb4l2IDhdiACEVM8aWKcLEaAlVQUYAO2wBnajaQ9kaTcKs6bN4UE5g1wdXlfACJfHqrSj\nNxtPUUevez+8EsbWFX/0BoWfvDkJsHsV+K6afbSgtM/PrX7y4fFRTY2Xtr7vSzt6g+kIu8KOXvc7\nGxy9b01PkDLffjlp4rhXSqm0z6L1fV/yG311eOUcxpIP3aDko9cVfwALP3ptptQGM7B7ze7mOKa0\nw1jOSN6ZjrDAo/f0ELbMhlfa0WNKgN3IO/4SRR3GWeGrQLMRFjjUwnfJWa1wlmAY0vjfosb2FCXE\nuxR7fsZS4GEs/8Kx2v9SgtLGM7NcHn5qJKv6ly7CmzAPMX6jkhuECm9kKnN474pL098+aznCoobX\nhe1CLGR4gzLPjkcIMABCUkIEICQBBkBIAgyAkAQYACEJMHhAaS3aEJEAAyAkAQZPMhWDwwQYPMaG\nQHCGAINnmHvBSQIMnjHsayfG4DABBk+SYXCYEjwAIZmBARCSAAMgJAEGQEgCDICQBBgAIQkwAEIS\nYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQ\nBBgAIQkwAEL6D4qEe3yAy4MnAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo3()"]}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Restricted Isometry Criteria\n", "----------------------------\n", "The restricted isometry constants $\\de_k^1,\\de_k^2$ of a matrix $\\Phi$ are the\n", "smallest $\\de^1,\\de^2$ that satisfy\n", "$$ \\forall x \\in \\RR^N, \\quad \\norm{x}_0 \\leq k \\qarrq\n", "      (1-\\de^1)\\norm{x}^2 \\leq \\norm{\\Phi x}^2 \\leq (1+\\de^2)\\norm{x}^2.  $$\n", "\n", "\n", "E. Candes shows in\n", "\n", "\n", "E. J. Cand s. _The restricted isometry property and its implications for\n", "compressed sensing_. Compte Rendus de l'Academie des Sciences, Paris, Serie I, 346 589-592\n", "\n", "\n", "that if\n", "$$ \\de_{2k} \\leq \\sqrt{2}-1 ,$$\n", "then $\\norm{x_0} \\leq k$ implies identifiability as well as robustness to small and bounded noise.\n", "\n", "\n", "The stability constant $\\la^1(A), \\la^2(A)$ of a matrix\n", "$A = \\Phi_I$ extracted from $\\Phi$ is the smallest $\\tilde \\la_1,\\tilde \\la_2$ such that\n", "$$ \\forall \\al \\in \\RR^{\\abs{I}}, \\quad\n", "      (1-\\tilde\\la_1)\\norm{\\al}^2 \\leq \\norm{A \\al}^2 \\leq (1+\\tilde \\la_2)\\norm{\\al}^2.  $$\n", "\n", "\n", "These constants $\\la^1(A), \\la^2(A)$ are easily computed from the\n", "largest and smallest eigenvalues of $A^* A \\in \\RR^{\\abs{I} \\times \\abs{I}}$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["minmax = @(v)deal(1-min(v),max(v)-1);\n", "ric = @(A)minmax(eig(A'*A));"]}, {"source": ["The restricted isometry constant of $\\Phi$ are computed as the largest\n", "stability constants of extracted matrices\n", "$$ \\de^\\ell_k = \\umax{ \\abs{I}=k } \\la^\\ell( \\Phi_I ).  $$\n", "\n", "\n", "The eigenvalues of $\\Phi$ are essentially contained in the\n", "interval $ [a,b] $ where $a=(1-\\sqrt{\\be})^2$ and $b=(1+\\sqrt{\\be})^2$\n", "with $\\beta = k/P$\n", "More precisely, as $k=\\be P$ tends to infinity, the distribution of the\n", "eigenvalues tends to the Marcenko-Pastur law\n", "$ f_\\be(\\la) = \\frac{1}{2\\pi \\be \\la}\\sqrt{ (\\la-b)^+ (a-\\la)^+ }. $"], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 4__\n", "\n", "Display, for an increasing value of $k$ the histogram of repartition\n", "of the eigenvalues $A^* A$ where $A$ is a Gaussian matrix of size $(P,k)$ and\n", "variance $1/P$. For this, accumulate the eigenvalues for many\n", "realization of $A$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAYMUlEQVR4nO3d7Xaq\nuhYG4GSPc8HrUnrHnB9UyreoCJnJ84yOvVt1VaHI60xCkruuSwAQzX93vwAAeIcAAyAkAQZASAIM\ngJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIA\nAyAkAQZASAIMgJAEGAAh/e/uFwBh5JyH77uu23/M8IDZLcsHLH/D1l1bD3v6O6FKAgxeMM6MrbQY\ngqrrumXMHPkN+8Y5mhZJJsNohyZEONOJ+TELqvFTSClIKjB4yVZ74PjGdEYlNPyGracABBi8YJYf\nyziZtRO+50hnGyDA4H2r5dHnkTPrRVveBSQBBp+YxUkfNuNmxnECbY1C3GpvnI0B2XkNRiHSJmOW\n4GaGDsJ7vHMACMkwegBCEmAAhCTAAAhJgAEQkgADIKTKA2xrNjkAoqs8wAColQADICQBBkBI5kKk\nFDvLX138SoAQBBhF+Vnc8u+GVwFEoAkRgJAEGAAhCTAAQoraB7a6hJKF2AHaES/A9ifXkFsAjYjX\nhNgv0751b87Z9FEVGv9N/X2BlFLcFZm3mhD7G8ffzB4TdHvrsPPZouu6nPM5f5vHn9iFZVC3eE2I\nO1ZPTM5WhZlf6dX1V3qdlV7pr0Trf2GeP6MLy6AS9QTYak1GybrtLBmlzr9Df9ad2m70LIswAwKr\nIcD66OrboPpbJFnJtnIrp7Q2E8cTm+2EO0+dc3KEQHxRA2wcUcP3cqt8y/TK6ee3zHp/dMbaBFS/\nH2h+Vp/077kcMxBW1AAjnmkv10uteR8Oxxg/1zzMHr85v/ILgRIIMC4xSqB3O6LOmec3b4x4HEYu\nqskgCgHGl02Lp3KGUQyvZKsmk2RQOAHGN43TazTKpij7STa8Yq2LUJp4M3EQxjS97nsdR+WUUtfN\nMrZ7fAGlEWB8R7T0GsvpZ6WpM2ezWEFRNCHyBZHTa7DetNhvWtiNgpqowDhbFek1ltPPvPLqqzEF\nGdxKBcapqkuvP/3mzELrUZDtz1P81dcFzRJgfEetZ+1hu6YrvHSbVwiYOxi+RRMi52nqCqqum21m\nl/7tTE8MnE6AcZI2O4TEGNxHEyIny6m9MFvMTdVnWDnTjkCVBBhnGDcePuaAX6i8LsmLme+VYvBV\nmhD5WGv11q6Vi6ANuIfvUIFxnhbGbhyzrMZWl4Qxwh4+oQLjM02NPHzR7Arox5yKP48v4CMqMDjq\nzdn0p0M8jO+As1QYYDlnLTMXabH8enN8Sk4/xnfAuaoKsDKXm6pSv6u76Y/sW+0Yayn74WRVBVhf\neDmZXuMvvSZFicLiiXmMmd4e3lX/II48dffLgZSWfWCG2sPrqqrAVukP+yqDEd6W+4NzOinw8IPj\nFp6qvwLjG5xcz5KnHwI6+xYOE2Bwv1khK8PgiAoDTNvL1z0aurQfnmg+B5VeMXimwgCDuFYGdwAb\nBBhvUn59yWwCKqUYbBFgvMjJ9BqzlnC7HRYEGJRqttyzDIMpAcY7nEqvM8swMQYPAoxXOHveQnMi\nrBFgEIHmRFgQYBCHDIOR+udC5HyzGfy40njn55ym/ZGu4qcpKjAOE1qF6LpJaKWU0s/GSptQMwEG\nIU2nALYMGy0SYBCVDKNxAowX6WUpyWwpFmiKAOMYHWCXyxvmD5vNYQ/NMAoRSrYcmrHSVJjTz18T\nYs6qZBqhAoMaqMNokADjFT7aF2y+CAvUToBBRUzVQUvi9YENndjLSQfG/dumJKBRs6k6vBGoV7AA\nyzkPyTT+fiC3TvT3WWF6i51cOhlGG2prQlwdZ8wHZnMUma+odCtj7b0pqFSwCuypvjiYFWrLB0DV\nflJKeTQ9h4OeKlVVga2GUzd1/auCu2QVM1WrJ8A0knybs2FEf381bxCqE6wJseu65SjEvsFw9S7g\njwEd1CVeBbZsDBy+104ISybpoFbxAoyLWaejAlKLKgVrQgTeNFwctugu1mhBUAIMGjQej6PCJipN\niNCKofDSLEwdBBgAIQkwDnERWB2Gv6MijAoIMABCEmDQFkUY1RBgAIQkwKA5ijDqIMAACEmAQYsM\nK6UCZuJo3dYyNKYXaoQ/M3EJMFJa+TA+rOSrjwQolCZEAEISYNCocTdY3nDjy4OnNCE2pF+6+u5X\nQbE2W5Lb4T0yCLErqgqw4QNj+fsdgA/VE2DjzwshPjtAKTQVElM9AcaOoTad9WqIeZ7a6Qlz/HCv\neiqVWQU23F7J5t3K5/NaeXewY/zGLzMpqq3A/na35pGPlXnscg5vEDaU/8avNsD+FP83gDt5gxBW\nPU2IyShEgJZUFWAAtMNMHACEVHMfWOMtijubPxml2eTOSc1fLLjc/MaPisZPF2l7D5R8YFQbYI1f\n1/x081vbIWONT/Hnuq4t/eY3eLoYbO2BYneIJsRGtTxVa9d1xb4hL7Cz+Y0fFXe/hJvt7IFiDwwB\n1qj+LFbmQcldHBUtl1+9rQabMg+MapsQ2dH4W5RVjR8V/dm55Z2wtQdK3icqsOYU+DGK2zkqUtln\n6mvsD/gqUM31cuPDipabPzQONL5neo03Fo0PBkeFSa6XeyDEgdH0exiAuDQhAhCSAAMgJAEGQEgC\nDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQrAcGR41n7N5Z1Hj2gNktT+f2fjpN\n/tOngEYIMHjBcm2arcf0Dxg/bLZg4Cfruew8hQyjHZoQ4Uwn5sfWWoIiCnoqMHjBVnvg+MZ0RiW0\nXE5w53mhTQIMXjCLpa0l2D9Pr52nSGIMUkoCDD6xWh593sQ36+JavRfQ5QtHHR8f2Fsddrg6hnD8\na58OyjAKEXoCDG5m6CC8xzsHgJAMowcgJAEGQEgCDICQBBgAIQkwAEISYACEVHmAmWsHoFYVBpjQ\nAmhBVXMhii6AdlRVgXVdZ2IRgEZUVYGt2lqxid7pqybu1MF2PnCi+gPMSfOAn8Ut/87+hZ//ToCJ\n+gOMj2yUUzmlbvojwMUEGH+6cZG0OyKm2/1xXGzl9WoM4FOVL6dipaXnXhy6mVfi6hX+HMBJVGAN\n246unNJjMfu17rHNu1JK/54E1PJJF8sKL+6XecAKAdaktaiYtvW9P+Aip7TMtm7nFz5eTLfe3mjo\nB7BOgDVmEV2reXP+06aflP791VJbxdY0rvSfATsEWDO22u4unL5ks5Fw4/G/edb/q42GRA2P0CwB\n1obpWT4/brzjLL9SVOX0b9zltt7eON6E+cvW8AgtEmANGJ36z+ro+qrxi3weZkCrqpoLkRWb6RVD\nTil1Xd6+VrpLqUv/9gaJAJVSgVUteHpN/YwzbJlY3eTqaaB+KrB6TdKrNrnPs65b3bQupZTz7xdQ\nKRVYpWZDHqo+jw/F5ZMOM+MSoS4CrEZ7A/ZqNmomXZsQpNXdArUSYNVxmu7tj8tXlkF8AqxeTs0P\ne22M8h7CMoijLgqLXcPQj/X98xj3serq1wo8owKriJPsS7YnZpyu1dlXb64zg+IIsArlNDkpmxXw\nifH+WQkz0QWFEmC1eJx5FxcsO/+u2JkCuB/6obcMyifAaNbeFMBP5mN8sa/RlPnwDQKsCpvlFyf4\nW8xsmUMvlGWmzIeTCTA4bLe3zBBQuJgAi++v/OJCO6tLT/8i2gnhS1wHBp/puseCL/NGwu4xrbAE\ng29QgQU3brZyHdiXPb2ceWfKj+EW/ZRwFgEGL1mNn3lcHUsy4COaEKugl6VUoxWlVxoYdyavuuG1\nQjQqsMic5kLZLctmjzHCHp4TYHC1IaWWhfNftuWssIZ9Aiw+p7m4LFoGHxBgYWk/rMuhRcskGYwY\nxAFl+Vu0bOW+30EfV78mKJIKLDgfyes1jqn5n9nU+CDAwukHWHfTH5P5iqr126445JUGRhhoQozo\nZ/q9mR0aMlxYtnafBkbaogKDkAz6AAEWmFn1SJPDYHFhmSSjavECbKfXZzwBT8V9QiufuKH3e2HZ\n3qCPZQtjxW8W6hYswHLOw5tt/P3AWxFSPxA/pWTaKqoWLMCe+v34OYqx2byoEo6mHJq2CmKqbRRi\n13XdoxVlfMvgxtd2Lh1gvGZjUvyUDF8kqqoqsJryaUv9W8iXGb5INeqpwCyhBC/J6Wdz2io1GREE\nq8DGzYOz0RyrdwGHDG+ZWWipyShYsABLa+E03NJObukA41skGXHECzDgCgeSbKfdvp1Pk9xIgIWi\nT4LLjS4pm92Ru98HLNsDDNDnCgIMeOonbc+IP/yoZZuL1TMKEbhGP3xxeXuX/rk4misJMOAd+fdr\nnmRdP+uH5m6+TxNiPBpqON3WcIwjYzE2r4w2cJEvE2BA79OxGE+STIxxNgEGnCz3pdvr4+8Nvucl\nAiwOnQqEklPqq7plQdYZfM8ZBBjwXatNiwbf8zmjEIGL5NFKmwOD73mbAAOutbYyWbe25Cbs04QY\njH4wqpGXPWRG3vMKFRhwp/V5PaxGxgECDLjf+tKaOaf8656XRdkEGFCQte6xxUT4kFISYECBVtoV\nNSqyIMCC8NalPWKMfQIMKNrK1WOjvjHdYy0TYEAEXTdLqi6ltQmIaYjrwIAwZpeOmcKjcSowIJiV\njjGaJMCAeObjOwzuaJIAK1FeuPsVQYnm4zu8UxqjD6xYs95pbf2wYbx4Zs7mUWyHAAPCy+PZOnJO\nj2mvLfFcN02IQA1mvWKCqwUCDKiHDGuKAAtFewg8s7I4C5USYEBt/jJsOaLXyN6KGMQBVCinn36e\njm69JjOstwYqMKByZpyqlQAD6vRGK6GmxUGIXSHAgHoZ91S1qgJM9yywSitileoZxJFzHq66H38P\nQJXqCbBNAQuyLiWjpOBc8yJs7czQbdzeoBAVQD2VyqwCG26vZPNSSm91SgM1nQSuND7hlJkU1VZg\nZe7uD1W4SUCpyj/hVDWIA4B21NOEmEYthzVtFACrqgowANpRbR9Yar4g29n8ySCXJndOav5ai+Xm\nN35UNH66SNt7oOQDo9oAa/yysKeb39oOGWv8UvedzW/5qEiPzW/wdDHY2gPF7hCDOBrV8pQlXdcV\n+4a8wM7mN35U3P0SbrazB4o9MARYo/qzWJkHJXdxVLRcfvW2GmzKPDCqbUJkR+NvUVY1flT0Z+eW\nd8LWHih5n6jAmlPgxyhu56hIZZ+pr7E/4KtANdfLjQ8rWm7+0DjQ+J7pNd5YND4YHBWz07Q90DcY\nln9gNP0eBiAuTYgAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgAD\nICQBBkBIFrSEo8ZLTmwt47C6is34lqeLUzxd52X5Mkpe8AK+R4DBC5aLq209pn/A+GGzFW8/WZBs\n/A9nTyHDaIcmRDjTifmxsxhuzrnwpXLhAioweMFWe+D4xnRGJbRcD3f8FLPFlKFNAgxeMAuMZX7M\n2gnfs9/ZJrSgJ8DgfTvl0SdmvWir9wKaIOCol8YHpo0hgqvDFF8alPF0oCM0QoDBzXRlwXu8cwAI\nyTB6AEISYACEJMAACEmAARCSAAMgpMoDzHxxALWqPMAAqJUAAyAkAQZASCbzDWarV8+MKkBrBFhE\nP4tb/iXZBjRGgFVmPdsA6iPAQuqWsZRzl9J6XC0qs+FnxRkQV+Wz0de2UMX3Lmtb20vaJIGSxavA\n9tfuO2VB9xJtZEn+bTP891i6d6UJsb/ryR6Z/f6/HahNEihUsADbX6x2uKWqwmulAXAZKgd+zeJf\nrbRDTp+0b5N87+kAvi1YgLVlWXX9llkn/fqUxnXbVp6NbxdmQDlqC7BlE+LsjB+mMtts0/vaE07C\n6d/q8/2GWc5pNBJk/pgoexgIrp4A22pdDHk+HafXXa//WXE2flmj8NNDBlykngCrxyi98vTHu8J4\nXJxthJncAq4WLMC6USfQrN5avSueSXrNOpyKCIlxz9l6bpVQPgINCBZgaS2cYrcWju2lV4meVmbD\nFvX/C/8HAkoSL8BaEHEVzuGKtDTtHusNxXJKyjLgHAKsGEP51XVfnHHjEk8qs402RhN/AC8RYGUY\np1ddxgMUV6dOmf5s4g/gKAHGdXJKz8bl/3s8MkAXIHAvAVaAesuvLcfH5cduSwW+6b+7XwCty+mn\nnz1lNav6qS37r6VrXylQFhXY3dorv3bsV2aLuT90j0HTVGAUKqefoThb3tv1Q0L64gxokgrsHvlv\nvZKUplNGGTW+tD/6w9wf0CYBdqNxI5jJcA85PvdHSsIMKifAbtMZMv6ZPK5WVxsShRlUTYAR23Is\n4npSCTOoTrwA25lyPuJs9MqvM0z2YU4ppX/9lPlPw2yWfoGOHCBYgG2tWpmmazHP7ipQ0S+uIgeu\nmB4/GIgkWIA9NY6x8S2DwoON7xnmy+82pkv+HZf/+4PjBEpXVYCtFmclJ5b2w7sM4/LTkdGMswEj\nQBmqCrAwXHtbkifz5aeUFGdQJAEGI12Xc36hOBNmcJ9gAdb9nl9+v++/6RsMV++CTzwvzoQZ3CdY\ngKW1cArR3bWkAyyYUXF2pNtMnsG3xQuw8HSAxTf78KGlEW4hwOBTj6RSnMGlBNg9VGG1ej7dcFKc\nwTmsBwbfMl5serPL87HetLXN4FUqsGs5QzXs5eIsqc9gjwoMbnCoOEtJcQY7VGBwv5eHNSbFGajA\nbuHUwy7FGRyhAoOivVGcrQZarMv84QgV2IV8UuZjfXH297Wmm34lc75QKRUYxJPXPgxtT6X/r/83\nj5+VYlQiXoA9nbG3/OWY4Qzzuiqnf93TCRuTpc6oR7AAG4fTalCtfjKFBk07z44tdZbUZ0QSLMD2\n9ZFWeoY5QXCLI8VZei3Mdt5rqjouUFWArZq9x7yv4GBxduzKs9URItsBCeepJ8D6oBr+G3SRMLjB\nkXWo06HB+nClegJsv2/sfoU3bMLDZB3qrts6dB/vsX+LfwUXCRZg4y6u0hMLqpBTSgemId7rV4Pv\nCBZgaa1JcHaLMIPvGddnafviMyMbuUC8AAPKYY0YbiTAvmg8ALJ/y/YjTG56OfBdBzvP9vNsa2i+\nlhWWBNi3TdpbHj/qLaB+BzvP1nJuOR7EW4YVAgy4wqM++53vam9mECMbOUaAATc4skzM5MacdZ4x\nYzkV4H45/QzLeG4/KE++aJ4K7AoukYHjDq3hmYzURwUGlG1cnG32iinOmqQCAyJ5+cozxVm9VGBA\nVC8XZ9RFBQZU4nlxdnhCfddNhxAvwJaT+R65C2hKTj/HrjkbP35g1FUMwQJsPPH86iT0/S3mpwcG\nr15zpqkxiqr6wAoPLdMKQAmeXnPWfwrWbVa+YBXYEbPyazY3aOEhB1zphW4zp47yVBVgfVZZHgx4\nw2Q2/ZW7Dc0vTlUBlsQVcIrxmWTZkKgsK0OwAOvHFA3f99/0bYb97QYiAicbTiY7SZZSdtq5XLxB\nHN3D+Jbx7bN7Ac7RdTsXTf8N/eAqwSowgAtsLQz994BHhq0M/dDAeBUBBrBq9bqXeVwdSTKti18S\nrwkRoEA7szL+ti5yNhUYwMn6DJvXZJoWzybAvs5qltCmzdbFPsnE2Mc0IQJ8V9+6uLjVkMVPCTCA\nK+TH1/RWGfY+TYgAl/lJjwz7a1fUovguFRjADXL6mYxXVIq9TgUGcJs8Xleznw9vdK+rx/YJMIBb\nzZeN/psU/45X8yvEssCaEAHutzejBxuqCrD8cPcLAXjZaoblDTe9xrLU04Q4LnhDFL8AxyxnZVSl\npZRSPSf6WYANtxeyeT4vAU8Vcr7qlT+cpJ4KbKbbWYPuDmX++YHilHHKShHOWtUG2J/i/wYAf5yy\nDqtqEAcA7ainDyyNur5q2igAVlUVYAC0o+Y+sMYLsp3Nn4zSbHLnpOavtVhufuNHReOni7S9B0o+\nMKoNsMYvC3u6+a3tkLHGLwLd2fyWj4r02PwGTxeDrT1Q7A4xiKNRLV/M33VdsW/IC+xsfuNHxd0v\n4WY7e6DYA0OANao/i5V5UHIXR0XL5Vdvq8GmzAOj2iZEdjT+FmVV40dFf3ZueSds7YGS94kKrDkF\nfozido6KVPaZ+hr7A74KVHO93PiwouXmD40Dje+ZXuONReODwVExO03bA32DYfkHRtPvYQDi0oQI\nQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmA\nARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEIS\nYACEJMAACEmAARCSAAMgJAEGQEgCDICQ/g9bmTmvBnYdpAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo4()"]}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["__Exercise 5__\n", "\n", "Estimate numerically lower bound on $\\de_k^1,\\de_k^2$ by Monte-Carlo\n", "sampling of sub-matrices."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARK0lEQVR4nO3d0Xqb\nOhaAUeh33nvaJ/dc0FIKGAMGpL211kUntTOpjuPwR0LG/ev16gAgmh+lBwAAZwgYACEJGAAhCRgA\nIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACH9V3oA\nULW+78d3LZ9+vPy04YPpJ09vWX7Cu6+w8Tl7vuzHfwjSEDD4YKNb4yd0k6K8Xq9Z9rp/A/Puq+35\nnO7fdM3+LztzCzlYQoQPhjB8/LS+7x9ohibByAwMjpnFbCjKJfOed8uD0xuXt0OzBAw+m07C7psD\nzb7y9kqjjIGAwTHvJkYP/ENOa8GUnwfYct8uxNlX2xknuxBhJGBQhukUfMmPEAAh2UYPQEgCBkBI\nAgZASAIGQEheBwZQXv/z35f9/bS97rPkuxDtVAbqN6vXHgrXRQzYu9dprl4fQcCAyp2o1x4tFC7Y\n8X3/ZREEDKjftF4bybkjcgkKF+z4vidg3hIJCGFnvY5+qavUX7j8mzgeu/QqwH4X1mv5FVYuknnP\nQmVZ2QK2nHIpFvCkE+900/+6qy7HBvO/m0Zxl2ArbB+XB89d4RvgKvUcdg6NpJ5h7xdsBrb6voK2\nbACNW32fnfTHw2AB69a+JeMt6b9bAO9Mf7/v2jgeupQUQHjTN+weMtb3/YmzcbHEm4EBsGpoWAtz\nr4EZGEB4Y7rSz7qmBAwgtqaiNSVgALG1s2Y4I2AAGQy7NppaRbSJAyC8NidhAgZwmdVLDt40I9pu\nVgsXdrCECHCNei6Y20K9OjMwgEts1Ov5lExf15yYgAF8pZ6J11QLDcv+n5f9+weUtXxbr6CHnYjD\ndg4M4KRr35SSowQM4Az1Ks45MIDDKq/X8n0TUzIDAzim/no1ck16MzCAvWYbDiusV5d91jUlYAC7\n1Dzxmq0ZDn9NXzJLiACf1VyvwSxX6evVCRjARyHqNb6Qa5h+DRenLz2ue1lCBNhSf71G7WzfGJiB\nAbwVpV5jutLPuqYEDGBdoHqVHkIZlhAB5kJslx+1NvEaCRjAP6JMvGbGrfMRL8t7joAB/BW0Xo0U\na0bAAH67qV43vWHY9ghbmIfZxAHQdWHnXqtaqFfnDS0BukvrVfywMwzg6DCKD/sEMzCgdZnmXqMW\ntiY6Bwa0K2u6Zh9kZQYGNCplvZpiBgY0J9brlHlHwIBWrG5nV6+4BAzIT7pSEjAgs2W6dCsNAQNy\nkq70BAxIxWphOwQMSEK6WiNgQHhWC9skYEBg0jUYrxr17uobs0sdRrzy4ZKAASFJ12hao9UyTS+K\nmOkCifECtvGLxsffQYDonOg6arw4/fDX4fCYI2PBArbxi8b4dtrLu4AEAqer77vxiDT9+JK7GhYs\nYB9NMza9ZSRsEI7VwtOGA+D4Z7IDYKqArU7Okn3DoClJ0jU9Cs2OSN/f9eFf3jo3Fl2qgAE5BF4t\nLGF5iqtLWqwZAQMqIl3nLFs1u2X7r0EFC9jGLxqrdwFRJFkt5EHBAtZt/qKhWxCRdHFOvIABOVgt\n5EsCBjwtfbpyvEy4fgIGXGy1TxsypatzLuNBAgYcTs4lknWL5wkYJFckTtuki0sIGGR2bb2Eh6oI\nGKQ1rZf2kM+P0gMAbqFepCdgkJB60QIBg2zUi0YIGKSiXrRDwCAP9aIpAgZJqBetETDIQL1okIBB\neOpFmwQMYlMvmiVgEJh60TIBg6jUi8YJGISkXiBgEI96QSdgEI56wUDAIBL1gpH3A4MwxnpJF3Rm\nYBCFesGMgEEA6gVLAga1Uy9Y5RwY1MuWDdhgBgaVUi/YJmBQI/WCjwQMqqNesIeAQV3UC3YSMKiI\nesF+Aga1UC84RMCgCuoFRwkYlKdecIKAQWHqBecIGJSkXnCagEEx6gXfEDAoQ73gSwIGBagXfM/V\n6OFp6kXN+j9Pz1f1z814Aev/PLqvxaM73rV6L9RAvajQ5NgZSbCA9X0/lmn68Ui3qJl6UZWNboU4\nlAYL2EfDJEzGqJB6UYno3RplC9iQrtlEbfkJ8DD1oqw00ZpKFbDVOCkWxakXpaTs1ihPwFZPiUFx\n6sXzcndrFCxgr9druQtxSNfqXVCWevGkRro1Sj5rMS2jIPXiAa1FayrYDAyiUC9u1XK3RgIG11Mv\nbqJbUwIGF1MvLveuWw1Ga0rA4ErqxVVMtj4SMLiMevE93dpPwOAa6hVLoMvX6tY7AgYXGOslXfUL\nkS7R2kPA4FvqFcUyXQU7sdpR3TpEwOAr6lW/OlNRfAAJCBicp16VqzNdXEXA4AxbNipX1WohNxEw\nOEy9aiZd7RAwOEa9qiVdrREwOEC96iRdbRIw2Eu9KiRdLRMw2EW9aiNdCBh8pl5VkS4GAgYfqFcl\nvKiLGQGDLepVA+lilYDBW+pVnNVCNggYrFOvsqSLjwSMwqad+OixkKhXQdLFTv0r9VOj75P/B0Z3\nqF7fOBQh9SpFujgk+fFdwGq22onHkraHej1Gujgh+fFdwKp1+Szn8vKp1zOki9OSH98FrE71rNG9\ny5563c3OeL6X/PguYLWZBUMnGiRdXMUuRJ5Tz8SrcasJKUK3+IaA8RD1KqKeVs1IF98TMJ6gXg97\nrFs6REECxr2c9HqYM0y0Q8C4kYnXk6SL1ggYd1GvZ+gWzRIwbqFeD5AuGidgXMxJrwe4egV0Asa1\nTLxuZcoFUwLGZdTrPtIFSz9KD+B+0x/92WHAXdfd1f/sX79+3/a7XpWNMOhdff/7tlf3+67Xq3u9\n/v61+Ajdlf+uWpmB8S0nve7wJ1r/MOWCqeTXunUx37tZNryc1ULYyQyM89TrWvYWwiECxknqdRVT\nLjgnXsD6Pz/u79YGLRvezUmvq0gXfCNYwKZxWg1VH2TzTFwmXt/TLbhEsIBtG5KmYfdRry9JF1wo\nVcBWzXpmdfE09fqGdMHl8gRsCNX45xgqxfqek16n6RbcJ0/Ats+NcZqJ1znSBXcLFrDpKS7FeoB6\nneDlXPCM5Id+bfuGeh1iygUPCzYD4xlOeh0iXVCEgDGnXvtZLYSCBIx/WDbcw5QLaiBg/KVeH0kX\n1EPA6DrLhjtYLYTaCBjqtcWUC6olYK2zbPiOdEHlBKxp6rWkWxCFgLUrd72uek8C6YJqCVij6q9X\n8XfFkS6onIC1aKxXDekqEipxggQErDn11OvLdIkQNE7A2jLbMV+K11QB3xOwRpWafkkXcBUBa0jZ\nxUPpAq4lYK0ouHgoXcAdBKw5T06/pAu4j4A14fnFQ+kC7iZg+T1cL+kCniFgyT156ku6gCcJWCtu\nnX5JF/A8AcvsgcVD6QJKEbC07l48lC6gLAHL6daLzUsXUAMBS+7aekkXUA8BS+iOU1/SBdRGwLK5\n/NSXdAF1ErC0vp9+SRdQMwFL5arFQ+kC6idgeVxSL+kCohCwJL4/9SVdQCwCls2J6Zd0AREJWAan\nFw+lC4hLwMI7t3goXUB0ApbH/unXrF7SBUQkYLGdWDyc1ku6gLh+lB4A56kX0DIBi+rEqS/1AjKx\nhBjenumXk15APgIW0qHFQxMvICVLiPEcWjxULyArAQvs4/RLvYDE4i0h9n+Oyq/FIXnjrjT2Lx6q\nF5BbsID1fT/GafrxaLhl9a4E1AtglGoJMWW0RvtPfakX0IJgM7CP+sU1/ma3JIjcxvTLdnmgHdkC\nNi4hzm6Jbs/ioYkX0JQ8S4jLuVcaexYP1QtoTbAZ2Ov1Wm41HLZsrN6VzLvpl3oBDcq5W2+UYDvi\nx8VD9QLalGcJMSX1AnhHwOr18dSXegEtC3YOrE3L6Zft8gACVqmNxUMTL4DOEmKdNhYP1QtgIGDV\nmdZrNv1SL4CRgNVLvQA2CFhd3p36Ui+AGQGryLtTX+oFsCRgNZpOv9QLYJVt9LVYLh5KF8AGM7Aq\nqBfAUQJW3vLUl3oBfCRgFRmmX+oFsIeAFTZbPFQvgJ0ErKTZ4qF6AewnYFV4/XypF8AhAlbM3+nX\nL/UCOMzrwMqY1mu8UboA9jMDK2D1klHqBXCIgBX1Z/qlXgBHCdjTlouH6gVwgoA9arl4qF4A5whY\nIb9enXoBfEHAnjNbPFQvgG/YRv+Qab2kC+B7ZmBPUy+ASwjYs37JF8A1BOwJfy85r18AFxEwAEIS\nsNutXvYQgC8J2HOsHwJcSMCeYvoFcCkBu5ftGwA3ETAAQhKwG9m+AXAfAXuC9UOAywkYACEJ2F3+\nbt/4af4FcD0BAyAkAbuF7RsAdxOwe9m+AXCTeG9o2ffjS4Pncdi4C4BkggWs7/sxTtOPR8Mtq3c9\nxvYNgAekWkI08QJoR6qADWbTr/5fz43D9g2AOwVbQtw29Gk2D3t4WubqvQDPyDYDs4oI0IhgM7DX\n67XcajisGQ63l92IaPsGwGOCBaxbK9Nwi7kXQFOyLSFWwfYNgPsJ2GVs3wB4koABEJKAXcP2DYCH\nCRgAIQnYpWzfAHiKgF3A9g2A5wkYACEJ2Lds3wAoQsAACEnALmL7BsCzBOwrtm8AlCJgAIQkYOfZ\nvgFQkIABEJKAARCSgJ1k/RCgLAEDICQBO2Ocfnn5F0ApAvYVL/8CKEXAAAhJwA6zfQOgBgIGQEgC\ndpbtGwBFCdgxrt4LUAkBAyAkATvA9g2AeggYACEJ2HG2bwBUQMD2sn0DoCoCBkBIAraL7RsAtREw\nAEISsCNs3wCohoB9ZvsGQIUEDICQBOwD2zcA6iRgAIQkYPvYvgFQGQHbYvsGQLUEDICQBOwt2zcA\nahYvYP0fG5/w5HgAKCJYwPq+f/2xDNV22E6yfQOgSsECtm0I2yVfyvYNgMr9V3oAt5vNya4qHABl\n5Q/YiWLZvgFQv1RLiAC0I9gMbLp3Y5xaDTs7yg0KgAKSH/pPtM36IUAIlhABCEnA1pl+AVROwP4x\nrh8CUDkBAyAkAfvL9g2AQAQMgJAEbM70CyAEAfvN9g2AWAQMgJAErOts3wAISMAACEnAJrz5MkAc\nAubNlwFCEjAAQmo9YLZvAATVesAACErAuq6zfQMgnqYDZvsGQFxNBwyAuNoNmO0bAKG1GzAAQms+\nYLZvAMTUaMBs3wCIrtGAARBdiwGzfQMggRYDBkACAgZASM0FzPohQA7NBQyAHNoKmOkXQBptBQyA\nNAQMgJAaCpj1Q4BMGgoYAJk0FzDTL4AcWgnYuH4IQA6tBAyAZJoImO0bAPk0ETAA8mkoYKZfAJnk\nD9gl2zf6vqI9IAbzjsG8YzDvGMw7VQ3mnfwBAyCl7AH73+//tX4IkEz2gAGQVP96ZZ6a/D4B9qv0\nOACiqb8OyQMGQFaWEAEIScAACEnAAAhJwAAIScA+6/8oPZC/KhlMbY9MJYMZx1DD41PnYFb/+rBq\nH5mCg1k+FDU8OBv+Kz2AGIa9mn1fxabNSp5MwzDqeWTGMRQczOwnfxxGkSEtnycFH5/lYAo+jWff\npq6mR6b403j6rxd/Dn9kBvZZVd+22p5Gw29nVQ2poNfrVc9DMRtM2YHNBlP2ObP8NhV8Gtf2nCk9\nhGPMwPZymF6q7bez6W/TrKrkO1UVT+PlAGp4HPYQsM+KP59mIxn/rGFIlah/raM4T+P6FX8a1/Mk\n2UnAdqnkO+oYzTcqec54Gtcs1ndEwD6Y/rbYRfvu3mo4xzt+bDAzVQ3J0/idqr5NZQezfJJU9eCs\n8hsQACHZhQhASAIGQEgCBkBIAgZASAIGBVRyPTAITcAACEnAoCRTMThNwKAYl6KAbwgYlGHuBV8S\nMChjdqke4CgBg5I0DE6zBA9ASGZgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRg\nAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIT0fzgRde1M\npjFnAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo5()"]}, {"collapsed": false, "outputs": [], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Sparse Spikes Deconvolution\n", "---------------------------\n", "We now consider a convolution dictionary $ \\Phi $.\n", "Such a dictionary is used with sparse regulariz\n", "\n", "\n", "Second derivative of Gaussian kernel $g$ with a given variance $\\si^2$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = 6;\n", "g = @(x)(1-x.^2/sigma^2).*exp(-x.^2/(2*sigma^2));"]}, {"source": ["Create a matrix $\\Phi$ so that $\\Phi x = g \\star x$ with periodic\n", "boundary conditions."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["P = 1024;\n", "[Y,X] = meshgrid(1:P,1:P);\n", "Phi = normalize(g(mod(X-Y+P/2, P)-P/2));"]}, {"source": ["To improve the conditionning of the dictionary, we sub-sample its atoms,\n", "so that $ P = \\eta N > N $."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["eta = 2;\n", "N = P/eta;\n", "Phi = Phi(:,1:eta:end);"]}, {"source": ["Plot the correlation function associated to the filter.\n", "Can you determine the value of the coherence $\\mu(\\Phi)$?"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAPHklEQVR4nO3d7XKj\nOAKGUbSV+79l7Q8mNLGJ4w9AesU5tTWVMV0bDbZ5LKGkS611AoA0/2s9AAB4h4ABEEnAAIgkYABE\nEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAA\niCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEGjBgpZTWQwDgcF+tB7An6QK4jqFm\nYLXWWmvrUcBTSpl84oJPDDUDu2dORsfq5CVKx/qfDwwesKndc1BKaf709zAGw/htDHO5mgyph1PR\nyTB6GEOfw4j4aDXUEiIA1yFgAEQaMGA9zMQBOFoXC6/H6WRlGe593wNrPQ7YEnHxHHAGBkES7pRD\npwQMgEgCBkAkAQMgkoBBA259wecEDIBIAgbNdL9LGbomYABEEjAAIgkYAJEEDNpwAww+JGDQmC31\n8B4BAyCSgAEQScDgbNYMYRcCBkAkAQMgkoBBSzbTw9sEDIBIAgZAJAGDBqwcwucEDNqzsR7eIGAA\nRBIwACIJGACRBAxO5XYX7OWr9QBeVr4vAPVuI9eDQwAMJixgpZQlTuuvHx+CntVqWgbvsIQIQKSw\nGdgzys9Pszf/amYGsKmkLQWME7Ba63z2ly+Wx9sNCiDGzU2ZhiN50jgBc9+LFF6nsIuwgK1nVzdb\nNjYPQYpShA1eExawaStOyyO6BXAddiECEEnA4DwJ98UhhoABEEnAAIgkYNCe7UfwBgEDIJKAARBJ\nwOAkD7YgzofsUYSXCBic5MGNLvfA4A0CBqeSMdiLgEEvLCHCSwQMgEgCBkAkAQMgkoABEEnAAIgk\nYABEEjAAIgkYAJEEDE7i55RhXwIGQCQBAyCSgAEQScAAiCRg0AV/lwq8SsAAiCRgAEQSMAAifbUe\nwMvK94+D1q2bBvPRzUMAjCQsYKWUJU7rr28euT8EwGDCAvaG8vMX+AgbwKaS9uvORgvY/RKiYgE8\nY321jIjZOAF7vLoIwGDsQgQgUtgMrNZ6vwtxnm9tHoKueGHCjsICNm3FaXlEtwCuwxIiAJEEDIBI\nAgYdSdi6DL0QMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYnMGv2IDdCRgAkQQM\ngEgCBkAkAQMgkoBBL/yN4vASAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAGDk/g5ZdjX\nV+sBvKx8/70U9ef1oPz8+yqqqwXA0MICVkpZyrT+evpZrOIvXwIY3YBLiDdhA2BIYTOwN1haJEsp\n7pbRRtza1WgBu59+KRbAM+JuxAy4hAjAFYTNwGqt97sQl1mXu18A1xEWsGlrSXB5RL0ArsMSIgCR\nBAyASAIGQCQBAyCSgAEQScDgcAk/Egp5BAyASAIGQCQBAyCSgAEQScCgI34bGjxPwACIJGAARBIw\nACIJGACRBAyASAIGQCQBAyCSgAEQScDgDH5CGXYnYABEEjAAIgkYAJEEDIBIAgbPKmW0bwTRBAye\nMkel57SU0vXwYHcCBi+wGx768dV6AC8r3x8y69215MEh2EUpGga9CAtYKWWJ0/rr6bte8yM3h2B4\nFg+5oLCA/WmdMQAGNlTANidn5edHU20D2FTSJvJDBWyTYrGXWv/b6dfna2oeHrzt/qZM5+xChHjr\nS03CZQf2ETYDq7XebzWcFww3D8G+3pjlKAocJCxg01aclkd0i3N0uIrY23jgBJYQ4W89z6J6Hhsc\nSsAAiCRg8JqeF+t6HhvsTsDgHcct3L0RId3imgQMgrkBxpUJGDxlPcvpfMajalyEgMGbmnei+QCg\nLQEDIJKAwR+6nejMK5mdr2fCcQQMgEgCBkMxIeM6BAyASAIGhztoVtTtzTk4h4DB36zLQYcEDAZk\ncsYVCBi8w5wMmhMwACIJGDzSai3OGiD8ScAAiCRgMBr357gIAQMgkoBBMJMtrkzA4A+hkbANhOEJ\nGACRBAzeZ5YDDQkY/OrzPikcHEfAAIj01XoALyvfn2nr3b31svq4e38UgJGEBayUspRp/fVCt2BR\nSuoWSnjGaEuIpZTitgOXN3dLvRhb2AzsT/MM7Gaidv8HIJ3Paewu7tP/UAHbjJNi8QkvH65jfbWM\niNk4S4gRpxueoZrwjLAZWK31fhfivGC4eQiuzCYOxhYWsGkrTssjusWZanUjCloaZwkRgEsRMAAi\nCRgAkQQMtrm/BZ0TMAAiCRiksuuWixMwGJO8MTwBg4+4VQatCBgAkQQMjmUpDw4iYABEEjD4lckT\n9EzAAIgkYDAymyQZmIBBHlmCScCgWyoFjwkYAJEEDDY8OfuxTREaEjAAIgkYAJEEDIZlhZOxCRgA\nkQQMgEgCBpEsD4KAARBJwGCbKQ50TsDgU7/91LPfBQWHEjAYnI4yqryAlW8P/sCZ4wGgia/WA3hN\nKaV+35pYf73+A6cPCoAGwgL22Jy0m4bd/Ot98wCYAicAQwVsk2IBPGN9tYyIWd49sN/Mp3v9T3hP\nDy+fB5+7ehge9GCcGdjje2MADCYsYOtbXIoFcGVhAZu27mndPCJmsKjVkiPDGuceGDTh8xK0ImAA\nRBIwACIJGACRBAw2uLMF/RMwACIJGOR5dYJoJz1DEjAAIgkYHMi9NDiOgMEOrNHB+QQMgEgCBj+Y\nS0EKAQMgkoDB4GwkYVQCBv2yngkPCBgkkTRYCBgAkQQMgEgCBrde3fVglwQ0IWDwqfm+lLtTcDIB\ng310Pg/TV8YjYPCpztMFoxIw+MdiIAQRMLhlRgURBAz+kS4IImDwz46Lh8etQ6oszAQMgEhfrQfw\nsvL9ybbefRB9cAiAwYTNwEop9VvZWqN5cAgO1fOLzic6hhQWsMdMvACuI28J8bH7udfNIyIHsClu\n7Wq0gM19Wj8NikWoWrtelmQ866tlRMzGWUKMON0A7CVsBrbeoLF8WFh2dtiFyOe8diBFWMCmrTgt\nj+gWwHWMs4QIw/twmdwqO4MRMAAiCRgAkQQMgEgCBjuwfwjOJ2DwH3scIIuAARBJwACIJGBwlK5u\njHU1GNiFgAEQScAAiCRgkMRKICwEDIBIAgb/mN9AEAGD3Rzxo9B+vBp+I2BwIXLISAQMgEgCBkAk\nAQMgkoABEEnAYJrsboBAAgYZJBZuCBgAkQQMrsLvGWEwAgZAJAGDfaznN+5XwQkEDIBIX60H8LLy\n/eG23q3oPzgEwGDCAlZKWeK0/noxP7J5CB7zkoEsQy0hihbAdYTNwJ5xM/0qP++nixwXV4q5JttK\n2u6joQI2n/2bRCkWuWq1oZHzPPjo36ehlhAnuWJoXt2wFjYDq7XebzWc1wznx21EBLiIsIBNW2Wa\nH1EsgEsZbQkR3rDjan/CjQMYhIABEEnAAIgkYHAh7hQzEgEDIJKAARBJwACIJGAQwO58uCdgME0H\n7G6wXQKOJmAARBIwuBwLkoxBwGA3lg3hTAIGQCQBAyCSgEHv3LKCTQLG1ckDhBIwACIJGACRBAwy\n2KMPNwQMrsidPwYgYHAt80zOfI4BCBgAkQQMdp6OWJ2DcwgYXJHKMoCv1gOAQywXaDd73jafQyeQ\nbpmBMSDTix05mXRLwBjNSxdcV+ffrM+Ms0Sf8pYQy/ebqf6ytFFK+e0QY1tfZ+eXQClTKRbBXnaz\neDifxslyIp0JC9g6TvehKj4o4iL7hFpfm1S9+ufhHEMtIdZazb1YW+ZhPG9zsuVM0qGwGdgbbqZl\nCjcqa1zwobhFrPEDpljXsflUz8tfD+6Edf4COXn57sHngD/PJOnWV8uImA21hMhlPX6vueBuSrhA\nwSNhAau1lm/r3RxtR0UP/qzUOS+T6Fj+uQzrThhdyVtCvF8SvHnEmuHVuJ7CNYXNwGDTnx9afKrZ\ni0kY/RAwsn1yJXUV/oSzR3MCRjDX0N2ZqhJEwIj3/DVX8BafhErk6ISAAS9zJ4weCBipLvWrN074\nz1Qj4ggYV/Hbb5fgPSZhNCdgRLrO9Gv+Lz2oE/JDNAED3mQSRlt5v4kDPpl+nfC7aPf9/z9zlnmF\nGS0jMQPjoswbdmESRkMCRpjr3P06kwKRSMDgovb9ECCBnE/ASPLh9Mu87QjOKq0IGPApd8JoQsCI\nsdfdr+U6a+qw5mwQR8DIYO/GQfaaNpmEcT4BA/bh4wUnEzACmH4dat9pUynmYZxEwOjdvldDFYRh\nCBgZ/NBShOVpcoY5gYDRNYuH59jxDGsYpxEw+qVeRzvo3HrKOIeA0aNlI4BLYahlV715GMcRMPql\nXtGsJXI0AaMvp31mV8czaRhHELCjlA7esj2MYXp6GOt01bp/YG6G0eTcZD0jn5ufx9+WE3s4Gz2M\nYTKMdwkYXTg0XTS3Xk5Mu0jSr6/WA+C6Wl3IBHLttLNR679n/PsLzwQfMQM7w4Mr9e6Hfvr1AnHo\nkO7+WF0fWv7373D9sdbEqO6f5fXrYZ23hy+npw49/cfavEeeHwYPlDr0NcNiRcfm56Z+f8FljXwJ\nStd5Hy6xhLj+ix5utvYeeujPYdz8xVQnDOmTQ++Ndv7XrbdBvfsCpqmPV/uDQ3++ER4cev77rjX8\nD+nc4DMwAEblHhgAkQQMgEgCBkCkS2ziOFkp/91ZXH6s/cwbjffftIdhNBnDejA9PCNT07PRwzNy\n84seaq2tXhjz9/UeefDWaPuefZKA7Wn9/lxeGTdfn2B5Rc4XiObDmH6+K5pcqqYOnpGGw1g/Cw2f\nkZvz0OoZWV+1W71Hbr7pdPoz8vhi1fDN8hJLiHuqtTZ/ppsPYNbJMDp5780X69ajuG1G25H0MIwr\n6+Fi9TkzsDH1cIHo4ZLdifXn/bZjmDwv0zTdTXqaj4T3CNho+nln3qytn2/+1ut/NtHDc9GVtiHv\nZHFsuf+3vhHIqywhDqj5FbOTN2T9NrU7J52cCnozt7P5WzVd+4Wm8TTc83a/y6vJMO6/qV2IU+v9\nXZ08IzeTnh5ORSfDsAvxDQIGQCRLiABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIw\nACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgk\nYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACL9H73SW/8/I2QzAAAAAElFTkSu\nQmCC\n", "output_type": "display_data"}], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["c = Phi'*Phi;\n", "c = abs(c(:,end/2));\n", "clf;\n", "h = plot(c(end/2-50:end/2+50), '.-'); set(h, 'LineWidth', 1);\n", "axis tight;"]}, {"source": ["Create a data a sparse $x_0$ with two diracs of opposite signes with spacing $d$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["twosparse = @(d)circshift([1; zeros(d,1); -1; zeros(N-d-2,1)], round(N/2-d/2));"]}, {"source": ["Display $x_0$ and $\\Phi x_0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAN10lEQVR4nO3d0Xaj\ntgKGUdQ17//K6gVrKDXEcRyD9Et7X/SkzpyxijGfJUgotdYFANL803oAAPAOAQMgkoABEEnAAIgk\nYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQ\nScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIg0YMBKKa2HAMDl/rQewCdJF8A8hpqB\n1Vprra1HAcAdhpqBHZmTcZ26LHYvBtb/fGDwgC2tX4NSSvOdwBiuGsAP/87mG6GHMTQfgDG8OIaI\nT/9DLSECMA8BAyDSgAFrPjEH4Abt12Ev1cNCM8MqZbF3MaiIg+eAMzAAZiBgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAwS+U0noEMC8BAyCSgAEQ\nScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABE+tN6AD9W/v761FrrV986/S4A\nIwkLWCllK9P+641uAUxitCXEUkpxhwuACYTNwL61zsAeJmrHPwDAg7hP/0MF7DROigXwiv3RMiJm\n4ywhRmxuAD4lbAZWaz1ehbguGJ5+C4BRhQVsOYvT9ohuAcxjnCVEAKYiYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgA\nkQQMgEgCBkAkAQMgkoABEEnAAIgkYPCWUpZaWw8CpiZgAEQSMAAiCRgAkf60HsAnlVLWL6qTEwCj\nGydgpZStW/uvARiSJUQAIo0zAztVl2X5u64Il2i+g9XafgwMJ2IJa/CAFefDuFTzvWtLV/ORMJZS\nSv+7lCVEACKNMwOrtboKEWAe4wRs0S2AmVhCBCCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAk\nAQMgkoABEEnAIFwpfhU9cxIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJH+tB7Aj5VS1i/q4Tbq27dOvwvASMICVkrZyrT/\neqNbAJMYbQmxlLKfhwEwqrAZ2LfWGdjDRO34BwB4EPfpv+uA/bQ9p39AsQBesT9aRsSs64D9qD2n\np8QAGFXXATuqtR6vQlzTdfotAEYVFrDlLE7bI7oFMI/RrkIEYBICBkAkAQMgkoABEEnAAIgkYJDM\nlbdMTMAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZA\nJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAJ2rVJK6yEYQxcDMIZOBmAMXY3hlwQMgEgDBmyA\njxUAfOtP6wF8knQBzKPUWluP4cNK+e8/StIYXl0WezlX6L8OQ83Ajvp/AeD37OXMKTVgD1MroQKY\nTWrAFAtgcgNehQjADAa8iAOAGaQuIfZsf35u/XywPXLPx4XT6zBvHsk2hvu3xvHvf+WR28awPnjD\nC/HtS99wDE9GddEwZn5HHJ+0yc5wBQG7xH4neHjzXP1GffjXhx30hpEcf3Th/q2xHRfW9+rDM94/\nhv0jy427xP742GQjLIftcP/OsN8hm7wjlsObosnx4cnud+cx6rOcA7tEKaXJj6DVWpvvfMcx3Lw1\nmm+Br8ZgOyy3b4QeDsfHMTQ5PrQ6KF3KDOwSx0/fM2uyNXrY+A9jOK4g3TCA257rxTF4ayyNNsL9\nu98NzMA+b+Z35tH9W2P9pNn2VTiOocl4epuR3zyY9WC9/+f9jmNotSfc/6Q3ELAPG+wDzi+12ho9\nvF0fznPc/Ow97IfHM7I3D6D+tbTbJR7G0OR16WFnuMjUE/mLNL/Cp89rru7ZGl8tWPU2hsmvQrz5\nrTHzO+L0KZofoz5FwACIZAkRgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACR\nBAyASAIGQCQBAyCSgMGbSllKaT0ImNif1gO4VnGA4Sp1/R/7GKOqtbYewjcGD9jSwWtQSmk7huYD\nGHIMW7Ze/zvH2wihY2g+gIgxRHwys4QIQCQBAyCSgAEQqf067KV6WGhmSLtzYE3HAdeIOHiagQEQ\nScAAiJR3Gf12cedxevvkW3CRUqwiQhthAdsvy54u0a6PRKzeAvAbQy0hihb3SPgRTxhf2AzsW8ef\nHn94ROQATkX89o290QK2LSE+PALAc/ujZUTMxllCjNjcAHxK2Ays1nq81HC9ZOP0WwCMKixgy1mc\ntkd0C2Ae4ywhAjAVAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAGDN/nFL9CWgAEQScAA\niCRgAETK+230T+6Z4nYqAPMIC9h666/j15vtjswaBjC2oZYQRQtgHmEzsFc8TL+2dcWVyAGcejha\n9m+ogK1b/yFRigXwiicf/fs01BLiIlcA0wibgdVaj5carmuG6+MuRASYRFjAlrMyrY8oFsBURltC\nBGASAga/lXC2GwYkYABEyjsHBpNbJ3zHc77bRNDpYCZhBgZJLFfCRsAAiCRgMBqzNCYhYPAz8gCd\nEDCIsW+njoKAARBJwACIJGAwAiuKTEjAAIiUF7Dy15M/cOd4AGgi7FdJrbf+On69PdJiUAA0kDcD\ne6LW6q5gAJMIm4G94WFapnDMoBS/0pcfi1vEGj9gisUY0o4t5NkfLSNiNtQSIszARzJYhc3Aaq3b\n54InV3MAMLywgC1nS4IPj4gZwAwsIcL7Ovyw1OGQ4CICBuNQL6YiYJBEomAjYABEEjAYU8KP8cCv\nCBgAkQQM3nHzuSjTKTgSMAAiCRgAkQQMUn21ruhSeyYhYBBDmWBPwKB3b1/B4dIPxiZgAETK+230\nTxzvtAJzqtX0i/GNE7D9XcHcIYxhbB16skc/aVUpzpwxrHEC9oXqcyi/8fbRv58dr5+RECXgg49z\nYPDMe0f/jzfjjY6aeDG84Wdg1hJ50xahn67C7evVdu+z7/O2Ukr/kzAzMDj3+6O/fsClxpmB1Vpd\nhchnvX4t33GWZh+Eq40TsEW3uIxr+aBDlhDhB76dkLnkD24jYABEEjB4xsohdEvA4PNkD24gYABE\nEjB4SVdXZ5jhwSJg8EFdRQ6GJ2AARBIwACIJGACRBAy+4YoJ6JOAwYe1Cp7QMhsBg1+RDWglL2Dl\nryd/4M7xANBE2O1USvnvDsv7r7dHWgwKgAbCAvbc2rOHjD38q3uG8bZtVzrdiXx8Il3cHGCogJ1S\nLIBX7I+WETHrOmAmTwB8peuAKRYAX+k6YEe11m1a9uRqDgCGFxaw5Wxa9vCImPFxtbpGA7qT93Ng\n0DMfn+A2AgZAJAEDIJKAQTBn5piZgAEQScAAiCRg8FuuPIQmBAxeslXqSa6UDO6U94PM0Ip0QVfM\nwACIJGCQzZX0TEvAAIiUdw7s+NvoX/kWAIMJC9j+zimnd1FZH3GDFYDhDbWEKFoA8wibgb3iYfpV\n/n+OW+QATpW0K4K6DthP27P+efe3BHjDk4/+feo6YG+0R64AJtF1wI5qrcdLDdc1w/VxFyICTCIs\nYMtZmdZHFAtgKkNdhQjAPAQMgEgCBkAkAYNIzvmCgAEQScAAiCRgAEQSMBiBU2JMSMAAiCRgAEQS\nMAAiCRgAkQQMgEh5v43+yQ1T3EsFYB5hAVtv/XX8erM+cvotAEYy1BKiaAHMI2wG9q1tFfGrR0QO\n4NTx+Nm5rgP2Rnu2JcQf/b8A2B8tI2LWdcB+1B7nvQCm0nXAjmqtx0sN13SdfguAUYUFbDmL0/aI\nbgHMY6irEAGYh4ABEEnALtf8Yp7mAzCGTgZgDJ0MwBg+RcAAiCRgAEQSMEjlqlsmN/gP/w6wyAtP\nre9f+zmf138dBg8YjG39hOZNzJzyfpAZ2EgXM3MODIBIAgZAJAEDIJJzYJ+33djl9BfnPzzy8af+\n9hmbj+G2mwb08EIsX/9XD/9CHO/nd/9G2J6i4d7Ywwvx+k6YdU8PAfuk/Tt2f3+y9evjI1eMYdtT\nT5+x+RiW/79VbjhqtN0Ircaw39StXoiHLdDkhdgfu5tshG//qy8dw0+PSPe8NT7IEuInrbclazuA\nhs/ezxh6eO+th+wextDDpmg+hjk1PyJdzQxsQM2PF80P3D3Yf/BvOIC2Y+jBwyyn7TD4LAEbSifv\n1YcF95utz7v/5/2avwT9aJjPHhbEtjN/+1OAfIolxNGYe9W/lkZbo4eNQCfWcPpAc5Gp1xYu0uri\nt+NFX/eP4ZVnnOoqxIaXe/XwQjzMe6bdCG3HMPBViAIGQCRLiABEEjAAIgkYAJEEDIBIAgZAJAED\nIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACL9C7NE\nXLx/e/7PAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x0 = twosparse(50);\n", "clf; \n", "subplot(2,1,1);\n", "h = plot(x0, 'r'); axis tight;\n", "subplot(2,1,2);\n", "h = plot(Phi*x0, 'b'); axis tight;\n", "set(h, 'LineWidth', 2);"]}, {"source": ["__Exercise 6__\n", "\n", "Plot the evolution of the criteria F, ERC and Coh as a function of $d$.\n", "Do the same plot for other signs patterns for $x_0$.\n", "Do the same plot for a Dirac comb with a varying spacing $d$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAATXklEQVR4nO3d63Kb\nyBqGUZjKfSe5cuYHsYyFDiBO/XavVbt2uSaOhbHMk69py/0wDB0ApPnv6gMAgE8IGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEE\nDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScCS9X3X91cfBMA1\nKg9Y38L1vYXPEWCm8oDVTLeAtgkYAJEELJPxC2iegFVBz4D2/Lr6AFjvlqthkC44QRPbwbpuGIar\nD2EdAQN4L+7ivlZipC0hppmOX7f/B2iPgAEQScCi3I1fD/8IoA3ugQEU6u6+VPX34dYSsBwvxi+g\nUqL1goDls5keTnfc95xgLSdgIYxf0KTbKqJRbE7AKtL38gbnOO1bTbdesAsxgfELYEbAiqdeAI8I\nWBW0DWpk/fA1ASub8QvgiQYC1tQW86Y+WaBtlQcse2wxfgE8V3nAAKhV3s+BvfixvunrhsXf/Fw7\nfnk9DqjO/LUQn706Yps/75wXsO7rK9T3/fxL1dQX7zE/zgwVeXGVu0VrejF8eGGsVd4S4uuvTd/3\nD36vaNxo4u4XsMxdsdqpVxc6gb341dfPhrM252vgIP2fo/5ZPPx5fO3qJte3gx46TmTAnn0VX8Qp\nqVsfj19ug0F1liwh3mlnFTHs83yx1Pv4jxLX4rYc8/h3gz5ZSHBVEh4uJs0vdLvcA0vMXtgENt2E\nc/cFe/hHeXYprn0cUIslm6srufqtFxaw7uVA3dRXDqje643Wze7duMnbhVizxAVPgIsIWF2UD2hG\n5QFL2pNn/AJYo/KAtct+eqB2eZs46mT8Aj61/CUTp+9cwb4PAQOIt/DnnSt71URLiAXYd/wKf0YC\no+lr9d698fEHrGznvQkMYL3jbjM/6srbaanNl0xsJmANvjhFg58yVOR2K2t8Y2zYs5tbH7xkYpe/\nithMwAB2dNF1f0tvxv7VdA9MwIoR/kwCLrTw99FX9qqJAlYjv1cFqrDwhV6Xv2Ti2w+VxS7EqskY\nUC8BAyCSgF3NkATwkQbugbV5Q6jNzxoO08KPVcVpIGAA29S08aEmlhBr55+NQKUEDIBIAlYGCxQA\nKwlYvUQRqJqAARCppYAVuJ3hnEMq8BMH2KylgAFQEQEDIJKAVc0+DqBeAtYGt8GA6ghYAcxJAOsJ\nGACR2ghYyyNOy587ULU2AgZAdQTsOidvrLCPA6iLgAEQScAAiCRgDbCPA6iRgAEQqbGAFbiR4czx\nqMBPH+BTjQUMgFoIGACRBKwN9nEA1RGwi7gdBbCNgDVGOIFaCBgAkZoJmJtAAHVpJmBlOjOrEg7U\nRcDa4zYYUAUBAyCSgAEQScCucNUinttgQEUEDIBI7QXMFobOSQBq0F7AAKiCgAEQScCuc8mWCvs4\ngFr8uvoA9tR/3doZXKYBalfbBDYMwzAMvU0KbzlFQLiqAmbwAmhHVUuI3WQV8eF/KaJvRh+APdQW\nsHEI+xGt6VgmHl3XDYPzAFSgniVE970AmlLPBDbdu+Fm2CJ9b1c9kKuegHVZ3Qo6VIAi1bOECEBT\nmgyYu2UA+ZoMGBYwgXwCdi7DH8BOBAyASALWNhMhEEvAAIgkYABEErAr2AQIsFljAVOOG6cCCNdY\nwACohYA1z0ZEIJOAnUgqAPYjYABEEjAAIrUaMKt5nY2IQLZWA3Yh2QDYg4ABEEnAsKAKRBIwACIJ\n2FlMOQC7EjAAIglY22yJBGK1FzCXbIAqtBcwHnKLDkgjYOcy/wHsRMAAiCRgAEQSsOZZ1QQyNRyw\nM7ct2CIBsLeGAwZAMgHjizERiCJgAEQSMAAiCdiJ7PcD2I+AoaxAJAEDIFKTATt54Aja3Rd0qEDz\nmgwYAPkEDIBIAgZAJAGj6zobEYE8AnYWhQDYVdsBs+kOIFbbAWNO1IEQAgZAJAE7mIEG4BgCBkAk\nAeOLfZJAFAEDIFKrATNtvOC+HZCg1YCdTC8B9iZgAEQSMAAi/br6AFbrv+7QDLN1uX5y82b+pxeI\nu5k0DHnHDLQqL2DdV5z6vp9XqohuAXC8vCXE14nq+743QwA0IHIC656MX92j4eyuZ/d/q+9tEXzA\naQGKlxewMUgv6rXkPwKQLm8JsXvSpHJXDhUU4ABhE9hYqbuNiOOC4TAMLzYoAlCZsIC9XiTUrR3Y\nSQ+EiFxCBAABO0z6HJN+/EDtGg6Y9UaAZA0HDIBkAgZAJAFjxuIqkEDADiYGAMcQMAAiCRjP2UkP\nFEzAjrlMu/QDHEzAAIgkYABEEjAesXkSKJ6AARBJwI5UwRxjNwpQKgEDIFLbAatgQgJoVdsBAyCW\ngB2gjvtGxlOgbAIGQCQBAyCSgPFOHSuiQHUEDIBIAtZ13TFDhk0QAEcSMAAiCRjPGSKBggnY3mx5\nADiFgLGAKgPlETAAIgkYAJGaD5h9CgCZmg/YQarpYjWfCFAdAQMgkoABEEnAWMZOeqAwArYrV3mA\nswjYF+0BiCJgAEQSMN6xkx4okoABEEnADmBkATiegLGYfS5ASQRsP67vACcSMAAiCZhbVgs4RUB5\nBGxvrvUApxCwnbgBBnAuAWMNnQaKIWAARBKwXbkBBnAWAZuwPgaQQ8D20EL5DJdAYQQMgEgCBkAk\nAdtPI4tsLayXAgl+XX0Ae+q/rq3DmS1xQQe4Qm0T2DAMwzD0ogJQu6omsIeD113MTh3OKjMMxk2g\nHFUFbNT3/bRSi4q1/dLcVBf7vq3PFyhSVQEbhy0zFkALKrwHdurjWVIDuEg9E9g4fl2zEbEdboMB\nxagnYIoF0JTalhCv0WA7zWHA1QTsp1XXZRdxgOsIGACRBIyVGlwvBYokYJs1e0G3ggpcSsAAiCRg\nnzJ/AFxKwFiv2VVToCQCxgbGUOA6Avbls6nCLAJwEQH7iMkD4GoCxkeMnsDVBAyASAK2gSmks5oK\nXEbA1nPJBiiAgM3oE0ACAeNTVlCBSwnYp1y+b8yswBUEDIBIAraSaQOgDALGBtZRgesIGHswmAKn\nE7CJ5fOEyQPgagK2hjkDoBgCxjaGUeAiAgZAJAFbz8zxkPVV4FwCBkAkAVvMhAFQktoD9rvr/6wP\nj1atYk0VuELtAeu6rvuoYXxA+IETNRGwPZk2AMpQe8D+7vRxzBYAhak9YF+sIh7LYAqcrpWAcRKj\nKnCW+gM2/Pk3HOwwhJkzAIpRf8DWkSiAEAK2gGWxJbQfOFcTAdtzFRGAMjQRME5lYAVOIWCLWSID\nKEkrAft8FdE8AVCkVgLGGQypwIlaCtjfNZdXg9cWzh5wvIYCdhsPPtmLaLYAKExDAQOgJo0FbNUq\nYmcpbD2jKnCWtgK2aRURgJK0FbAPmSo+YHgFDtZewNauIgJQpOYC9n4V8fs9zBAfMbACp2guYADU\nocmAHbaK2P/pbQ/5ZoQFjtRiwNbtRVy8IHb7aBoGcILUgPWP/nXfT5x/SD+ORMMADpYXsNd9Gr68\n+Sh7ryIq1g/2cQDHywvY6z4tHL+O+4lmv/35h6tHYaBieQF7bczbtGH9T2s/3JL3urXqVq+7/w7A\n7qoK2MPJbPjp+w/+7jMqzes1zZiGARyknoCtna4OvU1zN4q1yG0w4GA1BGxM17hyOHq/iWOvh36y\neNi5GQZwsNSATRN1e3vpFsSbZ6uI40d493GWl6nphtnHARwjNWC7eFWoNTPcswVDN8MAjtN0wLZ4\nsXg41fTNMLfBgCM1H7CP9iKuemc3wwCO0HrANg4JawesRhvmNhhwgNYD1nWrX1Zq4eLhlJthALsT\nsG+HpqXpm2EABxCwdauIH4xf3w/U4M0w+ziAwwhY13VLVxG31OvuLzbUsJHbYMDeBOyHM7vSXMMK\n1/f//geEELCuW7bQtX382uWvs7+7bmkYhBCwLy9XEfedltpaSCz8NtjDXBnFIIGA3Xsdlb3mp7Ya\nNiotCXeVGoZ//5u+A1AwAfvnxZyw1+Lh24/PSeYD1vTLr2EQQsAmHq0iHleXhm6G3ZJQQg/m6Zr/\n40XDIIGAPfAwWkf0pqGFxBIa9nDN8BkNg+IJ2Lf51ezoxcPuZ8PG/x30QE17vWb4jIZB2QTsp5Wv\ni3iEakt21RC2ZM3wGQ2DgvUrfn9xoL5f9wn2fdf9vr9OnXOz6kWxqrpbdsvACU+8D6autx+n6u8X\nyCJgs7/yMyTnx6Pykp0TsPm0tPHh9mohsB8Bm/+VH0PYtc2oM2ZHN+y42BjFoCQCNv8r3wErJxK1\nley4hh3dGA2DYgjYw7/Vdb/77m+J5+a0Lf7HOigD59RFw6AMAvbwb/17o/BzM41ZcMP2OstndkXD\noAC20T9Qwg/dLjGNVt62+33P8slFsb0eCiBgbxR+dcobvB7aeJYvmYc0DK4mYI8FLQsFvx7VLmf5\nwtU8DYNLCdhTKQuJXR0N++wsX34vSsPgOgL2SlDDbvIadrP2LF9er/lDBz1RIJ+ALVX4pSl4Q8dn\n7SmkXvMDKPyJAhURsDcuvzYuV0PDFl79i6rXSMPgdAL2XtBCYg2bEt+e5QLrNdIwOJeALZLYsNQh\n7LVi6zXSMDiRgK1W/nUpvmHPTnHh9RppGJxFwJYq9oL5WljDbuaX/oh6jTQMTiFgKyQuJHZZDXtW\npqB6jTQMjidgHyr/opS6oWP+z4S4eo00DA4mYOtkXZRSb4bd9H1qvUZZTxdII2CrZV1FIxv28BRn\nnfcbDYPDCNgngm6GTQU3LLReIw2DYwjYVuVfkeI3dETXa6RhcAAB+1DWFSm4YRXUa5T1jIEEv64+\ngGDDkHQhGv4MSemahTZ1U+XU9BnT9/W0GS7SD1V/F/X94Z/g7YoUcSJvVSi5B29DW/LBvxe9rxJK\nImDbH+L77YhzWfJks3ZGLOrgV4h70kCRBGyXR/l+u/zTuSQS54dhflTzY3h95MuPuYgJ7271ufzn\nDZRHwPZ6oH9vpJzOcmadJela8rfmH2HjPb/DM6ZhsI2A7fhY/94IPaPnJ+2zdC35OKu82KJpFIOS\nCdiOj/X9dh0ndXkY1l7o90rX2w/72UOcnTENg48I2L4P9/12red1yy20h3+32I0YMgaFE7DdH/H7\n7apPbdctHtEe/ghasd2aOvuwNQzWELAjHvT77arP7rd/n/LvpT079GB2d+Uo1sgTCD4iYMc99Pfb\ndZzj1S87Mu/Z3x8nIuu0nJoxoxgsIGCHPvr327mneePLZY2f+OsPEnRyLstY0DmCswjY0Qfw/Xbc\nmX5Yne2fRXrMTr0xZhSD5yoM2DRalwesC2zYQd1a/lhHP+gujGJwueuv7zvqv77PiwpYF3L9edaS\n0w74RczKPWmnZcwoBjNFXN/3dTeBzd9h/inP323395n8yeGPtfF9um7o7n991fXHMz2ka4/n7n1O\nfY49PDuFPX8S32f+bmV93S96n/Lr4BdanqTAZ0LfP10tLPBou+cHXKCjfvXawy9MykmBvdU/gRX1\nCRaylnjmXa5dzA+4zKM9b3/H5Qu+UICyru+7KDlgowszllKCh1IO/sptij8etcizA/sp8fq+UfkB\n6y5qWB37AFLGx4tfhurHA5d3dmAPhV7f91JswLpzG1ZHuqZk7PlDKhmtKPf6vouSA9ad0rD60nUn\nYl3xmtcyVjJqV/T1fbvCA9YdHJjq63UjYy8fOPbHxeGl0q/vG5UfsO6Y1bB20jUVsa545S9Fe7vh\nvrSTBS8FXN+3iAjYaMeLb5v1usnN2Oi8XzejZ4SLub5/Jihgo40/3tN4uqYiMta9+5FnMYMXwq7v\na8UF7GZtyaTrmZSSdUExG5V5EmlJ6vV9odyA3Sy5/qrXW0EZG4XF7Kbkc0p14q/vr1UQsJF/E6/y\n4usesV/xzr6vrLgpfp+97mLZp7iaq8S+Ik5LwCFuEfE1WMXP9izx9ut+7UD28RfxqNcIXuZp+ba/\nmvClz936rhK7iDgtAYe4RcTX4DOFvC5wmRZ+3at9GfffF39iw98THmO3J33FV4ktIk5LwCFuEfE1\n2KL6T/Aza09LtSVb65TynZG3koV8w0ZcWwIOcYuHv3QOnvj4e6G0p9kV39T7xa/1wpWj+DpUHjDg\nTIf8i3FbGuXwc8XXQcAAiPTf1QcAAJ8QMAAiCRgAkX5dfQB8brrH0r3M0W3v7+3kODPdo9PStX1m\n7p4eni2j+Xko/AkjYNkKfEpdZfqdNv0RlogfZznO/CdJWj4bU9N0ebbc3E7L7TyUfEIsIWbr+97P\nuo2GYSj5O+0q89PiOdOVfVG+0MPTUvITRsCyjZenYp9eFMhz5sa89dDdaSn5CWMJMZjvPdbynBnd\nrRwymp+Wwk+RCSxVmf8gomSeM1OFX5qvMj0t5T9hTNDBbJ2aswvxIadl6u66PF0fc1puIp4wAgZA\nJEuIAEQSMAAiCRgAkQQMgEgCBlcqf6cyFEvAAIjklTjgAgYv2E7A4Gx3L5Z/7cFALkuIAEQSMAAi\nWUKEsxX7yykgi9dCBCCSJUQAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEOl/yAk2NzHMix4AAAAA\nSUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 30, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo6()"]}, {"collapsed": false, "outputs": [], "prompt_number": 31, "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"}]}}}