Optimización numérica en presencia de un algoritmo rápido para algunos ejes

Question

Supongamos que buscamos $(x,y)$ para minimizar $f(x,y)$ ( $f$ continua, diferenciable, pero no necesariamente convexa) donde $x\in \mathbb{R}^m$ y $y\in \mathbb{R}^n$ . Podemos utilizar cualquier método estándar que funcione en $\mathbb{R}^{m+n}$ (por ejemplo, estoy trabajando en un contexto en el que el algoritmo de Levenberg-Marquardt funciona bien).

Pero ahora supongamos que tenemos un algoritmo especialmente rápido para encontrar el mínimo de $f(x,y)$ para cualquier $x$ (tal vez una solución explícita, o una iteración que se sabe que se comporta bien). ¿Cómo podemos aprovechar esto para encontrar el mínimo sobre todos los $x$ y $y$ ?

Una forma (normalmente mala) podría ser alternar la minimización $x$ y $y$ .

Una cosa obvia que se me ocurre es resolver explícitamente $g(x)=\underset{y}{\operatorname{argmin}} f(x,y)$ y minimizar $f(x,g(x))$ . Una de las razones para no hacerlo es que $g$ podría tener derivadas particularmente complicadas (por ejemplo, podría implicar un cálculo de SVD, aunque posiblemente la diferenciación automática seguiría funcionando bien).

¿Hay alguna otra forma de aprovechar el hecho de que podemos minimizar la rapidez a lo largo de $m$ de los ejes? Esta situación me ha surgido tantas veces que estoy convencido de que debe haber alguna teoría publicada.

Answer 1

Creo que el Complemento de Schur es el camino correcto:

Supongamos que queremos minimizar una función $F$ donde son aplicables los métodos de tipo Gauss-Newton como el LM (Levenberg Marquardt). Así, tenemos una función $F(\mathbf x) = f(\mathbf x)^\top f(\mathbf x)$ con $f( \bar{\mathbf{x}} ) \approx \mathbf 0$ como mínimo $\bar{\mathbf{x}}$ .

En la LM estándar, resolveríamos repetidamente la ecuación normal $$\mathtt H\delta = -\frac{\partial f(\mathbf x)}{\partial \mathbf x}^\top f(\mathbf x)$$ con $\mathtt H=\left(\frac{\partial f(\mathbf x)}{\partial \mathbf x}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x} + \lambda \mathtt I\right)$ y actualizar nuestra estimación: $$\mathbf{x}' = \mathbf{x}+\delta.$$

Supongamos además que $\mathbf x = \begin{bmatrix}\mathbf x_1\\\mathbf x_2\end{bmatrix}$ y podemos resolver $\mathbf{x_2'} = \arg\min_{\mathbf x_2}F(\mathbf x_1,\mathbf x_2)$ para cualquier $\mathbf x_1$ muy eficazmente.

Ahora, el jacobiano $\frac{\partial f(\mathbf x)}{\partial \mathbf x}= \left(\frac{\partial f(\mathbf x)}{\partial \mathbf x_1}, \frac{\partial f(\mathbf x)}{\partial \mathbf x_2} \right)$ y el hessiano aumentado $\mathtt H = \begin{bmatrix}\mathtt A & \mathtt B\\ \mathtt B^\top & \mathtt C\end{bmatrix}$ con $\mathtt A = \frac{\partial f(\mathbf x)}{\partial \mathbf x_1}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_1}+\lambda \mathtt I$ , $\mathtt B = \frac{\partial f(\mathbf x)}{\partial \mathbf x_1}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}$ y $\mathtt C = \frac{\partial f(\mathbf x)}{\partial \mathbf x_2}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}+\lambda \mathtt I$ . Así, obtenemos la siguiente ecuación normal:

$$ \begin{bmatrix}\mathtt A & \mathtt B\\ \mathtt B^\top & \mathtt C\end{bmatrix}\cdot \begin{bmatrix}\delta_1\\\delta_2\end{bmatrix} = \begin{bmatrix}\mathbf b_1\\\mathbf b_2 \end{bmatrix}$$ con $\mathbf b_1 = -\frac{\partial f(\mathbf x)}{\partial \mathbf x_1}f(\mathbf x)$ y $\mathbf b_2 = -\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}f(\mathbf x)$ .

