{
"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",
"\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",
"
\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",
"
\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",
"
\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",
"
\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",
"
\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
}