{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cb1e03ec-c339-4bd3-b492-a402b85413f6",
   "metadata": {},
   "source": [
    "# Optimisation Framework\n",
    "\n",
    "The $\\textbf{cil.optimisation}$ framework contains three structures, namely $\\texttt{Function}$\\, $\\texttt{Operator}$ and $\\texttt{Algorithm}$ that formalise a generic optimisation problem for imaging applications as\n",
    "$$\n",
    "u^{*} =\\argmin_{u\\in\\mathbb{X}} f(Ku) + g(u) \\equiv \\argmin_{u\\in \\mathbb{X}} \\sum_{i=0}^{n-1}f_{i}K_{i}(u) + g(u).\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "We let $\\mathbb{X}, \\mathbb{Y}$ denote finite-dimensional vector spaces, $K:\\mathbb{X}\\rightarrow \\mathbb{Y}$ a linear operator with operator norm $\\|K\\| = \\max \\{ \\|Ku\\|_{Y}: \\|u\\|_{\\mathbb{X}}\\leq 1 \\}$ and proper, convex functions \n",
    "- $f: \\mathbb{Y}\\rightarrow\\overline{\\mathbb{R}}$\n",
    "- $g: \\mathbb{X}\\rightarrow\\overline{\\mathbb{R}}$\n",
    "\n",
    "---\n",
    "<a id='BlockFunction'></a>\n",
    "**Note:** In certain cases, it is convenient to decompose $\\mathbb{Y}=Y_{0}\\times\\dots \\times Y_{n-1}$, $n\\geq1$ and consider a separable function, e.g., $\\texttt{BlockFunction}$, \n",
    "\n",
    "- $f(y) := f(y_{0},\\dots,y_{n-1}) = \\sum_{i=0}^{n-1}f_{i}(y_{i})$\n",
    "\n",
    "which results to the right-side formulation in $eqn.(1)$.\n",
    "\n",
    "With the above form, $f(y)$ is a __separable sum__ of decoupled functions, which is useful when we need to compute the proximal operator of $f$, i.e., \n",
    "\n",
    "$$\\mathrm{prox}_{\\tau f}(z) = \n",
    "\\begin{bmatrix}\n",
    "\\mathrm{prox}_{\\tau f_0}(z_0)\\\\\n",
    "\\mathrm{prox}_{\\tau f_1}(z_1)\\\\\n",
    "\\vdots\\\\\n",
    "\\mathrm{prox}_{\\tau f_{n-1}}(z_{n-1})\n",
    "\\end{bmatrix} = \n",
    "\\begin{bmatrix}\n",
    "\\underset{y_{0}\\in \\mathbb{Y_{0}}}{\\operatorname{argmin}} \\tau f_{0}(y_{0}) + \\frac{1}{2}\\|y_{0} - z_{0}\\|^{2}\\\\\n",
    "\\underset{y_{1}\\in \\mathbb{Y_{1}}}{\\operatorname{argmin}} \\tau f_{1}(y_{1}) + \\frac{1}{2}\\|y_{1} - z_{1}\\|^{2}\\\\\n",
    "\\vdots\\\\\n",
    "\\underset{y_{{n-1}}\\in \\mathbb{Y_{n-1}}}{\\operatorname{argmin}} \\tau f_{{n-1}}(y_{{n-1}}) + \\frac{1}{2}\\|y_{{n-1}} - z_{{n-1}}\\|^{2}\\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "# Gallery Algorithms\n",
    "\n",
    "1. [Simultaneous Iterative Reconstruction Technique](#SIRT)\n",
    "2. [A Fast Iterative Shrinkage-Thresholding Algorithm](#FISTA)\n",
    "3. [Primal Dual Hybrid Gradient](#PDHG)\n",
    "4. [Stochastic Primal Dual Hybrid Gradient](#SPDHG)\n",
    "\n",
    "<hr style=\"border:2px solid gray\"> </hr>\n",
    "\n",
    "## SIRT\n",
    "\n",
    "**Initialize:** $x_{0}\\in \\mathbb{X}$\\\n",
    "**Update:** $x_{k}$\n",
    "\n",
    "**for $k\\geq0$**\n",
    "\n",
    "   $x_{k+1} = \\mathcal{P}_{C}(x_{k} + D \\mathcal{A}^{T} M ( g - \\mathcal{A} x_{k}))$\n",
    "\n",
    "where,\n",
    " $$\n",
    "    \\begin{cases}\n",
    "    M = \\frac{\\mathbb{1}}{\\mathcal{A}\\mathbb{1}}, \\quad m_{ii} = \\frac{1}{\\sum_{j} a_{ij}}, \\quad\\text{ sum over columns }\\\\\n",
    "    D = \\frac{\\mathbb{1}}{\\mathcal{A}^{T}\\mathbb{1}}, \\quad d_{jj} = \\frac{1}{\\sum_{i} a_{ij}},\\quad\\text{ sum over rows }\\\\\n",
    "    \\end{cases}\n",
    "    $$\n",
    "    \n",
    "<hr style=\"border:2px solid gray\"> </hr>    \n",
    "\n",
    "## FISTA\n",
    "\n",
    "**Inputs:** $\\tau = \\frac{1}{L}$, $L$ Lipschitz constant of $\\nabla \\mathcal{F}$, $t_{0}=1$\n",
    "\n",
    "\n",
    "**Initialize:** $x_{0}\\in \\mathbb{X}$\\\n",
    "**Update:** $x_{k}, t_{k}, y_{k}$\n",
    "\n",
    "**for $k\\geq0$**\n",
    "\n",
    "1. $x_{k} = \\mathrm{prox}_{\\tau \\mathcal{G}}(y_{k} - \\tau\\nabla\\mathcal{F}(y^{k}))$\n",
    "\n",
    "2. $t_{k+1} = 1 + \\frac{\\sqrt{1+4 t_{k}^{2}}}{2}$\n",
    "\n",
    "3. $y_{k+1} = x_{k+1} + \\frac{t_{k}-1}{t_{k+1}}(x_{k}-x_{k-1})$\n",
    "\n",
    "\n",
    "    \n",
    "<hr style=\"border:2px solid gray\"> </hr>    \n",
    "\n",
    "\n",
    "## PDHG\n",
    "\n",
    "**Inputs:** $\\tau,\\sigma>0$, $\\tau\\sigma\\|K\\|^{2}<1$, $\\theta\\in [0,1]$\n",
    "\n",
    "\n",
    "**Initialize:** $x_{0}\\in \\mathbb{X}$, $y_{0}\\in \\mathbb{Y}$,$\\overline{x}_{0}=x_{0}$\n",
    "**Update:** $x_{k}, y_{k}, \\overline{x_{k}}$\n",
    "\n",
    "**for $k\\geq0$**\n",
    "\n",
    "1. $y_{k+1} = \\mathrm{prox}_{\\sigma \\mathcal{F}^{*}}(y_{k} + \\sigma K \\overline{x}_{k})$\n",
    "\n",
    "2. $x_{k+1} = \\mathrm{prox}_{\\tau \\mathcal{G}}(x_{k} - \\tau K y_{k+1})$\n",
    "\n",
    "3. $\\overline{x}_{k+1} = \\overline{x}_{k+1} + \\theta (x_{k+1} - x_{k})$\n",
    "    \n",
    "<hr style=\"border:2px solid gray\"> </hr>  \n",
    "\n",
    "\n",
    "# SPDHG\n",
    "\n",
    "**Inputs:** $n = \\# \\mathrm{subsets}$, $\\gamma=1.0$,\n",
    "$\\rho=0.99$,$\\mathrm{probability}$ $p_{i}$, $i=0,\\dots,n-1$\n",
    "$\\\\\\mathrm{Stepsizes}: S_{i} = \\gamma\\frac{\\rho}{\\|A_{i}\\|},\\, T = \\min\\limits_{i} \\frac{1}{\\gamma}\\frac{\\rho p_{i}}{\\|A_{i}\\|},\n",
    "i=0,\\dots,n-1$\n",
    "    \n",
    "    \n",
    "**Initialize:** $x^{0}, z^{0}=\\overline{z}^{0}, y^{0}$\\\n",
    "**Update:** $x_{k}, y_{k}, \\overline{z_{k}}$\n",
    "\n",
    "**for $k\\geq0$**\n",
    "\n",
    "1. $x^{k+1} = \\mathrm{prox}^{T}_{\\mathcal{G}}(x^{k} - T\\overline{z}^{k})$ \n",
    "\n",
    "2. Select $i\\in [n]$ at random with probability $p_{i}$\n",
    "\n",
    "3. $y_{i}^{k+1} = \\mathrm{prox}^{S_{i}}_{\\mathcal{F}_{i}^{*}}(y_{i}^{k} + S_{i}A_{i}x^{k+1})$\n",
    "\n",
    "4. $\\Delta z = A_{i}^{T}(y_{i}^{k+1} - y_{i}^{k})$\n",
    "\n",
    "5. $z^{k+1} = z^{k} + \\Delta z$\n",
    "\n",
    "6. $\\overline{z}^{k+1} = z^{k+1} + \\frac{\\Delta z}{p_{i}}$\n",
    "    \n",
    "<hr style=\"border:2px solid gray\"> </hr>\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "43cbf82c2f716cd564b762322e13d4dbd881fd8a341d231fe608abc3118da208"
  },
  "kernelspec": {
   "display_name": "Python 3.9.13 ('cil_22.0.0')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}