Utilizando el Complemento de Schur podemos resolver para $\delta_1$ :

$$(\mathtt A - \mathtt B\mathtt C^{-1} \mathtt B^\top)\delta_1 = \mathbf b_1-\mathtt B\mathtt C^{-1} \mathbf b_2$$

Después, podemos minimizar $\mathbf x_2$ para un fijo $\mathbf x_1' =\mathbf x_1+\delta_1$ .

---

Así, podemos aplicar el siguiente esquema LM:

$\text{repeat until convergence}$

$\quad\quad\mathbf b_1 := -\frac{\partial f(\mathbf x)}{\partial \mathbf x_1}f(\mathbf x)$

$\quad\quad\mathbf b_2 := -\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}f(\mathbf x)$

$\quad\quad\mathtt A := \frac{\partial f(\mathbf x)}{\partial \mathbf x_1}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_1}+\lambda \mathtt I$

$\quad\quad\mathtt B := \frac{\partial f(\mathbf x)}{\partial \mathbf x_1}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}$

$\quad\quad\mathtt C := \frac{\partial f(\mathbf x)}{\partial \mathbf x_2}^\top\frac{\partial f(\mathbf x)}{\partial \mathbf x_2}+\lambda \mathtt I$

$\quad\quad\delta_1 := (\mathtt A - \mathtt B\mathtt C^{-1} \mathtt B^\top)^{-1}( \mathbf b_1-\mathtt B\mathtt C^{-1} \mathbf b_2)$

$\quad\quad \mathbf x'_1 := \mathbf x_1 + \delta_1 $

$\quad\quad\mathbf x_2' := \arg\min_{\mathbf x_2}F(\mathbf x_1',\mathbf x_2)$

$\quad\quad \text{if } F(\mathbf x_1',\mathbf x_2') \ll F(\mathbf x_1,\mathbf x_2)$

$\quad\quad\quad\quad \mathbf x_1 := \mathbf x_1'$

$\quad\quad\quad\quad \mathbf x_2 := \mathbf x_2'$

$\quad\quad\quad\quad \text{decrease } \lambda$

$\quad\quad \text{else}$

$\quad\quad\quad\quad \text{increase } \lambda$

$\quad\quad \text{end}$

$\text{end}$

---

Actualización: Como alternativa/equivalente, podemos simplemente resolver toda la ecuación normal, (ignorar $\delta_2$ ) y minimizar $\mathbf x_2$ después:

$\text{repeat until convergence}$

$\quad\quad$ resolver $\begin{bmatrix}\mathtt A & \mathtt B\\ \mathtt B^\top & \mathtt C\end{bmatrix}\cdot \begin{bmatrix}\delta_1\\\delta_2\end{bmatrix} = \begin{bmatrix}\mathbf b_1\\\mathbf b_2 \end{bmatrix}$

$\quad\quad \mathbf x'_1 := \mathbf x_1 + \delta_1 $

$\quad\quad\mathbf x_2' := \arg\min_{\mathbf x_2}F(\mathbf x_1',\mathbf x_2)$

$\quad\quad \text{if } F(\mathbf x_1',\mathbf x_2') \ll F(\mathbf x_1,\mathbf x_2)$

$\quad\quad\quad\quad \mathbf x_1 := \mathbf x_1'$

$\quad\quad\quad\quad \mathbf x_2 := \mathbf x_2'$

$\quad\quad\quad\quad \text{decrease } \lambda$

$\quad\quad \text{else}$

$\quad\quad\quad\quad \text{increase } \lambda$

$\quad\quad \text{end}$

$\text{end}$

Cuál de las dos formulaciones es más eficiente depende del problema específico (estructura de escasez y dimensionalidad de $\mathtt{ A, B, C}$ etc.).

---

Tenga en cuenta que este esquema es mucho mejor que hacer la alternancia ya que la actualización $\delta_1$ toma la siguiente actualización sobre $\mathbf x_2$ en cuenta.

Actualización2: Por supuesto, no he inventado este esquema. No conozco ningún libro de texto que lo mencione, pero he visto algo parecido en el contexto del ajuste de los paquetes: http://research.microsoft.com/pubs/131806/Jeong-CVPR10.pdf (iteración de puntos incrustados).