Gradiente de la pendiente con restricciones

Question

Con el fin de encontrar los mínimos locales de una función escalar $p(x), x\in \mathbb{R}^3$, sé que podemos utilizar el método de gradiente de la pendiente: $$x_{k+1}=x_k-\alpha_k \nabla_xp(x)$$ donde $\alpha_k$ es el tamaño del paso y $\nabla_xp(x)$ es el gradiente de $p(x)$.

Mi pregunta es: ¿qué pasa si $x$ debe ser limitada en una esfera, es decir, $\|x_k\|=1$? A continuación, debemos encontrar los mínimos locales de $p(x)$ sobre una esfera. Puede gradiente de la pendiente de ser aplicada a la limitada optimizaciones? ¿Alguien puede dar alguna sugerencia? Gracias.

Answer 1

No hay necesidad de pena de los métodos en este caso. El cálculo del gradiente, $g_k=\nabla_xp(x)$, proyecto en el plano tangente, $h_k=g_k-(g_k\cdot x_k)x_k$, y normalizar, $n_k=h_k/|h_k|$. Ahora usted puede utilizar $x_{k+1}=x_k\cos\phi_k + n_k\sin\phi_k $ y unidimensional de la búsqueda para $\phi_k$, al igual que en una sin restricciones de gradiente de búsqueda, y permanece en la esfera local y sigue la dirección de máximo cambio en el nivel métrico sobre la esfera.

Por cierto, esto se puede generalizar para el caso de que usted es la optimización de un conjunto de $n$ vectores bajo la restricción de que son ortonormales. Después de calcular todos los gradientes, proyecto el resultado de la búsqueda de vectores en la tangente de la superficie por orthogonalizing todos los gradientes para todos los vectores, y, a continuación, diagonalize la matriz de productos escalares entre pares de los gradientes de encontrar un sistema de coordenadas en el que los gradientes de pareja con los vectores de la forma $n$ hyperplanes en el que puede girar, mientras que exactamente la satisfacción de las restricciones y seguirá viajando en la dirección de máximo cambio de primer orden; los ángulos por la que gira en el hyperplanes son múltiplos de un único parámetro (con los multiplicadores determinado por los valores propios), así que esta vez es reducido a un uno-dimensional de búsqueda.

[Modificar:] (En respuesta a la sugerencia en los comentarios de uso $x_{k+1}=\frac{x_k-\alpha_kg_k}{\|x_k-\alpha_kg_k\|}$ lugar)

Esto es un poco de desventaja en comparación con el gran círculo de la solución. Que efectivamente está buscando en la misma gran círculo, excepto que en este enfoque sólo puede generar la mitad de ello. Si el óptimo se encuentra en la otra mitad, entonces el valor óptimo del parámetro $\alpha_k$$\pm\infty$, mientras que en mi formulación compacto, el espacio de búsqueda $\phi_k\in[0,2\pi]$ mapas para toda la gran círculo. Además, esta formulación no generalizar para el caso de $n$ ortonormales de vectores, ya que tendrías que realizar una orthonormalization en cada paso.

Cómo determinar $\alpha_k$ no es un (nuevo) problema; en todo el gradiente de búsqueda que usted necesita para encontrar una manera de determinar el óptimo en un uno-dimensional de búsqueda.

Answer 2

La esfera es un ejemplo particular de un (muy agradable) de Riemann colector.

La mayoría de los clásicos de optimización no lineal métodos diseñados para sin restricciones de la optimización de las funciones lisas (como el gradiente de la pendiente que se han mencionado, conjugado lineal de los gradientes, BFGS, Newton, confianza, regiones, etc.) funcionan igual de bien cuando el espacio de búsqueda es una de Riemann colector (suave, un colector con una métrica) en lugar de (clásico) un espacio Euclidiano.

La rama de la optimización de la que se trate el tema se llama de Riemann optimization o optimización en los colectores. Hay un gran libro de referencia sobre el tema que se puede acceder en línea de forma gratuita:

Los algoritmos de optimización en la matriz de colectores, P.-A. Absil, R. Mahony, R. Sepulcro, Princeton university press, 2008.

http://press.princeton.edu/chapters/absil/

Parte de la teoría presentada en ese libro (y algunos más) es implementado en Matlab toolbox llamado Manopt, también disponible gratuitamente en línea:

http://www.manopt.org

La divulgación completa: yo soy un autor de la caja de herramientas.

El tutorial que sucede a empezar con un ejemplo sobre la esfera, que podría ser un útil punto de partida para la cuestión planteada por la OP (suponiendo que Matlab es una opción).

Mejor,

Nicolas

Answer 3

Lo que suelen hacer en teoría (en la práctica) es que para garantizar la convergencia de los algoritmos, para una restricción $g : \mathbb R^n \to \mathbb R$ $f : \mathbb R^n \to \mathbb R$ que deseamos minimizar bajo la condición de $g(x) = 0$, podemos definir $$ \overline f : \mathbb R^n \\overline{\mathbb R} \cup \{ + \infty\}, \qquad \overline f(x) = \begin{cases} + \infty & \text{ if } g(x) = 0 \\ f(x) & \text{ if } g(x) \neq 0. \\ \end{casos} $$ y podemos minimizar $\overline f$ bajo sin restricciones, y claramente $f$ tiene el mismo mínimo de $\overline f$. Ahora lo que quiero hacer en sus métodos numéricos es la definición de un similar $\overline f$ por darle un gran límite superior de los valores de $f$, para asegurarse de que su algoritmo no ir en esas direcciones. Usted podría perder la diferenciabilidad con este método, que puede ser bastante molesto puesto que usted está utilizando la diferenciación en su steepest descent método.

