Algoritmo de Chambolle-Pock

Optimización del algoritmo Primal-Dual para problemas convexos

En matemáticas , el algoritmo Chambolle-Pock es un algoritmo utilizado para resolver problemas de optimización convexa . Fue introducido por Antonin Chambolle y Thomas Pock ^[1] en 2011 y desde entonces se ha convertido en un método ampliamente utilizado en varios campos, incluido el procesamiento de imágenes , ^{[2] [3] [4]} la visión artificial , ^[5] y el procesamiento de señales . ^[6]

El algoritmo Chambolle-Pock está diseñado específicamente para resolver de manera eficiente problemas de optimización convexa que involucran la minimización de una función de costo no suave compuesta por un término de fidelidad de datos y un término de regularización. ^{[1] Esta es una configuración típica que surge comúnmente en} problemas inversos de imágenes mal planteadas , como reconstrucción de imágenes , ^[2] eliminación de ruido ^[3] y restauración de imagen . ^[4]

El algoritmo se basa en una formulación primal-dual, que permite actualizaciones simultáneas de variables primales y duales . Al emplear el operador proximal , el algoritmo Chambolle-Pock maneja de manera eficiente términos de regularización no suaves y no convexos, como la variación total , específica del marco de imágenes. ^[1]

Planteamiento del problema

Sean dos espacios vectoriales reales dotados de un producto interior y una norma . Hasta ahora, una función se denomina simple si su operador proximal tiene una representación en forma cerrada o puede calcularse con precisión, para , ^[1] donde se hace referencia a ${\displaystyle {\mathcal {X}},{\mathcal {Y}}}$ ${\displaystyle \langle \cdot ,\cdot \rangle }$ ${\displaystyle \lVert \,\cdot \,\rVert =\langle \cdot ,\cdot \rangle ^{\frac {1}{2}}}$ ${\displaystyle F}$ ${\displaystyle {\text{prox}}_{\tau F}}$ ${\displaystyle \tau >0}$ ${\displaystyle {\text{prox}}_{\tau F}}$

