matlab/sparsity_8_sparsespikes_measures.ipynb

{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Sparse Spikes Deconvolution over the Space of Measures\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 explores the resolution of the sparse spikes\n", "deconvolution problem over the infinite dimensional Banach space of Radon measures.\n", "I would like to thank Jalal Fadili for his comments and corrections."], "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_8_sparsespikes_measures')"]}, {"source": ["Sparse Spikes Deconvolution\n", "---------------------------\n", "We consider $N=2f_c+1$ low-frequency noisy observations $y = \\Phi \\mu_0 + w \\in \\CC^N$ of a input unknown\n", "Radon measure $\\mu_0 \\in \\Mm(\\mathbb{T})$ where $\\mathbb{T}=\\RR/\\ZZ$\n", "$$ \\forall \\om \\in \\{-f_c,\\ldots,f_c\\}, \\quad\n", "      (\\Phi \\mu_0)(\\om) = \\int_{\\mathbb{T}} e^{-2\\imath\\pi \\om t} d\\mu_0(t). $$\n", "\n", "\n", "We follow here a recent trend initiated by several papers (see for\n", "instance\n", "[deCastroGamboa12, CandesFernandezGranda13, BrediesPikkarainen13](#biblio)) that\n", "performs the recovery of sparse measures (i.e. sum of Diracs) using a convex sparse regularization over the\n", "the space of Radon measures (which is in some sense the dual of the Banach space of continuous functions).\n", "\n", "\n", "We consider, for $\\la>0$, the following regularization\n", "$$ \\umin{\\mu \\in \\Mm(\\mathbb{T})} \\frac{1}{2} \\norm{y-\\Phi \\mu}^2 + \\la \\norm{\\mu}_{\\text{TV}},\n", "      \\quad (\\Pp_\\la(y)) $$\n", "where $\\norm{\\mu}_{\\text{TV}}$ is the so-called total variation of the\n", "measure $\\mu$, which should not be\n", "confused with the total variation seminorm of a function. For an arbitrary measure, it reads\n", "$$ \\norm{\\mu}_{\\text{TV}} = \\inf\\enscond{ \\int_{\\mathbb{T}} f d\\mu }{ f \\in C^0(\\mathbb{T}), \\normi{f} \\leq 1 }. $$\n", "When there is no noise, it makes sense to consider the limit $\\la\n", "\\rightarrow 0^+$ of $\\Pp_\\la(y)$, which reads\n", "$$ \\umin{\\mu} \\enscond{ \\norm{\\mu}_{\\text{TV}} }{ y=\\Phi \\mu }. $$\n", "\n", "\n", "We focus our attention on discrete measures, which are finite sum of Diracs, and that we write, for $x \\in\n", "\\mathbb{T}^n$ and $a \\in \\RR^n$\n", "$$ \\mu_{x,a} = \\sum_{i=1}^n a_i \\de_{x_i}. $$\n", "In this case, one has\n", "$$ \\norm{\\mu_{x,a}}_{\\text{TV}} = \\norm{a}_1 = \\sum_{i=1}^n \\abs{a_i},  $$\n", "which shows that $\\norm{\\cdot}_{\\text{TV}}$ is a natural generalization\n", "of the $\\ell^1$ norm to measures.\n", "\n", "\n", "Several recent papers (such as [deCastroGamboa12, CandesFernandezGranda13, BrediesPikkarainen13, DuvalPeyre13](#biblio))\n", "have studied conditions under which $\\mu_0$ is the\n", "unique solution of $\\Pp_0(y)$ when $\\la=0$ and the robustness to\n", "noise when using $\\Pp_\\la(y)$ with $\\la \\sim \\norm{w}$.\n", "We refer the interested reader to these papers and references\n", "therein.\n", "\n", "\n", "The goal of this numerical tour is to detail how to numerically solve\n", "this problem using an approach proposed by [CandesFernandezGranda13](#biblio)\n", "and also to get some intuition about the primal-dual relationships\n", "that underly several theoretical analyses of the performance of the\n", "recovery.\n", "\n", "\n", "In the following, we denote, for $p\\in \\CC^N$, the adjoint operator\n", "$$ \\Phi^* p(t) = \\sum_{\\om=-f_c}^{f_c} p(\\om) e^{2\\imath \\om t},  $$\n", "so that $\\Phi^* : \\CC^N \\rightarrow L^2(\\mathbb{T}).$.\n", "\n", "\n", "We denote $ \\Phi_x(a) = \\Phi(\\mu_{a,x}) $ the action of $\\Phi$\n", "on a discrete measure, i.e.\n", "$$ \\Phi_x(a)(\\om) = \\sum_{k=1}^n a_k e^{-2\\imath\\pi x_k \\om }.  $$\n", "\n", "\n", "We also denote\n", "$ \\Phi'_x(a) = \\Phi(\\mu_{a,x})' $, which defines two linear operators\n", "$ \\Phi_x, \\Phi'_x : \\RR^n \\rightarrow \\CC^N $, which are represented and stored as matrices\n", "in $ \\RR^{N \\times n} $.\n", "\n", "\n", "\n", "Some display options."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ms = 20;\n", "lw = 1;"]}, {"source": ["Sampling grid for the display of functions."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["P = 2048*8;\n", "options.P = P;\n", "u = (0:P-1)'/P;"]}, {"source": ["Set the cutoff pulsation $f_c$ and number of measurements $N=2f_c+1$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["fc = 6;\n", "N = 2*fc+1;"]}, {"source": ["Fourier transform operator."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Fourier = @(fc,x)exp(-2i*pi*(-fc:fc)'*x(:)');"]}, {"source": ["Operators $\\Phi$ and $\\Phi^*$. Note that we assume here implicitely that we use real measures."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Phi  = @(fc,x,a)Fourier(fc,x)*a;\n", "PhiS = @(fc,u,p)real( Fourier(fc,u)'* p );"]}, {"source": ["Set the spacing $\\de$ between the Diracs."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["delta = .7/fc;"]}, {"source": ["Position $x_0$ and amplitude $a_0$ of the input measure $\\mu_0=\\mu_{x_0,a_0}$ to recover."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x0 = [.5-delta .5 .5+delta]';\n", "a0 = [1 1 -1]';\n", "n = length(x0);"]}, {"source": ["Measurements $y_0 = \\Phi \\mu_0 $ (noiseless)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["y0 = Phi(fc,x0,a0);"]}, {"source": ["Add some noise to obtain $y = y_0 + w$.\n", "We make sure that the noise has hermitian symmetry to corresponds to the\n", "Fourier coefficients of a real measure."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["sigma = .12 * norm(y0);\n", "w = fftshift( fft(randn(N,1)) ); \n", "w = w/norm(w) * sigma;\n", "y = y0 + w;"]}, {"source": ["Display the observed data on the continuous grid $\\Phi^* y_0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAbKUlEQVR4nO3d63aj\nOhIGUJh13rvdT878IKEJdmyuUpXYe/Wa6UtOItuCTyWE6Idh6AAgm//VbgAA7CHAAEhJgAGQkgAD\nICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTA\nAEhJgAGQkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlAQYACn9\nV7sB1+r7vnYT4L1h9nvdlUCGYfj8RVU1HmBdjM+g73vNiNOGSM3ouj8/Qmt4lG5VmLdCMwK1oUsy\n+g/xTl0nTlfQjDhtCNSMR9/NQmvxx0JtCPJWaEakNsRpxnuugUEd4wB3Hlfj7zMMfCGE9qcQIag/\nL5JqeAx91/+8MAabpaifjhNgJQTpSRGaEaENXYBm9H3X/alwxetZ9bdilL0Z514xOuW77fgmQT6F\n9QQY1PCq/Jr0jz5CtrFJurP/QopVGwuugUFpz1e/5kQXrCTAoLi35ddoXJEIa/TfajekNAEGFbwv\nsxRhrDfPrbtlmACDom52huFaY2KNl9/G/32ZYa2GnACDslbMH8J688Ujvy0kGYZhzK3GltcLMCht\n5Qyhy2DwngCDctZP3rgMxhpjaTVVV+8LrMbKr06AQVHmDznbminE539qY9WiG5mhKKUVp9taV02l\nWPaaTAUGhWwe7/4dXAbjLNmz6iUBBqVsnD9s7mxzL31/4a8d2kuvzhQilGT+8D5azItwVGBQQv7r\n5bRjWruYvSwTYBCay2BcYRiG7OnVCTAoZNcCelOO8IYAg0KkEZxLgMHlXACDK1iFSAumhAg6q28D\nDriACozcpttixugKW+scmT+0jgNeEmCkNwxf6TX+JmyG7fQ3ZlEJ9QkwEvstq0Jl2MHGBJ0UhQAE\nGFnNZw7nws0lugBGcQ3sNL+GACOx36qTaFWLBfRwBQFGSmvGl/cYg8ILDeyysYZl9GT1NVX4tEJv\nLHeCrObo+677c/i7/B36rlfGwYIAI59FMs3P7P2j7x//zvV9X3s68YwLYMPQ9Y/j34ZmTVe8blJ4\nTUwhktJUfi3qkvGPY1kW5FhWOXGp+XqNm6zdmAgwsvrt9t5FYNzsiOZexsQaC6/xf19mWKshZwqR\nZOZH32/FzfAYxuKs7pWwcy6AkdOl+6fMe/582nB80NeLrx+G8elfDTwDbC59gE2fx21ngW9oGFad\nHaYu8X3HWKCOsTjLhGobpzB1XEDiKcTxiaLT74dvLRXIvPHhBPG3VDve+PN66eC+LmpHRJ5Nz1bu\nvk+JbwZDjZVfXeoAa+OJomwynvn3ncojj2w+t82OiPxiMYW45suakTjAVup/qt0cjvratzf8/My5\nfa3Fkw+nmeafDn6fdKfK9NfAPmpy3HFnu2fS9ATubM384fwLUmRY+xUY7YlffnXdtluYhStXa7KP\nNVKBzdduNPk50W2clHu6J2YovyvHr6v8rZuFM6SvwKaD/6xZYELbuDPTMAzdn677Uzq61mTtjh5r\nISLMpQ8w7mbr/GGO+cY1LESMrU+u9vu3RyNTiLBGuVlED7G8GXM/VajAyGH3AHEcWvaPPsgFsLlN\nw15nSFgQYOTxy8YWH7UziwjMCDDuYlwBUWCqP+fVBMhHgHEL02OaS3ABDIoQYCTQ98lSwaQlFCDA\nSON4KhSbRbyOW8FgIsBo3L9b3YvMIu7YLmQDt4LBjACDNKykhzkBRgbnXQC7fApu71p/YCsBRnTj\npNwpqTB9k9SXwYCRAON2rpuIk4tQkgDjdizkgzYIMMI7dgFssd/gtbOIG5uadAtwCEKAEdqJF8Dm\nrptFvHwFx99BBQkjAQbnUE1BYQKM2xkeQ94tOdwKBhMBxk2dnwSpdmuEBggwYkuSChddqwPeEGDE\ndUoq/LbfYIRZRM+hhyMEGHdU9PFgwDUEGJzBFohQnADjvs6aRSw9D+lWMOi6ToAR2pVljVlEyE6A\nwWFll0oKXRgJMII6a17u/X6Dx2cRjyyVtBciHCHAuC+ziJCaAAMgJQFGVKUuLB2ZRez7NHuFQHsE\nGHEVuLPqlFlEd4BBFQIMUnIrGAgwIjpxdd7H/QZ3P13l+Pzh/r0Q/6r5QIBB13UHZhHNH0ItAoyQ\niq+MiLA5/XqW/kMnwAirZGWzYymH9YdQnQCDHzYVYeYPoSIBRji15vH6R7+yCMsy0whtE2A0buV+\ng/NaatV/ccb8ob0Q4QgBRjz1ri2tKcKm0DF/CHUJMCKqkg0birAYyzfcy8zNCTD4YSrCXmZYoPLL\nvczcngCDf6ZYeplhX3+MUX4BAoxYIixrGKfmFhn29eDKoesilF+AAKN5W/cbHMNpkWF93w1DNwwn\nX3bavxeizThAgBFOgAm6eYE15tZXko2ppvyCGAQY4QRJiEWxJb0gmv9qNwAiGp+xssgw6QWhCDB4\nTVxBcKYQAUhJgBHIFWvoI+83eLRtfwebcXBnAoxIAixBBLJo6hrYfDx75A4bKnLlCViptQps+Fa7\nIXA53Zybay3A+r6PfM0DgLM0NYXYfc8c9n0/FWGLPFOcAbyUbvTfVAX2MpyGn8q3ipUuOnYif+iR\n28YNpTtVtlOBzasuUrIEkW/z2wOs6+E37QTYMAxT/SvJknKqonvadrJ/9DoGL7U2hZio+IVTNHYv\n8/OmyeO+lPVaRFxNBRjczt8Gh2sv6y0ZxrN2phCB7JaPsJldCHOJlGcqMEK4bvlu5JXBkdtWy79L\nX33Xdf+eJtopwngiwIjB+Pr25vk0pddoGNqcLOUgAUYUVprt0NiKpR9rN36+tMZeKacQYEAsv06s\nenwMPwkwoL5FMr2stxRhLFiFCIQwzh9+XNeS9L7mZUInfAkBqcBoXOQb2yO3raI370rS8/50d/b4\nq7Oi8iQCjPoczGyS6+6Dl3uLdLr9GQQYISQdWYeQf2nD1P61yZTnpovn9Brp8KcQYEB90wn946xq\nulP/bw22x+NxAgzIJ8Us4pr1JjLsCAEGhJAik9Zbk0zpqsloBBiNi7zf4Clty76ScV6mbHgtGS6D\nrcwnRdhuAozKHL1stfKOsRQUYUe4kbnrPL+8Nu85jdl6t3XSu7OrU4H9uMewUxBADXvKqQyziGuI\nrt3uHmCLuzTGGOsfvRiDAn5Mfmw5jUeeRdx39nDO2eHWAdbqPYZ9H/TAhmfZD7eXtr6oJt+EAu4b\nYL+l1+ILcplHlxgbRd5v8MS2ZeyucNBNA+xjemW8HjZ/BPv0IPbgGZbrHY4r+dOK9/fSeJfBjizH\ncDhsddMA61bU7LmK+sUj2EcpMiz7yZdT7KhFI18G2yHXCSeIOwbYpmFOijHRy/Qaxc+wwDN8XCvF\nwbVJe68ouDsGWLd6sJNoTPTuEUppXgS3k+gQW+m3VzRek55+/UYEbnK7ANvRP4J3qb5fFVGRizAO\nuvUY5U/kzcK+PF+ffhlj7cX51W4XYN3GXhK8S608dO98got8eovctmKOvAehDs/fRrqJr0+Hd6+t\npBY97MUI6OVlpMcQeaOX5zYvXubU8pW1WjHBS1uKCdUtj3g+S7y/Pj3WYfN/DX62ieZeAdb9PJt3\nPzvWVNRnOZxeDt8WdwhMu4oMwxBxuPd36B6120AN9xm+vL8+HfGozOOOU4jd70X9m7o+5vG2eAnP\n97ct9niMdrRkGStwhRPrjGgde+T69NVuFGBTYf6+zPptFvG6hu3z3Onf3J09/qW0aFvMMdbVhscQ\n4Xbm5zd/0/XpxRff86Pc4UYBNvf+VP5bXR+tV/2Y//y0t0gXtQjjHDnvB2+pN77YVXXdZ7KcCoo3\nXA7rLgE2Zc/6hQyLQyt+r3rfwq9/DTBWLewmeyHm1eR7sCOYW8ryYu4SYN24vGd1Fwm+yHURwytL\nw1AZHK2cJbuKR+vLzrzt6TCf1hLz0o0CbLS+V738yrC9akM4xbnxM+esF8edfxzVnlqYH4C7j6/p\nPww11ozsFgF25NrP/L8K0quOXO8N8hJGTc4dsdKJXTFUrx7t2Zs43ItI4BYB1u1dhhe2Sx296lt7\nuMrpwvZVNpkPT8PO98RxlwA7MmkWeYXrjsY09hAK8jq9E1bp1c/7++weTCy25DjQqLtoP8Cm27/2\n9apoK1xP2wC0ahFWchAQ5orfC5HbVsaZhWO9Ln3uaeH2nWKD9gOsO6NDhOpS0zG/Owaqx3DXWcHB\nmUJ06TNOFKaCN2k9wP58/f+RbhFnheu5z18IlcrQhnMTKNQFi4BaD7Du/JF+3bHe8fLrn6qziEaa\nV0hxvrt0t/XCw7KLXsvXdnd1y8o/n7+kutYD7G/XnXGujLlp9JH+HWTKhTPdfla27qjorFNEhLFd\nipFQ136AdZcN9Ip/wFckaJZuSnua7M8nZs/0/lR7URnGQ80HWH9Fl6pVvpw5f1jvVRQ+ICPvNxi5\nbQWc++rL9+d5Tz43j6d3xkzJe80H2GminWpy9+wMgzuSKjnb/+PZe6106uqF7HoCbJta98lf93MT\ndVb4qJkU6equE/47dF2CM4MA26DuffJX/PTcZRw/pTh3Xz1gqjIguyJpdjxu4oYEGBWkONVykeuG\nTSUHZItQua5LFx5lfj0aN8kRKsA2Kz+LePVPNL6DHQpES51ZxDyXqAXYNrVmES/8uWU7a43bD+LG\nc+S2XerIjrcrv38x1/0ss4gfCbC7qzBXkGd8R0blu/TVYVx4FjHL/GEnwHZY7Mpx+UXpIj/L+K4Z\nPsoC70CxN7lwlqTrPE0FWP+t2E8sMzK6fN5SSdSM23+UxYqVVp+rl2tlcjsB1vf98K1AhpXvuOkG\nR79JNEHBuWrdOnmpAv15fC3NnAFO1Dezmc0YYIvf3/YiOcBB8dPhv9oNuNxFn8F8GdV1j4co81PG\nH9T9ufD7f/2UR9/9jX9QJHbi5zgfEZ5l6sOXrkL8qlQu62lfd0pd/0L+/cS+G4ZrzwDd0xkmxei/\nnSnE8oou1TV7QCuuPukXuIpTMr1GV88iZkirFwTYTouOe8n9xT+/5aWH5fhyxGR2qtu5pCflZ+PH\nem0wV3287W7tBNi4dmNUeJbqwq1xSr6OIgvYnGEpQ0/bKtf6w1E7AdZ13bQKsdSPKzTEa6MwauNV\nkMVF/W36toXLu+lsU2CyJ5GmAqyuk59yUnD+sJzb36VEGdfOt39/8wpbflz0unLOH3YC7Lgxaa7o\nWIUPj2Hour+DOonRdUsQS8pbW7x04b6LOYfIAuyQEitoJQob6TOjC1fq13iDr7orIHNnEWBnOuvE\n0eb8IWXEnqdNfbqcF5EVF4mcPEBJO3/YCbBTXDGL2M7hMfueFobRFewGbVaipw5QrrsCUoYAO+ra\nJymUPQLHy2BXfffYlQGNmU7Kp9d8FYtIQ8AFAXaOqU8fjxzzh3CWE8/480O7epCcMrQdNx47/n0q\nEmAnmHrzWWFTfV+MNude2CLFVniFRRhNnj5NEuFF7SbAYqlefl0xiygOS6peHLxUpQ8Mj6/bQk7M\nYrEeigA7zVmziDFPQEfZhP72ao30T+94QXrywfNMA/OHnQA7y+mziNWrluoNoD15y5dQh8OJ0ySp\n5w87AXaug0VY9fnDf6wY5Brly5ezsmd4DNECePdLa6P86gTYiU4pwqov3+jCzJDAcacvpg9ydJxS\nhGUvvzoBdrrdx0mg8uvbOUt1H716rrzjn11LFy0PvpRQ84fHRasjjxBgZ5ofJzs6fagzRjMbi9yR\nEcPZAs4fHpqtaWL+sBNgV+j7zWf/+bERZLgX7XCF3c5aTB9uHLZ9pDK9CRHmeI4TYCfbXYTN/8Pq\nfWtqTJA0JbXK13QfQ3cgeyIfAjuLsFbKr06AXWRTERaw/Bodz1EXwBhVH5MdEXD+8J8tx1dj5Vcn\nwK6wowgLVX6dLtzECzXUzYDjs4gxu/H6Vn299obKr06AXWdlERa3/Bq+2haqVRRzxV6ItTLg6Cxi\n+CNg00C5pSGyALvEpiIscvl1pD2Sr5aYtUJG0xMsI7+lYwvfp2yrg1EBdqGPRdii/IqWXse194rI\n6/SNfWP582vN/PXUyubKr06AXedrgdDvA595rwo7Mjo4cGv2ZEE2OwqpRIH3FUuvMmw6z4Q9yRwh\nwC70dbT8HbqnDHhOr4AjoyNDtun1Rp574Z7WZ1L8+cPJfNOs+a9h+JdeAU8yBwmwa73MsHl6fX1Z\ncx2ra/RFJRJkxB2kGaP+0W8rwgK1fZ0//ZhY069Jk8ejALvcIsMWtVf8S1+tXv5t3PGdXs+rO4L0\n8E3NyDiF8LWaYzHZk+Eks9t/tRtwC9/HwNA/+vE+jP7x/U+xO9a4mH56su1KiS4ecENfvfrjkZfw\nHvzxUF0crcFPMkcIsKLmQ6R0vWrTOG54DP3fTKNXLrUqMIqYuvGbJmUsvybpTixHmEKsYHgMuTrZ\nju2JIaCpG68pvzKm190IMD7YujOW+UOCWyyn+u1fiU+Asdam+cMu5/QLzVsUYa9HWsqvJAQYq0zH\n+fvxqdFrEMfPv1fshRjH13Xopwyb/z3xCTA+m47nNUVY6GdPwGIt1XeG9X3W1VV3JsDY7Lcya/73\nxrB0gXf4XGTYMHw9ZyRma/mNAGOtr1uwX90s2aW9N6BtZnTfmHrydOOU3ptO3/bFyr5v/AWWNL9v\nZnHAz//4vFEWVfR91/3ZXwCdcuyMFZgukVGKk6cbmdlgyrDnG/7nJ8rw3Z7SdAmuIMBYa3pG89cf\nzbc0Lf7oG1wDY5s1D34FKECAscGaQbmBO1CGAGOz908uByhDgLHNrxt4W2nGTxbxczUBxmaL1Rzz\nvyeOCB+HlT5cSoCx07j7zvSbCKdLTtT2Xoi0QYCxx9fuO2YO+cTghusIMPYbY8zpKTIXomiYAIN2\n/TW4oGUCDICUBBgAKQkw4IWDeyG69kYBAgy4hJvAuJoAAyClph6nMr/10sMgoDo3Q3OppgKsk1sQ\njCOS67Q2hdj3vS1wYCQ8aFubFVjf91MptsgzJRqsMT+IuIl0o/+sFVj/0/iXL4+34aeyzYT6rGhn\npXSnyqwB9hxL6cYOUILdpGhXO1OIwzC8L8WAMtR8lNFOgHVyC8JwFzMFZJ1CBODmBBjwwvH5DFel\nuZoAA65iUp9LCTAAUhJgAKQkwKBlJvFomACD9rkxiyYJMGjdrs04dm9t0z96N4FRhgADICUBBkBK\nAgw4n7uYKUCAAZewAJKrCTAAUhJgwAue7UB8AgyAlAQYNK5kKeWOaUoSYHALxaLFXcwUI8DgBnZt\nxgHBCTAAUhJgwAu790Ls3MVMKQIMOJ9F+BQgwABISYABkJIAAyAlAQacw6MsKUyAAS/YC5H4BBi0\nbwwj+zzRGAEG92AzDpojwIAz9b2bwChEgAGQkgADICUBBrxwZC9EKEOAASdwExjlCTAAUhJgAKQk\nwOBG3MtMSwQY3MIwlLiX2coPShJgwAu790J0FzPFCDAAUhJgAKQkwICjrA2hCgEGnMBdzJQnwABI\nSYABL9gLkfgEGNyL61U0Q4DBnVx5L7OajcIEGNxFgVuM3cVMSQIMgJQEGHCIi2rUIsCAFzbthegm\nMKoQYACkJMAASOm/2g04qu/7ca5juu9y92MggO7nLcyOJiJLXIH1fT8daWOMjewgAO+9WXaxOHzW\nH00OO8pLHGBjXNVuBSSzdcHF+gxzOFJY4gBbqf+pdnOgsnMPAmvoW5LuVJnjGtji3dy2wNewEL4N\nw/lzfdbQN2N+tkyRYTkCTAhBLY4+wsoRYB/N12443mC3xTIoRxORpQ+w6QBzpMEpHEpk0f4iDgCa\nJMDgdobHYPUgDRBgwE5TCva9m8CoQIAB+1lDT0UCDG4nwx0+8JkAg3sx10czBBgAKQkwAFISYHBT\nVtKTnQCDOzq+enC+hh6qEGDATlMKWhhCFQIMgJQEGAApCTC4nSseawnlCTAAUhJgcF9W0pOaAIOb\nOrKSvn/0tvGlOgEG7OdaGhUJMOAQN4FRiwADICUBBkBKAgxuyuUrshNgcEfjhavhMexYSW8JIkEI\nMABSEmDATiYhqUuAAftZQ09FAgyAlAQY3N2mdRy2TyQOAQb31fd7dkS0BJEgBBjclMtXZCfAAEhJ\ngAF7WENPdQIMWLs0Y/FlJiGpS4DB3W1alGEFB3EIMABSEmBway5lkZcAg/vadBHLLcxEI8CArluX\nT9MFMHUbEQgwYNd+HBZzUJsAAyAlAQZ85gIYAQkwuLvxgtbwGN6nlDvAiEaAwa25lEVeAgz4YFGZ\nWYJIEAIM+Oe3WcTF/KG6jQgEGPDFVS5yEWDAD8sJw0cv2IhJgAH/LmstssrqeSITYHB3zxe05rm1\njDSJRhj90PTV2L5v/AXCKfr+R4xNAfY8ebj4SlqV4uSZoIlHpPgMoLr1sSTAbiLFydMUIgApCTCg\n61zcIqH/ajfgqKnO7WfHX/zKF0IZhlUBJuQIJXGA9U8Hk9yCqznIiCNxgD0XXuPvFzG2yDkhB/DS\nc1UQXOIAezZF2jylJBasZIXhzc3PlinCLEeAramiBBUc8fEyWIYTGveSI8A+hlOKWxYgOwcZoeQI\nsI+GYZiqNEkGcAfpA2yKK7kFx/12Gcz8IQG5kRn48n4QaIhINAIMgJQEWAlBFqRGaEaENnSa8bYN\nz40q0MwIb0UXoxkR2pCFAAP++W2e0PwhAQkwYGleA6gHCEuAAT88F1vKL2Jq/P5fs8mw13RmcBDd\nVPx0aDzAAGiVKUQAUhJgAKQkwABISYABkFL6zXzfqLg//fsfXebhL2/aUPKd+e1nFf50snwiL/+1\nTDMifCKFn5/+sVcUaMObZnjCxkfNBtj8lFT4aWFvfnSxZf0fX/7Lp1cXbkaZNnxsRpkPZeUnUqsZ\n85N1xU+k5MHysVfUPUYqnsESMYVY1DAMETqiNkzinBr6vq9+2+LYhghvSJBmEFyzFRgfVT9HVD9f\nh1Ks+nnfgLptCKXYFOKaZvCSALujIEfmdMqu1YDxR0//W/ENqf5ZhFI9QSNM301PmZ8/bp4FU4g3\nVf0EUfGnT4ZvXdU3JMi7QShjcBrZvNfyXEGcVYjPSzmqLG0af27dhV7za+NlGvC+GdO/VvxEXv5T\nsWZE+0RKVjw+kexaDjAAGmYKEYCUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgA\nKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEG\nQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKf0fS4XOQpwt6IcAAAAA\nSUVORK5CYII=\n", "output_type": "display_data"}], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["f0 = PhiS(fc,u,y0);\n", "f = PhiS(fc,u,y);\n", "clf; hold on;\n", "plot(u, [f0 f]);\n", "stem(x0, 10*sign(a0), 'k.--', 'MarkerSize', ms, 'LineWidth', 1);\n", "axis tight; box on;\n", "legend('\\Phi^* y_0', '\\Phi^* y');"]}, {"source": ["Dual Problem\n", "------------\n", "The Fenchel-Rockafellar dual problem associated to $\\Pp_\\la(y)$ reads\n", "$$\n", "      \\umax{p \\in \\CC^N}$\n", "      \\enscond{ \\dotp{y}{p} - \\frac{\\la}{2} \\norm{p}^2 }{\n", "     \\normi{\\Phi^*p} \\leq 1 }$\n", " $$\n", "which has a unique solution $p_\\la$, which is a projection on a closed\n", "convex set of $y/\\la$, since\n", "$$\n", "      p_\\la = \\uargmin{\\normi{\\Phi^*p} \\leq 1} \\norm{y/\\la-p}^2. \\quad\n", "      (\\Dd_\\la(y))\n", " $$\n", "\n", "\n", "We denote $\\eta_\\la = \\Phi^*\n", "p_\\la \\in L^2(\\mathbb{T})$, which is a trigonometric polynomial.\n", "\n", "\n", "The Lagrangian dual problem associated to $\\Pp_0(y)$ reads\n", "$$\n", "      \\umax{p \\in \\CC^N} \\enscond{ \\dotp{y}{p} }{\n", "     \\normi{\\Phi^*p} \\leq 1 }.\n", " $$\n", "It does not have in general an unique solution. In [DuvalPeyre13](#biblio), is it proved that\n", "$ p_\\la \\rightarrow p_0$ when $\\la \\rightarrow 0$ where $p_0$ is\n", "the solution of $\\Dd_\\la(y)$ having a minimal $\\ell^2$ norm. We\n", "denote $\\eta_0=\\Phi^* p_0$ the _minimal norm dual certificate._\n", "\n", "\n", "Strong duality holds between the primal problem $\\Pp_\\lambda(y)$\n", "and its dual $\\Dd_\\lambda(y)$, and thus the\n", "values of these problems are equal.\n", "For any solution\n", "$\\mu_\\la$ of $\\Pp_\\la(y)$, one has\n", "$$ \\eta_\\la = \\frac{1}{\\la} \\Phi^*( y-\\Phi \\mu_\\la ). $$\n", "\n", "\n", "Furthermore, the\n", "support of any solution $\\mu^\\star = \\mu_{x^\\star,a^\\star}$ of $ \\Pp_\\la(y) $\n", "satisfies\n", "$$ \\forall i, \\quad \\abs{\\eta_\\la(x_i^\\star)}=1.  $$\n", "This property is at the heart of the numerical scheme proposed by\n", "[CandesFernandezGranda13](#biblio) to solve the primal problem.\n", "\n", "\n", "In [DuvalPeyre13](#biblio), it is shown that this certificate $\\eta_0$ is\n", "important, since it some how governs the stability of the recovered\n", "spikes locations (in particular how close they are to $x_0$) when the noise $w$ is small.\n", "\n", "\n", "Douglas-Rachford Algorithm\n", "--------------------------\n", "Note that this section is independent from the other one, and in\n", "particular this section has its own notations.\n", "The Douglas-Rachford (DR) algorithm is an iterative scheme to minimize\n", "functionals of the form\n", "$$ \\umin{x} F(x) + G(x) $$\n", "over a Hibert space (assumed for ismplicity to be finite dimensional)\n", "endowed with a norm $\\norm{\\cdot}$, where $F$ and $G$\n", "are proper closed convex functions with intersecting domains. We assume that the set of minimizers is non-empty.\n", "We also assume that one is able to\n", "compute the proximal mappings $ \\text{prox}_{\\gamma F} $ and\n", "$ \\text{prox}_{\\gamma G} $ which are defined as\n", "$$ \\text{prox}_{\\gamma F}(x) = \\uargmin{y} \\frac{1}{2}\\norm{x-y}^2 + \\ga F(y) $$\n", "(the same definition applies also for $G$).\n", "\n", "\n", "The important point is that $F$ and $G$ do not need to be smooth.\n", "One onely needs them to be \"simple\" in the sense that one can compute in closed form their respective proximal mappings.\n", "\n", "\n", "This algorithm was introduced in [LionsMercier79](#biblio).\n", "as a generalization of an algorithm introduced by Douglas and Rachford in\n", "the case of quadratic minimization (which corresponds to the solution of\n", "a linear system).\n", "\n", "\n", "To learn more about this algorithm, you can read [CombettesPesquet10](#biblio).\n", "\n", "\n", "A Douglas-Rachford (DR) iteration reads\n", "$$ \\tilde x_{k+1} = \\pa{1-\\frac{\\mu}{2}} \\tilde x_k +\n", "  \\frac{\\mu}{2} \\text{rprox}_{\\gamma G}( \\text{rprox}_{\\gamma F}(\\tilde x_k)  )\n", "  \\qandq x_{k+1} = \\text{prox}_{\\gamma F}(\\tilde x_{k+1},) $$\n", "\n", "\n", "\n", "We have used the following shorthand notation:\n", "$$   \\text{rprox}_{\\gamma F}(x) = 2\\text{prox}_{\\gamma F}(x)-x $$\n", "\n", "\n", "It is of course possible to inter-change the roles of $F$ and $G$,\n", "which defines another iterative scheme.\n", "\n", "\n", "One can show that for any value of $\\gamma>0$, any $ 0 < \\mu < 2 $,\n", "and any $\\tilde x_0$, $x_k \\rightarrow x^\\star$\n", "which is a minimizer of the minimization of $F+G$.\n", "\n", "Solving the Dual Problem with Douglas-Rachford\n", "----------------------------------------------\n", "It is in general impossible to solve numerically $\\Pp_\\la(y)$ because\n", "it is an infinite dimensional problem. In contrast, $\\Dd_\\la(y)$ is a\n", "finite dimensional problem, so that there is some hope to be able to\n", "solve it with an algorithm.\n", "\n", "\n", "Recall that the dual problem reads\n", "$$ \\umin{p \\in \\CC^N} \\norm{p - y/\\lambda}^2   \\quad\\text{s.t.}\\quad  \\normi{\\Phi^* p} \\leq 1. $$\n", "\n", "\n", "As detailed in [CandesFernandezGranda13](#biblio), the constraint\n", "$\\normi{\\Phi^* p} \\leq 1$ can be re-cast as imposing that the\n", "trigonometric polynomials $1-\\Phi^* p$ and $1+\\Phi^* p$ are sums of\n", "square polynomials. A classical result (see for instance [Dumitrescu07](#biblio))\n", "ensures that\n", "the convex set of sum of square polynomials (SOS) can\n", "be described as the intersection between the set of positive\n", "hermitian semi-definite (SDP) matrices $\\Ss^+$ of size $ (N+1)\\times(N+1) $ and an\n", "affine constraint.\n", "\n", "\n", "The problem is thus encoded using a SDP matrix $X \\in\n", "\\CC^{(N+1)\\times(N+1)}$.\n", "We define two linear mappings\n", "$X \\in \\CC^{(N+1)\\times(N+1)} \\mapsto Q(X) \\in \\CC^{N\\times N}$\n", "and\n", "$X \\in \\CC^{(N+1)\\times(N+1)} \\mapsto p(X) \\in \\CC^{N}$\n", "extracting sub-matrices as follow\n", "  $$\n", "      X =\n", "          \\begin{pmatrix}$\n", "              Q(X) & p(X) \\\\\n", "              p(X)^* & 1\n", "          \\end{pmatrix}.\n", "  $$\n", "The additional affine constraint is denoted as\n", "$$ {\\Cc} =\n", "      \\enscond{ Q \\in \\CC^{N \\times N} }{\n", "      \\forall c \\neq 0, \\sum_{i} Q_{i,i+c} = 0,\n", "      \\text{trace}(Q)=1\n", "      }. $$\n", "\n", "\n", "We thus re-write the initial dual problem as\n", "$$\n", "      \\umin{X \\in \\CC^{(N+1) \\times (N+1)}} F(X) + G(X)\n", "$$\n", "where\n", "$$\n", "  \\choice{\n", "      F(X) = f(p(X)) + \\iota_{\\Cc}(Q(X)), \\\\\n", "      G(X) = \\iota_{\\Ss^+}(X), \\\\\n", "      f(p)=\\frac{1}{2}\\norm{y/\\la-p}^2.\n", "  }$\n", "$$\n", "\n", "\n", "Regularization parameter $\\la$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 13, "cell_type": "code", "language": "python", "metadata": {}, "input": ["lambda = 1;"]}, {"source": ["Helper functions."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 14, "cell_type": "code", "language": "python", "metadata": {}, "input": ["dotp = @(x,y)real(x'*y);\n", "Xmat = @(p,Q)[Q, p; p', 1];\n", "Qmat = @(X)X(1:end-1,1:end-1);\n", "pVec = @(X)X(1:end-1,end);"]}, {"source": ["Fonction $f$ which is minimized in the original dual problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 15, "cell_type": "code", "language": "python", "metadata": {}, "input": ["f = @(p)1/2*norm( y/lambda-p )^2;"]}, {"source": ["Its proximal\n", "operator is\n", "$$ \\text{prox}_{\\ga f}(p) = \\frac{p+\\ga y/\\la }{1+\\ga}. $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 16, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Proxf = @(p,gamma)( p + gamma*y/lambda )/(1+gamma);"]}, {"source": ["The proximal operator of $G = \\iota_{\\Ss^+}$ is the orthogonal projection\n", "on $\\Ss^+$. It is  computed, for $X=U\\diag(\\si_i)_i V^*$ any SVD of\n", "$X$ as\n", "$$ \\text{prox}_{\\ga G}(X) = \\text{Proj}_{\\Ss^+}(X) = U\\diag(\\max(\\si_i,0))_i V^*. $$\n", "For convenience, this is implemented in the function\n", "|perform_sdp_projection|."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 17, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ProxG = @(X,gamma)perform_sdp_projection(X);"]}, {"source": ["The proximal operator of $F$ is  the concatenation of the\n", "proximal operator of $f$ and the proximal operator of $\\iota_{\\Cc}$\n", "(which in turn is the orthogonal projector on $\\Cc$). Note that $\\Cc$\n", "is an affine set for which the projection is easy to compute as\n", "$$\n", "      \\text{Proj}_{\\Cc}(Q)_{i,j} = Q_{i,j} - \\rho_{j-i}  + \\de_{i-j}$\n", "      \\qwhereq\n", "      \\rho_{c} = \\frac{1}{N-\\abs{c}} \\sum_{i} Q_{i,i+c}$\n", "$$\n", "where $\\de_{c}$ is the Kronecker delta. For convenience, this is\n", "implemented in the function |perform_sos_projection|."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 18, "cell_type": "code", "language": "python", "metadata": {}, "input": ["ProxF = @(X,gamma)Xmat( Proxf(pVec(X),gamma/2), perform_sos_projection(Qmat(X)) );"]}, {"source": ["Define the reflexive prox operators shortcuts."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 19, "cell_type": "code", "language": "python", "metadata": {}, "input": ["rProxF = @(x,tau)2*ProxF(x,tau)-x;\n", "rProxG = @(x,tau)2*ProxG(x,tau)-x;"]}, {"source": ["Initial point of the DR iterations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 20, "cell_type": "code", "language": "python", "metadata": {}, "input": ["X = zeros(2*fc+2);"]}, {"source": ["Parameters $\\ga>0$ and $0 < \\mu < 2$ of the DR algorithm."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 21, "cell_type": "code", "language": "python", "metadata": {}, "input": ["gamma = 1/10;\n", "mu = 1;"]}, {"source": ["DR iterations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 22, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Y = X;\n", "ObjVal = []; ConstrSDP = []; ConstrSOS = [];\n", "niter = 300;\n", "for i=1:niter \n", "    % record energies\n", "    ObjVal(i) = f(pVec(X));\n", "    ConstrSDP(i) = min(real(eig(X)));\n", "    ConstrSOS(i) = norm(perform_sos_projection(Qmat(X))-Qmat(X), 'fro'); \n", "    % iterate\n", "\tY = (1-mu/2)*Y + mu/2*rProxF( rProxG(Y,gamma),gamma );        \n", "\tX = ProxG(Y,gamma);\n", "end\n", "p = pVec(X);"]}, {"source": ["Display $\\eta_\\la = \\Phi^* p $ where $p=p_\\la$ is the solution of\n", "$\\Dd_\\la(y)$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAASWklEQVR4nO3dbZPa\nuLYGULiV//+XuR+YQ9NAg7Fla7+sVVNTOTk9iYyl/WzJHed8uVxOAJDN/80eAACsIcAASEmAAZCS\nAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACk\nJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgA\nKf2bPYB9nc/n2UMASOlyucwewgfFA+yU4R4ARJOi+3eECEBKAgyAlAQYACkJMABSEmAApCTAAEhJ\ngAGQkgADICUBBkBKAgyAlAQYACkVfBfi+Xz2/kNOv1/mFm1KRB4bZFFqB3Y+n1O8gJIDPMyEUBMj\n8tggkVIBdrlcNLP8JXJORB4bhFXwCPHBc2l4DrlX5ePna65fvu7X8TVTvuZ8Pt/fwdtPzxrPy695\nHmGoz9DX7PA1j1Vl9ngevyZdI1XwcdH9M7DVz8Ou9/Fy+fkBidwtw+uP/rt/Qe7j73z9GWGQ4bGH\n+0qSoqqk+GaCUkeIo9xPr/sJRwpP6XX9wfkU6D5e68L5NkLTDFYQYK+F7zz44OEO3h6PTg+J/w3g\n8bsQTbnCHrZcQaZiAQUDbOO293lWqSyJ/N49X/6XDb+OECMUjuexXUUYG2RRMMB2orLE9/LRwkNC\nzM2wp078xdgo6Wlank6qymYC7AV1hD0sfHSvqBXjhu5HgP3y11QTafEt/86uyM2vmVbSy9vqXm8n\nwL4Qs+SRwlffOW2mwRIC7JG2KK/l9y7yJoxKPs4xk3ALAfbj/UwSbJG9uXd/vVzgyBv61/br5djM\nNFhIgFHEuroftv8NOzC+9WZmalY2EmDfUVYCWn1TjikfKd4bxESqymoC7BdVpqGA5cM8rCHg1CpG\ngH1BWQlr9a055p6aObAHAfYfvVJSQ27cfnd/y69sThbwsXfR3GwhwL6mrBSzdwVRoWAnAuyHQpPU\n+xu38OXOe/QlH3/NN2MzG7P7akZpi9cRYCQ2atlLC8hIgJ1O39RBla6wsV3wkO+e15jDGwKM3Ea1\nFFoThls4qcy91QTYGvriCKY8tfrq11GYYFcCjOL+ehdiBO/HJv/yWjHpAs/TuATYfxSLjMbetVGv\nqB+7/VLX4C8CDICUBNjXHa69WgQ77Uu2b8I8/WIFE2YdAQa/hColoQbDV9y7AwiwlTyZmG7XArHu\n/tp+wZEE2Omk4vBbtPmgW8pl9f1yo78lwMjnq3W+8F2IG3+X06rt1+qxAScBBi+N+pZ6YD/dA2xd\nhdI3TxftFnj6xUYmzwrdAwz+snwTtutGTV3LyF07hgAjmSOP9b46SFSz4GACbD0PSFIY8i7EN7/G\nlsPD5WMz2bLYeKfc6K8IMI1zPkfestvv9bKyePQFEwkw+OBlhp3P0mu866dqF8JC/2YPAL4wq7Td\nHoY9DEB6DfTQH/hs+ah1gG2phpeLPrGdWSW1w2R72M5ee4WkGbZ62B1u9FiOEEkmaVEDhhNgkEP5\n3vy+NfEmFJYQYJtYYEda+9qUuFu2yGM7Upl1NORCynwaB+geYAoIxGRt8lH3AAMi+CuubEd4Q4CR\nSduuvPCFiyhW6xtg25dN4ZoCEWRcYhvHnPGSJ+obYOSy4W+5jdvhRx7bkVRt1hFgkEa9vJv+t9WQ\nmgAjDX06cK91gA0piNpD2E+irmVgKVBVFmodYGRhPReWKKKIRoBBDp0LvQ6GlwQYMIdYYqOmATZq\n5XRuig+24a+oiHuTIo8tjkQf0pChJrre6ZoGGCRl1wI3AozolOzC7DbYQoABCXTrY7pd7zp9A2xg\n62eqwbeKrZpil5NF3wAjkS3dRuT3DX47trYHbm0vnPcEGAApCTBCC7x9Yqtv91UmAw86BphlAAw3\n8JzTkelCHQNsLFNtbz7hBwU6sAKXQAQCDICUBBhx6dO5sRHnmQCjuMjvG4w8NpbbqdPSwH3UNMCG\n1w1TjWOUibx1F2Khca9pgJFFmXoNDCfAAEhJgBGUw6KqSt7Z4UcFzh6WaBdgJRcPb1R6F+Ldfzh2\nIDmo6TxoF2B7sK4AjifAiEtnwLOeu09eEmDA0Sq1JrsGqrR+r2OAVVo8VVm3b5jAcNUxwAAoQIAN\nY9Mwln1GSRuXiVnBPQFGcZHfNxh5bBCfACMce9klfEpx7NSHaG8+6hVg1jwUYCFz1SvA9qNXgoUs\nFkYRYESkxgEfCTCKK/kuRII44AaaI2+0CzCtfXCW6xKdp3Hna+dBuwADZtGdMJYAG8n6BDiMACMc\nZ0R8FKdZNF0nEmCQVZwizk6k43uNAsxqj889Kk9FZqBGAUZPkd83GHlsEJ8AG0YtGsLHyEcmCVcC\nDGCNw068Ha3/pVeAadwis0q/km4yu78M1yvAAChDgBHIHruKyO8bjDy2+Hx4CLDBLCroI91BbjEC\njBAE/zq5PjflfgUf2htdAizXOgfgoy4BBlRiX8JJgBGHkgR8RYCNpASv44CXdA6etNbIS40CTLr0\nFPl9gxvHFvjKHu1Uf5X15hoFGJElqsVAEAIMgJQEGJM5BQLWEWDAEYafEk8/dj5sANOvNKwWAeb7\nhYLbdX1Gft/gkLEFvj7YV4sAA6AeAcZMdg/AagIMSGxKD6TxCkKAMZkH1OUp90P4GJ91CTDfLxSQ\nBbmd+UZnXQIMgGIEGDMdsIEo/C5EfH7NCTDmcH5IarIzAgEG7E6538gH+NK/2QPY3azvsi0w4Z4/\nurEXVeAjCqLGfINv2YHx2jW9Lpf//rn/yVG/OAxhOrUlwHjhll43txhLVyzKvwuRg7lpcQgwHj2n\n183ADHPk1YRyP5AP80GLADu4VqYuzW/S62p7hlmEwBAtAoyxtmdY6oyPpvmH2fzymxNg/Pi4/bpZ\nnWG2X7COqH4mwFgp6fd0wHZTssRaeybA+GXFyly+rpbv8KjEHR/Cx/isfIC550ut6O+G//mwPUR+\n32DkseUSeQayn/IBxhdWlNPl/4nt165UcBoSYHvJVVC2f0vh+19BelFDrnVdXv0AUzQPcMuwl8vb\nmodRDlxNCUpn+QBTO5famPR/PQ+7/U+dREN6F3ZVPsD4bFSVuX9f4u2f2/81S+T3DUYeWyJ6o+Gy\nTMz6f50KB1NNjne5pKk4MJAd2C7SFfF0A4ZZLJY4BBhADrLzgQDrztETuzqs5prJDQkwAJ4l6AgE\nGM4lirAF2ZtPOBoB1lqHBRn5fYORxwbxCbAddYgHiKBVJ6Cw3AgwYBfqLHsTYN21al2BjxJ1HgKs\nr0TTFILQ8IUiwCgu8vsGB45NYT316Mnc6HsCbC/mGcCuBFhrUhZ4lqUyCDBgL1nq4BIdzifTEWBN\nWY0UUyksWUiAQR36kibc6CsBti/zDGAnAqyvJkcukd83GHlsPOtwu3L13AIMighVXmfVwVz1l40E\nWEcWeUluaxOhOpW5BBgUoa7RjQAD4EeiTkiA7SjyPIg8trGavAuR056z2o2KSYABu+jTJDGLAGtH\nLwkFWMin0+nf7AGMdDuQ8cdraOt8tvWhizoBdj6fb7l1/2OglZ0ivENFSberc4S4u4BzosNSBMqr\ns1N5uQML8F1e1yFNH8a9S7DxMFCQ+TZxGDv91gFXzR5D+vXpxU+HOkeIf5l7D64BGmceRBsPYw28\nv1vO4SdOs51+6/M53KrZY0j3n16A7v8zR4hAHXukTIZKPkywmP6gzg7scrn4LkSau1x6VVuaqxNg\np5C5paAAO/FHJhwhttN8xtOBrrEJAQaMV6xPKnY5L2VMfQHWSMYJyiwBD+ThgQCDanQqNCHAgFKa\nbB2bXOZ7AuwIOmJIyuKNTID1omtjbyp+XunqgwADXkjxJiGaE2BdKEdNpGuid9Jkwje5zL8IMIB3\nOvQESYNQgO2uw+wHOJ4AAwbTtHEMAdaIskITTaZ6k8t8Q4ABvJb0ydA6GeNQgLXQah1yGnHHvQuR\n+ATYQUQIHKzJomtymS8JMABSEmDAMPV2Ax1OUvPeNQHWRYd1yJV7TRMCDHgh+7sQm6R4k8v8iwAD\neCF5grcgwI4wt0uyDoH3ku7kBBgwUqhSqHurTYAB5NY2pwVYC6GaYo7RtqgN1GHhpJ4nAgyoqUP8\nNCfAjpO606Gb5u9CtFpTEGDFWYccxmSbonOnIcCAysTqR3kjUIABkJIAg4Ly9tRx5PoMe240BdhB\nJi6GXOuQILK/C/HK5H8v+00WYACkJMCAYWrseLLvS/oQYJVZh3BqsBBq9A0rCDAoq3zhZrvU4SfA\nDqWgQAoZy3rD8iLAgMoyRhELCbDirF7W+fZdiGXa/zIX8lGBKxVgQH0FijXPBBjUZPPdSs/bLcCO\nc/AM03LCOn3CIPuVCjCguOxlmr8IMOCFFe9CLJAT2c8tso//WwIMaKFbcX+vxqchwCor0BGzUY06\nxRIN17sAO5qCApE1jIG8BBhQn1h6VuAzEWA12edxpCzz7f04s1zFezWuYiEBBtBLmZATYMAL374L\nsYaWF52YADuU5cGRzLd7HT6NDtd4T4CV1W0qwxJ/nZ6VOVU7LbuWGvVBgAE0UimqBRgwQIqO/q9B\nXmt6ikvgngArqFKHxSwr3oWYRd0rO52axbAAm6D2+iEa822JYnX//U0vc7ECDGjkWrvv63urgC92\nsQIM6Ohayks+/Sp2OW8IsJr6zGCmS9fU31ZHupEPUak4CLCjVZo9kNRtGV4uZZfkczzXC2wBBpVV\nrc7bFY6uU5v7LsCqqddkMUXPdyEW8/y9KsXuqgADICUBBmxVrK+v4f4PDJTcfp0E2Cy7HvTVm6bA\nCvcZVrIs/Js9AGB353PN+sVHte+7HRjwQuF3IVKGACtFzQH6EGAT1N7U04qeiYkEGAApCTAozo6f\nqgRYNaoV0IQAAyAlATaNp99EtvxdiDb9zCLA6pCIvGF6UI8AAyAlAQZASgJsDo8NADYSYKXIRUZZ\n8i5Ez9WYS4BBfTobShJgAKQkwGYaeALjMAfoRoAB6zmcZCIBBkBKAqwOvTDvOWemGAE2jbwhsuXv\nQoRZBBgAKQmwChwNcTyzjukEGLTgRJB6BNhk2liAdQQY8MKSdyHCXAKsCAdEHM+sYy4BBkBKAmym\nIQ2skx6WM1uoRIABkJIAAyAlATbf9lMdz9I5mKNIIhBg0MVXjY53IRKfAAMgJQGWm5McoC0BNplz\nGpIydZlOgEEvdu2UIcBC2FJTNMLswbsQiU+AAZCSAAO+Y29GEAIsMXUE6EyAzechFocx2ahEgAGQ\nkgDLTUMNtCXAovBAi8MsmWzv34WocyICAZaVwAOaE2Ah6GcBviXAEhN7HM/WnzgEGPSi76EMARbI\n8t5WF8zevAuR+AQYACkJMOA7DiEJQoBF8W1RUETYwgEhBQgwAFISYLEs6Yv1zgAnAQYsp3kiFAEW\nyPLHWh6AscWS+fP+XYgQgQBLRgsMcCXAwhFRRGZjRhwCDJrSKpGdAItl2cOJ/ccBEJ4Ai+iv1ljL\nzGGe34Vo+hGNAIOO7OMpQICF876yqDsAVwIsqOfjGgc4APcEWCa2X4z1bVdkBhKKAIvoWibui4vt\nF8ADARbaNbeu/9b8MpEWioAEWFC3uJJe7OTTtwuZc0T3b/YA+JMCAvCGHRiwiI6KaAQYtObhFnkJ\nMABSEmDAC/fvQrRLIyYBBn15rEVqAgyAlAQYdLfkhNBejYAEGPCOB2CEJcCgNVsr8hJgAKQkwIAX\n54T370K0SyMmAQb8yQMwIhNg0J0NFkkJMOB0+nuzJd4IS4ABkJIAA15ss87nswdgBCfAgP/cEuuW\nXc4PiUyAAT/OPzuvy+n3O+khGgEGnE4/m63L3b9PJxlGYOdL6TOCl2vv+ZKfv8zX+JqeX/P7p+eP\nx9dM/Jr46VA/wGpfIIy1sOejvBTF0xEi8E78KkZbAgz4Ia5I5N/sAQCxyDCysAMDICUBBkBKAgyA\nlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlOq/C9FfxwdQUoK/8QUAnjlCBCAlAQZA\nSgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYAB\nkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJg\nAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQE\nGAApCTAAUhJgAKT0/0NOTK4ryKG1AAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 23, "cell_type": "code", "language": "python", "metadata": {}, "input": ["etaLambda = PhiS(fc,u,p);\n", "clf; hold on;\n", "stem(x0, sign(a0), 'k.--', 'MarkerSize', ms, 'LineWidth', lw);\n", "plot([0 1],  [1 1], 'k--', 'LineWidth', lw); \n", "plot([0 1], -[1 1], 'k--', 'LineWidth', lw);\n", "plot(u, etaLambda, 'b', 'LineWidth', lw);\n", "axis([0 1 -1.1 1.1]);\n", "set(gca, 'XTick', [], 'YTick', [0 1]); box on;"]}, {"source": ["Solving the Primal Problem with Root Finding\n", "--------------------------------------------\n", "Following [CandesFernandezGranda13](#biblio), one can hope to solve the\n", "primal problem by using the fact that the support of any solution is\n", "included in the saturation set of $\\eta_\\la$\n", "$$ S = \\enscond{t}{ \\abs{\\eta_\\la(t)}=1 }. $$\n", "A key remark is that this set can be computed\n", "by finding the roots of a triginometric polynomial on the unit circle\n", "\n", "\n", "\n", "We assume that $\\Phi_{x^\\star}$ is\n", "injective, where $x^\\star$ is a vector containing the points in $S$.\n", "This is obtained in the case that $\\eta_\\la$ is not a constant polynomial.\n", "\n", "\n", "Display the magnitude of $1-\\abs{\\eta_\\la(t)}^2$, which is a trigonometric\n", "polynomial which is zero at the locations $S$ of the Diracs of the primal\n", "problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAVxklEQVR4nO3d23Ii\nSbIFUDjW///LnAc0NC0Q5CUu7h5rWdlM23QNciI9fEciVdb1drtdACCb/5tdAAAcIcAASEmAAZCS\nAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACk\nJMAASEmAAZCSAAMgJQEGQEoCDICUBBgAKQkwAFISYACkJMAASEmAAZCSAAMgJQEGQEoCDICUBBgA\nKf0zu4C+rtfr7BIAUrrdbrNL+KJ4gF0yXAOAaFKc/n2ECEBKAgyAlAQYACkJMABSEmAApCTAAEhJ\ngAGQkgADICUBBkBKAgyAlAQYACkJMABSEmAApCTAAEhJgAGQkgADICUBBkBKAgyAlAQYACkJMABS\nEmAApCTAAEhJgAGQkgADIKV/ZhfQ3vV6vd1us6tgvuv1+vjngC3xKE9tcEypO7Dr9fo8s1jZr04I\n1Ri/GlVtcEypALvdbg6MZBQ5JyLXxuIKfoT4y+v2ew25fr/ncrldLpfn3xihnueS5tbT6fe8nbmv\nny3PvhZvKjz5tf73L9vU/PqbN77O/X+b1GM//930a/3eNa1rPvB73pTUr3/CKvjtouc5Nff7Yc+t\nEmGZf7VuhJL6eRsYQbq9U21N+q1JbY/XGLzePXZctF186VbSr5dN8cMEpT5C5IPHofj+6/KSZ+VF\n3o2Vavt1+zVYp/YOuGX6vc1EBFhfj+0UTcyqGvo1dkMlRL9izvdb29pGTvx+XytS7/xoXtKvm+ZQ\nUf1BwQALMqoCdsDrwgQssqHHD/UEaYlnz7X1KO/MlT1f2/T19vnhIgoGWExzo+LtV8911OIz1/Hy\n5ocs2n+JgOvcqqSMcSjAgGbq3R/UeBef5f3xLgHWUag+eFtMqAo7Cfj54UOn2pq8auR1e/V6F3K+\n/IA3W3fPby3VVWpPgI0wt8m+7sOwG5VdYs6ySt0Vc4WbePfWErxbAba6wnuSKSb+AP2rSvF5+fvt\nFHub2wmwLqL109eZ0rbg6/XfX8xi8c9bORRTEGDjhG2Uhj9z/BpaYmwAK/yWTxcOy7J05QMsyXXo\nZvtoOzkEn/8g5K9fl6kxFvlZtMVq6/GTFMe+bszX7CRL2PRQPsAyNWJez8+pevX0UNdxJcFdsa7r\nEVd/vGaChasfYLPkOhYd/kPNj1urz+8312oUYMHPy/UnT4rl9EblAyzKVZ3Y9/2+9K6HjnvwR1dh\nB+sYi7/9w7Lvx/IBxlaHA2b77JBha5p4xY8FW8wWjVnVXAKMg+b+rRnEtEg/RMuSRZb9lQBrL05z\n761k+x3S4fRyE0Zb29uVepYIMO273ZaAyXXvFfmZfr1rO/mXqrQrJJ913n3qd7pEgMWRPUo3/szh\nZ6k3TEDBm2r65Z5ewDDBO6EHAVbcmU/5Xp+pcfg1Xy242aihzE8UF9iD/8wugIhut5/mnvVsBdJJ\n96OGfykw1tfhDoz3Xh+r8deDNo69OIxUPpbW3FP1A2zKdY3wZ/ib7NhfjzRsa8xMKfa8wbd69Fvk\ndettzTDIqH6AEZMZsZSF07CNTguYfRsKMCCrjWM97+cHe8Wsqh8Bxkyr7bfash/nLys1ZI13KsAq\nKzBQoJUaI5tnqwTYsN61SbaTr+d97TcNuY4FN9QqARaHgcJqIgzWvX9nQqtXayjCMkYjwJisd6JH\nfqaf2s5LUmZEBZZOgNXkPg9ebXlKdXY13sVGAoyZCpwBuSw2NAsoc70EGAApLRFgg4/57ir2KnMe\nnOJDv2nFZ2F/OqOh7PXvtUSAxTGyvVZr5b9EfqbfUrU1f6+BF6+xdd7pXgKMyQQtI/0VBgdCInWu\n1Nh3AgxooMZArCF1su4iwAhhnS1XVYor+DllF8ngFFdqIwEGNLDI9CeUhQKs0rnjs3Xe6eJc6GOp\n+bpulVZyqZPEQgE2wMZtUGm3EF+EfltqqvbQdgHLXA4B1lKZthis67pFfqaf2sb7623tfbvBlyfC\nqWUAAdbSIk0D2T1v1aW2bbE3K8CIotjWAnoTYDUF/3wDTjpz3Lnvjvsr3P+z2H4p9nY+WCXA4lzR\nOJVAE3lvnfNWfkalEbRKgBFcv01V+3mDudYtzqV4XrdKA/1ZnNXu55/ZBVRTdTOQ1O22xCA7YMGt\nWq8T3IEBxHUsdRaJZwFWTepDVuriV7bIuHwVvGNfyyt2pQQYcNaUOV5sFvcWPGuPEWAABa0Q8GsF\nWMkzCEy3wqxM6nno1btMawUYkXXaXZGf6Xemtt6nsebrFvg61Ff17C7A5ujaTyYFDDBsox3+Qs/P\nHCk5FgQYpLRrHlU9gPPV7fbzqyQB1owZ0YRl3GjXQhWbX5qEOwEG9Zn4lCTAKK72sxD7iVwb3C0U\nYL0/RSn2KQ30IBZpaKEA6y3CzoxQwxkOAbtYLhYnwCboN3dMNAbbe2bSortkP5L2JsBamr45tTvB\ntWrR6XvtwZ/pnEiAtSQ/mrCMbRmCedkLnwmwlkwKWI1dP5EAo7iqz0LsLXJtobhJmkiAVWPsEJbm\npC0BRixmXDouGbMsF2Du9ykgVBsPLibUexfecy0XYHGE2odwmCHeleX9QIBRXORn+qkNzhBgRGR4\nAl8JMFiCMwH1CLA2TAeCi/OtlPObJc57Ya61Aqx23wtRPtAe1LNWgAFQhgBrY++9Xe17wZMszha7\nVsmSUpIAa8PnM2FFfqbfsNoO9GfkdYM7AVaKmcNIzm13ndbB8n4lwKA+JxtKEmAE5fjZUKfFPJCL\nopSGBFgzdiZhdWrO8YeMpY41RspXAoziIj/TT2184Ap8JcCA49wlMJEAm6nhCcthjSk0XtcIdz74\nbMUAs+Xis28/CNXArlTXyxHqWge0YoABrZiw7sAmWjHA9AS0Mn432b88rBhgzoykVmCC24M0sWKA\nBVFgDPXWZMxFfqbfsdrGTP/I6wZ3ywVYnF3pEMpgWo5ilguwHns4yFyIk830duxa6xCKWS7AejBN\nGOzYmSnISQtaWS7AesSGucBgqU8/qYsnlOUCjCxajbnIz/RT2y7xKmIyAdaAEyUsrnm4SustBFgD\nWg1oy7F4CwFWgQTlq+YDUdd1ZXm3EGCwBAORegQYoRm7rTS/A/MZ14OlmEWAtXG4gxs9LanBi8Au\nU84WMVvdMWsWATZZkw1p/3wQ+Zl+K9emaTlv0QCLs3niVEIKZRomcHZHYYm+WjTAdAbMVSaJL93m\nSaUl6mTRAAvSGQ37vmQkl3xTlbhAXVnerxYNsCCC5Cgco4GZS4BRXMBn+j1kr80tAnMJMCCBwFnP\nNCsGWJxjY5xKYCSdTxMrBpijHEABKwaY018uDhy/nGlgi0klKwZY2z08fSJMLwBgihUDzB0Yy9L8\nVLJigLGUlZ83+Gr7/XrAdYtX0X80/CzExyobCTCAWIJHdRwC7Cyt1pXlLcx9xl+szEYC7CytBjCF\nAJvvfAS6TWGLhn3i3NaVHb2RAKO47M8bbPrldv3muOtWnrXfSIBBGifnWqi/vmfXK8Qf6O6ZphBg\nDehdxojTafETZTALMsW6AabhEnGx7qzDOuIcViJbN8AAWmmeNw4rWwgwSMOpfB2u9RbrBliQ/jhZ\nhmPaUkJd7sE7KMiGJZR1AyzULOCDk5Mr4DP9Hg7UNuzdfK1t8A6yYXm1boDVEHg400WcOe4OjOkW\nDbAymyHOOAMYbNEAM/dJqszZ685O5IxFAwyA7BYNsGLHWD6I/Ey/KbVt/JqR1w3uFg2wViLscWFM\nbRF2GTEJMHIwxYhPlw4mwELQ91BAk49DTIPtBBgsJNQHzqGKaUL2DCbATqm3A6nNhI3PVNlOgJ1i\nHABtmSrbCbCzHJeCK/YsxGE+1DZ+wgZeJ2YSYGdNPC6tc1Izvy4rXe6kdOl4AuwsXUsi2rUfJ4zx\nBNhZ57vWTGFl5v4rM2GjpQOsyc7RaowRp9PiVFKVUN9o6QCzD1cQ+Zl+u2ob/hdI/vn1RlYS+Or9\n1nCeGE0bLR1gifYGAL8sHWAA5CXAclvqowZ3zHG0arylGpjmBBikEWfcO0wQwboBFmcWQDojt4+t\nyl/WDbBoR8ho9QAEt26AkcvhY3jS5w1OF7k2uFs3wGxPlhXtdj9aPWSxboBlZ89DQDbmSALsOJ0K\nXS24xRZ8y2cIMIAofGtjFwF2nFZLocyzEFvZ2LeR1y2yk2PBqu8iwI5r2GqycCPb+7xoa1is+aMt\nb20CLLFiO58Plj0t5aqWwQQYrMUtAmUIMFiLb9IE56ZzOwF2ilYbaeXVbvjeYyZQzKr2anKZaizF\nGKsH2Mle0Wojrbzacb4HNuwYkfFyZ6w5tdUD7KSV7wnGO7bakZ/pl7E2M/qDJtczcFOEI8AgB3MN\nflk9wEINBWdbPijcHqG2IYmsHmBJh0LSsjljwSm/4Ftml9UDDNhLrhCEAKO4yM/0S1pbp6qfXzbw\nwhCIACMZo206d2AEIcAASEmAReFUy8pe+9+O4KulAyz1DkldPNPF/CQ2ZlWEtXSAkc702L5ef35l\nN30lP8i+vNnrT0SAHaRHF1Tpx+QC1h85U3c5/EYCXpTgBNhBZTZbec2fN3i7/Vz98+NmY2095trX\nrzzrOY33N5t6f8mhYQTYQXp0Nb8Ga6sMmyVm5fdVTZ1eZyz7xg8TYLCV+TLAyosc81QRmQCLRQfH\n9Pa6rDxqIQIBlo+Qm+KvuHI5aMipaBcBdpxWS+H88wb7RdSsZyFuad23tUnrr06OBSu8iwA7TquR\n9BBzuHWTvt+RTo4FK7yLALOZ8xl/dPh8rdMdZfxBpX7cgY0kwGxmEnBaWoQLvYsAC3QHtv0Fdfkw\nX9vDtYBZBBjJxAwMt+MwngCDL2JG5ixWgzgEGMWdeaZf7/uqWc8b3OKv2txrEsfqARZ4gJCGLoIp\nVg+wdMfJdAUDdLJ6gAUkokLZfne1yIVzu0kcqweY3QiQ1OoBRnmHnze46/937CS0pbaeT2L8/G/X\nuKPsw+KNIcBIyYA4z8cPoWjpAwTYEXNbzdwZpvZSm5j91O6cOARYLPo+tVyRoNn6OdAJLscBAox8\nbPUmcsVtIke/Idq6jgUIMPIZsNWdoDlMFA0jwI4wquZabf17vN/V1pCSBNgRs05YTnYHHH7e4IAP\ngjbW1uO6b/hrYkTcQVZuGAF2UNceFVQ8CzIQteVGFmoYAXa5HGo4PUpqQUKxJGs7jAC7XDQc/+V0\n8pZt0psV3kuAXS6pfuRMi0dW++rIdaIRYBQX+Zl+uWqrHc9kJMDCMSY26jr8z1yFwKl0StX3RV4C\nDICUBFgazr9jWGfIQoAFZYxy0QbwkQCDlmp/C7P2uyMdAWZPEpfmzMvd8wACLBPj7GH7Uhx4pt+w\ndf5aW+8h+OH1PQvxpF3rJ+2OEWARGR3ZZZlHOq2fnY917lZHaQJstymzKctAJBd9FYQLcYwA281Z\nqTBzhCZMiTEE2G7DZpxhmlTb4dVvFBqy/di8YwiwNIybY3I9b3DsV//8b+OuG9wJMBLrMWPXOSis\n806pSoAF9TxcHIUzin/V4leY14HDgfPEAQLsR8zNHLMqFqQVe7PCBwiwH7uOP+PPSk5nA6w2QTRV\nKC7HAQLsx67hNWbS3Rt6takKsJEAO2LwWcnR7K3IyxK5NihDgF0ugcfN7Ra3tiwiP9Pvc21zb74j\nrxvcCbDLZfakIA5DGxIRYNDLyYORNIXPBNjlYlIQkg8G4DMBdrmYFJm1unZ64CvnPKIRYJeLnVla\n5Gf6qa22jUtopQ8TYNBFllOR6TldllYJSIDBj4BzJGBJNOcMcZgA20erAQQhwABmcp99mADbR6uF\n0upy9LuxdsvOV5rkMAG2j1ajEgcyUhNgFBf5mX7Ta/twIJteW3bB/4amGgTYv9xd0VaZqWRr9GaF\njxFgcLnEC5uRf+cc07kQxwgw0nN67c14JSYBtpvNzGocEYhJgFHc12f69Z7OH14/8vMGI9cGdwIM\ngJQE2A8fDALkIsB++LwkqcInj8JvbR0GS1cC7IdhsbJ+V79MX5V5I9FIuDMEGEAvX4PfyeAMAQar\ncxPQz9e1tfhnCLAdtFpGn5/pN+aa/vVVIj9vMHJtcCfAqMDZgoCcAXoTYDtoR4A4BNgOjvnU41jW\nj4nRmwBjdb0n+IHXN/hq2HjpnSEOE2D/YXBk9Hn/R36mX4Ta/iohQm3wmQCDpX09/gsywhJg+7jZ\nr2TkaBYD0JwAg4gcleArAQZ8IU2JSYBRxLHP6IzmLXz+eYbV60eA/cssoxOtBT0IMIqL/Ey/LLUF\nLpOlCTAW5YMdyE6AbWXeReYW4STtPYVlP0mAwSAbp5WhBhsJMABSEmAU9+GZftM/eIz8vMHItSUy\nvcdqE2DUEXnkJh1kkZe0hqSNEYQA+82OXUHwqzx4qJmhJCXAdrDPAeIQYDBO8Ds/yEWAUYT7Y1iN\nAGNRMQPPLVpJby+ra32eAPuPmEONM7I8bzCaR22Ba2R1AoxStpxqZ518gyeBGwLSEWCb2NsA0Qgw\nGOrrYSj4jRrEIcBYkZBgGM3WjwB7wweGSb2dFJGf6Re8tsDVFSHbThJggElKSgJsKzs8i8/3DXPv\nKj53kTse2EWAAXTndNKDAPvNnVZ50y+xWYYeaEKAAX+aHvbwgQD7zlkpkezfZJobGPHXB54JMIqL\n9rzB53Ki1fYscm1wJ8DecxRNLePly1gzG709DDghnCfAWEuQqRE/ruJXCAKMaoJEVDrWjXQE2Cb2\ndgHxbym0WW3xOzAdAfaGOVLAY1gEfKLfo8EC1havooIscisC7AutVkmoo8lza8Vss1DLBa8EGAW9\nTt6YCcE6nAZ6KBVg1/9p9GpNXoZpIl/B13EWYcA9fbY5tY4FRLjcBfwzu4Bmrtfr449ePv/zeVot\nu/s4DnkdbxdpsRLXuq1Sd2AN3Yedbsvr6QpGDK6HgD32v5LiVQb/VecO7C8nPlH8936uVTEM90iv\nmJcxZnk6v5/72lraNuoHWIvPEkMf4dkm+EUMWF7AkspIsLYpMtZHiACkVOcO7Ha7PY4MHqQNUF6d\nALvILYCV+AgRgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQk\nwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASv/MLqC76/U6\nuwQA2rvebrfZNQDAbj5CBCAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIM\ngJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIA\nAyAlAQZASgIMgJQEGAApCTAAUhJgAKQkwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKQk\nwABISYABkJIAAyAlAQZASgIMgJQEGAApCTAAUhJgAKT0/5rHUaV+5/2YAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 24, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "stem(x0, sign(abs(a0)), 'k.--', 'MarkerSize', ms, 'LineWidth', lw);\n", "plot([0 1],  [1 1], 'k--', 'LineWidth', lw);\n", "plot(u, 1-abs(etaLambda).^2, 'b', 'LineWidth', lw);\n", "axis([0 1 -.1 1.1]);\n", "set(gca, 'XTick', [], 'YTick', [0 1]);\n", "box on;"]}, {"source": ["Compute the coefficients $c$ of the squared polynomial $P(z) =\n", "1-\\abs{\\eta(t)}^2 \\geq 0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 25, "cell_type": "code", "language": "python", "metadata": {}, "input": ["c = -conv(p,flipud(conj(p)));\n", "c(N)=1+c(N);"]}, {"source": ["Comput the roots $R$ of $P$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 26, "cell_type": "code", "language": "python", "metadata": {}, "input": ["R = roots(flipud(c));"]}, {"source": ["Display the localization of the roots of $P(z)$.\n", "Note that roots come in pairs of roots having the same argument."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAQZklEQVR4nO3d3XKj\nuBpAUXQq7//KOhc4agbbhD+BPrFWTfUknUwPwW62JQROOecBAKL5390bAAB7CBgAIQkYACEJGAAh\nCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJ\nGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAh/dy9AT1IKeWc33+zfPz+VQAO\nErBDppV6p1sA9ZhCPCTnvFCplNJy4QDYzQisorFt0wlGPQMWmLbZRMBq+fZEDPQE/Xhur2U2uKpY\nWzvE3OC7NyEYU4hVeCIC1GYEdrLxRV/OuTQs1mtAgCiCDbGjCzenAVzG8WErIzCeKKXBgeIC5tL/\npFhHCBjPUtKlYddwgF4g8AdZxMGz5DyMBw3HVYjOCIwHmb7eLR8rGQRlBMaDjK0qxcpZvSAwIzCe\nZVovWnbkJOXs3NLK83AWAYZjBAa05eNM71b51ymbRJuMwHgih7WWlYU2564Und1b4OOn+4Zu3EXA\nIJ6OrwGYDbl2Lxmdpej9ntrTT8vdc8qvO7ac6wkY0JByld70011/zs63gHAfuECcA4NgyvRax2br\nRa9UBmTTktEmIzAI46zptQd6H1d9HGaVWcTrt5AdBAzCOGt6LYQjP920SeXj2Xzg+/Tgt++kWaYQ\nIZKyfKMs1YPHEjCIZDo2ME7g4QQM4pEuGAQMgKAEDGhPSq9/9v8BLydu1PQPP/gNnMIqRKAx0xuN\n7LrpyOy+G1YV9soIDGjJrFiHV1vOrvoq9zycDtFmn5bvnP7O9JuH/46x3r/z/RsGw7IKjMCA3qy5\nHdTyrRE/juHKr99GeLPbKtb40ZgSMKBDTU0hSlolphCBlszmDPeeAztzkw5zh6pKjMCAxkwbtnfU\n8j6F+P4GK+UbpmewpvOEsz9h9p9v/U5OJ2BAe44d8Rfuc/jtdzZ9w6b/1i0W6zGFCEBIAgY8zp9D\nKEIQMABCcg4MqMjSO+oRMKAWU3NUZQoRgJAEDICQBAyAkAQMgJAEDICQBAyAkAQMgJAEDICQBAyA\nkAQMgJAEDICQBAyAkAQMgJAEDICQBAyAkAQMgJAEDICQBAyAkAQMgJAEDICQBAyAkASsrpTS3ZsA\n0KefuzegW9IFUJURWC0555zz3VsB0C0jsKuVkZm89aqMvWeP8Mcx+Y7voSemao4QsKvpVq++datY\n88ivfHaM/y9PpQ5MDwhitpWAwQYLlboyJ+P/a3q4EzMeSMDgb3+Orm7R1MbA9QSsLhOGfQj3MBqc\n8QQCBv+RUg9H/OmP0ObwEY4TMBiGro/y/f1EMBIw6GTUtdUzf2p64kJmeO5xPKXPV55BCEZgPIth\nR1H2g6vKCErAeAqH6W/KVWV2DrEIGJ3reHXGua7ZPzLJiZwDo3M5v/5hq3NPjznZxukEDPgs5zNX\neYx/2qBknMcUIvDVWafHZtFyPpJTGIHRD4vCKymDpyN/wjRXJnU5hREYPfCKvrbj+3YcxlnEwYmM\nwAhPvULwAHE6ASM881G32Dep6JHiRAIG7HHuGkXYQcCAncaxr4ZxFwEjGK/6W6Nh3MUqRCKxhq1N\nHhRuIWDEYKkhMCNgBGDgBbxzDowA1CscZ8W4gIAB57OygwsIGFCFhlGbgNEiB74+aBhVCRjNsWRj\n1MehX8OoR8Boi3oNvaSr8IBSiYDREPUaBXrz4va3kI65DoxWqNfImxfDSkZgtMIxehTlzYtLaLcO\nwgzaOIuAQXPKYLTNdA3DfIZz6x2WNYxTmEKE5jTbrWLcwn3Tm+W/bf/HpHFGYFxrfK1e/uG7xo/v\npUC735rZ489BAsblck5Dfv3qGBZW6evu0GoYBwkYF0ppyLmcPml8hLGG4+/BB7GD5wA3EjAuNT3i\nRz/6716GB5xCwLhQznlIk89u3JQTWIwA97IKkRu8MpaiRsy1xtACIzCulcdRWG73At0Vptve+AVb\ngZiJZSsB4wZ9HO6PL8NjyqJEthIwoBUaxiYCBvsZe8GNBIzqvKZmPYMw1hMw6rLKnK08YVhJwAAI\nScCoyPALqEfAqEW9gKoEjFrUC6hKwIB2WZHIAgED2mVVPQsE7Gr+NsImGsY3AsbJHGuAa3g7lUPS\n79E6v61YSGn6xlflq3noenlexz8aNxoHYZ5azAjYfimlUqbpx8X0d7yDFBzhLwvvBKyicRA2Zqy8\ne2/54nvwovMaGbZK5twPELCKxkSVwVlKw+tNiPOQkiM9MJunEbNtLOKo5dOM4ueP+2D4BVxMwKpY\neCXV61G+158LaJYpxP1yzu+rEMcJw49fAo4z1qcQsEMWVh7qFtRgST2FKUQAQhIwIBg3l2IkYBzi\nOALcRcDYz6kI7mIQxiBgQFBePCFg7GT4BdxLwAAIScAACEnA2Mn8IXAvAQNisxzxsQQMgJAEDICQ\nBAyIzUXNjyVgbONIATRCwAAIScCA8FzU8UwCxgZuHwW0Q8CATjhB+zQCBvSg1EvGnkPA2MD8Ic0a\nn5xmuR/l5+4NADhkNuQaP5WxJzACA2LLeZ4r9XoIAQN6UKKlXs8hYACEJGCsYmUX7TP2ehoBA/rh\nxr6PImD8rYUjQgvbQAjGYc8hYKxy40Fhmi4ZAwoBo3VlUsg1qsCUC5lpmmtUgW8EjKaV+wNNPwUY\nTCGyxu3ZGDfg9s0AmiJgQG8s9nkIASMGwy9gRsCA3ric+SEEDICQBIwlKXklCzRKwPiDk09AmwQM\n6JAXXk8gYACEJGAAhCRgAIQkYACEJGAscSacuFwB0j0BAyAkAQMgJAEDICQBAyAkAeMr58CBlgkY\nACEJWC3p190bste45XG3H+idgFWRUsq/QjYspdclYN4ZEGiVgPGm1GukYUCTfu7egMcpA7LsLhfw\neCFnaJohYFeL0q3y1yqlIcYWQ0DTA4KYbWUKkTc5DymNf61SGvKQ3BIRaJARWBXTtRtRhlxFSsMw\n5JzS7+kvCQNaJGC1hOtW8Tv2KqsQ1QtokSlEPivRUi+gTQIGQEgCBkBIAgb0yex39wSMr/z9B1om\nYACEJGAAhCRgAIQkYECH3FbwCQSMJY4CQLMEDICQBAyAkAQMgJAEjCXj+6kANEjAAAhJwPiDG0oR\nkeftEwgY0BBT1qwnYEArSr1kjDUEDGjFOO+XkglAVvm5ewMA5kOu8VMZY5kRGHC/nP/lavxgd71M\nPz6HgPE3RwSuMWsYLBMwAEISMP7mfhxcxtiL9QQM6IdXWo8iYEBXjOGeQ8AACEnAgH4Yfj2KgLGK\n4wIP5Ixa4wQMYG6aLhlrloABzJVLR9yYsWXuhQjwH27MGIURGNCDEyf6pjdmfP+UdggYGzgZwHMc\nvKcwFxAwNnBPKdrkTNUzCRjAZ6LYOAFjG4MwoBECBsRm/vCxBAyITb0eS8DYzPECaIGAARCSgAEQ\nkoABEJKAsZ/19NzLM/DhBAwIyep5BAyAkASM/dyVA7iRgAHxmD9kEDAOMggD7iJgHOWFMNfzrGMQ\nMACC+rl7A2JLv9Nn+e0FYZrMrL1/FYCDBGy/lFIp0/TjQrcA6jGFWFFKKVnhAFCHEVhF4whsNjhb\nmHUEFnS5dN5r3CMEbK3Z8+zP/Hz7hr671eUhBur5+OqWlQRsrU3h+XhKDNjNayPeCdh+Oef3+cAx\nXR+/9ATjdc1P+om5gicVHwnYIQsrDx/VLYDrWYXIydxcitN5NchHAsb5HG6ACwgYACEJGAAhCRgA\nIQkY0CJLgfiTgFGXwxA7uPCLNQSMuqyqZyv1YiUBozoNA2oQMK6gYaxk+MV6AsbVlIwF6sV6AsZF\npoMwDQOOczNfrvO6V/2QhmEYf/F6G9hNwLjIOOp6JWx4feyMB7CbKUQukvOQh3+5ytnwixdTyuwj\nYNzAokQKg3B2EzAuNAmXhjGoF8cIGFf7XcQxnhNz9Hqc6WJUjz9HCBjXyrkMvsZ/G4c9x/SxVi+O\nEzDuMGbMAOxhyryxenEKy+i5mQPZQ8yG2qaQOU7AgCuMrSoZky6OM4UIXGfslnpxCgGjLdZ09Moj\ny+kEjLa4Pqw/6b+3DDP84iwCRnM0bL32d9SYLtGiBgGjRRr2pxD7x3J5qhIwGqVhy6bXVDVLvajK\nMnra5fD3jWuqYBAwiMg1VTCYQoSgyukl9eKxBIww3Pl3qsFueXS4mIARxuRG9jdo8+jcSMZmV3rB\nNQSMYK5/E5bZm4Aw40ov7mIRB/FcfKz0JiAL7BNuJGCwxIL1ZXYFNxIwWGLBOjTLOTA6UfXslAXr\nI6cAaYqA0YnrF3c8jdNdtMYUIv2YTvedfqh98rHbmT/aJGD0pmTMAfcU9iTNMoVInxxzz2JP0iwB\nAyAkAeMpLPH4k11ELM6B8RRVl3iE5io3ghIwnkXGZqzRIC4B44kcsgu7gricA4OX7k//OMVFZ4zA\n4GX6ZmM9jUu6/KFgEDCYKof4bk6SOcVFxwQMPujmoN/NDwLvBKyulFJ2COnF+wmk2x/b2Sbdvj1w\nJQGrJTld3p33PNw7QWd6kIcTsFrGgZeM9e1j0tZ825/W/DnqxcMJ2NVK0kwtdmnNo7pmKtKz4yG8\nxj1CwE4wewoul0m38BSgmB4QxGwrATuBJgFcz504AAhJwOoyOAOoRMAACEnAAAhJwAAIScAACEnA\n4B/X4UAgAgYvpV4yBiEIGLyMlzy4Qy5E4U4cMB9ydfNultA3IzAYcv6Xq/ED9YL2CRi8zBoGNE7A\nAAhJwOAfYy8IRMAACEnAAAhJwAAIScAACEnAoHPujEWvBAy65e6O9E3AoFs5uy0WPXMvROjTdNRV\nPlYyemIEBn2a3dRxer9H6IOAQbfc3ZG+CRgAIQkYdM7wi14JGAAhCRgAIQkYACEJGAAhCRgAIQkY\nACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYACEJGAAhCRgA\nIQkYACEJGAAhCRgAIQkYACEJGAAhCRgAIQkYX6WU7t6EbWxwVbG2dgi4wWz1c/cG9CCllHN+/83y\n8ftXAThIwA5ZfomnWwD1fBg6sNXyCGz6JXMawAIH5E2MwCoan4vTvHl2ApxFwNaaDZ7+TJFWAVQl\nYGttCtLHSUUATiRgJxvTlXP+eA4MgLMYKAAQkhFYde/TiYEuEWt8LvTbMLf9PRxigB539xaNP4FH\noQ8R9xKwihYWzbf/pGx/xf/0r/37IaDlPby85Y2Iu3tH7T+Bh+CHiBa4lVRF48mwj19KKTX+F2xh\n40Nofw+H1v7uDfEEDn2IaIGA3WO20IPT2cNV2b212cNrmEI8x6arxFp7Ybj1ErdbhN7DnbF7a7OH\nVxKwc6x/wjV4zqO17fko9B7uid1bmz28noBdanxqhliB1r73PRnlIrz2t3D4spGewFVFeQK3Q+oB\nCMkiDgBCEjAAQhIwAEISMABCEjAAQhIwAEISMABCEjAAQhIwAEISMABCEjAAQhIwAEISMABCEjAA\nQhIwAEISMABCEjAAQhIwAEISMABCEjAAQhIwAEISMABCEjAAQhIwAEISMABCEjAAQhIwAEISMABC\nEjAAQhIwAEISMABCEjAAQhIwAEISMABC+j+nfz6FyJm2nAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 27, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "plot(real(R),imag(R),'*');\n", "hold on;\n", "plot(cos(2*pi*x0), sin(2*pi*x0),'ro');\n", "plot( exp(1i*linspace(0,2*pi,200)), '--' );\n", "hold off;\n", "legend('Roots','Support of x'); \n", "axis equal; axis([-1 1 -1 1]*1.5);"]}, {"source": ["We isolate the roots $R_0$ which are on the unit circle, those we are interested in."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 28, "cell_type": "code", "language": "python", "metadata": {}, "input": ["tol = 1e-2;\n", "R0 = R(abs(1-abs(R)) < tol);"]}, {"source": ["Note that roots located on the unit circle (those we are interested in) are\n", "actually double roots."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 29, "cell_type": "code", "language": "python", "metadata": {}, "input": ["[~,I]=sort(angle(R0));\n", "R0 = R0(I); R0 = R0(1:2:end);"]}, {"source": ["compute argument to obtain the position $x^\\star$ of the Diracs locations of\n", "the solution to the primal problem."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 30, "cell_type": "code", "language": "python", "metadata": {}, "input": ["x = angle(R0)/(2*pi);\n", "x = sort(mod(x,1));"]}, {"source": ["The first order condition of the primal problem $\\Pp_\\la(y)$ implies that the optimal\n", "coefficients $a^\\star$ satisfies\n", "$$ \\Phi_{x^\\star}^*( \\Phi_{x^\\star} a^\\star - y ) + \\la s^\\star = 0, $$\n", "where $ s^\\star = \\text{sign}(a^\\star) $ can be computed explicitely as\n", "$$ s^\\star = \\text{sign}(\\Phi_x^* p) = \\text{sign}(\\eta_\\la(x^\\star)). $$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 31, "cell_type": "code", "language": "python", "metadata": {}, "input": ["Phix = Fourier(fc,x);\n", "s = sign(real(Phix'*p));"]}, {"source": ["One thus obtains, since $\\Phi_{x^\\star}$ has full rank,\n", "$$\n", "      a^\\star = \\Phi_{x^\\star}^+ y - \\la (\\Phi_{x^\\star}^* \\Phi_{x^\\star})^{-1} s^\\star.\n", "$$"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 32, "cell_type": "code", "language": "python", "metadata": {}, "input": ["a = real(Phix\\y - lambda*pinv(Phix'*Phix)*s );"]}, {"source": ["Display the retrieved Dirac locations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAOcklEQVR4nO3d0ZKq\nPAKFUTJ13v+VMxdMMwiI2tiQHda6+KtLndMZFT6BmC611gEA0vzn6gEAwG8IGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiPTv6gH8rVLK1UMA\niFRrvXoIL3QesKGN16CUYhjtjMEwWhuDYbQ2hiHk079TiABEEjAAIgkYAJGaONn6dxo5mww0KOIy\nz5nme8uInWf/kzgAnml/H32axJw7hQhAJAEDIJKAARDJNTCAoxYXkF5eWpsev/nIZxMo3pxYETH/\n4isEDOCQT6c/zAOzGZtn+blJlt4nYABftnMMtLir1jrdMoZwumXexekBz25f3HITHQbsPofP7Gt5\nq255bHe2eSy1eIH+dLr55pHZ+jcuurXzyL51FbC7vXjsWLwZmvpY0/LYbu6dF2L9mMv3PJcP4Cpd\nzUKstdoR8EzLG3nLY+MXdnZEi3OAxz++jP/CPfd+XR2Bbfp0dhAd0APO9OyE3s7jXz540bl933rD\nx204HZ67eDnDhzt4dtmgBe9caOEEKfuHc8a5+C0RT07/R2DQvvb3FJzv/UO62+rqGhhMFtt8U7uA\nlsdGO257Zet9HQbM681o2vgbfEu0PDZI0WHAYK7lQrQ8Nmifa2AAh80n5ry9EOLPwy/7HBMxU2OH\ngAEcs5hWWsrLht127YzvEjCAb3ujYav/xdN1oTa/Z7b+vtDm4ze/c9ZNNQUMLrLYiSSfyenH5p59\n8dJ8Y++/Xk5sWK3Yu+jNIlfzNX+nSi3atnn7/iqLWQSMzjV6ln+94/j8Mzvf985LsH7M5xnYTMiz\nnLy/Ksezh6WH6hkBA/i2zz+L7BwVvZmfZ/9CN8dbawIGcEytH81C3PoHlleqnl3TGmZBmp8PfLa+\n4v6/nE7AAA77JFqLwKxv3Lxl80z45j/16b+cyxeZ4QrrPUgv+xQ4jYDBRebFUi/4nIDBZaZLE9cO\nA0K5BkbnWj7d3/LYbqKb6Qz3JGD069jEMLrnA0Q6pxDp1Hp5OqAvAsZtaBj0RcAAiCRgcJlSiuNC\n+DUB4zZcsYe+CBidWuRKvaA7Aka/av1ft9QLeiRgAEQSMAAiCRgAkfKWktr8o22Le60Qw/81/Gbw\nRoUjwgI2/fnRxc+LW9Z3AdAZpxABiBR2BPbS+hTi4s8lODID2BT3x2X6Cdizs4uKBfCOnY/+bXIK\nkd41vB1aCxGOCDsCq7WuZyGOx1ubdwHQq7CADVtxcrYQ4IacQgQgkoABEEnAAIgkYABEypvEAZ9p\neGqPaUdwhCMwACIJGACRBAyASAIGQCQBo3cNLzZoLUQ4QsAAiCRgAEQSMAAiCRgAkQQMgEgCBkAk\nayHSu4bXG7QWIhzhCAyASAIGQCQBAyCSgAEQScDoXcOLDVoLEY4QMAAiCRgAkfK+B1Z+Trmsv0Oz\ncxcAnQkLWCllitP85+GnXuMti7sA6E9YwF6aZwyAjnUVsM2Ds/I4y0vbADaVtDmxXQVsk2LdXcNv\nAG9OmrK+KNM4sxABiBR2BFZrXU81HE8Ybt4FQK/CAjZsxWm6RbcA7sMpRAAiCRi9a/hatLUQ4QgB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkfK+yAyfafjr7b56D0c4AgMgkoABEEnAAIgkYABEEjB61/Bi\ng9ZChCMEDIBIAgZAJAEDIJKAARBJwACIJGAARLIWIr1reL1BayHCEY7AAIgkYABEEjAAIuUFrPzY\necCZ4wHgEmEBK6XUH5uhUi+WGn5LWAsRjggL2L4xb1ePAoAz9D+NfnFMpnAAm+LOYPUTsPGpn/47\nhUqxAN4x31tGxKyfgE1PvROJAHcQFrD53A3FArizsIANW6cEF7eIGcAd5AUMPtPwBxoftuCIrqbR\nA3AfAgZAJAEDIJKAARBJwOhdw9/HtBYiHCFgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkrUQ6V3D6w1a\nCxGOcAQGQCQBAyCSgAEQScAAiCRg9K7hxQathQhHCBgAkQQMgEgCBkAkAQMgkoABEEnAAIiUtxZi\n+Zl2vF5Hbucu7qvhN4M3KhwRFrBSyrTNz3+ejLds3gVAT7o6hShaAPcRdgT2jsXhV3lc6UDkADaV\ntHVhugrY+OwvEqVYAO/Y+ejfpq5OIQ5yxVrD26G1EOGIsCOwWut6quF4znC83UREgJsIC9iwVabx\nFsUCuJXeTiECcBMCBkAkAQMgkoABEClvEgd8puHZPWYewRGOwACIJGAARBIwACIJGACRBIzeNbzY\noLUQ4QgBAyCSgAEQScAAiCRgAEQSMAAiCRgAkayFSO8aXm/QWohwhIABnGL+nT+fXb7BKUSAv7f4\nxrovsH+DgAEQScAA/tjm8ZaDsMMEjN41vJuwFuJdbF7xchnsMJM4nljsVrzVABrjCGzL+kOxj8nA\nEYsPwT4Tf0NXR2DlJzO+XgM0Z9wvlaJe39JPwEopU7fmPwPQJacQAYjUzxHYM+VXl68Wh29lcBks\nVf3te+A0rQ2vjm94/kb7b8gg/Qfs9+cSpzeZ05HhGn/5mhteKbVW5+H/TsTTGlHZ/gP2e7ZegIb1\nE7DxM+P087WDAeCv9ROwQbcA7sQsRAAiCRgAkQQMeFTr4IQ8CQQMgEgCBkAkAQMgkoABEEnAgEdl\nXPszYCUhbk7AAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDHlkLkRACBkAkAQMgkoABEEnAAIgkYMAj\nayESQsAAiCRgAEQSMAAiCRgAkf5dPYCPTdeW1ysF7NwFQGfCAlZKmeI0/3ky3rJ5FwA96eoUomjB\nF1gLkRBhR2DvWBx+Lb7OYrME2BT35b+mA/Zpe8bHLx6mWADv2Pno36amA/aL9sgVwE00HbC1Wut6\nquF4zrA8rn+jZAB9CwvYsFWm6pozfFEpQ62m8tK+rmYhAnAfAgZAJAEDIJKAARBJwACIJGAARBIw\n4JHvpRBCwACIJGAARBIwACIJGACRBAx49LguNjRLwACIJGAARBIwACIJGACRBAyASAIGQCQBAx5Z\nC5EQAgZAJAEDIJKAARBJwACIJGDAI2shEiIvYOXHzgPOHA8Al/h39QA+U0qZZvfOf54/4PRBAXCB\nvCOwHZtJA6BLYUdgv7A4JlM4gE1xZ7CaDthH7SmzK8/zQzHFAnjHfG8ZEbOmA/ZRe/avjQHQmaYD\ntlZrnT4XKBb8CWshEiIsYMPWdrW4xYYHcAddzUIE4D4EDIBIAgZAJAEDHlkLkRACBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkY8MhaiIQQMAAiCRgAkQQMgEgCBkAkAQMeWQuREAIGQCQBAyCSgAEQScAA\niCRgAEQSMAAiCRjwyFqIhBAwACL9u3oAH5u+X7n+hLhzFwCdCQtYKWWK0/zn4ade4y2LuwDoT1jA\nXppnDICOdRWwzYOzxZJu2gYvlDLU6jTGDcUtgNl0wL7SHhshwDvWF2Ua13TAtAeAZ5oO2FqtdT3V\ncDzXsXkXAL0KC9iwFafpFt0CuA9fZAYgkoABEEnAgEfWQiSEgAEQScAAiCRgAEQSMAAiCRjwqJQh\nZCUhbk7AAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDHlkLkRACBkAkAQMgkoABEEnAAIgkYMAjayES\nQsAAiCRgAEQSMAAiCRgAkf5dPYCPTdeWN1cKGO+1iABA98ICVkqZ4jT/eXHL+i4AOuMUIvDIWoiE\nCDsCe2l9CnHxdRabJbwwbjK2lPuJ+/Jf0wH7qD3Pzi4qFrxrvsWVomF3s/PRv01NB0x74DzrHZaG\n0bamA7ZWa13PQhyPtzbvAqBXYQEbtuLkbCHADZmFCEAkAQOGYdiaduiUBm3LO4UI/JXHSWjyReMc\ngQEQScAAiCRgAEQSMGCDL6XQPgEDIJKAARBJwACIJGAARBIwYEPEX9Pg5gQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjBgg7UQaZ+AARBJwACIJGAARBIwACIJGLDBWoi0r8OA2fDgoHEjsinRuH9XD+Cb\nbG9w0GIjKqWYT0+zujoCq7Xa2OC7fC6kWV0dgW1abH4KB7Ap7sNKasDez5JiAbxjvreMiFlqwGQJ\nzmFbo1ldXQMDDpIrgqQege2wBcIRtiBSOAIDIJKAARBJwACIJGAARBIwACIJGACRBAyASAJ2hkYW\nZWlhGC2MYTCMxsYwGEZjY0ghYABEEjAAIgkYAJE6/3OrziYD/E77deg8YAD0yilEACIJGACRBAyA\nSAIGQKQO/yLzZJqCeP5Elf1fXcoZc2d2xnDmM/Psd5386qS8Ipv3njOMFl6RxbThq56K+b332UbW\nztkujug2YPOn/uSXYedXnzat/+X//fGWv35m9odxzhheDuOcF+XNV+SqYcx31he+ImduLC/fFddu\nIxfuwYacLyA5hXiqWmsLn2iMYdLOZ8xSyuV7jXEMLTwhjQzjthrZU73U7REYL12+j7h8f92U045+\n9gdw7RiactopxHeGwSYBu6NGtsxpl33VAMZfPf33wifk8teiKZcX9NrTd6Na67SdatgzTiHe1OU7\niAt/+6T+GC59Qhp5NmjKGE6fbPb1fK6gnVmI66kcl0xtGn/vtRO95tfGzxnA/jCmey98RTbvOm0Y\nrb0iZx7xeEVeDqzxQLQ+PgDY5BQiAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyA\nSAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkY\nAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASP8FafkIp4TA38MAAAAASUVORK5C\nYII=\n", "output_type": "display_data"}], "prompt_number": 33, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "stem(x0, a0, 'k.--', 'MarkerSize', ms, 'LineWidth', 1);\n", "stem(x, a, 'r.--', 'MarkerSize', ms, 'LineWidth', 1);\n", "axis([0 1 -1.1 1.1]); box on;\n", "legend('Original', 'Recovered');"]}, {"source": ["__Exercise 1__\n", "\n", "You can see that the dual certificate $\\abs{\\eta_\\la}$ saturate to\n", "$+1$ or $-1$ in region that are far away from the initial support\n", "$x_0$. This means either that the noise is too large for the method to\n", "successfully estimate robustly this support, or that $\\la$ was not\n", "chosen large enough. Explore different value of noise level $\\norm{w}$ and $\\la$\n", "to determine empirically the signal/noise/$\\la$ regime where the support is\n", "sucessfully estimated."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 34, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo1()"]}, {"collapsed": false, "outputs": [], "prompt_number": 35, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Certificates and Pre-Certificates\n", "---------------------------------\n", "The minimal norm certificate $\\eta_0$ is important because it describes the\n", "evolution of $\\mu_\\la$ of $\\Pp_\\la(y)$ as a function of $\\la$ when $\\la$ is small enough. In\n", "particular, as shown in [DuvalPeyre13](#biblio),\n", "if the saturation set\n", "$$\n", "      S_0 = \\enscond{t \\in \\mathbb{T}}{ \\abs{\\eta_0(t)}=1 }$\n", "$$\n", "is equal to $x_0$ and if $\\eta_0''(t) \\neq 0$ for $t \\in S_0$,\n", "then for small noise and\n", "small $\\la$, $\\mu_\\la$ have the same number of Diracs as $\\mu_0$, and\n", "both Diracs' locations and amplitude are close to those of $\\mu_0$.\n", "\n", "\n", "A major difficulty is that in general $\\eta_0$ is the solution of a\n", "convex progam, and is hence difficult to compute and to analyze.\n", "\n", "\n", "It is however possible to compute a reasonnable guess by computing the minimal\n", "norm vector which interpolates the Diracs with vanishing derivatives. It\n", "corresponds to computing $\\eta_V = \\Phi^* p_V$ where\n", "$$\n", "      p_V = \\uargmin{p} \\norm{p}^2\n", "      \\quad\\text{s.t.}\\quad\n", "          \\Ga_x^* p = b_0\n", "  $$\n", "where we define\n", "$$\n", "  \\Ga_x = [\\Phi_x, \\Phi_x'] \\in \\RR^{ N \\times (2n) }$\n", "  \\qandq\n", "  b_0 = (\\text{sign}(a_0),0  )^* \\in \\RR^{2n}.\n", "$$\n", "\n", "\n", "It is easy to see that the solution of this problem can be computed in\n", "closed form by solving a linear system\n", "$$ p_V = (\\Ga_x)^{+,*} b_0.  $$\n", "\n", "\n", "Compute the $\\Gamma_x$ matrix."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 36, "cell_type": "code", "language": "python", "metadata": {}, "input": ["d = 1; % number of derivatives\n", "w = ones(N,1);\n", "Gamma = [];\n", "for i=0:d\n", "    Gamma = [Gamma, diag(w) * Fourier(fc,x0)];\n", "    % derivate the filter\n", "    w = w .* 2i*pi .* (-fc:fc)';\n", "end"]}, {"source": ["Compute $\\eta_V$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 37, "cell_type": "code", "language": "python", "metadata": {}, "input": ["pV = pinv(Gamma') * [sign(a0); zeros(n,1)];\n", "etaV = PhiS(fc, u, pV);"]}, {"source": ["Display the pre-certificate. In this case, it is a certificate and hence\n", "$\\eta_V = \\eta_0$."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARhklEQVR4nO3d27ba\nxraGUdjN7//K7AscFuYodBz/qN5bLpIsr8wSVOlTgSyfL5fLCQDS/N/RAwCAOQQMgEgCBkAkAQMg\nkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIG\nQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABE+nP0ALZ1Pp+P\nHgJApMvlcvQQvmgesFPCewBQTcTVv48QAYgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACI\nJGAARBIwACI1fBbi+Xz2/ENO/z7MrdqUqDw2SNFqB3Y+nyMeQMkOHmZCqYlReWwQpFXALpeLi1ne\nqdyJymODshp+hPjg+dTwHDm/ptmvOZ/Pp9Pzpczj/+u4Mb8eW6nX0K8Z8NfEXUg1/Lro/jsw34eN\n6W4ZXv/u7xyoMBcqjw1uIk6erT5ChNOLQlz/5vzv/3SMN2N7+J+ASQSMnh6uHW9fj1boxNPYDhoH\nhGsYsPrbXrZz7dN1Ctxu6rlNiWMb9mFsdeIKQRoGjME9770OGsg/7ut1U2RsEErA6GPKDqbsXqfs\nwKAsAaOVKVua/VPxcvsFLCRgNFF87zK9rMBEAsaI9tyE/fpTipcY6hAw+ni5gzn24QKfPzyMe/AB\nlCJgdDAjBDVvmvApIkwnYDRR8NQ/+96NalmFmgSMcdXchAETCRjxahZo9var4FYSahIwhlZ2E1Zw\nSFCNgNFBtV2L37kMOxAwsn3dqXx93uCBm7B3Y1M+mELAYGW2X7APASPe8lTU/Cas2nigGgGD9dl+\nwQ4EjGAr7lHW2oStNSQJhK8EjOb2fN7gr99+eRYiLCFgZFtxp1Lwm7BSg4FqBAz+Z0kO3XwIOxMw\nUm23O7HvgQgCBv+Y90HiFtsvmzn4TMAI5hQPIxMwePTrJmzTb798ngnvCBiRpp/Wvz4L8c3/a+pP\nWVKveWMDrgQMPvncMNsjOJCAwWtfd0e3em23j7JDgw8EjFQ7nNxvHyQ+77R2qBfwmYDBJ7c+3TJ2\n3zP1ggMJGHl++uZp+fMGL5d/Mvb8L2ebODbftMFLf44eAGSw2YJq7MAAiCRgRBpnPzTOkcKvBAyA\nSAJGmDHvaBjzqOEzAQMgkoDRXOXnDVYeG9QnYORx2gdOAgb1CTa8JGAARBIwkox8M97Ixw4vCRjN\nLX8W4nYqjw3qEzDC+EIIuBIwACIJGASw74RnAkYMXxgB9wQMgEgCBjHsQeGegNFc5ecNVh4b1Cdg\nJHHCB24EDIBIAkYGX//YfcIDAQMgkoDRXOXnDVYeG9QnYJBE8uBGwIjhSyDgnoABEEnAAIgkYATw\nxc+VD1HhnoABEEnAaK7y8wYrjw3qEzAyONUDDwQMwvhGEK4EDIBIAgZAJAGjuoWfmFV+3mDlsUF9\nAgZJ3MwCNwIGQCQBI4BtB/BMwACIJGCQx80fcBIwinOmBt4RMJqr/LzBymOD+gQMgEgCRnV2KQ+8\nIHAlYABEEjAAIgkYzVV+3mDlsUF9AkZdTu8feHFAwACIJGAARBIwSnPLOPCOgAEQScAgj40pnASM\n9io/b7Dy2KA+AaMot4kDnwkYAJEEDFLZpDI4AaMu3xABHwgYzVV+3mDlsR3rfP77F3wgYEAt993S\nMD4QMIjU9fPVa7Eul79/nTSM9wSMipyzxnSr142G8YGAARBJwIASnrdfVzZhvCNgAEQSMIpa6yaF\nys8bXD62ZvuSd6+HTRgvCRhwPHFiBgEDSvi8HS28keYwAkY5LsZ5x9zgnoBBqjabElliHgGjucrP\nG6w8toLaBJu1CBhwvOlxEn1uBIyKXGuPQ5CYTcCAGK5suCdgkK3BDubXLDU4ZFYhYNTi3ARMJGDA\nYWZcr/gUkRsBo7nez0KEkQkYkMdHzZwEjIJsS6Zr8Fo1OASO8ufoAbCJh+tT5wgKmr2LulzswDid\n7MBauq3ty8UfpAS0JWDd3P5c9tuuK6hhWwyy8vMGK4+tPi8eAtbKrV4PghrGDLnvrA+3WULAgDzK\nx0nAOnm3/bqyCaOUAafi+fy/v1iFgLXS4A9ljxgkRaSU4DlaKSMvTsCamL4erBw4xPXWKvcGr0jA\n+piyd7G/6Sf3Pc0d+U9efravYasQsH8+mA79ePrXMSce42yVnzdYeWz1Rbx4H76Z1rDlRg/Y/e+a\nuk0yU+oQXnZGE9HgyoYO2PPFUWLDPt98+Mya4XArrq/KS3Xi2qx8CMWNG7DPv+f31H1W9T46ONzE\nevkgcYlBAzblt0wF+XXAZQ+w7MDqG/MMaMIMbtCAfZVyWbRkhPWPbhWVnzdYeWxbG6E9E48x5WxT\n0IgB+/VLI6Cygqf+gkNqacSAnX75LVP1J+K8Eus3bOqnJWY9zjNcwH4KUvFZVT+u03U6Fj5b970u\nuEhnH6BV8KvhAnb6fcY3nlWND20oBU/isIOxAtbvfO3MBQXNWJgp31mUMlbATnPvOC84q5YPSfw4\nyupzr84KrTOSEQwUMBPrpTovy0ZBrfy8wcpjY3+mw68GCthsZWdV2YHBnqothIXjqXNNWd8oAVv+\ne79Kzaq1BlNt5bNEqSn6TsQgZ+t9dAWNErCFep/oD191hw8Aiuh9qlndQAFbPjNKnWdNdLhXZHmu\nsjCLHEt9QwSs2Wxodjhbq/y8wcpj207Xa68h38yDDRGwVXRddV2Pi6F0msadjmVrowRsrTlR5CKr\n3xTvd0Q78wJWsOK7UORUU1z/gJkHU3iV2MHW0+zAaWwFHaJ/wFZU4SLXOoH2KpxqIgwRsHVnQ7+E\nHLta+r2ejKnx87HKah+whlcyrs7g2YHrQmmO0j5gK0/rxpuVrouw8vMGVx9b/Tex8LtRS4EX6vgR\nfNU+YJss6PqnCRhTs7V51OGkvIztA9bNZo9s3+Q/G/HTYbmOczggYgL2s44zlQ6Kz8yUi/pfbXdc\nxd/QCgRspv1X4w4/sespZhDevpZn/JanmrUIGEfapcp1l+O6Y2t5+o7glT+KgM1x1Hzd9OdahPRQ\n+IrlZy1PNSsSsPn2XCed1iTDSjktTtdvYWYdkYDxj/2nb7+TGgdqOZ2yorInAZtp/3XScmVCun4L\nM+iIBGyRZldGQROXd5rNycr2ean3XJVxk+dc+UE7y218B9r1pdvhPd/tB11/1m6zeM/jGkTZl3Tn\nObzPz2r+g+rXoX/ANj3A83mP66Nrhfd5o7r+rEGs+JKuu3Zazqt9lv/poCPa+uS5Ch8hZig/kaCK\nfotlz0xmEbClEt91YCEPyKhAwKrbedbufPXa72L5WJVfz8pj4yruPRKwReLebxjEDld+ey7/rX9W\n6PZOwFaw9Xu/fyZDZzMwFAErrXFIGh9aD6vfgrinrvuVTX9u4udJArZU4rv+Wb8jAj7IvZoUsHX4\nQ+2oI/d8tK5m+xWngmcCxmvNvgOHHjZamKGLUcBWsNF77zqaNkLPjy81W5jRhyNgq4meB/Cg8p9k\nPUWnZN60PKglBKy0rn8ea/i5kRgtP3Bb/aByuyhg61h9BjjFM0/uyShCs4WZfjgCtqb02fCg2eFw\nCLNoXS5Q7glYXb1nau+jo4KN5liFqbviZUGFw5lNwFaz4jyocNEaPa2BryqcZxYSsJU1mBPQyYpL\nssjqdnF5I2BFFZmjW6zYImeB3pa/yCs+C7HIZG5mlXWU/tYI2JpW+qPcV/iPAKcNTtBFzvjLh9Hj\nPCNg6+sxM4B7Bdd1wSHtTMAqanOVt/9/GZ450T+4viANlqGArWzhnCi40goOiRQmz3Ya5Gc5AduE\ndfuOV2YHK30X2+StWvFEX7AZs9+lgscyg4Ctb/bMaLOvh2a61Px06nUsAraVHrNETemhx3p8cF2e\nMw6tzboWsFraTKwPRjhGbtq83T0OpFnIBWwTM66MKk+symODr/rdWnUzexPWg4BBT8Oe1N7xgvT7\nll3AtvLTlVHlidXsIcUwQ821eTXyJkzAgBeWPwux1Cl14b3BDVS+Sp5NwDY08cooYmK1WcaMrOs0\n/nqqiTjJzCBg7KrfEqK9iEn7oWFds30SsK31uDIqPjyeecuezXhNsk79L882t39sOSUEbHNfr4xa\nTqxnWecCVlFwbveeh7ezze2v+3/fj4Dt592VUYq4AbNEm2chLhR36n8Y8OWSdwjT/Tl6AEO4XP6e\n/Z/PCY3n1rOhDpaCbitxityCj7PQ7MB28nwdlHVlFDRU+Gx6mUz74uzAdjXsesi9mI12Ph825dLf\n8fTxD8IOjB9Y1aQb+bkV/QgYOxl290mcoW4PjiZgTGU904NNWBsCxm9mLHtnikTLn4UYesVj+xVE\nwKAnp+AP3m3C1CuLgPEDC5s2nhvmo4I4AsbPZqxz5aOghwcv3f9LIggY23JVO6CgN/0+V1nPFuAk\nYPzKCh/EOM9CvHbLxE4kYGzOqeFAw2SIEQkYc0w8LTp7AtsRMH5mRwVUIGDMNHF3pXZj8r6zAwFj\nKz4/BDYlYMzh+jqCt4neBIz5vu6xnEBzLX8WImxNwJjp8/nN54fD8tazGwFjkZdnK09EBXYgYMzn\nz1UCDiRgrOD5kd62X3W4wqArAWORW6iuZ0n1amPJsxBNAPYhYCz10DAnL2Aff44eAB2IFrA/OzAA\nIgkYdLbz5tgNI+xJwACIJGAARBIw4AXPQqQ+AYP+9vxqSvjYjYABEEnAAIgkYABEEjDghRnPQvSb\nwNiZgEFz7qqgKwEDIJKAARBJwIDV+LiSPQkYDMEdFvQjYABEEjDgBc9CpD4BA1bgI0r2J2AARBIw\n6M/HgbQkYABEEjDghRnPQrTPY2cCBkAkAYNRuFGQZgQMgEgCBixlb8chBAyASAIGQCQBgyH8eo+7\nZyFSn4ABK9A79idgAEQSMAAiCRgMxP3udCJgwAvTn4UoihxFwACIJGAARBIwGIU73WlGwIClpJFD\nCBgAkQQMgEgCBmOZeNf7xGchuoeeAwkYAJEEDIBIAgZAJAGDgWxxv7t76DmKgAEvTH8WIhxFwACI\nJGDATDZpHEvAYDjCQw8CBkAkAQMgkoAB87mHngMJGIxlYnImPgsRDiRgAEQSMGAOtzJyOAGDEckP\nDQgYAJEEDHjBsxCpT8CAmdyoyLEEDIYjPPQgYABEEjDgZ74gowIBg0GJEOkEDIBIAga88PVZiO4E\n4XACBiOSHxoQMAAiCRjwG3d/UISAwbikiGgCBrzgWYjUJ2DAz9wDQgUCBoMSIdIJGPADnyxSh4AB\nEEnAYGh2VOQSMOA3vjyjCAEDXvj6LEQ4nIDBuESKaAIGTOULM0oRMBidLBFKwIAf+NSROgQMeMGz\nEKlPwGBo03dUikY1AgZAJAED7K6IJGDAVO7goBQBg9FNyZItGgUJGHA6SRSBBAx44flZiD4/pBoB\nAyCSgAFfdlc+XaQmAQP++hAqnx9SkIABn9h+UZaAAafTf3usW648C5H6GgbMwoOFrovofP67mHx+\nSE1/jh7AmqQLlrhcTuezzwyJ0WoHdrlcnn/zCjDXdTUJGkW12oG99Lwte46cX+PX+DX//Zrz6XT5\nL13Xf/zfL646Zr9mnV8T9yHWOXTL8vBCP7wHt3+8/3vgq5enMItoQBEnz9QdWP1XFnqw1iir1Xdg\nwEJyRZDUHdgHViAsYQWRwg4MgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAED\nIFLDZyE+iPsTbgCYIuBPfAGAZz5CBCCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoAB\nEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQB\nAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBI\nAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACR/h+VRjW5H/EZBAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 38, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf; hold on;\n", "stem(x0, sign(a0), 'k.--', 'MarkerSize', ms, 'LineWidth', lw);\n", "plot([0 1],  [1 1], 'k--', 'LineWidth', lw); \n", "plot([0 1], -[1 1], 'k--', 'LineWidth', lw);\n", "plot(u, etaV, 'b', 'LineWidth', lw);\n", "axis([0 1 -1.1 1.1]);\n", "set(gca, 'XTick', [], 'YTick', [-1 1]);\n", "box on;"]}, {"source": ["__Exercise 2__\n", "\n", "Study the evolution of the pre-certificate as the separation between the spikes diminishes.\n", "When is it the case that $\\eta_0=\\eta_V$? When is it the case that $\\mu_0$ is the solution of $\\Pp_0(\\Phi \\mu_0)$ ?"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAaAUlEQVR4nO3d25ai\nyhIFUDhj//8vcx7osihQTG5JRjDn2A/dvVG5RORKULEfhqEDgGj+d/cKAMAeAgyAkAQYACEJMABC\nEmAAhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhPTf3SvAyfq+H//gdwZ4FJX/QM7A\nElpv4P7Hp/+7/OunhaEp+yr/7b+v9AjtcAaWzTAMfd9/6uTp/1oupmOJa1/ljzU//nXl4bTJGVhO\nO6LobZ6Ng8K+J4T69hXqeL41S7Xdz0Y1zsCyefXhsvf2zS7XJ7bQiN2V//bMTNmHIMBSmbVueftN\nT7P0LeHsrnxCE2CpvK74jcq7ejkJFWMEsrvyl1R+IAIsoen70jPTPl/5NAdEtKnyx7J/2xFEYeQC\nICSfQqSUb8bwTCq/WS4hUsRlRp5J5bfMGRgAIZlcUGrlHXJITOU3S4ABEJJLiACEJMAACEmAARCS\nAAMgJAEGQEgCDICQkt+Jo+/7rhu6rvA2MOVLHnnUvldpzelbsfsJNz1wz6tE/KqJWx+14VO9xRgE\n2q/85AHWdUNXfBj6fs8B2/qofa/SmtO/2bl7t2x64I5XiZsECcosuk/1FmIQCFH56S8hBjgGQTXf\ngNCuCOkQQPoAA7iBiKpAgPEUBhRIRoCxkzyAfVx+P4sAY6dwTRhuhYlOyV1NgAEQUvoAMwUCyCl9\ngDmLz8z7cPBk+QOMxMxO4MkEGK1wOgVsIsAATlYyGzNjO06A0QQXA4GtBBgAIQkwgNpccjiFAAM4\nn4iqIH2AeZ80sypvgxuHoFEJf9Cy7/v2fyzuUrNfovu0NwoXI4Snlf2R6lX5aaQ6A+v7PsSviF5q\nuQfe7pPCxRpn5OmUfdd1W6rXvsokVYANw2AyxdM8sOzPDSGRFlfCS4gzy+pcdvvPMsNr4c/LlDzP\n2jLLxa55rdef++nCH5b5s1jJay3fGTpr/2xfZnh79ez4a72qoovp6x54O3CHWmZewMXPU1T5bW/7\nhcsESvSE182nw9mmNwb6fs8lqa2P2vcqW57/Nez2P436aXz/939/Ft7wZsDpWzGu9dX7f8dqX328\nzjIr9dxvif2t3m78c+Hm/gzOmx+4SUnZ7K75OkKUUKpLiHT/QmhY/MubBbuum8xbH3cZinP1/b//\nKvip1T8vtuXuTdeuZZwTmPAE2B9bK6/5Sh0vj9y9FgWkJ+VmJb21eCZztTANwlsJ3wM7ciaxvRMa\nrf5h6MZzrLerN7l2MUz/JZygq32Fe0+gXxU1noRVWJdp9Ra+9PSS3eRdhgtXkqs5A0tu0/tDkNWn\nN5ySnfpXu4rbCAH2CMuaztG3ObYih58zm8tf6NMAPb70vuE7x6D/2oocm1NCgGXztXafU9zUsayo\nu2rsU4bV/7xfyWtdsT5Pm9IJsPymXf2QCynktvWLK+sPyVH8r9jOsTmFBNgjHLm6suQcjpnpoHnv\nAPoq9emM7YEesuEC7I8cH6P/6eHZ/T5+/5va8a37R03xaMqsWj/cVGK+8NuKnT22zV7e5IGNmfBj\n9OwW5a4T7PCoS0xP2MYVzX6953TOwP7Y8T2wB7qiNx7Sb3eZvg907q7+9GzXHdBnNt26x7aPAGOz\ndkaQx/Yt5zqxpHO8DRGFAOOfFt575zrVvqfFLZZH9gnRKMBSqVOyTTWGEfmrq7+nVe1bGU0VHi0Q\nYKn8/TRz0SjiJvTEdezGp8vfGDq2Ni15SFsLsD+iV3CtH7Oo8SpAieij1hEC7JAnlw6BXPRF4/r1\nf93k6cRnLn8qc8GDBNgvxVRIbFNCnXA1AXZIvswrHHTybXhW9b+nRX2f+jH9URZg2YgWaqpWb+nH\nYnYQYGkV3uRwuljEX79sZ02adcuPD9QJth138lx/bKZyesJcVoDRhN0DxxO6FHhLgLGH2KA+VbeU\n6ZRxBwH26+GlsMnp+8rYxO0UYTgC7JfyLWdfhbA+zzg4Cyl5+NMmhU/b3tsJsEPUK3DQwWFkfTaZ\ne4wSYGm5FyLd59Et1mH/NAqfey9EwhFgR7XTBWdNtXJP2TjRevG30xqbqP9ABBixGW5oSlOx3dTK\nXEGA5ZG+WAGmBNgvc/mIxPYnyT4lWOdA1y+nI68Y6PBdRIDloZqBRxFgvxLM5aebsONeiFCokao5\n/V6IKSXeUAF2VFPFcXxlEqQ4Ly18SrCpBiEZAQYQ1cNnnALsqKYKqKmVaZMTAr5SJFEIsBtojxtt\nynhHaoWfspxpcz1zT2oF2K826w8udWnZ5x49uZ0Au0Gtb7S4F+JRkfdN4FUvd8W9EG856GbP+wiw\nX5FHK8hMby7JvE6A8ZbeeIh9waA8aIQAO0oz0yZnLfU1u8+zDlMCjGfJ2snE1WrsNbpaUwLsqFaL\n7yqPCoDIGxt41VsQ4tA/bfBZEmBpXXQvxOg9E339IzoSBuuPdS/EEolrXoDdoPHGSVzuW9kVn9gz\nuTU+Rr0IsCSiFBwNUjwEJcBucMXs1YyYNh2vzFg/MklNAiyJapNos3V4hgCtLsCSqDZnNDlNZusB\nNYNpgaMw6nPfCm/7B42GjfOOrcvve0jh03bnPfPKSp77QoUvetZDdqx8yBbZWPmbduOefXigYC7q\nlyte9+BmdlcehR3LB0iHAKt4RN9v28C+3/xzGzsmsBft8hOfeflUrz05DoxXbEKFnbl15a/b2Ktt\nqvxNu3HHPjlSmdf1y+7XfbtvD5bK1ZW5a/kA6eASYhI1Lyk0X9UfbV3zuFvKcS7TtU+A5XHuaLvS\nvRr74XZU2r6aUWlnyfo7rgIsj2pl57wkn0BjVrMq94U27ATYLS4aLBQ0nEUahSDA8rj4t+F1GHTd\nM3ohyiYKsNquq4woNXcjF8reUjkEJcB4w4g2Engz9XeIUjxRvnoWYDyI0fAWdvu58uXQbgIM1hkt\nnqswKiTKXQQYQH6F58GxwliAUeQ5P1/7WIVH2PVAvdAOAQa06IEx8cBNPkiA1aZGgYMuPQ8OdJIt\nwPhI1nZd9/M7FMlVGLPCldOmfRJo0M9EgNWm0G8UbgyFcyVrAQEG63J1/AG+xUxrBBhFctz/LcVG\nhJRpz1/UCy3sonDnZwIMIIwjGdNCRp5LgPHe21oPN0Fjk6/HN98I+JWab5kAY4MHjl+wLlnCxepx\nAZZBshZqTKiGhicRYLVdETbVJk2SMrGmpt6NVFpT++SlzbW6hQCjyHj/t+ids2tYbGMovdspidJI\nLB20vBfiiX1RYRflOAojAVZb6AyIXvqhdz5cKmJ3CzC2kQHUodL4SoBlUHPqFHGaRrmV4ytREkh2\nEAVYBtcV5XI4S9YA8FW1SdvX5rp6TcJ1twDLwFkR7Qs3OJbQevcSYBTJcS9E1jnIJRL0QprcFWBJ\nxO8p2lV5vGtqeH1IZzW1z8sJsCRCfz8aKKQrpwRYbRfNdJQ1Z3lbogrsk3B7JtwKrxBgQD1BL1XV\ndNcuihhsAqy2i6rEuMApGhnFGlmNl3z9lWOLBFgSVzf88v5vPIQjP5OvF+JukAAD5uKOaNGtzERP\nPCitneDuJsAoZVB7rLPGu6Dj5qfVztQRQQ+NAOOLTF1KiaBjGVuNrR26wQXYXOjD+SiOVAV2Mi0T\nYHOmn9BNoqtyRzQbmS18Pe7ElxufatyouIOeAKNIg/d/a2+N8piObsxU6IU6e34Y/v0X1393rwBN\nG4Y/vRS61tnEsaZ9zsCAqr6eXjSYnctVevINjtshwAAa0mB+N0uAARSZnQZJmtsJsLnyU/XnnNRX\n29Ln7NLHMujvZtctCTCK5Lv/G2zyyo++1w2tSBhgB4urfJrzsAnR0F0fY5t2qUFkKseY2v5GVNvP\n7e+KFqT6GH2OHm7Kzy59VlbHouzr6Pu+64auG7quf/1Dne9HOsKfpDoDG4ahwe/bJtJ3LQ2XDvVI\n2VdkPteWVGdgby0H3GW3T5YZxj+vLlPyPOvLDNMlz3itoe/769Z5+ZCL90/5Mhcer58Zd1Rf98Db\nuUiVZYa+77puXkU/O/zPss2s80x//Wt9L9fzXmtZ+TG8GYlCmO3i2d5//fXtULv6tJvfhtmx//Y9\nquYTTp75ez+c91qbt+Lqh7wdZ29UWPbLvzZl5f5711XycTV7ofu7K27ZLS2X0EvUM7A6dcNS+zWd\nmJ1/o2EYViYQVxjHojinQzeIGmAtaKSwLl2NsUU/XaaDHeIOyq8M0wuNSPUhjtGR2tr00EZquMJq\n+JhA+xIcoBCpVq0XcvzcydUSBhhxbRrCQox3wHVcQgzPOE5Eb99sdqoxZW985QwsvMdW+WM3HBgJ\nMBoik4ByAiw8lxAJ5/bfhyQHAfZHxA8ROGsBnkmA/SEMAKIQYPs1knaNnAjCVn7gmIMEWAY6n9BM\nwthHgBGVUS80sy6OE2AZGMoJ6lW68owdBFgGj23+x254JqZf7CbA/tBLUM1r/mEiwj4CLIM0uZtm\nQyg0DNKL/QTYH1t/jrkRhgDggQL8aPQRb38F/Juh6wofNe66S1+i/rPda9O27NvwzS8RsUd2VT78\n0X7lh2zORrz9PQgA6nAJcT/pBXAjAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZA\nSAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYAB\nEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJg\nAIT0390rwMn6vh//MAzDvWsCNan8B3IGltB6A/c/Pv2v2b90k6EBWlYSXctiftsRn3qEpjgDy2YY\nhr7vP3Xy9H/NFnv9deXh0Kz1yu9W52EqPyhnYDmdMnkc+9lJGIGsp9TbfFr+4+tJlH3jnIFlMz2R\nmv2vr7PL8SGzxb5ObKEFRyq/W5x+KfsQBFgqs9Ytb7+VS4vQvt2V332YtxGCAEvldcVvdLCrtTRR\nHKz8t28G0z4BltDKjHLa57NTLp9CJrpN51Jv3+JV/LGYawAQkk8hUso3Y3gmld8slxAp4o0Bnknl\nt8wZGAAhmVxQyqeNeSaV3ywBBkBILiECEJIAAyAkAQZASAIMgJAEGAAhCTAAQkp+Jw43gOG4iF81\nUfkc137lJw+wLsIxoGVxk0Dlc0SIyncJEYCQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgC\nDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEWz7m/0xPiV3/ITRGyT8IftOz7PutP+b36fPzD\nwc0899m4V9Cyr1+Ey7CMuN8YpToD6/s+8VRuuWlHNvbcZ+NGccteEXJQqgAbhuEZk6mh687azCfs\nruRSlH399f9tIqkZV8JLiDMlVwziLDNdchgn3wdea/ac/afLOG1se+1loo9rX/fA2w2sucy44N8/\nzM8mr1mf6Yv2r4e0tn/uWiZQ5Ye8br5u+mZA0DcG3nq9WfDzD0PXdbs37u+z/XuWLLvqTFFKaLae\nIVa77/tXbr3+cPVqLyu/+f10jxAllOoSYmKL9Op+Zo6nPFuYCRe5vNKrm82lLvan8uOcbzAnwCKZ\nvtvx+uPW9huXH4bu+FNBLD/F/1v5d64NhyV8DyxfUb4iZzTdwGHo+r7r+9LLICtPRWhBD+UwdNOz\nrgpTqHE/TTJsWwfRFGdg/KOBqelTVl2XYS4w5CPAYlhJl/Krf7PTr5Vl4DlM3eISYK0rSZSSDJNM\ntGaZHFdnycrza5CIBNizrA8QpqLUUT8t1l9R5QclwAIo6a71k7CSi4ezheEud1Wgyg9HgDXtbUd9\n+p78pwxbT69A37onk2lB1inCWQuo/AQEWCqvDHv15qZzr01LwhXuqkCVH5EAa93Wvnot/4oxnUlT\nWnsDjLgSfpE5jd1ddzyxfK+Te91VgSo/FmdgQH5iKSUB1rRbuk6rc7U2v86h8sMRYPGce9e7oPfQ\nI5NLi/DTpfhPL+oNs0AEWKNu76LbV4CU1BUnEmBAi26JOtcjYhFg7dJLpFS/sLe+otPEKAQYb8hO\noH0CrEUmgDzcWVOoHa1k9haIAIvn3Hu4rTybHOVcn282fUOpuRdiAgKsUaaBcGPESLcQBBiQmblg\nYgKsOY1M/bQ9VwhRVyFWkk6AAW26PUUamUqyQoC16PbWhdP5FRVOJ8DiqXkvREMAFdxyQ871FzWJ\nDEGAtUVgwNTBjpBDuQkwoJJwcWJC2TgB1px2mrydNQFYEmBAo45MoY6fPJnAtU+ANcT1CrKKW9tx\n1/wJBFg81e6F+LPAia8Gb7gXIvsIMKCG3VfkdgfN8WuAriI2ToC1RcNAaxo5Vev7f//xIsBYI1Ch\nBdPckmEvAqwVipKsjtT2vinUid3Uwhxu3Jxh+F0Zw8VIgAF8cfsvk72iq4VAbYcAa0hhada8F+LI\ndI9LXXEvxK9PecsNGE8xrriu7AQY8MmJnxo4GBb3Dta3J93tK9AsAdYEkyla41MDLfi05889CXvN\nVMIdaAHGF2Z/DzT91MDBsTLNz4DdNbhf2oCv0Ar68RAB1go5QSNmI1p395suLbRGC+swc+IqjU91\nfLJSnwDbY3rGHehgw1fL9Bo1OILXV7nZS17uyCp9OtYHn7YmAbbZ69BOJyzHy2jL8lXvhQgv+4rl\na/idXoQlcZug8o/frT/6TEWAbTC9Xrz8Wkb8dliTe+sYrUzJV/69jtsr8JbNL0vii140QJT9d/cK\nXG69J7c+z9unGgY/fH6Ot7vxxp3T912INn7nqtXu+w1H5PbUOd2mzT/4Qlc/eYJhJ/0Z2L9DdNb7\nVeuT03ztWtPs2uz01Lbye4053t3cfmm661od1J4cmesuOl7D0HVdgF2ZPsD65Wi451kKentfhoXo\ntwqD2vLabFc9yabPP3npCEfojUtW++ETtcrpflFstzxT2Sr/JcSX2QcuthZHyfLjtcRqFxnS+LqH\nZ/OPknbddAimT5jp2JWX4nWDWv39efUravB29HFvCFbiwweNXptcMm8ZipfcvXyICe2l67l1p60/\nz9LKM08f8n6xiD3yU/nlO3bTISgvhhPLpnANr26os2r13Fe55EC3X/n5A+ztBq58ImO52NY9VP6o\nQOfyhav6aYcff+Z9Ci+trLz6ji1qwbjaxUetaLGty597ZMtftNkzsPJauu4q0ZYlA1T+gy4hTn29\nnFiYcJziup08e+bX0POQq0AXXdM+/rHbi1Rbqwr1s+n5C49IoBlzofQf4lgz+3TA8s/7jvSmN7oz\nFdMOlcfB5bf32D2onf5OZLWnyqrNWcWlHh1g3d+Pvc2i60jD3PL1w0tVO0niXHd9bjDxjRODVmy+\n069OgI1en5mefW77oK89nKyYtooV4SntHtTuKt12aubSr3Ps8MxvOAiwq6zXU9Y6a/9d3wf6VIrH\ni3D9GSrXQvs9dWN3pDz96gQYd8naUW1aZtjxjynl+4JXuaY+wTF7VPtBfiIBdqH1mW873VjuUb2R\nzLQaT/yQ7duSSPwG2FSUdog74HwlwK61MvMlZUe1bFaNQfd/Ix107zcXP/k04AQ91l8JsMtdNPMN\nrZEx6IHO/aTS22sMt4yYd1VUg5U8/ZJr7vTqBFgdOWa+QVebTBRhieHdDxamJMAqOf0z+m0q/5Xb\n3PvhOWaTsxun/O1/qKSwO87akPSjTSfAgIO+3pjtdG8vWt6lwauIz/HQeyFyIw2fT/qZfh1aYytn\nYGyjx7jX20+O3JWgkvteAowbaHvOYkb1ZAKMDS763itsUv9dt3UnVvXt2xKLAONM7oVIHRELTXec\nzoc4qE0Xc4pGCumsn/d0cWIHZ2AAR4mfWwgwNtvdq5ocOJEAAzjkxHtnsIkAoyotSlYuMNQnwNhm\nPYHK74UIT7PSHfpmHwFGJVqUxFxauIUAAziHWVplAow99jWqWSp8ojt2EGDUYGZKehKoPgHGZhoV\nTmR6t5sA40zu9sbDraSR7jidAKMSzQuf6I59BBg7lV/3cIWEh5BDlQkwgDOZsVUjwKjBzBTeknZH\nCDD2EEjwltaoSYBxpuXd3kwweaC3Zf/pXogybzcBxn7CCbiRAONyJpg8SnnBmwIeJMDYSSwB9xJg\nXMgEk8cqLH4TwSMEGIeIKFgSS3UIMM60vNubTuaxZtO7WXeY/B0nwNhvPZz0J3ApAQZwvpJrD65P\nHCTAOGr19yMqrge051N3uD5xCgHGJfQncDUBxiHOseCTsTuWk7nxX/TOcQKME4wN2ff9eLc3/QlL\nn+6FyG4JA0yVVPYzzfy32+3/W9jtbXqdhP3M7XrTuxP9d/cKnEkP32rouv7nD/qzHmXPY6U6AxuG\nYflFWq7W9/0kuoaffzOqVqLsmzftjmHyLxyV6gzsreVIuux2yxxfpuv6r83Z2jqvLxM9g7/ugbcb\naJnTlxkX/OmO7m2DNLXOgSq/Dzp3m+3i2d5//XX6Z65TlnAhNVVChWW//Cs3Kg65toQooahnYO3v\nWTidso9oGIaVmQdHRA0wmjI25NilmhNmXg2iO86VMMCUyF3s+RvZ+e1zjE6X6lOIADyHAAMgJAEG\nQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACElPD3wGbe\n/p43pKfySc8vhAIQkkuIAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYAB\nEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJg\nAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAE\nGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJD+D1wjsnR5DRtmAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 39, "cell_type": "code", "language": "python", "metadata": {}, "input": ["exo2()"]}, {"collapsed": false, "outputs": [], "prompt_number": 40, "cell_type": "code", "language": "python", "metadata": {}, "input": ["%% Insert your code here."]}, {"source": ["Bibliography\n", "------------\n", "<html><a name=\"biblio\"></a></html>\n", "\n", "\n", "* [LionsMercier79] P. L. Lions and B. Mercier, _Splitting Algorithms for the Sum of Two Nonlinear Operators_, SIAM Journal on Numerical Analysis, Vol. 16, No. 6 (Dec., 1979), pp. 964-979.\n", "* [CombettesPesquet10] P.L. Combettes and J-C. Pesquet, _Proximal Splitting Methods in Signal Processing_, in: Fixed-Point Algorithms for Inverse Problems in Science and Engineering, New York: Springer-Verlag, 2010.\n", "* [DuvalPeyre13] V. Duval and G. Peyre, _Exact Support Recovery for Sparse Spikes Deconvolution_, preprint hal-00839635, 2013\n", "* [deCastroGamboa12] Y. de Castro and F. Gamboa. _Exact reconstruction using beurling minimal extrapolation_. Journal of Mathematical Analysis and Applications, 395(1):336-354, 2012.\n", "* [CandesFernandezGranda13]  E. J. Candes and C. Fernandez-Granda. _Towards a mathematical theory of super-resolution_. Communications on Pure and Applied Mathematics. To appear., 2013.\n", "* [BrediesPikkarainen13] K. Bredies and H.K. Pikkarainen. Inverse problems in spaces of measures. ESAIM: Control, Optimisation and Calculus of Variations, 19:190-218, 2013.\n", "* [Dumitrescu07] B. Dumitrescu, Positive Trigonometric Polynomials and Signal Processing Applications, Springer, 2007."], "metadata": {}, "cell_type": "markdown"}], "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"}]}}}