Lo que podrías intentar es usar algo como un mollifier : la idea es muy simple. En lugar de hacer una gran brecha entre la superficie de la esfera y el exterior de la superficie (me refiero a fuera por "no en la superficie, lo que puede significar que dentro o fuera de la bola... de todos modos), que hacen que la brecha suave mediante un golpe de función (consulte la http://en.wikipedia.org/wiki/Bump_function para las imágenes) que es $0$ sobre la esfera y tiene un valor muy alto cuando está a una distancia $\delta$ a partir de la esfera. Si usted llama a un golpe de función $B(x)$ que tiene la propiedad de que $B(x) = 0$ siempre $\| x \| = 1$, y luego definir $$ \overline f(x) = f(x) + B(x). $$ Desde $B = 0$ sobre la esfera, $\overline f$ $f$ tendrá el mismo mínimo con respecto a su restricción. Lo que estamos esperando es que el mínimo que encontrará no será fuera de su restricción. Si esto sucede, realice $\delta$ más pequeño. Esta técnica podría no funcionar si su función se comporta mal, cerca de la frontera de la esfera (es decir, va a valores muy bajos cuando esté cerca de la esfera en una extraña manera, incluso a pesar de que hay un mínimo en la esfera). Si, tiene buen comportamiento a pesar de que no esperamos ningún problema.

Un ejemplo de un golpe función sería la de definir $$ B(x) = \begin{cases} \sin^2 \left( \frac{\| x \| - 1}{\delta} \right) & \text{ if } \left| \| x \| - 1 \right| < \delta \frac {\pi}2 \\ 1 & \text{ if not. } \\ \end{casos} $$ Hacer $\delta$ lo suficientemente pequeña como para multiplicar la protuberancia de la función por cualquier factor muy importante que usted necesita para que sea práctico. No sé si es la mejor idea, pero uno que tengo. Tenga en cuenta que mi "bump función" no es un golpe de la función en el sentido matemático (no de forma compacta compatible), pero es sólo una práctica nombre en este post.

EDIT : las Mejores opciones de $B$, se espera, como uno de mis comentarios mencionan $B(x) = c(\|x\|-1)^2$, $c$ lo suficientemente grande. Supongo J. M. sugerencia (la lectura no-lineal de las referencias de programación) podría dar más indicaciones.

Espero que ayude,

Answer 4

Gracias joriki y Patrick. Francamente, yo prefiero joriki la respuesta de porque considera que la restricción específica. Y Patricio respuesta puede ser aplicado a más problemas genéricos. Para el registro, sucede que tengo otra solución, que creo que es muy similar a la joriki. Agradezco cualquier comentario sobre el método.

La idea es utilizar una rotación de satisfacer la restricción $\|x_{k+1}\|=\|x_k\|$: $$x_{k+1}=R_kx_k$$ donde $R_k$ es una matriz de rotación ($R_k^TR_k=I, \det R_k=1$). Nota de rotación acaba de cambiar la dirección de un vector, pero sin cambiar la longitud. Con el fin de determinar $R_k$, yo uso el eje angular de la representación de una matriz de rotación (aka Rodrigues' fórmula): $$R_k=I+[\omega_k]_{\times}\sin\theta_k+[\omega_k]_{\times}^2(1-\cos\theta_k)$$ donde $\omega_k$ es el eje de rotación, $\theta_k$ es el ángulo de rotación y $[\omega]_{\times}$ es el sesgo de simetría de la matriz de asociación con $\omega$. El eje de rotación puede ser elegido como $$\omega_k=x_k \times n_k$$ y el tamaño del paso, $\theta_k$ podría ser determinado por un unidimensional de búsqueda.

---

EDITAR: $$R_kx_k=x_k+\sin\theta_k (\omega_k\times x_k)+(1-\cos\theta_k)(\omega_k\omega_k^T-I)x_k$$ $$=x_k+\sin\theta_k (\omega_k\times x_k)-(1-\cos\theta_k)x_k =\sin\theta_k (\omega_k\times x_k)+\cos\theta_k x_k$$ En la ecuación anterior, utilizamos los datos: $[\omega]{\times}^2=\omega\omega^T-\|\omega\|I$ $\omega^Tx=0$ (nota:$\omega\times x\ne 0$). Elija $\omega_k\times x_k=n_k$, es decir, $\omega_k=x_k\times n_k$, la ecuación anterior es exactamente el mismo que joriki la respuesta.

---

EDITAR de nuevo: demostrar que es de gradiente de la pendiente. Para $x_{k+1}=x_k+v$, mientras $v^T(-g)>0$, es el gradiente de la pendiente. Y $p(x)$ está garantizado para disminuir dado lo suficientemente pequeño paso. Para nuestro problema, podemos utilizar esto para probar? En nuestra ecuación, $$x_{k+1}=x_k+\underset{v_k}{\underbrace{\sin\theta_k (\omega_k\times x_k)-(1-\cos\theta_k)x_k}}$$ así que necesitamos comprobar si $-v^Tg>0$. Es muy interesante que a diferencia del común de los casos, la evolución de la dirección depende del tamaño de paso de $\theta$. Debemos saber que el gradiente de la pendiente es sólo válido (matemáticamente) para suficientemente pequeño tamaño de paso. Al $\theta$ es muy pequeña, $\sin\theta=\theta$ $1-\cos\theta\approx0$ (de segundo orden). Así tenemos $$v_k\approx\theta_k n_k$$ Obviamente $-g^T\theta_k n_k>0$. Así es el gradiente de la pendiente.

Resumen: con el fin de satisfacer la restricción, tenemos $v_k=\sin\theta_k (\omega_k\times x_k)-(1-\cos\theta_k)x_k$ (complejo). Pero con el fin de demostrar el descenso de la propiedad, tenemos $v_k\approx\theta_k n_k$ (simple y esencial).