$${\displaystyle x={\text{prox}}_{\tau F}({\tilde {x}})={\text{arg }}\min _{x'\in {\mathcal {X}}}\left\{{\frac {\lVert x'-{\tilde {x}}\rVert ^{2}}{2\tau }}+F(x')\right\}}$$

Consideremos el siguiente problema primal restringido: ^[1]

$${\displaystyle \min _{x\in {\mathcal {X}}}F(Kx)+G(x)}$$

donde es un operador lineal acotado , son convexos , semicontinuos inferiores y simples. ^[1] ${\displaystyle K:{\mathcal {X}}\rightarrow {\mathcal {Y}}}$ ${\displaystyle F:{\mathcal {Y}}\rightarrow [0,+\infty ),G:{\mathcal {X}}\rightarrow [0,+\infty )}$

El problema de minimización tiene su problema correspondiente dual como ^[1]

$${\displaystyle \max _{y\in {\mathcal {Y}}}-\left(G^{*}(-K^{*}y)+F^{*}(y)\right)}$$

donde y son el mapa dual de y , respectivamente. ^[1] ${\displaystyle F^{*},G^{*}}$ ${\displaystyle K^{*}}$ ${\displaystyle F,G}$ ${\displaystyle K}$

Supongamos que los problemas primal y dual tienen al menos una solución , lo que significa que satisfacen ^[7] ${\displaystyle ({\hat {x}},{\hat {y}})\in {\mathcal {X}}\times {\mathcal {Y}}}$

$${\displaystyle {\begin{aligned}K{\hat {x}}&\in \partial F^{*}({\hat {y}})\\-(K^{*}{\hat {y}})&\in \partial G({\hat {x}})\end{aligned}}}$$

donde y son el subgradiente de las funciones convexas y , respectivamente. ^[7] ${\displaystyle \partial F^{*}}$ ${\displaystyle \partial G}$ ${\displaystyle F^{*}}$ ${\displaystyle G}$

El algoritmo de Chambolle-Pock resuelve el llamado problema del punto de silla ^[1]

$${\displaystyle \min _{x\in {\mathcal {X}}}\max _{y\in {\mathcal {Y}}}\langle Kx,y\rangle +G(x)-F^{*}(y)}$$

que es una formulación primal-dual de los problemas primal y dual no lineales planteados anteriormente. ^[1]

Algoritmo

El algoritmo de Chambolle-Pock implica principalmente alternar iterativamente entre ascendente en la variable dual y descendente en la variable primaria utilizando un enfoque de tipo gradiente , con tamaños de paso y respectivamente, para resolver simultáneamente el problema primario y el dual. ^[2] Además, se emplea una técnica de sobre- relajación para la variable primaria con el parámetro . ^[1] ${\displaystyle y}$ ${\displaystyle x}$ ${\displaystyle \sigma }$ ${\displaystyle \tau }$ ${\displaystyle \theta }$

Algoritmo Algoritmo Chambolle-Pock Entrada: y conjunto ,  ${\displaystyle F,G,K,\tau ,\sigma >0,\,\theta \in [0,1],\,(x^{0},y^{0})\in {\mathcal {X}}\times {\mathcal {Y}}}$ ${\displaystyle {\overline {x}}^{0}=x^{0}}$ criterio de detención .

 ${\displaystyle k\leftarrow 0}$hacer mientras  se detiene criterio  no satisfecho  fin hacer ${\displaystyle y^{n+1}\leftarrow {\text{prox}}_{\sigma F^{*}}\left(y^{n}+\sigma K{\overline {x}}^{n}\right)}$  ${\displaystyle x^{n+1}\leftarrow {\text{prox}}_{\tau G}\left(x^{n}-\tau K^{*}y^{n+1}\right)}$  ${\displaystyle {\overline {x}}^{n+1}\leftarrow x^{n+1}+\theta \left(x^{n+1}-x^{n}\right)}$  ${\displaystyle k\leftarrow k+1}$

• "←" denota asignación . Por ejemplo, " el elemento más grande ← " significa que el valor del elemento más grande cambia al valor del elemento .
• " return " finaliza el algoritmo y genera el siguiente valor.

Chambolle y Pock demostraron ^[1] que el algoritmo converge si y , secuencialmente y con una tasa de convergencia para la brecha primal-dual. Esto ha sido ampliado por S. Banert et al. ^[8] para que se cumpla siempre que y . ${\displaystyle \theta =1}$ ${\displaystyle \tau \sigma \lVert K\rVert ^{2}\leq 1}$ ${\displaystyle {\mathcal {O}}(1/N)}$ ${\displaystyle \theta >1/2}$ ${\displaystyle \tau \sigma \lVert K\rVert ^{2}<4/(1+2\theta )}$

El método semi-implícito de Arrow-Hurwicz ^[9] coincide con la elección particular de en el algoritmo de Chambolle-Pock. ^[1] ${\displaystyle \theta =0}$

Aceleración

Existen casos especiales en los que la tasa de convergencia tiene una aceleración teórica. ^[1] De hecho, si , respectivamente , es uniformemente convexo entonces , respectivamente , tiene un gradiente continuo de Lipschitz . Entonces, la tasa de convergencia se puede mejorar a , proporcionando un ligero cambio en el algoritmo de Chambolle-Pock. Esto conduce a una versión acelerada del método y consiste en elegir iterativamente , y también , en lugar de fijar estos valores. ^[1] ${\displaystyle G}$ ${\displaystyle F^{*}}$ ${\displaystyle G^{*}}$ ${\displaystyle F}$ ${\displaystyle {\mathcal {O}}(1/N^{2})}$ ${\displaystyle \tau _{n},\sigma _{n}}$ ${\displaystyle \theta _{n}}$

En caso de convexidad uniforme, con convexidad uniforme constante, el algoritmo modificado se convierte en ^[1] ${\displaystyle G}$ ${\displaystyle \gamma >0}$

Algoritmo acelerado de Chambolle-Pock Entrada: tal que y conjunto ,  ${\displaystyle F,G,\tau _{0},\sigma _{0}>0}$ ${\displaystyle \tau _{0}\sigma _{0}L^{2}\leq 1,\,(x^{0},y^{0})\in {\mathcal {X}}\times {\mathcal {Y}}}$ ${\displaystyle {\overline {x}}^{0}=x^{0}.}$ criterio de detención .

 ${\displaystyle k\leftarrow 0}$hacer mientras  se detiene criterio  no satisfecho  fin hacer ${\displaystyle y^{n+1}\leftarrow {\text{prox}}_{\sigma _{n}F^{*}}\left(y^{n}+\sigma _{n}K{\overline {x}}^{n}\right)}$  ${\displaystyle x^{n+1}\leftarrow {\text{prox}}_{\tau _{n}G}\left(x^{n}-\tau _{n}K^{*}y^{n+1}\right)}$  ${\displaystyle \theta _{n}\leftarrow {\frac {1}{\sqrt {1+2\gamma \tau _{n}}}}}$  ${\displaystyle \tau _{n+1}\leftarrow \theta _{n}\tau _{n}}$  ${\displaystyle \sigma _{n+1}\leftarrow {\frac {\sigma _{n}}{\theta _{n}}}}$  ${\displaystyle {\overline {x}}^{n+1}\leftarrow x^{n+1}+\theta _{n}\left(x^{n+1}-x^{n}\right)}$  ${\displaystyle k\leftarrow k+1}$

• "←" denota asignación . Por ejemplo, " el elemento más grande ← " significa que el valor del elemento más grande cambia al valor del elemento .
• " return " finaliza el algoritmo y genera el siguiente valor.

Además, la convergencia del algoritmo se ralentiza cuando , la norma del operador , no se puede estimar fácilmente o puede ser muy grande. Al elegir los precondicionadores adecuados y , modificando el operador proximal con la introducción de la norma inducida a través de los operadores y , se garantizará la convergencia del algoritmo precondicionado propuesto. ^[10] ${\displaystyle L}$ ${\displaystyle K}$ ${\displaystyle T}$ ${\displaystyle \Sigma }$ ${\displaystyle T}$ ${\displaystyle \Sigma }$

Solicitud

Ejemplo de eliminación de ruido

Una aplicación típica de este algoritmo es en el marco de eliminación de ruido de imágenes , basado en la variación total. ^[3] Opera sobre el concepto de que las señales que contienen detalles excesivos y potencialmente erróneos exhiben una alta variación total, que representa la integral del gradiente de valor absoluto de la imagen. ^[3] Al adherirse a este principio, el proceso apunta a disminuir la variación total de la señal mientras mantiene su similitud con la señal original, eliminando efectivamente los detalles no deseados mientras se preservan características cruciales como los bordes. En el entorno discreto bidimensional clásico, ^[11] considere , donde un elemento representa una imagen con los valores de píxeles colocados en una cuadrícula cartesiana . ^[1] ${\displaystyle {\mathcal {X}}=\mathbb {R} ^{NM}}$ ${\displaystyle u\in {\mathcal {X}}}$ ${\displaystyle N\times M}$

Defina el producto interno como ^[1] ${\displaystyle {\mathcal {X}}}$

${\displaystyle \langle u,v\rangle _{\mathcal {X}}=\sum _{i,j}u_{i,j}v_{i,j},\quad u,v\in {\mathcal {X}}}$

que induce una norma sobre , denotada como . ^[1] ${\displaystyle L^{2}}$ ${\displaystyle {\mathcal {X}}}$ ${\displaystyle \lVert \,\cdot \,\rVert _{2}}$

Por lo tanto, el gradiente de se calcula con las diferencias finitas estándar , ${\displaystyle u}$

${\displaystyle \left(\nabla u\right)_{i,j}=\left({\begin{aligned}\left(\nabla u\right)_{i,j}^{1}\\\left(\nabla u\right)_{i,j}^{2}\end{aligned}}\right)}$

que es un elemento del espacio , donde ^[1] ${\displaystyle {\mathcal {Y}}={\mathcal {X}}\times {\mathcal {X}}}$

${\displaystyle {\begin{aligned}&\left(\nabla u\right)_{i,j}^{1}=\left\{{\begin{aligned}&{\frac {u_{i+1,j}-u_{i,j}}{h}}&{\text{ if }}i<M\\&0&{\text{ if }}i=M\end{aligned}}\right.,\\&\left(\nabla u\right)_{i,j}^{2}=\left\{{\begin{aligned}&{\frac {u_{i,j+1}-u_{i,j}}{h}}&{\text{ if }}j<N\\&0&{\text{ if }}j=N\end{aligned}}\right.\end{aligned}}}$

Se define una norma basada como ^[1] ${\displaystyle {\mathcal {Y}}}$ ${\displaystyle L^{1}-}$

${\displaystyle \lVert p\rVert _{1}=\sum _{i,j}{\sqrt {\left(p_{i,j}^{1}\right)^{2}+\left(p_{i,j}^{2}\right)^{2}}},\quad p\in {\mathcal {Y}}.}$

Entonces, el problema primal del modelo ROF , propuesto por Rudin, Osher y Fatemi, ^[12] está dado por ^[1]

$${\displaystyle h^{2}\min _{u\in {\mathcal {X}}}\lVert \nabla u\rVert _{1}+{\frac {\lambda }{2}}\lVert u-g\rVert _{2}^{2}}$$

donde es la solución desconocida y los datos ruidosos dados, en cambio describe el equilibrio entre la regularización y el ajuste de los datos. ^[1] ${\displaystyle u\in {\mathcal {X}}}$ ${\displaystyle g\in {\mathcal {X}}}$ ${\displaystyle \lambda }$

La formulación primal-dual del problema ROF se formula de la siguiente manera ^[1]

$${\displaystyle \min _{u\in {\mathcal {X}}}\max _{p\in {\mathcal {Y}}}-\langle u,{\text{div}}\,p\rangle _{\mathcal {X}}+{\frac {\lambda }{2}}\lVert u-g\rVert _{2}^{2}-\delta _{P}(p)}$$

donde la función indicadora se define como ^[1]

$${\displaystyle \delta _{P}(p)=\left\{{\begin{aligned}&0,&{\text{if }}p\in P\\&+\infty ,&{\text{if }}p\notin P\end{aligned}}\right.}$$

sobre el conjunto convexo que puede verse como bolas unitarias con respecto a la norma definida en . ^[1] ${\displaystyle P=\left\{p\in {\mathcal {Y}}\,:\,\max _{i,j}{\sqrt {\left(p_{i,j}^{1}\right)^{2}+\left(p_{i,j}^{2}\right)^{2}}}\leq 1\right\},}$ ${\displaystyle L^{\infty }}$ ${\displaystyle {\mathcal {Y}}}$

Obsérvese que las funciones involucradas en la formulación primal-dual establecida son simples, ya que su operador proximal se puede calcular fácilmente ^[1]. El problema de eliminación de ruido de variación total de la imagen también se puede tratar con otros algoritmos ^[13] como el método de dirección alternada de multiplicadores (ADMM), ^[14] (sub)gradiente proyectado ^[15] o umbralización de contracción iterativa rápida. ^[16] $${\displaystyle {\begin{aligned}p&={\text{prox}}_{\sigma F^{*}}({\tilde {p}})&\iff p_{i,j}&={\frac {{\tilde {p}}_{i,j}}{\max\{1,|{\tilde {p}}_{i,j}|\}}}\\u&={\text{prox}}_{\tau G}({\tilde {u}})&\iff u_{i,j}&={\frac {{\tilde {u}}_{i,j}+\tau \lambda g_{i,j}}{1+\tau \lambda }}\end{aligned}}}$$

Implementación

• El paquete Manopt.jl ^[17] implementa el algoritmo en Julia
• Gabriel Peyré implementa el algoritmo en MATLAB , ^{[nota 1]} Julia, R y Python ^[18]
• En la Biblioteca de Discretización de Operadores (ODL), ^[19] una biblioteca de Python para problemas inversos, chambolle_pock_solverimplementa el método.

Véase también

• Método de dirección alternada de multiplicadores
• Optimización convexa
• Operador proximal
• Eliminación de ruido por variación total

Notas

1. Estos códigos se utilizaron para obtener las imágenes del artículo.

Referencias

1. Chambolle, Antonin; Pock, Thomas (1 de mayo de 2011). "Un algoritmo primal-dual de primer orden para problemas convexos con aplicaciones en imágenes". Revista de imágenes y visión matemática . 40 (1): 120–145. doi :10.1007/s10851-010-0251-1. ISSN  1573-7683. S2CID  207175707.
2. Sidky, Emil Y; Jørgensen, Jakob H; Pan, Xiaochuan (21 de mayo de 2012). "Prototipado de problemas de optimización convexa para reconstrucción de imágenes en tomografía computarizada con el algoritmo Chambolle–Pock". Física en Medicina y Biología . 57 (10): 3065–3091. arXiv : 1111.5632 . Bibcode :2012PMB....57.3065S. doi :10.1088/0031-9155/57/10/3065. ISSN  0031-9155. PMC 3370658 . PMID  22538474. 
3. Fang, Faming; Li, Fang; Zeng, Tieyong (13 de marzo de 2014). "Eliminación de neblina y ruido de una sola imagen: un enfoque variacional rápido". Revista SIAM sobre ciencias de la imagen . 7 (2): 969–996. doi :10.1137/130919696. ISSN  1936-4954.
4. Allag, A.; Benammar, A.; Drai, R.; Boutkedjirt, T. (1 de julio de 2019). "Reconstrucción de imágenes tomográficas en el caso de un número limitado de proyecciones de rayos X mediante la técnica de pintura de sinogramas". Revista rusa de pruebas no destructivas . 55 (7): 542–548. doi :10.1134/S1061830919070027. ISSN  1608-3385. S2CID  203437503.
5. Pock, Thomas; Cremers, Daniel; Bischof, Horst; Chambolle, Antonin (2009). "Un algoritmo para minimizar la función Mumford-Shah". 2009 IEEE 12th International Conference on Computer Vision . págs. 1133–1140. doi :10.1109/ICCV.2009.5459348. ISBN 978-1-4244-4420-5.S2CID15991312  .​
6. "Un algoritmo proximal genérico para optimización convexa: aplicación a la minimización de variación total". IEEE Signal Processing Letters . 21 (8): 985–989. 2014. Bibcode :2014ISPL...21..985.. doi :10.1109/LSP.2014.2322123. ISSN  1070-9908. S2CID  8976837.
7. Ekeland, Ivar; Témam, Roger (1999). Análisis convexo y problemas variacionales. Sociedad de Matemáticas Industriales y Aplicadas. p. 61. doi :10.1137/1.9781611971088. ISBN 978-0-89871-450-0.
8. Banert, Sebastian; Upadhyaya, Manu; Giselsson, Pontus (2023). "El método de Chambolle-Pock converge débilmente con y ". arXiv : 2309.03998 [math.OC]. ${\displaystyle \theta >1/2}$ ${\displaystyle \tau \sigma \lVert L\rVert ^{2}<4/(1+2\theta )}$
9. Uzawa, H. (1958). "Métodos iterativos para programación cóncava". En Arrow, KJ; Hurwicz, L.; Uzawa, H. (eds.). Estudios en programación lineal y no lineal . Stanford University Press.
10. Pock, Thomas; Chambolle, Antonin (6 de noviembre de 2011). "Preacondicionamiento diagonal para algoritmos primal-dual de primer orden en optimización convexa". 2011 International Conference on Computer Vision . págs. 1762–1769. doi :10.1109/ICCV.2011.6126441. ISBN 978-1-4577-1102-2.S2CID 17485166  .
11. Chambolle, Antonin (1 de enero de 2004). "Un algoritmo para la minimización de la variación total y sus aplicaciones". Revista de imágenes y visión matemática . 20 (1): 89–97. doi :10.1023/B:JMIV.0000011325.36760.1e. ISSN  1573-7683. S2CID  207622122.
12. Getreuer, Pascal (2012). "Eliminación de ruido de la variación total Rudin-Osher-Fatemi mediante Split Bregman" (PDF) .
13. Esser, Ernie; Zhang, Xiaoqun; Chan, Tony F. (2010). "Un marco general para una clase de algoritmos primal-dual de primer orden para la optimización convexa en la ciencia de la imagen". Revista SIAM sobre ciencias de la imagen . 3 (4): 1015–1046. doi :10.1137/09076934X. ISSN  1936-4954.
14. Lions, PL; Mercier, B. (1979). "Algoritmos de división para la suma de dos operadores no lineales". Revista SIAM de análisis numérico . 16 (6): 964–979. Bibcode :1979SJNA...16..964L. doi :10.1137/0716071. ISSN  0036-1429. JSTOR  2156649.
15. Beck, Amir; Teboulle, Marc (2009). "Un algoritmo rápido iterativo de umbralización de contracción para problemas inversos lineales". Revista SIAM sobre ciencias de la imagen . 2 (1): 183–202. doi :10.1137/080716542. ISSN  1936-4954. S2CID  3072879.
16. Nestorov, Yu.E. "Un método para resolver un problema de programación convexa con tasa de convergencia O ( 1 k 2 ) {\displaystyle O{\bigl (}{\frac {1}{k^{2}}}{\bigr )}}". Dokl. Akad. Nauk SSSR . 269 (3): 543–547.
17. "Chambolle-Pock · Manopt.jl". docs.juliahub.com . Consultado el 7 de julio de 2023 .
18. "Numerical Tours - Un recorrido numérico por la ciencia de datos". www.numerical-tours.com . Consultado el 7 de julio de 2023 .
19. "Solucionador de Chambolle-Pock: documentación de odl 0.6.1.dev0". odl.readthedocs.io . Consultado el 7 de julio de 2023 .

Lectura adicional

• Boyd, Stephen; Vandenberghe, Lieven (2004). Optimización convexa (PDF) . Cambridge University Press.
• Wright, Stephen (1997). Métodos de puntos interiores primal-dual . Filadelfia, PA: SIAM. ISBN 978-0-89871-382-4.
• Nocedal, Jorge; Stephen Wright (1999). Optimización numérica . Nueva York, NY: Springer. ISBN 978-0-387-98793-4.

Enlaces externos

• EE364b, página de inicio de un curso de Stanford.