{
 "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": [],
   "source": [
    "using PyPlot"
   ]
  },
  {
   "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": [
    "<h3><font color='magenta'>\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",
    "</font></h3>"
   ]
  },
  {
   "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": [
    "<h3><font color='magenta'>\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",
    "</font></h3>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}