{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Chambolle-Pock Primal-Dual Splitting Algorithm\n",
"==============================\n",
"$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have seen in the lab 3 that total variation denoising can be performed using the dual forward-backward algorithm. But the setting is restrictive: this algorithm cannot be applied to general inverse problems. \n",
"\n",
"This tour explores the primal-dual proximal splitting algorithm proposed in \n",
"\n",
"A. Chambolle and T. Pock, \"A First-order primal-dual algorithm for convex problems with application to imaging,\"\n",
"_Journal of Mathematical Imaging and Vision_,\n",
"vol. 40, no. 1, 2011\n",
"\n",
"and further analyzed and extended in \n",
"\n",
"L. Condat, \"A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms,\" _J. Optimization Theory and Applications_, vol. 158, no. 2, 2013."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"from __future__ import division\n",
"%pylab inline\n",
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Convex Optimization with a Primal-Dual Scheme\n",
"---------------------------------------------\n",
"\n",
"We consider a (primal) optimization problem of the form\n",
"$$ \\umin{x} f(x) + g(Lx) $$\n",
"where $f$ and $g$ are convex functions, whose proximity operators can be computed, and $L$\n",
"is a linear operator.\n",
"\n",
"The dual problem is \n",
"\n",
"$$ \\umin{u} f^*(-L^*u) + g^*(u) $$\n",
"\n",
"The (relaxed) Chambolle-Pock algorithm takes initial estimates $x^{(0)}$ and $u^{(0)}$ of the primal and dual solutions, a parameter $\\tau>0$, a second parameter $0<\\sigma\\leq 1/(\\tau\\|L\\|^2)$, a relaxation parameter $0<\\rho<2$, and iterates, for $k=1,2,\\ldots$\n",
"$$ \\left|\\begin{array}{l}\n",
"\\tilde{x}^{(k)} = \\mathrm{prox}_{\\tau f}( x^{(k-1)}-\\tau L^*(u^{(k-1)}) ) \\\\\n",
" \\tilde{u}^{(k)} = \\mathrm{prox}_{\\sigma g^*}( u^{(k-1)}+ \\sigma L(2\\tilde{x}^{(k)}-x^{(k-1)}) \\\\\n",
" x^{(k)}= x^{(k-1)} + \\rho (\\tilde{x}^{(k)}-x^{(k-1)})\\\\\n",
" u^{(k)}= u^{(k-1)} + \\rho (\\tilde{u}^{(k)}-u^{(k-1)})\n",
" \\end{array}\\right.$$\n",
" \n",
" Then, $x^{(k)}$ converges to a primal solution $x^\\star$ and $u^{(k)}$ converges to a dual solution $u^\\star$.\n",
" \n",
" In practice, like for the Douglas-Rachford algorithm, it is always interesting to take $\\rho$ close to $2$, e.g. $\\rho=1.9$, instead of $\\rho=1$ like in the paper of Chambolle & Pock. Also, for fixed $\\tau$, the higher $\\sigma$, the better; so, one can set $\\sigma=1/(\\tau\\|L\\|^2)$, which leaves only the parameter $\\tau$ to tune.\n",
"\n",
"With this choice of $\\sigma$, the algorithm exactly reverts to the Douglas-Rachford algorithm when $L=\\mathrm{Id}$ (replacing $\\sigma$ by $1/\\tau$ in the algorithm). So, it is a natural extension of the latter.\n",
"\n",
"\n",
"We recall that being able to compute the proximity operator of $f^*$ is\n",
"equivalent to being able to compute the proximity operator of $f$, thanks to the Moreau identity\n",
"$$ x = \\mathrm{prox}_{\\gamma f^*}(x) + \\gamma \\mathrm{prox}_{f/\\gamma}(x/\\gamma) $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Image Inpainting\n",
"---------------------------------------------\n",
"\n",
"Like in the lab 1, we want to reconstruct an estimate of the Lena image from a random subset of its pixels. So, we want to solve \n",
"$$\\umin{x} \\mathrm{TV}(x)\\quad\\mbox{s.t.}\\quad Ax=b,$$\n",
"where we keep the notations of the labs 1 and 3: $A$ is the degradation operator which multiplies the image by a binary mask and $\\mathrm{TV}$ is the total variation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### \n",
"__Exercise: write the code of the Chambolle-Pock algorithm to solve this inpainting problem and apply it to the Lena image.__\n",
"\n",
"Try different values of $\\tau$ and $\\rho$ and observe the convergence speed by monitoring the decay of $\\mathrm{TV}(x^{(k)})$.\n",
"\n",
"Compare the inpainted image with the one obtained in the lab 1, with Tikhonov instead of total variation regularization.\n",
"

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Image Denoising\n",
"---------------------------------------------\n",
"Like in the lab 3, we want to denoise the noisy Lena image by solving \n",
"$$\\umin{x} \\frac{1}{2}\\|x-y\\|^2+\\lambda\\mathrm{TV}(x).$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### \n",
"__Exercise: write the code of the Chambolle-Pock algorithm to solve this denoising problem and apply it to the Lena image.__\n",
"\n",
"Try different values of $\\tau$ and $\\rho$ and observe the convergence speed by monitoring the sum of the primal and dual energies, like in the lab 3. Compare with the accelerated forward-backward algorithm on the dual problem developed in the lab 3.\n",
"

"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}