{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Forward-Backward 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 Forward-Backward (FB) algorithm to\n", "minimize the sum of a smooth and a simple function. It shows an\n", "application to sparse-spikes deconvolution."], "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_4_fb')"]}, {"source": ["Bibliography\n", "------------\n", "Excellent review papers on proximal splitting algorithms include:\n", "\n", "\n", "Amir Beck and Marc Teboulle,\n", "_Gradient-Based Algorithms with Applications to Signal Recovery Problems_,\n", "in \"Convex Optimization in Signal Processing and Communications\".\n", "Editors: Yonina Eldar and Daniel Palomar.  Cambridge university press.\n", "\n", "\n", "P. L. Combettes and J.-C. Pesquet,\n", "_Proximal splitting methods in signal processing_,\n", "in: Fixed-Point Algorithms for Inverse Problems in Science and Engineering,\n", "(H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and H. Wolkowicz, Editors),\n", "pp. 185-212. Springer, New York, 2011.\n", "\n", "Forward-Backward Algorithm\n", "--------------------------\n", "We consider the problem of minimizing the sum of two functions\n", "$$ E^\\star \\umin{x \\in \\RR^N} E(x) = f(x) + g(x). $$\n", "\n", "\n", "We assume that $f$ is a $C^1$ function with $L$-Lipschitz gradient,\n", "which means that\n", "$$ \\forall (x,y) \\in \\RR^N \\times \\RR^N, \\quad\n", "      \\norm{\\nabla f(x) - \\nabla f(y)} \\leq L \\norm{x-y}.  $$\n", "\n", "\n", "We also assume that $g$ is \"simple\", in the sense that one can compute in closed form\n", "the so-called proximal mapping, which is defined as\n", "$$ \\text{prox}_{\\ga g}(x) = \\uargmin{y \\in \\RR^N} \\frac{1}{2}\\norm{x-y}^2 + \\ga g(y). $$\n", "for any $\\ga > 0$.\n", "\n", "\n", "The FB algorithm reads, after initilizing $x^{(0)} \\in \\RR^N$,\n", "$$ x^{(\\ell+1)} = \\text{prox}_{\\ga g}\\pa{ x^{(\\ell)} - \\ga \\nabla f( x^{(\\ell)} )  }. $$\n", "\n", "\n", "If $0 < \\ga < 2/L$, then this scheme converges to a minimizer of\n", "$f+g$.\n", "\n", "Sparse Regularization of Inverse Problems\n", "-----------------------------------------\n", "We consider a linear inverse problem\n", "$$ y = \\Phi x_0 + w  \\in \\RR^P$$\n", "where $x_0 \\in \\RR^N$ is the (unknown) signal to recover, $w \\in\n", "\\RR^P$ is a noise vector, and $\\Phi \\in \\RR^{P \\times N}$ models the\n", "acquisition device.\n", "\n", "\n", "To recover an approximation of the signal $x_0$, we use the Basis\n", "Pursuit denoising method, that use the $\\ell^1$ norm as a sparsity\n", "enforcing penalty\n", "$$ \\umin{x \\in \\RR^N} \\frac{1}{2} \\norm{y-\\Phi x}^2 + \\la \\norm{x}_1, $$\n", "where the $\\ell^1$ norm is defined as\n", "$$ \\norm{x}_1 = \\sum_i \\abs{x_i}. $$\n", "\n", "\n", "The parameter $\\la$ should be set in accordance to the noise level\n", "$\\norm{w}$.\n", "\n", "\n", "This minimization problem can be cast in the form of minimizing $f+g$\n", "where\n", "$$ f(x) = \\frac{1}{2} \\norm{y-\\Phi x}^2\n", "\\qandq g(x) = \\la \\norm{x}_1. $$\n", "\n", "\n", "$f$ is a smooth function, which satisfies\n", "$$ \\nabla f(x) = \\Phi^* (\\Phi x - y),  $$\n", "it has a $L$-Lipschitz gradient with\n", "$$ L = \\norm{ \\Phi^* \\Phi }. $$\n", "\n", "\n", "The $\\ell^1$-norm is \"simple\", because its proximal operator is a soft\n", "thresholding:\n", "$$ \\text{prox}_{\\ga g}(x)_k = \\max\\pa{ 0, 1 - \\frac{\\la \\ga}{\\abs{x_k}} } x_k. $$\n", "\n", "Signal Deconvolution Problem on Synthetic Sparse Data\n", "-----------------------------------------------------\n", "A simple linearized model of seismic acquisition consider a linear\n", "operator which is a filtering\n", "$$ \\Phi x = \\phi \\star x $$\n", "\n", "\n", "Dimension of the problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["N = 1024;"]}, {"source": ["We load a seismic filter $\\phi$, which is a second derivative of a\n", "Gaussian.\n", "\n", "\n", "Width of the filter."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["s = 5;"]}, {"source": ["Second derivative of Gaussian."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["t = (-N/2:N/2-1)';\n", "h = (1-t.^2/s^2).*exp( -(t.^2)/(2*s^2) );\n", "h = h-mean(h);"]}, {"source": ["Define the operator $\\Phi$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["h1 = fftshift(h); % Recenter the filter for fft use.\n", "Phi = @(u)real(ifft(fft(h1).*fft(u)));"]}, {"source": ["Display the filter and its Fourier transform.\n", "\n", "Fourier transform (normalized)"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["hf = real(fftshift(fft(h1))) / sqrt(N);"]}, {"source": ["display"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAASt0lEQVR4nO3d25Kj\nOJQFUJjo//9lzQNTDAkYY4xBW6wVHR2ZpC8q2dbmCBn6UkoHAGn+5+4GAMARAgyASAIMgEgCDIBI\nAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAIMgEgCDIBIAgyASAKMx+n/\nGraMf5r+euxhv2nYzu3fPMvvHgou9t/dDYAblFI2fj2g7/vxQaY/f9Oq041ZVUqZ/fzrp4ZfEGAw\nj59xy3SUH//0dqyf3njnlvHX6ZO+eq6xtavNe/Xr8omm/yIZRiIBxhOtDu7jlmmQzEb8V6P8Mu2m\nPwx3X83I6Q2W99r+V0xDaJnB31eEUDkBxhPtH9D3HyJ6GxgbDzW9/UcHpYQTTybA4KX9xdDbIFlW\nSKt/3b7Zp08KbRNgsGL1GNjsr8u7zG48WygxC6RXj7+zApstmHx1xG7ZeGiGPTj42J7S53fl0emP\nrJIjlO+BwcfuHe5nKxXhsex5ARBJBQZApAYDzOwKwBM0tQpRdAE8R1MV2PRkBwC0rakKbElNRvVK\n13mXUqP664HGA6xLeA26nC/iaOe5dp4d+HYR/RnRyC6qnXc34b2mphAhTcBABtUSYABEajDAIspz\nGJTSJUzVQI0aDLBEKaGrnc8U0Z8Rjexy2hlBgAEQSYABEEmAARBJgME9+r5zNAS+IcAAiCTAAIgk\nwACIJMAAiCTAAIgkwACIJMDgZk6HCMcIMAAi5V3QcrzM2uo5MVOuEAjAl8ICbHox0+WFTcctKdc8\nBeCwsAA7YHZhbMEGsKpPOxjbWoAtpxAlFsAe09EyIszaCbDt2UUAGmMVIgCRwiqwUspyFeJQb63+\nCYBWhQVYtxZO4xa5RQoXA4PvmUIEIJIAAyCSAAMgkgADIJIAAyCSAAMgkgADIJIAAyCSAIP7uSgz\nHCDAAIgkwACIJMAAiJR3Mt+NU847Gz3Ac4QF2MZVK6fXYnZBS4DmhQXYW9MYA6BhTQXYanHW/12e\nLNsAVvVpX+ZoKsBWSSyAPZYHZSpnFSIAkcIqsFLKcqnhMGG4+icAWhUWYN1aOI1b5BbAc5hChKv1\nfWdfC74nwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwACIJMAAiCTAAIgkwKAKpXQJp/+GiggwACLl\nBVj/z8YNrmwPALcIOxv97FLLy9PPSy+Ah8irwDasRhoATQqrwA6Y1WQSDmBV3AxWOwE2dP34f1e5\nBPjIdLSMCLN2Amz72BgAjQkLsFLKuF8gsQCeLCzAurUpwdkWYQbwBE2tQgTgOQQYAJEEGACRBBgA\nkQQYAJEEGACRBBhcqu87X/SAUwgwACIJMAAiCTAAIgkwACIJMAAiCTAAIuWdjX55OZU9fwKgMWEB\nNr301+plwIYtrhAG0LymphCFFsBzhFVge8zKr3FecSDkqFYpztPBnWajZf2aCrCh912gGeCAjV3/\nOjU1hdiJK4DHCKvASinLpYbDnOGw3UJEgIcIC7BuLZmGLRIL4FFam0IE4CEEGACRBBgAkQQYAJEE\nGACRBBgAkQQYAJEEGACRBBhcx7l64UQCDIBIAgyASAIMgEh5J/Pd4FT0AM/RToBNL8Q8uygzAO0x\nhQhApHYqsBdKwnWxSVVJne9Nzg/U8ebe1HyAmUvkV/r+/5LjxrfY7Q2gVX3f159hzQcY/MoYG3d9\nPdnXonm4dgKslGIVIrco5YYskV7QToB1cot8O7NQekFnFSKcYgge4EoCDIBIAgzCmD+EgQCDc5hF\nhIsJMAAiCTBIYv4QRgIMgEgCDE7jMBhcSYABEEmAwUW+P3zlABhMCTAAIgkwOJPDYHAZAQZApLyz\n0W9cM8XlVACeIyzA+v7/r7A8/Xk0bFn9EwAtaWoKUWjRMEsQYSasAttjVn71fw+pCzmAVX3aAqSq\nA+zT7BluP7uZxOJiO6+qDLXZ2PWvU9UBdiB7xBXp5B/sVHWALZVSlksNhznDYbuFiAAPERZg3Voy\nDVskFsCjNLUKEVplUhGWBBgAkQQYnM8ZEeECAgyASAIMgEgCDK5gFQacToABEEmAQe1Ub7BKgAEQ\nSYABEEmAwU9881UwXyODPQQYAJHyAqz/Z+MGV7YHgFuEnY1+erXl2ZWXx42XNwp+yBJEeCWvAtuw\nGmkANKmpAFvV/3V3c9ZV27AZ7TzmVRVVWztfiWhnRCO7uttZ/1A5U/UU4qwTt6ur6RWZp6WYmoy7\nDIsJvQFJMR0tIzKs6gD7KHu2j41BFuEHb1UdYEullHG/QGLRPBkGGxof+iOqYFpXuu7Y+3D4bHoP\nc4/606HxAAOgVe2vQgSgSQIMgEgCDIBIAgyASGHL6LfNVtivLrifbbnL7KSO4/bV7wncaLWdr7r3\nXjV341LlzetC+nN8f/qwf+Nt71XSzqWmAqz7292zL4q9PRHwNVZX9q827N7vt83aufwQVtLOqQq7\ncVXlzRvV3J/T9+eyeZU0OOXD3k1GzmXvdVV+2AdNTSFW1bOvlFJWT6Jf21fWVttZuQq7MVrN/Rnx\n/gz6sN/dhINaq8C6yU7E3Q35wKxO5xjdeC79+QvV9mpt1dUeqQG2ep7fWbVbgz3nI769wR+dNLkG\nqw2uv9lZ9Ocv1NmrFY6cO6UG2Ku+ru01eNueGvZ6bm/Ap1anZeL+FTXTn79Qc69W27Bt9Xbop5Z7\n5XELkypc8JPSzq765s1U3rwupD9zVyHW06sbs1mvfq1HOwEGwKM0tQoRgOcQYABEEmAARBJgAEQS\nYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJg\nAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmAARBJgAEQSYABEEmA8Uf/XsGX2/+8f86P7\n/ujGFz8aXOm/uxsA9yilbPx6ymNecMePjFlVSpn9fE0D4FwCDLqu68ZBfCzCZqP8eLNud97M7j59\niulzjRGy5+lm950++Op9Z3ffbsCefxTUQ4DxUKvje7eIk9mIv7z9q8ec3mUjHlazas/TLeN2+lw7\nnx2iCTAe6tMq6sTHPOvpJBMPJ8DgpeV0XD1PJ71AgMGKVwelZjfYeISNY1r7b7/RvG4xYzm970eP\nBqHsxMERdxVAv3hexRyhfA8MjrhrxJ+WVvBw9rwAiKQCAyCSAAMgkgADIJIAAyCSAAMgUuNfZLbg\n+PGWi2y9Jd6adZoee6j616g3HmBdwmvQ5XyTNKudfd/NGvvJqeSvUGF/rnXacGqPe9qzX4WduSqo\nnXc34T1TiLRpORB3XVdKwEB8o1edFjKa8TgCjCa92cc1HC+tpteoFJ1GdQQYrRnOc3t3K8Jsp9dA\nhlGbjNnYw1KmmznLnoH401s2T6exFDF4qsB4KPXEATqNqggw2qE+OODTTpNh1EOA0YgD6WUshmgC\nDJ7rWM0q+KmEAKMFhycPjcWQS4DBQ31zyFDwUwMBRrwv124YiyGUAAOOEPzcToCRzdL5Y/QbDRBg\noJg4SL9xLwFGMGXEMfqNNgRcD2y8kMPsxFyz7dPrPdR/Ci9qMxQT3jif0m/cqPYAm55QcnlyyTG6\nhh/kFsBz1B5gG5Zx1a9dcHd2JT4h1ww7/sfoN16Ju25pcIANpmXZrCCbboS3zIYdo9+aMR0tI8Is\nOMBm9ZagehQjJpC9CnF6eOzelkCEHwW/9fTcovYKrJSyXG04bpz+6dViRdjPbBgEqT3AukUgvVpw\nKLeeQ8YAXfoUIrDfT4PfLCLXE2CEUX4BAwEGf6gkIIUAg0e4oHKV/VxMgJHE/CEwEmAwp5I4TNdx\nJQEGQCQBRgzzh4fpOpokwIAzmUXkMgIMVhiFoX4CjAwmwYAZAQaNuz771a9cQ4ABEEmAEeCW+UNl\nBFROgAHnE/9cQIBByyx+oWECjNrdOAQrI6BmAVdknun/jSjLSzBv/Am42BD/Pov8TliA9X0/htP0\n59GwZfVPALSkqSlEodUeu/Df0Hu0LawC22NWfvV/D2IIObiMWcQsfdoh36YCbOj9WURJLL5hCOY5\nNnb969TUFGInrtoiOYANYRVYKWW51HCYMxy2W4gIg0riXwnL74QFWLeWTMMWicWPGIKhTq1NIdIM\nmQFsE2DAbzmhCT8iwKBB6leeQIBRI+Mv8JYAg/dMgn1JB/ILAgyASAKM6pg//JIO5CEEGHAFs4ic\nToDBLsZfqI0Aoy6mv4CdBBhwEVUs5xJg0BQlLM8hwKiIwRfYT4DBXmbAoCoCDLiOnQBOJMCohfnD\n7+lDHkWAARBJgMEHzIB9Tx9ylv/ubsB7/b83e1mbHOn7ftjeTz4Tq7ekZua+gE/VHmBjPs1+7v4m\n1kBuATxH8BRiKWWWWH3fL1MNHiKoijWLyClqr8A+Ms4lbhRqqrQKBY280LC4AqCdAHuVTBKLcw3V\ng7cV7dk+RlOh4CnEqYi+BkZmEfle7RVYKWW2CnE2Q/jqZqRQzZxCN/JAtQdYtwikV7/KLchiMpYv\nNTKFCFcy/QU1EGDcyQ44cJgAA26jluUbAgziKWR5JgHGbQy7wDcEGBxh7ussepLDBBgAkQQY9zB/\neBY9yWMJMOBmZhE5RoABEEmAcYM2Zr3UDXAvAQbcz94ABwgwrtZG+VUJncmTCTCgCoowPiXA4Dhj\nLtxIgHEpU17AWQQYUAsVLR/JC7D+n40bXNkeuItylof77+4GfKbv+/LvIzv9edxyR6PYy4ALnCiv\nAttQSikGSK5l1utc+pP9wiqwA2ZlmYS7i/ILKhc3idV+gEksmtTwDsFQhLX6r6vZdLSMCLOmphCp\nlvEIOF1YgJVSxlWI09Uc97aKh3PY5nS6lD3yphCXU4KzLeYMAZ4grAIjkfnD0+lS6AQYUCeziLwl\nwPith9QKRlu4ngADKmW3gG0CjB96SPkF3EKAQZhH7RYowtggwOAchlq4mADjVx5VKADXE2BA1ZS2\nvCLA+Anl14/oWBgJMKB2ijBWCTDO99gqwTj7O/qWJQEGQCQBxskeW35d4OF9qwhjRoBxpoePsJ1B\nFi4kwIAY9g+YEmCcRvkFXEmAQQb7BwNFGKP/7m7Ae/2/d2v5+/Gdbe8nb+rig345w+toGGH1Bvxa\n7QHW9/2YRm9/llvwBHYRGDQ1hdj3fW9y4Q5Gk1/TwzMmEunqr8A+Ms4lTkuxWaSp0k5nbIU2xBUA\n7QTYq2SSWFzPHNcFdPLpNnb969TIFGJEXzfJCMKNTCQ+XO0VWClludqwlLK6fforF5Be19DPsKr2\nAOsWgTT++mo71MAE1zX085M1MoXI9YwaVMJE4mMJMI6QXpfR1XvIsGcSYHzMkLqTURV+SoDxGelF\nnewuPJAA4wPS62I6/CMy7GkEGHsZTA8wpF5Mhz+KAGMX6UUKGfYcAoz3pNctdPthMuwhBBhvGEa/\nZDC9hW5/AgHGFulFriHDxFjDBBjrhk++9LqLzj/F0IcyrFUCjBXDB94AehbTWTcqRf83K+BkvlzM\nvv/tvASnGzNMx7ZEgPH/fMJ/x0nTbzedTvRCtEGA0XU+1TWRcz8lxloiwB5tPDDgk3wBRVg9xFgb\nLOKoQn/5IeZxkeHw3+57ZRwKz21nnQkX0Z8HGjmu77hytX1EZ6ZQgT3I9INT4Sj5BNtFWJ3p1byx\nz01IxGkqwMZdm+IN2HXd2tdfdMztTCRWa5lks+3Upp0A6/t+zK3pz0/wak7iSX2QZJlhDsZUZfZC\nbMz5ecnu1c5Avxpgj5ltfsq/sy2zj54XMU4jg+cr9YdDOxXYqvpfgJM85d/ZNC8ifMYqRAAiCTAA\nIrVzDKyzChHgSZoKMACeo6lFHLMKbFmQ1VOizdZMjttLKfU0snvRzlfde6+au3Gp8uZ1If05WW/s\nw37c296rpJ1LTQVY97e7Z6vqK/mi2OqJZFYbdu+32WbtXH4IK2nnVIXduKry5o1q7s/p+3PZvEoa\nnPJh7yYj57L3uio/7IOmFnFU1bOvlFKW7ez7vrYzpK22s3IVdmO0mvsz4v0Z9GG/uwkHtVaBdZOd\niLsb8oFZnc4xuvFc+vMXqu3V2qqrPVIDbPbyr84c1mC1nTO3N3hPI6uy2uD6m51Ff/5Cnb1a4ci5\nU2qAverr2l6Dt+2pYa/n9gZ8anVaJu5fUTP9+Qs192q1DdtWb4d+arlXHrcwqcIFPynt7Kpv3kzl\nzetC+jN3FWI9vboxm/Xq13q0E2AAPEpTqxABeA4BBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAk\nAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQB\nBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQCQBBkAkAQZAJAEGQKT/BW+NFMnG+s99\nAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["q = 200;\n", "clf;\n", "subplot(2,1,1);\n", "plot(-N/2+1:N/2, h);\n", "axis([-q q min(h) max(h)]);\n", "title('Filter, Spacial (zoom)');\n", "subplot(2,1,2);\n", "plot(-N/2+1:N/2, hf);\n", "axis([-q q 0 max(hf)]);\n", "title('Filter, Fourier (zoom)');"]}, {"source": ["We generate a synthetic sparse signal $x_0$, with only a small number of non zero\n", "coefficients.\n", "\n", "\n", "Number of Diracs of the signal."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["s = round(N*.03);"]}, {"source": ["Set the seed-number (for reproductibility)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["rand('state', 1);\n", "randn('state', 1);"]}, {"source": ["Location of the diracs."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sel = randperm(N); sel = sel(1:s);"]}, {"source": ["Signal $x_0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x0 = zeros(N,1); x0(sel) = 1;\n", "x0 = x0 .* sign(randn(N,1)) .* (1-.3*rand(N,1));"]}, {"source": ["Noise level."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = .06;"]}, {"source": ["Compute the measurements $y=\\Phi x_0 + w$ where $w$ is a realization\n", "of a Gaussian white noise of variance $\\si^2$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["y = Phi(x0) + sigma*randn(N,1);"]}, {"source": ["Display signals and measurements."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAWvElEQVR4nO3d3Zaj\nLLcG0LhH3/8tuw98289WY/xBYMGcBz3SVamEIPIAkmQYx/EDANH8X+kCAMAdAgyAkAQYACEJMABC\nEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmCQwPDX8ie3HypRoaBxf0oXAMIbhmH+\nVof5tu95gLcJMEhplWTL6dT0q9VPtndYmR9HIsKKAIOnljm0jJnVzGx5//knu3cAzhBgkMB2CfGk\nn7ll+gXf2MQBTz3ZrzH+dXwf8zPYMgODBHaXEM8Ez/EdtnOv1cIj9MzSBOSQZBlwfhCLivAxA4NX\n7c7MgCQEGLxIbsF7bOIAICQr6RCJTRwwE2AAhGQJEYCQBBgAIQkwAEKyjb6A+bMXXH8Ezlt+bIve\n42MGlt+yCfp8O+CkVXeh9/gIsMy2bU4rBLin8W309X2G925t11ZIWrJqchpbUN866hcPaP3p0P41\nsOLH4N9vNdy5w9slrOGDX5WhSAH22tujIiS5BtPhgUhShu3RHMfP92B7Wob6Rv87LCECp7gGk9Yw\nqMOnBFhW2+FO6YEgnJK/q5369yZ7+enjwObbJ1/gqq/QdXx6WEKszdTshkH7e4WKbcN2thfxsCbf\n9a73WDEDe13x1fbiBchThnkk+21IW7weihfgU8ew/V49JJyH5TkQl1Zci8wya2iQDwkwWuD9CbvG\n8Z/EetJfxe/rstL88hBg0Lgpe54nkGsw71GZ9wgw4Kw5C/N3uLp4tmzigP/xWXOVGMfYx2JV/u1v\nJ7ZjPGQGRjzLDdZTN5Hk/Qne51SVeZ4XtIu34pqBACOYb/s1lv2dzoIaaJBvE2Dw+XyZb5mEQc1c\nA4PP58tFi8YGzqGvKsGWGRh0wRU+2mMGRnj3JhPb6Uj0nW/QGzMwgkny6RLfpiPRd75BVwQYIYkZ\nQIDRHRsOoQ0CjO7sztuan8x5Xy3tEWDQC++rpTECjB6ZjkADBBidsuEQohNgAIQkwAAISYABEJIA\nAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAIKd43Mg9/v/di3HwE0LD4SoztbwFoSbAAG4ZhTqbl\n7ZncAuhEa0uIwzAMvpoQoAPBZmA/TTOw1URtewcAVsKN/psKsN1wklgAZyx7yxBh1s4SYojqBiCV\nYDOwcRy3uxCnBcPdXwHQqmAB9tkLp/kncgugH+0sIQLQFQEGQEgCDICQBBgAIQkwAEISYACEJMAA\nCEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkw\nAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACE9Kd0AVIahmG6MY5j2ZIA8LZ2\nAmwYhjm3lrcBaFI7Hf1ugA3D8PlkeoHj+Pk7A8ztzFMXLF6H1DYNqD8c2pmBfZepIynYYZ15av1p\nTmqbJtSeYO0H2NtTzGH4b5wy3zj/J8+f9OABz9ynTnmKOmVM2gPx8+f37pbc9LwFm0TZ+nn7hQc6\n1w4MEUZhdiECEFI7M7BxHO1CBOhHOwH2kVsAPbGECEBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQ\nkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAA\nhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAjpT+kCXDYMw3RjHMdvv9r9LQAt\nCRZgwzDMybS8PZNbAJ0IFmA/TZOwZYwtp2UfCQfwxaq3rF9rATbl03JyJrEAzjgY+tep6gC7OnmS\nVQD9qDrALgXS7iUxAFpVdYBtjeO43YU4RdfurwBoVbAA++yFk8tdAB3yRmYAQhJgAIQkwAAISYAB\nEJIAAyAkAQZASAIMgJAEGCQ2vZ9+GD4RPkwOAhNgkNIqtGRYVeZRxUvDi+Xjk4EAg2R0WzV7++gY\nu+QnwIBOZY40khNgAE/JqiIEGCTj06S7tXvotYe3xfs0+s7NA73phjOkNuP4z2C8qgO0LJgZw6ey\no8MNAiyBgqEyDE7C6jgidZqOy3y2Jj9MNY9dWiXA0nsvVG6PmuetvU4qljpsEq++3t4qszjXwNpn\n4WhrOWlWJxCUAIvkxvhu2zvrr9XAktogLgGW1dsfBLB8Is7rvLo6f/nEJcDySdJNLCdh47g/J1tt\nKoHJtwYDQdnEkcBy99GlDuLeJfTjPxFaQCcEWBqBBraBiprNk4uLHW7kg0oIsC7oYZeevx9o+7Gt\ngWp49XYliEuAlfS81/PBHE90W123F72hKgIsn+XAP0mvsbtFfju+1kOltTt9iTUJ+5RrFYZcJCTA\ncstwxuoUXrW7BKfOz/g25IJ7bKO/7+S3rw6lLzgUL8D5Mrz6hbbF66F4ASopwyVNtgRlSEWA3eTb\nV5PLVqXP306+mjQ8mUNoObuWh8bHffGNAEsm/zm27TcbW415Z/Sd5lmmNwXffmuwDvqA2uAkARbb\nsgNtLL0a1m0H3fyQi8yGsekW9OYi72699dozpZGnSms4cDWUoaD55X97yZ3XTy3qT4fGA+w99qEl\nl6dKd4c0mY9dDWWonPOLMywh3pTwMj6TPFVawyqW1nKVGmOXGRidKvsOJG82h+cEGAAhWUIEICQf\nJZXeMPw3r533QM7T3O1Pkj/1z2csXoa3C7AsSfED8fn+qps/EKs9wOM45q+E+SkKtsYaDsT5Rpjt\n9ExCgKW0PGPnRjPf3v7kjTLMLXX3GYuX4fPvqZKh1yhbCaXKsKzqUgdiVQNFDsSy7y5SCT9f9atl\nuNoj5Tk1ErKEmNI4jmUPeQ0NroYy1HDuTV12DWWooSqKl6FPxXukt5mBNah4f1G8467BcuBfsABl\ny1CD1SynbDFIS4A1pZJzdbXgntn0vMt/8yt+COpRMD5rWBCbr/wtLwGSiiXE1ph7jX99CtVGDZVA\nJabgNKB5SddrCy8ptfltu+krfxnOPGNXuxALbveq4UCs5j3dVkLZMjS8C1GAARCSJUQAQhJgAIQk\nwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAh\nCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQb/\nGIZhGIZv/61Q5cWD9/wpXQDgfyE0juPydrkSQQACDM5aRcty6rMMnuUdljemO6z+dk6p7Z/s3m33\nSbdPPd95+YcSkcYIMFibgmH57+dvSCxjZhsqPxNiGU6rxNr6drfVgyzvdvxCrtQBBCDA4ILV5bHl\nr66u/p28dnXmbhYe6ZMAg7OWU67V6tzn4nrd7pLg7bsdP7XpF60SYLBju+y2muVsdyeuLoNNdoPn\n5ITp0rzqW8LtlgraYGgGjTMDo1XeBwYtk140TOMGICQzMABCEmAAhCTAAAhJgAEQkgADIKTG38js\nLZxAUePnE7UXqn+PeuMB9qngGBR/I07xAihDJQVQhvwFGIb9Lqh4JfwsQ4jRvyVEAEISYAAvijCT\niarBAAsx8wXgofLrsAntfn1tSy8Q7hmGj/OgiKlPilj5ITrPpmZgyy/JBZYsTNCe9nchbr+0qVRJ\nAJ57bz4d7vpL+wEmsYAi3oiDVyNm2VuGCLOmlhABqnI8fo6QEVUTYACE1GCAWTPsmSEt9KPBAAN2\nSfeqOBzPCTAAQhJg0IVUK+vDYOpALQQYACEJMABCEmAAhCTAAAhJgNGah1sMbFKAKARYMXrJh1Rg\nfr6WhaoIMABCEmAAhCTAmtL5qlrZBa7OKx/yE2Ct0Y2y4sIVrRJgJQmbZjiUkJ8AIySBQYdMplcE\nGCRTeayOY+0lhEsEWBn6EeAeA5GZAKN3277gXgdhbYdd8uY9AqwY/V1++hEq4WpWEgIMoAyTs4cE\nGAAh/SldgMuGvyOWcTMDHxaDme1vacw0ej1/nM8MdS3sQCDBAmwYhjmZlrdncosb5JYaIKJgAfbT\nNAlbxtjw78BbwvXs6qQtp2oLRj+GaFfkWguwKZ+WkzOJxdturGRqlVToYOhfp6Y2ceTJqgiHlZDe\naFomdjSsnQDLOV74+VRCDpa+nRHOFJ4IFmDjOA5/LXdzfPvVe/KfeGdS8+TrDtRrVFjUCouUxPnX\n5d1LVCLeNbCDnYcudzUmQy85pb7uuCCLnNwWbAYWiD7xvKt1VbDLe3hY6+ys6ywV/CTAqIK8JzQN\nuAgBBmfppDjPvDYDAXbHcbvUzbGkPUy+1UOea5AHT+EAxSXAXlFw5LXqDqJsX7aTgiK0utAEWCRJ\nTjZn7EM3sjbhWyBCq63t1VYerhJg10Rs8bvfOFyPPFVaKiEiNpja1F+H9ZewVQLsRWVX9rtlNbJO\nBYdNB+2hqsEcVwkw4JTbs9hKxhOVFIOEBNhlV8/hYXDm3HG1u6zhMtLJyV+Hc8TdQ1PDISM0ARbD\ndKr32fFV9ZJvdLjzn/x8IUVebFXV+x5J2SQBltjcHZw5YWrrnXlPiK8Bq7Y1VluwmXO5CAF238Eb\nM5v03vlZvMbe7n2Kv8BZuH72+TLjwVshnx8XFwjKEmA31dMlXTKftJV0ZCc7kRqKSlxBz1Z+EmCB\n6dZzel7blXejlRcPtgTYBe0FRpFX1F41nlFzPNRQtpe+5TxPY3tSgZWshQQlwAq4tPg+LbLX0MUc\nOLgS0NLEJUNJfnZnoTu7b69uqtiDVvQwIX56ko71tM8OdRFgz7+EsHgZgj7y/AiudXeo5jnT0sH7\nHKpqtMuBrHnbpIsAO6na5nswxMsz+lvm0PKHBzV2ZkVodZ/3zsn6p7BnnPmSgeKmguUp6quVUG0N\nPxeiIZ3UfoClGgZm6AS3J39Ou887x8yZUkX56pZsUm3UPv55wuo9+V08Zyz/9saH11y6c/KN8slb\nbIUTpjZirKkAG/5a/XxuPavfbLuA5Ctd8ydonL//6sZzpU6eJ/OzIg76u/xzuGVTfGmDQ1Xmd3fc\nrufME7KH6f7zS3FfvfBRdqycUDsBNgzD+Nc2w/6954XmWEnX/4b5yvkndYNevtuMq5YDjm/T4vlG\nA33Q7PY4bzlCvVchu2Pc7Xg3T7//cJx3PC6vqqN7rp0Au+rgiB405Us98pPh5G7BkltmWJKH2j7a\nqxl2Y3LZ6uftJl85SPJXNZQq7bMc/LZ4m0mYfMPw+XwCDD+7CLBpQnYwmK12ovDqMsLxH27/9tsH\nileu/hKecfIaZJIX+3Ce8fZS3s/i3Z6EnX/A915j/sj/d6R+vHpVnS4CbFy0zeMBVG2eF+nJ9fNv\nifXtUuKZMtTZLYZQfIB/LOfM+6eEGVawx7hxxB9feK6vEzzURYCtnPyc+JP3fMO09pitDLtPcf6H\nL92HyW4zyNB1XnKwkpGqqD8X5H8+0aU82Fb4q68uoUtTyfN/W6d2AmzauzEZTx+HbU7kTI7d8hSx\nmiT9LMbx0uv25w/P/4fV8tLHmWfzs+NOMsc9+XQ5H+qgKX5L9CRN5XzDrsTw9/N6luVf9WlRAviS\ndgLs8/nMuxD//dnyP/u3Dx8zVdmS7ZVI25/+bOXJm3j0cya5b930cTf6UjXeftich/XniTydJveK\nlOEUyOmgrraj9nCaCrA9/+vstwcy6DF725mlmDN3O37w25V/Kb+LT55e9a0O33jVB3mw/dV7Z9aq\nq317j+vt3z5/3quPf2lE3kzX13yArd1ulNkOeVVt62dzX3YoJ5MvValy2nbQy2vjBSep7/XmP8+U\n9y6onFxmyNMMMl99TCXETPq5XgLsyVGJdURfcnwaN7C8XltpayvPeQVLfn56UfnUPMPMcjvujNjk\negmwh04e2nuD8Xtvu3lj4P/cpSJVWP7KnblWkXYLxo33iUdxpqjhdnOkEuUFCjCuidKyZ/eW4yC/\nSkalf8tQ9yz18/kIsEpU/oEgNUi+/bIx712cD1rtb2+yqEesD2dIS4A1opKxWyXqT7tXSzguPpr2\nYavQqMqqvyWXJcDo19VPijtDj0N+3Y4zBFgtyl5771blZ37lxXuVNsxPAqwKPfdTl+jUmmHRe0WF\n3CDAkqm/b23jDGngJYST/LtLIAkBdl+SlTrd8XvaCOyedXL4rPnfJsDgh0660dveqB91vvJtNNZ5\n8gkw2Gdc/JP6oSwBRo2szeahnkPzAQgCLLec3z0R2rZaVNQbqq3Vl1Ymm5k1VnvgchJgzbKFgYS0\npYKaCd3kBFhKRU7ylgaVaelzk1CNZan/AwIMjtzuPpofVehYKU6ApdF8b1WJb51mneulFRaJat1r\nLZ23MQEGbcoQ6p33nhT3p3QBLhv+TnbGzdkzLOZB29++WaTpGc/eM9Uz8jYd9EQ9vO24hl3n/iZY\ngA3DMCfT8vYsZ279+7yn7pOkFU6Pc+mFOgF6c3DEb7QfrspTww5ia0uIwzAMumrQu9GBYDOwn6YZ\n2Gqitr1DQuISaEO40X/VAXY1e3bvkGFR0VC3oApXRyssEtHl6WSWvWWIMKs6wC5lz+4lsW7pQ6lT\nt1fg5hfuxEyo6gDbGsdxuwtxiq7dX7XKW0YgIumVVrAA++yF0/yT5nPrQBsffdvt8By4obVdiAB0\nQoC1wJQFzntvln/yEpcTNhUB9kighhioqABnCDAoyTU/uE2AAR0pOFwwUklOgMF93+ZPuirIQIA9\npauKKMPbcTQMeJsAy02/VpxDQEGaX0ICjPCK9Ag+EwiKE2BZGXzBx95LEhFgAIQkwIC+2EnfDAFG\nCyL2C66iwUMCjOq4QAKcIcAACEmAAVmZXpOKAINOWaolOgEG9wkAKEiAQUl9RmCfr5rkBBgAIQkw\nuMk0AspqMMAG7w4F6MCf0gVISXQRizkcPNHUDGwcx1GXAJTjE8JyamoGtms1LZNwALvCLWJFDbDz\nsSSxwnHEoIhlbxkizKIGmFhqlQMLnNTUNTAA+tFggJmcAWXphPJoMMCAY3bK0QYBBkBIAoxO+TIR\niE6AAaRkYJSNAINO6WeJToBBj6QXDRBgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJg\nAIQkwAAISYC9bij9zUvFC1BhGYp8kFJtldBtGYoXQBlSEWAAhCTAAAhJgAEQ0jA2/bUKDSzy8prx\n89E84Kv606HxAAOgVZYQAQhJgAEQkgADICQBBkBIf0oXoEHD8N/WmHkP5LxTZvuT5E/98xmLl+Ht\nAixLUvxAfL6/6uYPxGoP8DiO+SthfoqCrbGGA3G+EWY7PZMQYCktz9i50cy3tz95owxzS919xuJl\n+Px7qmToNcpWQqkyLKu61IFY1UCRA7Hsu4tUws9X/WoZrvZIeU6NhCwhpjSOY9lDXkODq6EMNZx7\nU5ddQxlqqIriZehT8R7pbWZgDSreXxTvuGuwHPgXLEDZMtRgNcspWwzSEmBNqeRcXS24ZzY97/Lf\n/IofgnoUjM8aFsTmK3/LS4CkYgmxNeZe41+fQrVRQyVQiSk4DWhe0vXawktKbX7bbvrKX4Yzz9jV\nLsSC271qOBCreU+3lVC2DA3vQhRgAIRkCRGAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBI\nAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQ\nkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACH9P9Dohyo6\n2870AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; ms = 20;\n", "subplot(2,1,1);\n", "u = x0; u(x0==0) = NaN;\n", "stem(u, 'b.', 'MarkerSize', ms); axis('tight');\n", "title('Signal x_0');\n", "subplot(2,1,2);\n", "plot(y); axis('tight');\n", "title('Measurements y');"]}, {"source": ["Sparse-Spikes Deconvolution\n", "---------------------------\n", "We now implement the FB algorithm for the sparse spikes problem.\n", "\n", "\n", "Regularization strenght $\\la$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["lambda = 1;"]}, {"source": ["Define the proximity operator of $\\ga g$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["proxg = @(x,gamma)perform_thresholding(x, lambda*gamma, 'soft');"]}, {"source": ["Define the gradient operator of $f$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gradf = @(x)Phi(Phi(x)-y);"]}, {"source": ["Define the Lipschitz constant of $f$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["L = max(abs(fft(h)))^2;"]}, {"source": ["Gradient step size $\\ga$, should be smaller than $2/L$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = 1.95 / L;"]}, {"source": ["Initialization $x^{(0)}$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x = y;"]}, {"source": ["Perform one step of FB."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x = proxg( x - gamma*gradf(x), gamma );"]}, {"source": ["__Exercise 1__\n", "\n", "Compute the solution of L1 deconvolution.\n", "Keep track of the degay of the energy $E=f+g$.\n", "isplay energy decay"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAP0ElEQVR4nO3d3Xqj\nRhqFUciT+79l5oC0hgaEEOKndtVaRx7HmUh026+/olTqh2HoACDNP08/AAA4QsAAiCRgAEQSMAAi\nCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyCSgAEQ6d+nH8AJ+r5/fTwMw4OPBIDb9NE/8afpeol+RgDsVMME\n9irWas8AqFLwPbB3uZIxgBYEB2w0XTC0eAjQjviAAdCmGu6BbfizmmhREeA75a9pVR6wUfl/DDfo\n++wdp+dyNV5ciilX4yViM4ElRAAiCRgAkaoKWMTMC8ApggM2LlUvo2UJe5XLMuVqvLgUU65Glho2\ncRi8ABoUPIF1a78u+QUKoBHxE5hiAbQpewIDoFkCBkCkJgJmkwdAfZoIGAD1ETAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQqfqAeSswgDpVHzAA6iRg\nAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIrUSsN6ZiAB1aSVgAFRGwACIJGAA\nRBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAi1ROw3nm9AC2pJGDqBdCaf59+AL+S\nLoA2ZU9ge+o1DDc8EADulj2BDX/qZA4DaE32BAZAs7InsD36vu+64c8H/x/aAJiKW8qqP2DDMIx/\nKNIFsGH6QzIiZoUGbPvaSREA7oEBEKnQCcyMBcA2ExgAkQQMgEgCBkCkQu+Bfcs9M4DWmMAAiCRg\nAEQSMAAiCRgAkRoKWMLJXgDs1VDAAKiJgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkYAJEEDIBITQRsGJ5+BACcrYmAAVAfAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyBSWwHr+6cfAQAnaStg\nAFRDwACIJGAARPr36Qfwq/7v+1rDMDz1SAC4U/YE1i92ZSw/A0CVggP2atXwx+zzAFQsOGCj6Zqh\n9UOAdsQHDIA2BW/iMG8BtCw4YEvj3a9Z2P7cElM7gC1xGwjqCdhqvV6fSftzAbjb9OdnRMwKDdj2\ntVudsawoAjSl0IDtNN1J/+wjAeBmhQZsf5CkC6BNJwRsdbnvhq5ELNECcJHjAdvux22LezvzOQz2\ncQBU5UjAps3Y6NP4ZXZYAHCFrwO2P0izHZmnN0wUAVr2dcCOZUNsADiXsxABiCRgAEQSMAAiCRgA\nkQQMgEgCBkCkgwFzjBMAz/ppAuv7XskAeMTxgDkjCoAHHTxKaqReADzl6wlsGIaxW8MwWEIE4Ck/\nvR+YCQyApxy8ByZdADyrudeBWfIEqMPXAVu977X8pNtjAFyquQkMgDoIGACRBAyASAIGQCQBAyCS\ngAEQScAAiHTwKKnV13h54RcAtzGBARDp6wnMKYgAlMAEBkCkhgJmdASoyZkBc4AvALdpaAIDoCYC\nBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIB99OZZVzfgG4TYsTmOOuACrwdcBWDzxcftK5\niABcqsUJDIAKnHkP7BGzOc99OIBGZE9gq4uZjzwSAG4WPIGNrZqOXONn+r43hwFUL3sCm/nYLV0D\nqEbwBGbMAmhZVRPYclHx/Vde/2gAuNLBCWx1r8RTGyim/11jGUAjgpcQX17Rer16ermzY/rl9z0y\ngBxxu7i/Dtg9I872dXz3GIZhWP6Lsy9O+wMCuMnmr/4lquoe2B6vP6CEPx0A3ip0CXHfRgyv9wJo\nV/YEtjxBeM+/pXoAFSh0AttjvOP18abXhr4XM4BU2RPYslUWFQEaETyBjY4Vaxj+28RhCAMIlT2B\nAdCsdgNmPz1AtK8DduzVbYW/Jq7sRwfAiiMT2Orev+2vP/BfuYG7XwC5Dh4lNW3Yu20UEWfs2s0B\nEOrgLsTpqYPHzi0EgF8c30Y/PQN+45+WzxAGkOiE14EFteojDQNIcfkLmSPms9cQBkCKawP27sD4\nAg+St5AIkOXCFzJvVGr1nSfLUfBDA+A/7Z7EsWTwAggiYH9xvhRAigsDtrFOWOA9sBcNA4hw7SaO\ndw0rtl4jGzoAynf5NvrCW/WRhgGUyT2wddNoWUsEKJCAvaVhACUTsC0aBlAsAftAwwDKdOEmjo9n\nbaTs75ielDh+EPLAAWp2YcBee+hTQrVhfAYyBlCOa5cQK0jX1OzZ9L1FRYDHXH4PrL6GyRhACS5/\nIXOVZiuK04/r6jVAuexCPG45jXUGMoC7/DqB7Xlbr8pWEWeW01hnIAO43h0TWMnvXXmWcRozkAHc\n5tcJbHWv/PQz48clv3/KuQxkAPf4dQJbfaXXNF2NdGtmeyAzkwH8zi7Ea717e0wzGcCPBOwmSgZw\nLgG72+pNsk7JAL50TsBmezSm2w5bvhO2YeOQ++n/dNkA3jltF+Jyr7xo7fRudbEzlgG8d8IENt1z\nOP3M8mO27SlZJ2YAXdedeA9MqE60/S6aYgbQnXgPbPYZPTvLxljWiRnQsBMCtnpSVDtHb9xm+TYu\nM7PPuPxA3U47zHe5C1HDLrW9zNgZzoDanbOEuDxKqoUDfMshZkCDvJC5Nl/FrNMzIJaA1ezjbbPO\ncAbEuvwkjts48uMjK41ATZzE0SgxA9KdcxLH4/WyZ+QXYgYkumQX4s3GBUwNO8VXrzYTM+BBNnGw\nxaFWQLHiA+bl0rcRM6AoXwfsq5W6q9Oy58HMvkbtTiFmUJ+4GzHxE9jHICnW1Rw3DHV4/NVQ3/o6\nYPf0YPvavd6BTJyK4rhh4E7ZE9iyc17OXI79u/P9cQEHFBowBarM/ntmnZ4B+xQasD1mkTN7pfC6\naeAUwQGjApYZgcMEjFJYZgS+Uk/ALB7WxDIj8FE9AaNWlhmBVQJGEieAAC8CRioxg8YJGDUQM2iQ\ngFGbnWczKhmkEzCqZSyDugkYTdgZMyWDIAJGczbWGJUMgggY7doYy5QMyidg0HXvxzK3yqBYAgZ/\nscAIKQQM1ikZFE7A4AMlgzIJGOz1sWQyBncSMPjaxx0fSgY3EDA4TsngQQIGJ1AyuJ+AwZmUDG4j\nYHAJJYOrCRhca7tkMgaHCRjcZGyVgQzOImBwK0uLcBYBg2coGfxIwOBhbpLBMQIGpVgtmYEM3hEw\nKI7tHrCHgEGh3CSDbQIGpVMyWCVgEMN2D5gSMMjjJhl0Aga5DGQ0TsAgnv33tOmfpx8AcJphWMlV\n389HNKiDCQxqszGQmcaoiQkMqrUcyExj1ETAoHIyRq0EDJogY9RHwKAhMkZNBAyaI2PUQcCgUTJG\nOgGDpskYuQQMWM8YFE7AgP/MMmYUo3DZJ3H0a99eg8MG4Aezo+4d4UGxTGDACiuKlC97AuvMW3AZ\noxiFM4EBW4xiFCt+Auv+vhNmIIPTGcUoU/AENnZrto9jdVsH8DujGKWpYQJ7TV2vpE3nsFnSjGhw\nmFGsbnEDQHzApkEahmH5B6BYcK5hmL9Vpm+yOmz86l+mQgO2fe3Gq6xM8JTlKObbkfsF3wMDnuXY\nDp5V6ARmuoIIRjEeFDyB9X2/ugVR/OBmRjEeUegEtse4ZSPiTiNUz84O7hc8gXWLYWsYBuMXPGV5\nmD1cKngCGykWFGU6inmhGJfKnsCAAjmzg3sIGHA+y4ncQMCAq9idyKUEDLiQ5USuI2DAtZbLiTLG\nKQQMuINRjNMJGHAToxjnEjDgVstRTMY4RsCAu81Gsc6KIocIGPAMK4r8SMCAJ1lR5DABAx5mRZFj\nBAwoghVFviVgQEGMYuwnYEBZjGLsJGBAiYxifCRgQKGMYmwTMKBoRjHeETCgdN4ek1UCBmSwnMiM\ngAExLCcyJWBAEsuJvAgYkEfD6AQMCOWWGAIGpLKc2DgBA7IZxZolYEA8uxPbJGBADSwnNkjAgHpY\nTmyKgAFVWS4nylitBAyozWw5sZOxSgkYUCcZq96/Tz8AgAuNDZt26/XxLG/EMYEB9VtOY52BLJ8J\nDGjFq2EGsjoIGNCc5bpip2SBBAxo1GrGOiXLIWBA01bXFZefEbMCCRhA1316jzFjWYEEDGDOWBZB\nwADe2jmWdWL2BAED2GVjLOvE7AkCBvCdj+/b4obZPWoIWD/5GzT4+wLcyBrjg+ID1v/9V6bvew0D\nHiFmN8sO2FivV7F655oBZdgfs07PjgoO2Kxe48d93xvCgKLsv2e2/Ho2BAdslXQBJfsYs+Xn/VR7\np4aA2cQBJFq+3+aq5ef9nBsFr7a9u+M1fUbLr8l9vkBTvrqnf8oPttkPzPJ/WtYwgc02cczugZX/\nZ3APtwanXI0Xl2KqqKuxcz57908PPI/t3/4LVGjAtq/duz6NmzgufFgAD1l9R+kNq/+0mDqfo9CA\nAbBttUZNVa3QgJUzxQME+XZQ2/MFxSpowfeAdy9k9tJmgE2ff/KXH4dCJ7CvbOycic4zwJ3ifuHP\nnsC6wH2fAJwiPmAAtOmfpx8AABwhYABEEjAAIgkYAJFq2EbfuD37MD8e2F/fif7Lt4ub/aNR3Vfj\n92dax6XYuVe57quxcczjnudV5sWxCzHb6iu1P57I+dUXhFoNWFNXY/utGNq5FHueSPVX493vc+kX\nxwQW7N3JI8tftT4eTVLZ2SU7D4Ou+Gq8O6RmpvpLsf97ZPk1S4lXY89D3fO8yrw47oHF2zh55N3f\npNfnP35BondLJe1cjeWv2+PH+59pNZditH06T8VXY/sR7nlehV8cAavfx7O1mjp8q9mrMQzDu7m8\na+xSrKryagx/bH/N6sdffc1TF8cSYrCU76I7FfWGhI+L3nRwivE9Aqd/KzZ29xDHBFaVxr85y1/S\nudPsajR7cV7Lp6Ou4W+Q+ghYPXxzds0//ZnZClKbDVs+6zavQ5UsIdZAujqLhwuz2xJt/tRefmts\n7EIkjoBlW+4Sbtm737VdHF6abXmVBCzeVz+dP37rNvW97Wq8uBRTzV6NPc+rqIvjHliwj39RZq/+\nmX1+zxcEGf42/eTr466Nq7F8IqtnKbVwKZbePakGr8ae51X4xTGBxTt8GMxXX1CTdq7Garc2vuDj\n/0Oc1zb65eeXX1z91XgnbvB6MYFVbuO8nJ1fUJN2rsbvz7SdS7HnM9VcjZn0i2MrDgCRTGAARBIw\nACIJGACRBAyASAIGz1jd3g3sJ2AARLKNHoBIJjAAIgkYPMM9MPiRgAEQScAAiCRgAEQSMAAiCRgA\nkQQMgEgCBkAkAQMgkoABEEnAAIjkMF8AIpnAAIgkYABEEjAAIgkYAJEEDIBI/wNurNWU9Lh/7QAA\nAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Display the result."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAATEklEQVR4nO3d0Xak\nKrcG0OKM/f6vzLmw299WK7EsBBbMebFHdpJOUah8LkQr5ZxfABDN/7VuAADcIcAACEmAARCSAAMg\nJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMCg\ngPTX+r+3/871V4HJ/de6ARBeSinnvP16/d9HX6X4S0AsAgxKWnJlCZhtqbR+/9fvbK1/R1zBkQCD\nb22z6rRI+iG3Tn8HuEKAQQG3J/d+zi3lF/zAIg741jdLNvJfP/yC4gxOqcCggOOU4MXg+eF3drXX\n6UQlzMzsBDzONCA8QQUGT1EwwaOcGAIQkkUcAIRkCrGB9bK96he4brvix+jxUoHVt90FrY4GLtoN\nF0aPlwCr7LjP2QsB7hl8EUd/d4Ce9nZvjWQku13OzhbUu4H6wQ3afzqMfw2s+Tb495l4J7/wdAt7\nuAlJG5o04Gx/+6oJRa7BTLghirThuDVzfr0Ptm/b0N/Z/wlTiMAlrsGUlZI+/JYAq+p4utP6RBAu\nqT/ULuP7kKN8Sq+1crr+BndjhaHjNcMUYm+W3S4l+98jdOwYjtVexM1afNW70WNHBfa45rPtzRtQ\npw3rmey7U9rm/dC8Aa8+Ttvv9UPBOqzOhvhoxrVJldnDDvklAcYI3J9wKud/Euub8Sr+WFeV3a8O\nAQaDW7Ln+wRyDeY5OvMeAQZctWZh/QHXEM+RRRzwP54114mcY2+LXfuPP11YjvElFRjxbBdYL8NE\nkfsT3OfUlbXOCzrEm3GtQIARzLv1GtvxzmBBD+yQTxNg8Hq9qbcUYdAz18Dg9Xpz0WKwE+fQV5Xg\nSAUGU3CFj/GowAjvXjFxLEeir3yD2ajACKbI0yXelSPRV77BVAQYIYkZQIAxHQsOYQwCjOmc1m3D\nF3Puq2U8Agxm4b5aBiPAmJFyBAYgwJiUBYcQnQADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJ\ngAEQkgADICQBBkBI8T6ROf393It8eARQ2nwkxvGnAIwkWICllNZk2n69klsAkxhtCjGllHw0IcAE\nglVgv1oqsF2hdvwFAHbCnf0PFWCn4SSxAK7YjpYhwmycKcQQ3Q1AKcEqsJzzcRXiMmF4+iMARhUs\nwF5n4bR+R24BzGOcKUQApiLAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQ\nkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAA\nhCTAAAhJgAEQkgADICQBBkBIAgyAkP5r3YCSUkrLFznnti0B4GnjBFhKac2t7dcADGmcgf40wFJK\nr1elN5jz628FWNuVl27YvAnpbQbQfziMU4G9V2kgaThgXXlp42lNepsh9J5g4wfY0yVmSn/OU9Yv\nrv+T71/0hz945Xf6VKepS8aU3RC/fv/erxW3vG7DXaJt/zz9xgMdaz9IEc7CrEIEIKRxKrCcs1WI\nAPMYJ8BecgtgJqYQAQhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQB\nBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJ\ngAEQkgADICQBBkBIAgyAkAQYACH917oBH0spLV/knN/96PSnAIwkWICllNZk2n69klsAkwgWYL9a\nirBtjG3LspeEA3hjN1r2b7QAW/JpW5xJLIArfjj171PXAfZp8SSrAObRdYB9FEinl8QAGFXXAXaU\ncz6uQlyi6/RHAIwqWIC9zsLJ5S6ACbmRGYCQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGhS33\n06f0ivAwOQhMgEFJu9CSYV1ZzyoeOr3Y/n0qEGBQjGGrZ09vHecu9QkwYFKVI43iBBjAt2RVEwIM\nivE06Wmdbnr7w9PiPY1+cuuJ3vKFI6Q3Of9zMt7VBto2TMXw6mzrcIMAK6BhqKTkIOyOLdKnZbus\nR2vxzdTzucuoBFh5z4XK7bPmdWmvg4qtCXeJR9/vbJ3ZnGtg4zNxdLQtmvUJBCXAIrlxfnccnY3X\nemBLbxCXAKvq6QcBbF+I6ybvrsnfPnEJsHqKDBPbIizn85pst6gEFu92GAjKIo4CtquPPhog7l1C\n//mfCC1gEgKsjEAntoGaWs03FxcnXMgHnRBgUzDCbn1/P9Dxsa2Benh3uxLEJcBa+n7U82COb0zb\nXbcnvaErAqye7Yl/kVHjdIn88fzaCFXWafkSqwh7tdsrnHJRkACrrcIRa1B41OkUnD6/4t0pF9xj\nGf19Fz99NbW+4NC8Adfb8OgH2jbvh+YN6KQNHxlyT9CGUgTYTT59tbhqXfr97eS7ouGbGsKec2q7\naTzui3cEWDH1j7HjuDnYbMwzZ99lXmW5Kfj2rcEG6B/oDS4SYLFtB9DB0mtg0w7Qw59yUVnKQ+9B\nT07ynvbbrCNTGXW6tIcN10MbGlrf/ru3PHn/9KL/dBg8wJ5jHVpxdbr09JSm8rbroQ2dc3xxhSnE\nmwpexmdRp0t7mMWyt3xKj3FKBcak2t6B5GZz+J4AAyAkU4gAhORRUuWl9KeuXddArmXu8TvFX/rX\nV2zehqcbsG1J8w3xev+uh98QuzXAOef6nbC+RMO9sYcNcX0nrHZ4FiHAStoesetOs359/M4TbVj3\n1NNXbN6G17+HSoVRo20ntGrDtqtbbYhdDzTZENuxu0kn/PquH23DpyNSnUOjIFOIJeWc227yHna4\nHtrQw7G3DNk9tKGHrmjehjk1H5GepgIbUPPxovnA3YPtiX/DBrRtQw92VU7bZlCWABtKJ8fqbsK9\nsuV1t/+tr/km6EfD+OxhQmy98re9BEgpphBHo/bKf70a9UYPnUAnluB0QvOQqecWHtJq8dtx0Vf9\nNlx5xalWITZc7tXDhtjVPdN2Qts2DLwKUYABEJIpRABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyA\nkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgAD\nICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTA\nAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABC+q91AwC4JKX/fZ1zu3Z0QwUGEMA2vY7/OycBBkBI\ng08hJmcpwAhOZgxTer1eDw5xuftpysED7NXBNkgptW1D8wZoQycN0IZOGnCvDcez8Zxfp8FWpA0h\nzv5NIQIQ0vgVGFDE7oy8dQ0TXV7782JP5mwV4p4KDLgjwgxTp26vJ8z5T25Jr0X7ueBH9TDZDQM4\nHWQdW/d82Zkp1ej5EIOnCgygHpVrQQIMoJ7uq5pIBBjwO8Puo3TvPQIMuGQ7yK6rCbhh13V68jYB\nBlxlCVwp1hMW4T4wgEe4betpKjCA8jw8vgIBBlCYuKpDgAEQkgADKMwVrzoEGAAhCTCA8tzsVYEA\nI6SUXCend+udXtLrIe4DI541upYvltHBPTcwGxUYwRwLr2M1pjiDGQgwxiTDYHgCDP7h6hpEEe8a\nWPo7uhw/LTRtBp7+P0uU3mxza3t1bQB1PsMXKgsWYNtPuT79xGu5Nbyc9xXS8TuvgbLnS6cLXmAM\no00hppSSCaDRbT+K4vQzKW4M06d7TfRd6XTBCwwjWAX2q6UC2xVqx19gPOti+ntbWBkH4c7+hwqw\n03CSWABXbEfLEGE2zhRiiO4GoJRgFVjO+bgKcZkwPP0RXLfuNcOs2Ttd8ALDCBZgr7NwWr8jt2Bn\n+5wtxweDGWcKMRZ3y1KZ9GI8AqyB7a05YgzgHgFWm1tzAIoQYACEJMAACEmAARCSAKvtuBjM8jCA\nG+LdBzaAL5/aB8BLBQaTsNiV8ajAYGQDf0onqMAACEkFNo7tubazbN5x8ZVhqMAGsbvC4YIHMDwB\nNiwZxinlF8MQYACE5BoYjMzngTEwFdiwjFbA2ATYIHZxJb2A4QmwceT8J7ekFzADAQZASAIMgJCs\nQixvvQGr8lSeG7+AqajACjs+O7XJ684ZZsu7nvO9w4QEWEnHobPOYGrIfolwOrPuhPbG5wgwRmCA\noCutzmVnI8AACEmAfWyZEOjqfMqNX8CEBNgHdrl1zLBjkLSKltkirdT7dd0CArGM/ivHzwZcn51a\nOUIuPrN1ng+9vPHuTq9bjN1LEJoAC+BKqXfj7ww2Orc6dYCjnPeHm93yCaYQp2aWrCCdydb22aTS\n6yEqsKtOhyf7Ja/NvuFjt8oqNffAqFRgvTsNTif7T7i3BscdP9CKAOvd6Rha6jzU+ezOtkN0Toec\nHLA1VIClv5744wUXarfiQy+v8JlqEMU418BSSvnvqLP9uqDdyqJPX2F3peTTf1WEpXqE8G63t9+y\nNU6A1bFm2O30IoqLJysWTBc36sHS6oOWBvZIpdLEaQX20HQiwPD6T4fxK7Cnt8E6HffzvNztVfgX\n51J2r35sTKBpw0+beuOt/VxGX99Y7176YpMabpS2K/5/eOMV7ld5uttPj76jzo/HEGf/Qy3iGFI/\nz1cE6Mo4AZZzXlchdlj5fpND1nZXppMhhKGmEDvMra3bC0Du/ZOpBJog5WdDLoo5vimKGKcCi2KA\no7Ef2w89KfsBKO6Za2WSgd4eVcRQFRhTeXqkM8R0IlB5/cNTMaO8hVhUYEAvQpdfnopZnwADevHo\nkz8ZjwAjKuMaTE6AMQ6RNgDLZ7jOIg4CM7oNKehmHfIGgM6pwADKsOywMgFWwHoHUtlbkYBwfJ5c\nTQIMgJAEWHmKMIAKBBgAIQmwEaxX4ADmIcDKq3n9drdsZJtk8gwYmwAroKu1s7s8E2PAqNzIXEbz\n3AKYjQBjLscZVycfEJQAC2b3gUMA03INLBKh9SUdCCNRgYW3mwHzOFFgEgKsnu3q9udyRWL94Pi8\n8Am5CsgwTCFWshs3DaMAXxJgNZSKq+OZsnPnT+X8v07bfj0JZ06MRIAF09VN03HpOhiAAPvKxYc2\nlR0ul7rho7/p48pYSG5GIsDuO14M75PLb8CQBNhNxxj4ORh2Z75OhGll8quAjMQy+nqajBTqLU65\no4ABqMAACEmADc4EETAqAXZToFuyXH4DhuQa2H2BkiBQUwEuUoEBEJIAAyCkeFOI6e/i33yYF0ub\ndcHHnwIwkmABllJak2n79UpuAUxitCnElFJyfybABIJVYL9aKrBdoXb8BXhOnU8upUPrYJNSyE0f\n7uy/6wD7NHtOf0FiUdPx0cl2wEkMsOm3o2WIMOs6wD7KntNLYlBThEN+apUrpIgZFkvXAXaUcz6u\nQlyi6/RHAItHKyTnLk0EC7DXWTit35FbwHUFM+z06f4GpKeNtgoRGjJgQU0CDEry6ORp2fT1xZtC\nhM4ZuaIovqVs+spUYMAUVEjjUYEBsxBag1GBARCSCgyms72f96UuISwVGMzleLuSm3AJSoABEJIA\nAyAkAQZASAIM5nJcsmERB0EJMJjONrGkF3FZRg8zapVbVvBTkAoMqMR6fcoSYHCT4fh7+pBvmEKE\nj22HXVNh0IoKDICQBBh85nTWy1TYPSpXviHA6EhKf5Jg/aJDp2OugfiKnK3gpyTXwOhFt4lFWXKL\nUlRg9GvgSAtRa0LnVGDwseI1RErqEviYCgxqU3JBEQKMfilKgB8IMHphiRrwEdfA6MsMuZWzWUQo\nQIBBA0tOewwVfMMU4uNS65Pt5g3QhncN2M2aNmlDfc3b0LwB2lCKAGNG6z1YQFwCjOkcnyUPRCTA\nmMsxsWQYBJXy0FeQB5jkpbTTHd5+Anv9p8PgAQY7p6c0DgKIyBQiACEJMOZyLLaUXxCUG5mZjsSC\nMajAAAhJBVZeSn+WxqxrINeVMsfvFH/pX1+xeRuebsC2Jc03xOv9ux5+Q+zWAOec63fC+hIN98Ye\nNsT1nbDa4VmEACtpe8SuO8369fE7T7Rh3VNPX7F5G17/HioVRo22ndCqDduubrUhdj3QZENsx+4m\nnfDru360DZ+OSHUOjYJMIZaUc267yXvY4XpoQw/H3jJk99CGHrqieRvm1HxEepoKbEDNx4vmA3cP\ntif+DRvQtg092FU5bZtBWQJsKJ0cq7sJ98qW193+t77mm6AfDeOzhwmx9crf9hIgpZhCHI3aK//1\natQbPXQCnViC0wnNQ6aeW3hIq8Vvx0Vf9dtw5RWnWoXYcLlXDxtiV/dM2wlt2zDwKkQBBkBIphAB\nCEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkw\nAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgC\nDICQBBgAIQkwAEISYACEJMAACEmAARDS/wNrlCRHKvnaFAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "subplot(2,1,1);\n", "u = x0; u(x0==0) = NaN;\n", "stem(u, 'b.', 'MarkerSize', ms); axis('tight');\n", "title(['Signal x0']);\n", "subplot(2,1,2);\n", "u = x; u(x==0) = NaN;\n", "stem(u, 'b.', 'MarkerSize', ms); axis('tight');"]}, {"source": ["Over-relaxed Forward-Backward\n", "-----------------------------\n", "It is possible to introduce a relaxation parameter $-1 < \\mu < 1$ to average the current\n", "iterate of the FB with the previous one.\n", "In this case, one must set the descent paramter to\n", "$$ \\ga=\\frac{1}{L}. $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = 1/L;"]}, {"source": ["The algorithm initializes $z^{(1)}=x^{(0)} \\in \\RR^N$,\n", "and then it computes\n", "$$ x^{(\\ell)} = \\text{prox}_{\\ga g}\\pa{\n", "              z^{(\\ell)} - \\ga \\nabla f( z^{(\\ell)} )  }.\n", "$$\n", "$$ z^{(\\ell+1)} = x^{(\\ell)} +\n", "      \\mu \\pa{ x^{(\\ell)} - x^{(\\ell-1)} } $$\n", "\n", "\n", "The regime $-1<\\mu <0$ corresponds to under-relaxation. Setting $\\mu=0$,\n", "one recovers the unrelaxed Forward-Backward. Setting $0 < \\mu < 1$\n", "creates over-relaxation.\n", "\n", "\n", "Note that convergence of the iterates $x^{(\\ell)}$ is only guaranteed\n", "for $ -1 < \\mu < 1/2 $, and that one can only prove convergence of $ E(x^{(\\ell)}) $\n", "toward $E^\\star$ in case $ \\mu \\in [1/2,1[ $.\n", "\n", "\n", "For a in-depth analysis of this scheme, see the book:\n", "\n", "\n", "H.H. Bauschke and P.L. Combettes,\n", "_Convex Analysis and Monotone Operator Theory in Hilbert Spaces_\n", "Springer-Verlag (2011)"], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 2__\n", "\n", "Impement the relaxed FB algorithm and display its convergence rate for several values of $\\mu$.\n", "isplay"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3dfZKj\nOBLGYbEx92r3yabmZO0+GfsHVZQKfSAJgTJTvyc2NmpcNsZyobdTCLGs6+oAANDmf6N3AACAFgQY\nAEAlAgwAoBIBBgBQiQADAKhEgAEAVCLAAAAqEWAAAJUIMACASgQYAEAlAgwAoBIBBgBQiQADAKhE\ngAEAVCLAAAAqEWAAAJUIMACASgQYAEAlAgwAoBIBBgBQiQADAKhEgAEAVCLAAAAqEWAAAJUIMACA\nSv+M3oFLlmUJH1zX9fk9AQA8jAoMAKCS7grMUW8BwKyowAAAKqmvwNzPM2EUZAAwiUVvjx+dweHI\nMACYg4UKbE+sLdKWZTk8En0yAEA19QHmB9K6riQWAExCcYCRTAAwM2YhAgBUIsAAACopHkLcTnf5\nA4nhIwCmkpqcjBLqOk/FAbZN2eDvFYBPXS8shMa+VPcQ4uEvdV1X/nYBYBKKK7ANiQUAc9JdgQEA\npkWAAQBUIsAAACoRYAAAlQgwAIBK6mchAsCcqm6FGL3MS/ssbiowAFAmXMNB42XI11GBAYBKmRsf\n5p9vBhUYAGiSiqsJizACDAD6CIf17guVw718C1+1eO7Zr0cRYABg3BZX9k6bcQ4MgHH3ddTqTiod\nTpsty6L6xBgBBgAPyRc9d2fJYdSRCgwApFNUY9yUcKrLrAwCDAAeYjVIRiHAAOAhpwVWW8IZGAxs\nwyxEANBkC7kwtDLhl1q5Q3tFSAUGAN3s8/r8wNjTom9gZAqvQz5tUzbsFWoEGAB0sMeDHx5bbHQv\ndMI5hKdvcXiJ9tprQ4ABQB9hKtyXE/ktR39rI7R8nAMDAKhEgAFAH/ZOMglHgAEAVOIcGAB0YO8M\nk3xUYAAAlQgwAIBKBBgAQCUCDACgEgEGAFCJAAMAqESAAQBU4jowABDH3sK7d6ACAwBBUvfuQogK\nDADE2asu0iuDCgwApEjFFTEWRYABQB/h0F9b8PgnvTgBlkGAAQBUmuIc2LI4/hEDTGv5uGv8bf2g\nZxnJeID9+vXr9Xr9+uWce4/eFwCzy48oMlpYy3iAAQB1klUEGAA8hBqrLwIMAB7SMITIBPoMZiEC\ngBRbhoWhRekWRQUGAN0syxKG0Pbzuq7lOUThVYIAA4AO9sjZ48p9XctcVT+t63pIL8qvlFkCjEvB\nANwtTJq27CGxCnEODACgkvEA+/v37+hdADALTlw9zHiAffuXPywAMGWWc2AAcCtOXD1vlgrs1/sX\nxT0AWDJLgAEAjCHAAAAqTRRgv369Ru8CAKAb+wH28cFMegAwyH6Aud9/9gxjHgcAmGE/wH7/+f6Z\nUUQAMMN+gDn3M8QAACZMdyEzq/oCkM9fler0EunoElYzXFhtvwJbX6/ff5jKAUCH7Q4sh0dG7Yxw\n9gPsgNNgAORbv9Q+v/aFqk0TYJwGAyBeqtiiCIuaJsC8C8L4SwBwh3Dory14/PqpvJZaPA1vqtFE\nAbZjFBGAGVtczXnabIoAW3//ZgQRmNey3PU/SQ5nv2bIsCkC7JO3JAcAPG/JurLltlFH7WYKMM/r\n9Rq9CwCesq53/U+AeeYchqa7kBkARpk2aW4ySwX2eRrMOxU2wfgwAFkahhBnOJXVbJYA23EaDIBY\nW4kWhlamdEut3DFDtTdHgMW+SCbTA+huzxI/VPaMCdfLSK2dkSnLDo/vmTfbRWBu8nNgLOwLoJdD\ndPm5UlUMret6CKHTlx9eMkPttZmjAvtyOA3mOBMGoJ+wnGqbIphf1TDz4GwzEi0EWGHtvP7+vf1w\nOA0m73pEAMA53QGWOnsZcXYajBgDcNFU558ksHAOrGrplN9/3I8xxH+/XvLfNmC9b7Pf/gEAbqC4\nAqu+bGJdP0cRvTWlfr1/ff723+U7zH6unQYAp2Y7/ySB+grssALYSRG2ru79zj3hZ0G2OWySP1EA\nkKBufqco0Yv1Dg+GE1iX99s59+e3c39+f3x8ll/v93v5SCTff6cTWOt3HcA9auesY6ex6dRXYKfC\nKyqWsAjzvrljknnjitEwoz4DgCH0Re6urQJzzi3v95/fzjn38etjf/D9969zP/InWZNtziqzndoG\nBvTRWEYIobHp9O3xrjnA3J5h3iii2zPMHTPnJMmcWz/WqrkeapsckE5jLyyExqbTt8e7KwHmYhn2\nHWC7sHo7C7NPxfXZ2RsCqKCxFxZCY9Pp2+NdGGD76a6SAHPOvZe3+/PbOZfLsK8tRnagMMw2rZGW\n3QUAP2jshYXQ2HT69tgX1lsuiLTcbQi+ijDnZ9g2xSMzJpgq6aryzHc523brygrFmJrGXlgIjU2n\nb4990au+ygNsE6nD/GmKhWe3rpdop/rlXIjkgw0ae+Go2tXlaxewj25BXdPp2+OD/NdW+JVsGeZP\n6Pj78eGcW18vf1st+9d8Fu2iO9POtyef///AKBp74YPTf5d3eUl0I+qaTt8eVyn/SlIZtvsRZl9b\nv7RzmbHNZ0LO99/6ecXbU8mXsqWgC6IRKKGxFz44PTMSfX74kvyrottR13T69rhK1VfyXt7u58Vh\nhwzbRcLMe8vCtztXtucD0u5gyzwZ+Zfnx2H4AwzQ2Av7GtKoZD524Vurazp9e1yl4St5L+89wz4+\n/rrff37eAjMil2c/96ZqT0o1DHYPz7yY9WP93LGtHJSdhYX2jPR/3r8xsrO7sb3w4d3bUiR8SX47\nBJhZzV/J6yuTtgzbfj5Nsk1pnh08s+795a9bZvjl+KOjWy5aScdCmYrTXoJOGGDhExre141uujb6\n9rjKla/kR4a57/h6rS/3tShwicY8i3og5J466bR8LHvJ9V17KbHtsP//7r81zAk8H5mp1eNuervD\n0T0qwPbnfO8YsxAN6BJgLsgw9xVj3280JM/ynkk7YTTmYhiH60dD7xMfrpyP6ADL3+9pvyFUWwUW\nbi2775GNqIsDfXtc5eJXks8wF8TYj7cuO2aey7OMu/s55X9jfgoKD8UuWdggjM9x2Sl6CPGOACtZ\nk6hh51XQt8dVrn8lpxn2+bR0kn3uicASrRY514OBGnH07nxKFJ2iA6xkC+FLmMSRom+Pq3T5Sl4/\nE+Xj4+/rtW5z7o/PPIux7x0zkGdRU45bjqIuCx8oEIUHWEMF9uQ0eud+PF/+oUaAFYlm2PZzmGTl\nMbYzm2cp+yXK9zH9h32T48wUwVJZqD3AwlfVTkF0BZmX3g4BJknfv+bDcOKeYa5TjPlKIs1ImOXd\nF3Km//IlGBOE/7VMwOvFzxI/RarqoZIlXqOPZF5S+L4/t6ngECHA6mQyzN0QY7vTPJsizFLuCDnT\nx4U6FcOk4wIsOpG9+ZKscFOH30YfTL2k8E3VxYG+Pa5yx1dyuD7skGHuzhjbEGYtHl/iC0MM7IU1\nBoBP4/7r2+MqN30lfh32TsTJ3TG2y+cZYVaqS8KZPppUIMCaadx/fXtc5b6vxK/DwiJs91iMOZLs\nPmSbHmMDzDWN3QlBgIlz61dSUod9/vZnjN2XYTvC7CG9RiZNH4ZPIsCaLcvi/v38Wdo1fykE2CWv\nRBJE8+z5GNsQZmOQbSNoLCOE8ANsQ4AN9sBfcyrDNocke3JEMUSSCXI93kwfuc0IsGYEmDgP/DW/\n39890eu1RvNMVIxtMmFGkg12252+Z0CANdPYdPr2uMozX8meYftsjjDGwkHFUSOKoVSYkWSyXAk2\n04e5T2MvLITGptO3x1Ue+0rCDPv6z5f3nHfkhV6MDcywDUmmVVu2WTz2NfbCQmhsOn17XGV4gDlt\nGbaJJhkxpkxtqpnoCjT2wkJobDp9e1zlya+kV4Y5MTHmYklGjOlWlWoKOweNvbAQGptO3x5XGRJg\n7izDnJ5SbEOMWVYeaRr6Co29sBAam07fHld5+CsxnGGOGJtEYZ5J7Tc09sJCaGy68z3uslD/KM9/\nJfkMcwUxJjnDHDE2G215prEXjvI73pJPVPv86BbUNV1yj/M3XtsJ/8BDvpLMybCvx1+Fm/r4+yEw\nwxwxNq2SbmFon6CxFz5oqBnuuB+YCpE9Lkzy64H/gFFfiV+HuYLhxAyxGeaIMcjLM4298EHbHZmr\nXpLajrqmO+6xsU8uJMBcuhr7+u0r81vJGeaCGCPDJiUjzMT2RYVSaeTS3XLm/pZVTaGx6fTtcZWx\nX8np+bC8H3d/lp1hjhjDwWme3XNgjj3kD+/eliLhS/LbIcDMkvCVXIkx1RnmiDFsHizOood8uPpo\nL4dDkgB7mL49riLkK6kdUfz5zNfPTb177NGNiDGcuLM4mzbAXM2oY2o7EnrLKvr2uIqcr6RjhjmF\nMUaGIa53mAkfQszP7l7Xta2cSm2WANNN4FdyOkExhRiDcT3CbM4AO2w5tZ2qnVdB3x5XkfmVNFdj\n7+X98evDf0TdiTEyDKUyfX22FpEcYCVbCF/SazunLxHYW+bp2+Mqkr+SthjbRvP9GPv4+/nznmSn\nF5k9X7r5MUaGoVpxmAkPsIYKrPmEVsNLxPaWKfE91vhJouR/kDDGXFmSHSYofv/8s0TLvvW78Jld\nUIqhg/ww47pqD7DwVYUnwC5WfuHOq5ALsLZWEEXRV9JQkIUZlkmvPa7GTmtkjiK6iYXBMrTL8vtM\nv3Kq6ktP14UqycWGRlDUW+6SAbb9oO7zHKj7Snot4eHXZJ/P9E6VjZ0Pwn2f0ZnXfQ8MsEOKRMup\nK5s6/Db6YPT55W+qq7d0YYCpWOGwnMavxFXOVExVVJlrX7YwGz6tMZVkG/IMDQYe8kp7m53G/WcI\nUa7yGMvf8TmfZMNjzJ0l2YY8QwkCrJnG/WcSh3SFK1FtOZTPnlSSCYkxV5ZkG/IMUWMDzGn+F7/G\n3lLfHlfR+JWEmq99Tm4wlmRyYmxHnqEWAdZMY2+pb4+raPxKUi6ubR/fZpBk0Rhzo5NsQ57hlKVD\n/mEam67ofmC9rjN4nsavJOOODHPaYsxHpOHA2CH/pGVZPhtuXd2yjL2zdiECTJmbMswlYsz1uG7s\nsZHJ8jzzkW2W2DvkH/MdYDvxLUmA6XNfhrmyGCuMn+ErWrXlmY9sU8fkIf8MAkwcq3/Nt2aYKxhU\nPM2e0/Qq31R314NtQ7wJZPWQf4DGpiPAtLo7w1wsxvylqjLBcxp1Yu/S2SvbNiTc8wwf8nfT2HQE\nmGIPZJgLYuw0w/Zwqq3S5MRYSt942xByfdk+5G+lsekIMN2eybDP9/pKskyGladX+JKqVwl0R7w5\nEq6S+UP+PhqbjgBTb8+wuwPMnWVYQ3odXtjwWkVuSjhHyHnMHPK1y9Lmn3+6wr3T2XS5xXxPyf+0\nGr+SBk9mmIvdVDPYn3ftNifJsFOE3EUGDvmSsKl9CQEWIf/TavxKGjw5kPj9pokYuxI/N8XYvlDk\n6/UykI6EXIaBQ77thpb5l5Q0i8am07fHVTR+JW0eLsK+33d57xnWJRs6ZtjwC9EGmvaEnPZDPpVG\nLp1hJed9CDCVNH4lzQZmWPigf//MWhczrPwStOa3sMHkpMqxh/zh3RvmChTOQqh9iX+HrE30PJm6\n3lLfHlfR+JU0675ofcVbp2855prCLL/6VDgSmAktXReiidI94R6INwIsfDB1YsjgObDcU5XMPPRp\n/EquOGSYe7wac13DrLaW+rEbBZ0vM0cu6phwvbItesiHx0Uvh+NrYIC59Khj4bCkxt6SADNoyJyO\nyG5kk8yVhdnd44GUYvcZsmSX8ADLz5LbR/mqAiyzWT/ASi6OUtdbEmBmjTolFtV9jLEvSrHn3bTU\nsvAhxJsC7LDl1Hby29TYWxJglonKsJ3MMCPD5LiUbb9/Sw6wki2EL+m1nfwTNPaWBJhlQsYSMzJh\n9nyStd0yZru8rOqNCMg2RcEmO8AaKrDTafSpPcm/hABTQONX0peEaR0lhCRZ7Tr6F5FkXfwINuUB\nFr6qLYryW0jVeep6SwLMPi0ZtkuF2TNJ1iWiqsKPGOto+DkwF5v7V9V5Fq4Llc/FkiWWCDDpNH4l\n9xl4oVib7pdIFyrJsL6rhJBhvQw85A8hkZkB2Lapw28zERV9u5InqOstCbDpqIsxJybJuscMM0e6\nGxtgqnsbjftPgM1I3aDiJoyx4fPvryPD+iLAmmncf317XEXjV/IYpTHmgiQjxrAbPoSot8PR2Fvq\n2+MqGr+ShxFjQlRl2MwL7ecRYM009pbckRnO6Twx5ubLsLYZklOFGYd8M41NR4DhGzE2XGraSOEU\n/HzCzZBkHPLNNDYdAYYj+et3hAxnWKh8lZCLW9CIQ76ZxqYjwBChMcPczxhTnWEuET/NwdN3a5Jx\nyDfT2HQEGOLIMJPMX0DNId9MY9PpDrCSNVeG76ReMhezP7VnGAGWYjjGOOSbaWw6Agw5ZJhhJi8+\n45BvprHp1AdYfh80fiWiKB1IdGRYmb4ZJqGw45BvprHp4gFWaPinJcAeoLQIc2RYmS4ZJmfGI4d8\nM41NZyHAMqssa/xKBCLDbLuSYYXXVj8WY2YO+dPF42ufP/tq9NKkspZzYN0ZCDBHhmW1ZVj+VUNG\nFA0c8iX37qp9idX7gf1v9A50sH7Z/vPwVS0/jdhB9fbcChdOFM4PrcxNn9GQLqeZ936//cf73sza\nvEO3luHfOTP6ktMn6KUvcncl8000/ptCJr1FmKMOK7ZnTD7Pasu1J6c7aj/k/bBJPRJ9SaYnLJ+a\np67phFZgS9b2HGP/lBBObxHmqMOKnS696JrSaJ46LDr807Adv2ejl8sQGmBAX6/1tccYGVYimjTN\ntZSZ68wgir6aMY8hxFupHkjcMC/xVCqlro8EFg5RVm3q4O/fv+Ehf1/Zd/gghw6n4XrZwuG+8Lcu\nPepYOCypsbdUXIGF5bmQy6uhAnVYSjS0upzHKhmiPPV6vfSOQ5acHKnlx1J0O6dP0Etf5PpKJo+q\n/oACbUWY3grMMaejTOEdyK5s9vo1Z2HWRiuwx5xWYPnw2K9qrarAwi1Ht1PyBHW9pb49Pjh8K+Fv\ntX9AaQyMIjoyrEw+MK5vtu2Cs/wLxx7yQ4YQy7eTf4LG3jKyEkfDZxD7ycXumF42AsxxMmyo8gyr\nDVHhAdZQgZ1Oo0/tSf4lNgLsn/Ch2rS3NKKKebzW15Zh7+VNhkkzzx04D7YMC7OkNopcMGXjYmko\nUzxy8+Ny5c8ZTuO/KeQzU4Q56rBxUkXYxegaXoG52Ny/qsAoObUffST1/JInuNFN1ya5x4V1lfAP\nrPErkc9kgDky7HFhhl1fPnHgIX/oMzOrYLRt6vDbTES1lSUae8vIEOImtbTg4beAavtAomMscZyb\nZos8L+wY27rK/Kuivz19I5Od9vl1YGvMA3sGsQwUXj5Ca5RUSilNLzzvas0ovD7TWBSrYGkUccPJ\nsFH6LvU7fAhRb4ejsbe8tMepDyynIeTsiTH2AsyRYSYQYM009pbtS0llPu368y7JsMdSboVYZQoN\nOL3yPMVrIUICjXdXSeGuK4AuyVmIQLNUqsmv2/xJiQCE4xwYGqVOg5XUZMKTjJNhenHy4gp1vSWz\nENEuujJ9NNjU1WRk2CROb37WMDGy77zKEsvZu6xqb0CTZ7x/J8BuFWbV6ezEMMlkZhgrdMwjsyBI\nc/xcX1LkinyeWQoz4/07AXarQ1xVza0/JJnAGCPD5uEn1q23jb7pVtS5+8ukf2UgyYz37wTYra4E\n2GELVa96DBk2ifsWv8/cPLpq+yU3oT650YzFJDPevxNgd9tPgzVf2kyGQYK7B/2akywfrm1FXjTM\nNMaY8f6dALvb9QBzZBgE6DhyWPhGB9H3LYzVtvRN1WSKkqy9fz+driohOQiwu/U6lSV8bSoybAZb\nDDw54SJVXbXVVc2zH/UWZFevA3MygiqFALtbx7kYZBjmlD+/1XyqrDaJwxiTn2EdrgOTnBDCd8+G\nXsEjPMAcGYY79ZqyeP0qtD3J7AeYcATYA6KXM1/ZVK+t3YEMg3yPnc8bznj/ToCp0zEOb3JYLJEY\ng0CTZBir0UMiyYvcv9YX69ZDuMw8+xKv16vthQ87KVBKVsaUXOJQgWkkfyBxQykGyRpOho1dAatB\nhwqM5Z8xp0NiUYpBlPL4eX25cW/ucV6ghHPl/Uf29JJZ6FCBKaWlCNtQikGs05NhdyzY+JiiIcTw\nOWGGycwJAkwpXQHmyDAIlsowdQOGIQIMQqnLMEeMQaToyTAD6eUIMIilMcA2xBik6bUovjQEGOQy\nk2GOGMNo990yZqDSafThJI7tQcnp5Qgw5fQG2IYYgzTPL1h8q9JZiCH55ZcjwPSTvzDHKWIMuElp\n/+7HmKJIIMC0016E7YgxoLt/Cp9HDABXbHHlx9j2MzEGNCsqUMJRRC15RgVmgJkibEc1BnRxvpRU\n9BwYy0cBzQ7LATvn3sublaiAWldnId63Z11QgdlgYCpHCtUY0KzDdWCSEWA22BtFPCDGgAbcDwwY\nj0FFoAEVGHQwX4TtqMaAQkUV2GHKBjM4gPtQjQGFTq4D2xeL0juTHtCI68aAU1wHBjVSo4j74z5L\nI40MKgJRxvt3AsySMMCi0eUjxgDDjPfvBJgxfoYd0iuVapYyzBFjgMd4/06AGVM+Wug/01iGOWIM\ncM6FAVY1w1B+NhBg9pQXWLYzzBFjmB4BBsvMZ5gjxjAx4/07AYYZMswRY5iS8f6dAINjFQ/AKOP9\nOwGGzTwZ5ogxTMN4/06AYTdVhjliDBNgNXrAJtZUhHnGCxQqMPhmK8J2VGMwyXj/ToDhYNoMc8QY\nzGEIEZgFg4owxniBQgWG0MxF2I5qDAZQgQEzohqDAcYLFCowRFGE+ajGoBQVGDA7qjEoZbxAoQJD\nCkVYCgUZtKACA/ADBRm0MF6gUIEhgyLsFNUYJKMCA5BENQbJjBcoVGDIowgrRzUGaajAABShGoM0\nxgsUKjCcoghrQDUGCYz37wQYTu0B5siwSsQYxjLevxNgKEGGXUGMYRTj/TsBhkJ+hjlirB4xhucZ\n798JMJQ7ZFgUwZZHjOFJxvt3Agy1iLHriDE8w3j/ToChF8YYaxFjuJua/j0TRcvy3bMcnkOAoS9i\nrBYxhvvo6N+3iAp31Y+unf80Agx3YNZiLWIMd5Dev2eqKxcEW5hzBBhuQoY1IMbQl+j+/VBghcOD\nh8ejj0j+gFCNDGtDjKEXHf17dAix5EECDLciw5oRY7hOR/9OgEEsMuwKYgxX6OjfrwTYYVMqPi90\nYTngi4gxtLEfYCo+ILQjw66L3pmFJEPGmP49Ov19dzrbsPxBAgyPIcN6oSBDIQIM6IOTYX0RYzil\no38vDDCm0WMsirDuiDFk6Ojf8ytxcCEz5CDD7kCMIUpH/85SUtCCgcT7EGM40NG/pwLMBRnGYr4Y\njiLsVsQYdsb7dwIMQ5BhdyPG4Jz73+gdAIBqr/V1SKz38o5eSQbDjBcoVGAYhSLsMVRj06ICA6Ab\n1di0jBcoVGAYiCLseVRjUzHevxNgGIgp9aMQY5Mw3r8TYBiLImwgYsw84/07AYbhyLCxiDHDmMQB\nwDKmeBhmvEChAoMEFGFCUI0ZY7x/J8AgBBkmB3fONMN4/06AQQhmJApEQaad8f6dAIMcZJhMxJhe\nxvt3Agyi+BkWItUGIsY0Mt6/E2CQJp9hG5JsFGJMF+P9OwEGsU6TjBgbhRjTwnj/ToBBl0OqkWED\nEWPyGe/fCTBoxHQPOYgxyYz37wQYlCLDRCHGZDLevxNgUI3Ln0UhxqQx3r8TYNCODJOGGJPDeP9O\ngEE7xhJlIsYkMN6/E2AwgAwTixgby3j/ToDBBgYSJWN14FGM9+8EGMwgw+SjIHuY8f6dAIMZDCRq\nQYw9xnj/ToDBEoowRYixBxjv3wkwGEOG6UKM3ep/o3cAAMx6ra9DYr2Xd3TSBxoYL1CowGAPRZhS\nVGPd/TN6BwBgCltc+TG2/0yStTFeoFCBwSSKMO24dKwLKjAAeNqeVWFBRoyVM16gUIHBKoowSzg9\n1oZZiAAwGJMV2xgvUKjAYNhWhFGBGcPpsXLG+3cCDIYximgYMVaCIUQAECccVHSMKwaMFyhUYLCN\nImwSzPKIYho9MIX3eyHk9MpcBO0mDjPjBQoVGMxLFWH+7Vd8xJh2nB7bGe/fCTCYFwZYKrp8xJgB\nJJnx/p0Awwz8DAvTKxVsZJgNM58eM96/E2CYQbTkSuUTd3Y2KTU70XaYGe/fCTDMozyZyDDDpkoy\n4/07AQZEkWG2TRJjXMgMzMgPrZJJH9AlFVRbsJm5INp4gUIFBmRwHfQMMlmlvSAz3r8TYEAeY4mT\nOC25NIaZ8f6dAAPymFs/m3yS6Yox4/07AQbkVU3BhyUGksx4/06AAYUYS5zTe3m/1pfSWYvMQgQQ\n994mq72XzDTF/G/V2T9Lr8/lbyS/wSdb0n+vLaL2u7f4iSU8vRwVGIBdpgMtWddjf1p+feEh5V3J\nW2dW4brypu7nEl9VK6SkdruwXD59eXL9zN9/HAE2HAEGVLleBKTWY6xdKGR/SWHqhO8b9s7RTeU/\nsh/JqSQoTPdwP6NtldqH1JbDS/oyyzqXfDvRd5SJ+4EB6CnaI+cHIV1Quvk/h5117Vv4zym86Uzm\nOafvntnJ/JYLpfan4aNpZ7xAoQIDqljq9QrrG11KPlSvDy6/DjPevxNgQC17nT7ayA8wZiEC+EF+\nt6UdLdwLAQbgqrBH9h9p7q+n7ehv/eD5jeta5ZlJHACODjMADzMUMot3pGblVZ2VybxFaiOZbjf6\nEcKnXZS/F3bm7e4O6eYdU4EKDECdrfvbO8FUF3x4WtXG237b8Y2qnnlHiZlpvdrqtrlJ5VfAxuc4\nMIkD6OL6jVcOV3e5WHGQuU4rfwlX1bW6p5tqu+A6f6HbxQZMFcHyM+ZWVGAASg0pgO7eyK0bvHXL\nk6eXowIDYFJtBfbYLl0sYa9swR4qMABToN+3hwADYJC9uLL3ia4zPsLGECIAWEUFBgBQiQADAKhE\ngAEAVCLAAAAqEWAAAJUIMACASmpWo49OiF+WyDrKzJsHgB21RwMAAAMESURBVBnoqMCiQQUAmJn0\nCuw0uqi3ynFZ946m2NEUPlpjp6IpRFdgFF4AgBTRFdie//kk838r/58MAIAuRFdgeVtuHbKNog0A\nJqFglNN9xdJhV/esCgu1wtINAJAiPx1EDyGW8Jt4XddDYsn/AgAAbcYEWL4wKkwdwgkAZqb4HBgA\nYGZjKjCKJwDARYorsGVZolMQSUcAmIHiSRzblA3mGQLAnBRXYC4ottZ1pfwCgEnouA6szTwrdJRc\nPHDaGsaaKzOePE9TXP+kJpvCTdkambUNSz6XzMaxGWBT3Wbl9MNef4JG+YvffVabIvws0/5VcIyk\n/j1X8rkkN47ic2CnZliPI7X4SPivrdPWsNRchRcaGm6KQ4fV/EnNNIWb9Rgp2dWSzyWzcXSfA4tK\nNZ+iv7lah+VI/F+dtoa95koNlczTFOE/t7efyz+pmabYzXmMlC+DnnpceOMYDLBN5u91QqetMU9z\nTdsU4RSnaZsiymRrrF/yz4n+XPWcUY1jeQhxBloOpMeouAvfY1RPOuhlv97mMMY1bYNYYrYCm9bM\nB6f8IZ0ncaeh3T6Cul85OucBYg8BZgoH58yfPXQYQZo2w8IPPm1TGMMQohFEF4OHB/k7Dc0jPDRS\nsxChDgGmXjhLeFqpf2jTMvDNHOfGTBFg5v9Yqzro09Yw31w7mmJHU/imbY2SzyWqcQyeAztc8nJ4\n3JjTv5XT1jDTXOtP/oP7z26apnDZq3nmaYooWmNT8rmEN47lCszqv5JCzevBVD3BjHmaItpTZ55w\nugWN1sRtK+ZsjSh1hdfOYAXmYn+a2v+tdMVpa8zTXPM0xfVPaqYpHK2RVvK5JDcO83AAACrZrMAA\nAOYRYAAAlQgwAIBKBBgAQCUCDBgjOrcbQDkCDACgEtPoAQAqUYEBAFQiwIAxOAcGXESAAQBUIsAA\nACoRYAAAlQgwAIBKBBgAQCUCDACgEgEGAFCJAAMAqESAAQBUIsAAACqxmC8AQCUqMACASgQYAEAl\nAgwAoBIBBgBQiQADAKj0f7cdp7uXnslLAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["FISTA Accelerated Algorithm\n", "---------------------------\n", "It is possible to use an adaptive relaxation parameter $\\mu = \\mu^{(\\ell)}$\n", "that changes from iteration to iteration.\n", "\n", "\n", "This strategy is used in the Fast Iterative Soft Thresholding (FISTA)\n", "algorithm introduced in:\n", "\n", "\n", "Amir Beck and Marc Teboulle,\n", "_A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems_,\n", "SIAM Journal on Imaging Sciences 2 (2009), no. 1, 183--202.\n", "\n", "\n", "In this algorithm, the relaxation parameters $\\mu^{(\\ell)}$\n", "are automatically set and tend to the limit value\n", "$$ \\mu^{(\\ell)}  \\overset{\\ell \\rightarrow +\\infty}{\\longrightarrow} 1 . $$\n", "\n", "\n", "The step size should be set to:\n", "$$ \\ga = \\frac{1}{L}. $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = 1/L;"]}, {"source": ["The algorithm initializes $z^{(1)}=x^{(0)} \\in \\RR^N$,\n", "and $t^{(1)}=1$ and then it computes\n", "$$ x^{(\\ell)} = \\text{prox}_{\\ga g}\\pa{\n", "              z^{(\\ell)} - \\ga \\nabla f( z^{(\\ell)} )  }.\n", "$$\n", "$$ t^{(\\ell+1)} = \\frac{ 1 + \\sqrt{1+4(t^{(\\ell)})^2} }{2}$\n", "  \\qandq \\mu^{(\\ell+1)} = \\frac{t^{(\\ell)}-1}{ t^{(\\ell+1)} } $$\n", "$$ z^{(\\ell+1)} = x^{(\\ell+1)} +\n", "      \\mu^{(\\ell+1)}$\n", "  \\pa{ x^{(\\ell)} - x^{(\\ell-1)} } $$\n", "\n", "\n", "For this scheme, one cannot prove convergence of the iterates, but one\n", "can prove that it reaches an optimal convergence rate of the iterates,\n", "namely\n", "$$ E(f^{(\\ell)}) - E^\\star = O(1/\\ell^2) $$\n", "while the convergence rate for the usual FB is only\n", "$O(1/\\ell)$."], "metadata": {}, "cell_type": "markdown"}, {"source": ["__Exercise 3__\n", "\n", "Compute the solution of L1 deconvolution using FISTA.\n", "Keep track of the degay of the energy $E = f+g$.\n", "isplay energy decay"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAdN0lEQVR4nO3d65Kj\nuLIGUDix3/+VOT+YojFXcVdKa0VPRI+LwjbV5qsUidR2XdcAQDT/9/ULAIAzBBgAIQkwAEISYACE\nJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgA\nIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgpMID\nrG2btv36RQDwgMIDDIBSCTAAQhJgAIQkwAAI6X9fvwDgklafEjfpuu7rl3BMFQHWtk20nwscEO68\nQ4Yi/iZkCBGAkAQYACEJMABCEmAAhFRFEwfwkO0r/313ydo2ek+4SAUGfKNt24idb+RDBQZclVJL\nTbYRXVynAgM+YPyQ6wQYACEVH2CGKQDKVHyAATnqr4EZSOQKTRxQmkfbIxYTZ7EjQ9cGT1OBAZ+R\nalyhAoPSvD8sd6KNvvlLr7ZtDSRyjgoM+Ibc4iIBBkBItQSYkXbIjQtgXOQaGPAGU/pyu1oqMCBD\n0osrVGDAeef6D+EWKjAAQio8wLqm6UyHCFCiwgMMgFIJMABCEmAAhCTAAAhJgAEQUvj7wCa397vj\nBKASsSuw+eQ0ZlcDqETgABuyqvszeRyAgoUfQhyPGXZdJ73gTdufuP7j+e00vv2zu7hQpPABlqJr\n2rbxzxfyIlp2jbN/90ClbHxoh/kLHGAFHH0ow4kpfQ2WbFu8wL9xnCfbL/5ykLJNLIEDbG7751HA\nTwuKYcA/ReKl/XFDwNr2KduEE7iJY2I3n8aNHgDZWouWlCuO47/Pt0/ZJpASKrDE0qptG/kFcQ2f\n9PlH/tD9oBun9bU9H/rddzLWd3rsR4fartgBNi+KgfbJJYS6m/qhTp/WU+7+3LhctBgDw/Z9Tgz/\n67pD5mIHWOPfFmRgMRUe6tpY+8inX9optZ1kHNsbg5AlXQMLHGAFHH14wl1F0gu2O+sS99Asjbat\n7Xn36YZvT9z+qPTrWOmGAcaNnadsE07gAOul/OoHPOpEG33z9+Htk+biaf3oSXn36cINHk6O4eIh\nTdkmlvABBgT1yQl0twHkac+F4rlfI0ILHGCF/SSgZhc/zonfntyx3A5/yWQIMXqp9JBy7gMDYrnx\npDyfYyJx52s9jcNto58kR8pdXGPz97tYaJqJA+CMh6b0XetNSJ91aXhwfulrPsiZXsPNQ+jcLWXb\nX5p0/483mDxRyjbhFF6BtbF/OlC4W06g853sttofegGHirDx8GMfY0OSHZoM6OhLncfVuW1iuX+E\nNytt03Zt/5f+95SPXw/c7onrNGzYPuBxfxwRX3nhFdiTMxIA8KXSAwzgPillio7B1wgwgFThBtnK\nVk0XYtfq6ACeJuHepAIDICQBBkBIAgyAkGoJsK7t/wOgELUEGACFqaYLsWmapmlbk3HAnVLmVl+c\nA3B7ptqUW6nWdqgPsB4qMOBta9PpugWYQ+qqwIAnnJthfaMmW/xS+qJfirBKqMCAD6zNtn6uCCtg\naStOEGAAhFRZgOmkh4wdWjFrbQ+N6XSr4RoYFOfR0/flYbrJ0sDXx/2MH1ZLgAFXLVY8G4kyrpM2\nejdgmwCD4gSJgfmNX7c0EA4VnjgsngADrroldZq/W8HONeW77lWhCpo4/BYGOXnthmWRVrzqKrC2\nabtGpEFR5kWb9KpBBRUYkJ95wJyIHP2HlauuAgO+1TdZNMd7F0/QylE2FRjwtsUblq/fxTzZ2127\nIls1VWBtZyYOuFdKTqxtk54x5/Ygw4qnAgMgpBoDrG3UYQDh1RhgABRAgAEQkgADIKSKuhC7ZxeZ\ngM+YdYI6VRRgTaOTngJpFqdahhABCKnSANNJDxBdpQEGQHQCDICQBBgAIQkwAEKqI8DGfcatnmOA\nEtQRYAAUR4ABEFK9AeZWMIDQ6gqwTmgBlKKuAAOgGAIMgJCqCTCd9ABlqSbAAChLdQE27uPQiAgQ\nV3UBBkAZSliRebye+tbqtF3XDFtamhkguNgVWNu24/RqfsNsjbvBAApQQgU2VF0p6TXRNm3XaEoE\niCdwBbYWVydiDIBwAgdYb3zRa+sC2N8Wz74aAN4SPsDO6fqxQwDCKuEa2LbJiOKQWt1fL0fbKswA\n4ik/wBbGFScXybq2bfdHHwHISqVDiE0zjTGdHwCxVBxgTTNpoJdhAIEUFWCpDfTz4cK/WTnaVowB\nxBA4wPrLVvPQOnA562cOqsWHAchU4AAbtH9OfO8witg2rQwDCCR2gM2LrdTya9hs0mT/m2FiDCBb\n4dvor/e/d03XLwzWz4s4nrO+cZcYQK5iV2CXrBRh/VeUYgCZqzjARsZXwv49+Ft4iTGArNQdYOtF\n2PD1xRiTZACfqzvARjZWBZvHWCPJAL52wySAi/3rmcwtmDTL4d/rH6an3wiz7cTK400DVOF8gCXe\nd/VtkqVO0/ubYSlrNO++e2EG8KgzAdb+TGCxUawkbfaoowHWHMmwpe9eJc8A7nU4wPpYOvFdn2TY\ngef9zbBDAba+p2XCDOC6whfCOhac92XY0i63FP1DAHhE+Jk47jSahKNrmzvyK3WO4PmXRBrANgH2\nazyR1N2zSM1vKdswXTVangH8EmAzT2bY5HkmDpVoi3sAqIcAW9B2wwqXr87meyjS1r4q1YBKaOJY\n+cam7cbxkM1ROjf3RzYvH+A2KrBVX9Vh2xZfxW6qrW2Qx3sCOONsgfLRfV1HXXmd/cz0XTO6JBbh\nLQ9Oz9MY6l0C9bpUgZ27qTmWts/4Pg2yqcNSrL3Sc+VanPcN1OJ8gNWQXv/EzLBFNw5CBj8SQGyH\nAyyHGQ7f0TVdP4rY9vNyDOJn2NxdqVbcgQHydTjA+tDqry1VWoQ1ZWbY3IlUU6gBr7l0Daz46JoW\nYeMMq9XRm9XWNij93w7wuJMBVnx0rSroYthdXFQDPuE+sB3jIqyZrBMmw9ZJNeBpJ5s4JhXY/MGS\nL4/VdzHsLlINuJEKbF9fdbkY9hANkMA5AuywaYYpwh6gARLYJcBSDRfDFsiwV9ySan5QUAwBdoCu\n+gxdX4NGpEFQAuwaA4n5OVqoiTQISoAdszWQSMYmmXSoSpNnkCcBdtK/CRIVYQEdGniUZ5CnkwHW\nLn3cFx+EENIjTZ5BJlRghy3MUq8IK1HiqKM8g6+cnI0eaiPPIDdt2YHUL/vyyJ7/Wjn+zY44nLqK\nPqQsShw+908DbmQIEW5woj4TZnDRnQVKhhP4PleBNYowkqXUZ/7JwFEqsBv86+aAJeNwcvEM7iLA\nzhvf1PzfX7qm6x/QjsgKzSBwFwF2iYk5uMjFMzhNgF01Xi2saZpWEcYFKXmmOIOeALvHcA1MQcaN\nXDyDDQLsZl3TNb8XxvR3cAuDjTAhwB5gnTCed3SwUZhRHgH2oK5t2k6TPW/YHWx0jyLlEWBvkGG8\naTvMlGUU484Ay2oOjo8ZRSQPiWHms0tE//f1Cyhc93eC0J3I57ruvz9zbfvfHwjkcAW2OOHh/MEM\n50UEehtlmZqMQFRgj3OLGDlbK8vUZORPE8djXAYjlCHD1GREoQJ7XtsqwghkuyaDfKjA3qalnigW\nazIFGflQgb1kHFrqMGJRkJGn8BVY+/sZyqvv8fcymLVXCE1BRm5iV2Dt7DfA+SNZ+HtVLoZRAAUZ\nmThZgS3mxMvhMTzdUHX1j7Rtm1cdBiXaKMh8/nhH7Aqs+R0zDJFb4yKs//Pt64GL5gWZe8h4x+EK\nLERIZGS4DDb61XRyMWwtwzQrEkh/Yli8h8w5g4cEbuIIHaV9OG2XXxruCce4Im8q6nLR2pSMY5+9\n380P8TjJJvWZDCOuxVHEgk45fKycAFubZTiXN3iw3XjIMAFGdGKMh2Q6hLjd0Lg7F37W2jblszvU\nYQYSic4NZDwk0wBLNO+kz9f4pua0DIPCbDR6NJKM4zINsPRAChBdg4Pz0yvCKNL2tPeBPtB8LvB9\nYJlOurFt7bML9Vmc9t49ZKTLtAJLtxhjkcqyTYowamBokXMOV2Dn6p6Q1dJDzn4czdlB2axDxlFn\nKrCjXX8PpVcxZda28W1h6jBqoGuRRCdvkxpn0toeUrZ5Wkb3gY0d/yy6tZlquY2MNSevgXVdN+TT\noXu2OMdaYlTL9FSsuaFAybmNItMKrPn9LB6vwxRhVGt+vsnzI84Lcj2/3yRGgDXGEuEY44o0LwTY\nt/VZvgHW/M3H4XoYnCXGKvfsfWBr+ZF1rrzm7BFwPQx6Lo9V7sEg2U6pdzIsRlKeahBOvx42STtF\nG6VSkNUm/EwcRblpkt/ddTIbMUaJtmf0aIRZcQLPhViOU5+qIYEmcWV0kcotTufRM9FiYZ4dYfv8\nGliMIcTehW6OxTBbq7H04lOb7cSKcoZg7tkhxPH9zpPHH33eSky6ORK7E00QTG3G55v5CensLS18\nL06BckqkCqy5WoSda+uQYVQrcSwx0CmkNpo4cnJwxcvlfSQEkgmCodmrzDa+JNIyEapAOS5YBdZc\nnee3OVJRqcNgzYnfJGOdacqgAquXOgzWLKbRdqop1N4XrUA5KF4F1hy+b+ViIaUOg9MUat96diaO\n7Q3MxLHqbIadSyAZBnc5dxU74lkqB4/fB9Z82jQfNcB6L65BK8PgOad7s+Kevd7xxmz0Auw8GQYl\nums2kNCnt+uCn9/3lBNgjQyDKtySbf0tOaFPfimCn9/3hA+w5u15AmQYZOuJWRxDnyDjn983lRBg\nzdv9uZZfgYhuj7f8z51X7wPbbTVszHx43XyViCdHBxanWBRjkLndU0J50/C/cSNzIWXQ5yYx9uS6\ns31cTWJMhkFoiWeLQDl3Q7TMe+XHjwwl2icZVmZ2uioGcH1By8U7vcbRVWB+fG6yYN/Dvy8JLSBP\nVmQO64sMO7fc898quP/9ufWlAfUSYJG9mGHnLCaWGANuIcCCeyvDThRh21vKMOCiewJs0kw//l9X\nwh6X5bGdtH6M/yxuA3DU1QAbkqkdmXyJxw2HOo8ibLtxUYYBt7jhPrBJu3zzG11irGZrHYzW0gSu\nu+1GZkH1sX7yzqef5C94NlIncXGyyXwfAEfdE2DzCaXk2WfemoN6McMOLa2ZEocAa25o4licDjFl\njkRu9u6VsHupxoCj7pmJo2mabmTyJQpz8b7mxV0BHHVPG/3iVFJ84JUibM2h8cPBjXEIVMWNzOXK\no6X+/b0BlRBgxfmo/D1XfgGc9vhMHHwgs/ua398bUIOrbfRd1/VxpZM+L8NtYa901Su/gPfdUIEt\nBpX0yshjddgkri6mlyIMOKTEBYtHylyROd04uh47Dn3e3FV73V7MzeNQmQhl0MRRtFcWW5nMMZ+P\ntYXHLEgGZRBgpTuXYcOWr/fj3DWQeGjKfCCiwyNshzoMPx++q30IcbA2lpj+03zxMG6vxnLo2xf3\ncHH/QCZUYHWYhNbwJ93R7S+4smDYfCHNe/cP5KPwAkUF9mMjgTaO0uS73jqeR+uk3cLr4v6B3BR+\nfhdgt/kixk6XR4mB5PY1CM0QImkmifXKiOKJXDnUEunOMwit8AJFBXa/eW49fIQX0+XGmkkdBkHd\nsyIzFenjahxjD89WJVeARYYQOaXrFjobYzKQCEGVE2Bt25oF/23zGAN4SzkBxmfiX2VUhEFEhQSY\n2utjr6xA9igZBuGUEGBaDbMQP8NuN85CuQi304XIA15ZRfN2XdP1MdM27aHWx41wWswwfZVwi/AB\npvzKyLAMdOmuz5Qvw+C62AGWculrso20e9aQYYUWYbvRtZ1Mw7ffuwoo1Cl2gDUJgSSxAnh9do8U\nkwy7ZUKQISAXnwI4JNMA2y6t+kwyeJipo0XY4s+6f/CLn+84Y9bqrSup03+vng64LtMASzTPuf4R\nwRbGfKXNFyepWjOpkyZfuvcpFGFwWuAiZqNKG96UKu0zw08nZaWx+TZrS0i/6+k1w/QlwhWBK7BJ\nMqm9gtlOuHFD43f9IK/lijoMTijhRmZytH1fc0q3fQWzLAotuEKA8Zi1DEsfHqwpw7R1wFHlBFjX\ndcYPg0n5eVWQYcA55QQYOZoUYeNlw9J/2yhl4bE1ijA4R4DxsMWBxKO18mT7PsZKDDMZBukCdyES\n1bmR3vldYr0TJV1+Nu482zD/Fl0hVKXw26TcB1asxdor+M868bawlKiTZNRABUZMi80dMWcQPiSx\nUDNZMDUovEBRgdXilpk7Mpj+Y7sIS5kZxLgi9dDEQREu5s28H+Sj9pC1jsS+ZWW+2eIeJl+dfC8U\nQ4BRiu25Pzasbf91l+PipPjzfFq0GGP3vjz4XOEjbIYQq3O0KXFt2PDT4cR758Kf7M2IIsVQgVGx\ncdpNUurT33sWM+Z08EyqMaUYxSi8QFGB1SixCEupsYq4yWzw9Oow8DIVGNXbCKfT19WypA6jMAKM\n4qSkzonSSoZBZgQYbCpi8HBMhlEMAUaJEpfTTAynsgYSGxfAKIUAgxpZw4UCCDAKtbse9KGxweKK\nsEaGEZ8Ag4MKyjAIzWz0lKvr/q0EPVlO7ERrxrC3h3yxQMywDlnbtLdcGEtcEQZuoQKjGtfvSn5u\nIHFjPsYgJlMGm0GYFwgwipb/tavJlMH9pFZrszI+4PqVsI2skmE8SoBRunsnOdxt0J/82TbZYPza\nXsywu3T/4tetZrxBgFGBPgzmM/Ze2VvzmytrcbURY/PCa+2Jmmcz7EoRtnbRS4bxAgFGHR7qhlis\ntOZpNI+fxOVa8p4HZLtlQx8HTyt8snaz0fOUxKbB3cop5d/nK5PiH20gTJzbXl8iz1GBwSmTMmtt\nfHI7cor47Uoy8RUBBhfMmwbXtkl5cGMPvVeuhKVIr6vM98Fz3MgMr4hTbN11UzM8TQUGEeR0Q9vR\ny1qKMB6iAoMKzGNvqSJMmVlKUwb5EGBQtMRJqv7miuyapn0gmG6fdBEaQ4gQxolRxENb/m3ctU3X\n7nyvECIHKjAoVMq90tsJ13919L1XLmIpwridAIM40pd0Sbz3eemrbdN24yf5W4wmtxaM+euRi7Ux\nhAgBpVROZ3VN13a/V8J+d3g6J+4KmLX5763hUhsVGJTr2s1nbdcMpdgNzYdt2zTDN08HJw/sZmkK\nq8lSZEqxSqjAoCx3TJz4Lxj+9tFdKWxOTNW/tqeVCRjna7goxWogwCCU12f0GGfYmcrm6KJoGxvu\nTR88eVCGFU+AQUxrNU3vcs7dMwq3tNh0f4Ft4xrb8p7SJr+3nGZVBBiwbD6QeL49ZHeq/uQ9pySr\nDKuEAINodqur+4YZ/ybbP7XDlXLwZ29pGXaihUSG1UCAQViTM/6j8/w+cO2tnfcizt7CuB3jaI7q\nRSyeAAOOuPdO6vk27U9D/L9NTqWRifDL5j4wCGg+Jcd97Rs32Hsxw7RSo4dG76htJxMKX6mlumbY\n89L8xUSmAoPI3lwe7IE1yX5ibD3trj1H2nz8BCTAIL6syq8r/l5/f9/0+f6R3u+N0tPe/UaGhSfA\nIKZP1mhOCci0NF1dMPOWuT+aZnIL2r/o6tsqFzcjGgEGpXiz/HpmFHF6YeyWZ+m6Zt7NMY4xGRaW\nAIOwch4wPP7ahvS6YUnoo6vJyLCYBBhENpyCXwuz7Sc6kgTjqmjaMX8lWvZuoF5uqX8uw/pLccdn\nLmaXAIPgJhd1XvPM6fjftbF3yqODh26URQmvah5aYuxWAgy4yfFmyEkrxw1zZ2y+huUiLDkpJ6G1\nE2Mbe5NhNxFgwJeG6RYX0utoEXZ9ruGVPWxk1f6AZDd6i4sbcFYJAdaOfP1aoAJvjlieG0hcf4UX\nJ5dajNvV5sn56K4Mu1X4AJuElgyD9yx+3L7qjbx4N/d6Uq7NJrycYbsvQ4bdJ3aA9XHV/fn65UDF\nnjsX393NcagI295mmmGJr9DJ6iaBA2xIr+GR/u+KMHjc52XW2uNXXthmUq41mAyP/0wdkr5mm/PV\nBaXNRq8Og1e17c9E8g99AOez71/c3990+G3TbrQ+Ji5FNp1c/+hB6I/hjdaOVXGnx8AV2EATB5Rv\no2S5MT7PFkZD+XVgGpEn4mT7PrPizpDhA2y3iaP99eJLg6J9OAh2xzPu3nN2eiXoIy/i1mM42cm8\nd/+uJ8pGCUOIw7Bhn09t284vjAGPe/qztjiQeLn8mo4i/j1L1yaXU3+vod9+e1jyKe36GGb/v6P1\nQosZS8w0wLZLpbV86rpOjQXveeHq19ozzquNJ57lYPmVeGlt7bnO20ivyePD0Ssiw8IPIQJf6s+D\nn9zafH1Pe/30SWuS3RXhp2Ns477pubK6HzOtwIz7QRjvf1rfr/ySnSnCrriSQ/HrsMAV2PyuL+OH\nUIvFDoUTu9lbGHonIZZy9PBsVaeropSRw42niy9wgA0mHYaqN+CoMwtDJ1SB52ZcPOzoSa+UgcTY\nATbPKukF3OLKwtCHBw9PnLhyHUR9U+wAa0YTIZoOEbjRzsLQe/lxctr7pxeOGRRRhIUPMIArJtXS\njYN+4wy7cbc/o511/9IuwACaZi1mFiuV5OG7rWXD1p5o7xVOpq06n47xizABBjC1fBFr8R7qI7tq\n/5usMK2zcfLw5ne91C2SGQEG1G5jpcqmWV+CMvmi+2QF5+Z43vybmHH4vm5vYeikVxZ7CFKAAeyZ\nn+iPn/rneTONnJV9pi+qeV7MUUQBBvDPah4Mt05fu4E6qWwaz88w2uDfN3bLJWNtA4kCDOBtiUN/\n0/RKqJMOZ1jkVg4BBtAMtdWbzzj8/d9wYrccbLsv7IMFXDIgwAC+sTWc2M7S66FbpyMTYABf2h5O\nfKO0CtuLKMAAPvaTYd3y438PbYXN1SIs2mUwAQbwvWE4cQihn+U0o0XLOwQYQC52BgwTxvpOFmEx\nRxEFGEBmJnHyZvkVqtQTYAC5OhUn9bTUCzCAjF1YuLL4fnoBBpCfSVy9cI0q4GUwAQaQpSFRTk0c\n3P/lTBEW5zKYAAPI1bWJg0+LMvYowAAKdLGVI0SGCTCAkqVHURvtKtj/vn4BAHxsEnJdG6OnQwUG\nUKbEVo4Qo4WLBBhAvQ6tOpYbAQbAX3pFGDkcCDCAYm2PIg4Phqu9egKsFm2cmxNf4GgMHIqxqo7G\nTnpFOBQCDKBki0VY3MaNMQEGUJedxo3hMlj2RZgAAyjcuAg70HaYfUNH22X/Eq+oakQbYNX8TL93\ndsw/HQoPMAB6oW/5WmQqKYAqlBFaY66BARCSAAMgJAEGQEgCDICQBBgAIelCDG9yr9vifRHjbc5t\nEE7/jio/GtffaRmHIuUz0pR+NNp29aaplPeV58FxH1hsi3dqj3+m1zcIajHAqjoa8/dS5z+MlDdS\n/NFY+30u+sFRgQU2/LsZ/q30j8x/1ZpsMLe7QSzb76KGozE5YZ1+p9EPRfpnZL7NXMSjkfJSU95X\nngfHNbDwxp/DlF+dxo/vbhDR2lBJPUdj/ut2//f0d1rMoehtfEaaoo/G9itMeV+ZHxwBVr7tT2/K\nBiWp9mh0XbdWlzeVHYpFRR6N7s/2Not/P7TNVwfHEGJgUT5Fb9q4Ul2h0E0Ht+i6rm3b8b+Kje4e\nwlGBFaXyD2f+QzpvmhyNag/OMHzaayr+gJRHgJXDh7Op/u1PTEaQ6syw+buu8zgUyRBiCURXY/Bw\nZnJZos6z9vyjsdGFSDgCLLZ5l3DN1n7XdnAYVJvlRRJg4R06O+9+dKv6bDsaA4dirNqjkfK+sjo4\nroEFtvsPZXL3z+TxlA0C6X6NHxz+3tRxNOZvZHEupRoOxdzam6rwaKS8r8wPjgosvNOTwRzaoCT1\nHI3F3NrYYHcP4Qxt9PPH5xsXfzTWhCu8Biqwwm3Ml5O4QUnqORrX32k9hyLlkWKOxkT0g6MVB4CQ\nVGAAhCTAAAhJgAEQkgADICQBBt9YbO8G0gkwAELSRg9ASCowAEISYPAN18DgIgEGQEgCDICQBBgA\nIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQTOYLQEgqMABCEmAAhCTAAAhJgAEQkgADIKT/\nBwgTwV/btDVDAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 30, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo3()"]}, {"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"}]}}}