Optimización no lineal con restricciones

Question

Deje $w$ ser algunos de vectores. Quiero resolver: $$ w=\arg\min \sum f_i^2(w) \\s.t. w\geq 0, \sum w_i = 1$$ Que puede ser formulada de la siguiente manera: $$ w=\arg\min ||{\bf f}(w)||^2 \\ s.t. \quad H(w)\geq 0, G(w)= 0 $$ donde $H(w)=w,G(w)={\bf 1}^T w-1=\sum w_i-1$.

Me parece que este es un problema estándar. Las restricciones del problema puede ser resuelto usando el de Gauss-Newton (GN) método. Sin embargo, me parece que no puede encontrar la manera de resolver la limitación de problema. Que método debo utilizar? puede GN ser adaptado a estas limitaciones?

Gracias.

EDITAR:

una. Voy a añadir una fórmula para $f_i(w)$: $$ f_i(w) = w^T A_i w - z_i$$ where $A_i$ is a symmetric matrix and $z_i$ son conocidos escalares.

b. Tratando de encontrar mi respuesta, me encontré con el de Levenberg–Marquardt algoritmo, que aparece casi adecuado aquí, pero tengo una matriz de $A_i$ y no un escalar que varía con $i$.

c. He pensado en hacer lo siguiente: puesto que en cada iteración de la solución de Gauss-Newton método que resolver un lineal LS, yo podría resolver en lugar restringido LS con mis limitaciones. Es que un camino a seguir?

Gracias de nuevo.

Answer 1

Como se señaló en los comentarios, esta respuesta no dar una solución exacta.

El uso de una variable adicional $t_i$, que se utiliza para modelar $f_i(w)$, el problema puede ser representado como un cuadráticamente limitada problema de optimización cuadrática (QCQP): $$ \begin{align*} \min \quad & t_i^2 \\ \mathrm{s.t.} \quad & t_i \geq w^T A_i w - z_i \\ & \sum w_i = 1 \\ & w \geq 0 \end{align*} $$ Usted puede preguntarse por qué la primera restricción se utiliza "$\geq$" en lugar de "$=$". La razón es que asegura que la convexidad (dadas dos soluciones factibles, el segmento de línea entre esas soluciones es también posible). Sin embargo, como se señaló en los comentarios, la formulación no es equivalente.

La ventaja de la escritura de este problema como un QCQP es que hay una gran cantidad de algoritmos para resolverlo. Me centraré en los métodos de punto interior, que consideran que la (a nivel global) de la solución óptima en el 40-60 iteraciones independiente del tamaño de su problema. Matlab tiene quadprog. Software comercial como CPLEX, Mosek, Gurobi ha enlaces para muchos idiomas, incluyendo Matlab y Python. QCQP es un caso especial de optimización convexa, por lo que en Python se podría utilizar scipy.optimizar.

Answer 2

Su problema de optimización tiene la forma \begin{align} \tag{%#%#%}\operatorname{minimize}_w &\quad F(w) \\ \text{subject to} & \quad w \in C \end{align} donde $$ F(w) = \sum_i f_i(w)^2 $$ y $\spadesuit$ es la probabilidad de simplex (que es un sistema cerrado, el conjunto convexo).

Usted puede resolver sus problemas de la forma ($C$) con la proyección de el método del gradiente: $$ w^{k+1} = \Pi_C(w^k - t\nabla F(w^k)). $$ Aquí $\spadesuit$ es la proyección de $\Pi_C(y)$ a $y$. Proyectando la probabilidad simplex es un muy bien estudiado operación y se puede realizar de manera muy eficiente con los algoritmos existentes. (Se tarda sólo un par de líneas de Matlab. Ver algoritmo 1 aquí, por ejemplo: https://arxiv.org/abs/1309.1541 .)

El proyectado el método del gradiente tiene un tamaño de paso de restricción, así que asegúrese de seleccionar el tamaño de paso de $C$ lo suficientemente pequeño para que el tamaño del paso de restricción se satisface. También hay una línea de búsqueda de procedimiento que se puede utilizar para seleccionar tamaños de paso(que creo que es probablemente una buena idea).

También puede probar a usar una versión acelerada de la proyectada el método del gradiente, tales como la FISTA, aunque no estoy seguro de cómo van a realizar para nonconvex problemas.

Edit: El proyectado el método del gradiente es un caso especial de la parte proximal del método de gradiente, y una línea de búsqueda de la versión de la parte proximal del método de gradiente se describe en las ecuaciones de 1.26 y 1.27 (p. 18) de la ponencia "Gradiente Basado en Algoritmos con las Aplicaciones de la Señal de Problemas de Recuperación" de Beck y Teboulle. Esta línea de búsqueda de la versión de la parte proximal del método de gradiente también se describe, con un poco de notación diferente, en la diapositiva 6-20 de Vandenberghe del c 236 notas: http://www.seas.ucla.edu/~vandenbe/C 236/conferencias/proxgrad.pdf.

Para mayor comodidad, voy a escribir la parte proximal del método de gradiente con la línea de búsqueda en detalle, usando la notación que creo que es un poco más clara que la notación en las anteriores referencias. Este algoritmo minimiza $t$ donde $f(x) + g(x)$ se supone que para ser diferenciable (con un Lipschitz continua gradiente) y $f$ es un convexo, semicontinua inferior de la función. La porción proximal del operador de $g$ con el parámetro $g$ se denota por a $t > 0$. Si usted no sabe lo que la porción proximal del operador, que está bien, simplemente reemplace$\operatorname{prox}_{tg}$$\operatorname{prox}_{tg}$, y se han proyectado el método del gradiente. Aquí es el algoritmo: \begin{align} &\text{Initialize } x_0 \text{ and select } t_0 > 0, 0 < \alpha < 1, \beta > 1\\ &\textbf{for } k = 1,2,\ldots \textbf{do:}\\ & \quad t := \beta \, t_{k-1} \\ &\quad \textbf{repeat:} \\ &\qquad \quad x := \operatorname{prox}_{tg}(x_{k-1} - t \nabla f(x_{k-1})) \\ &\qquad \quad \textbf{break if } f(x) \leq f(x_{k-1}) + \langle \nabla f(x_{k-1}),x - x_{k-1} \rangle + \frac{1}{2t} \| x - x_{k-1} \|_2^2 \\ &\qquad \quad t := \alpha t \\ & \quad t_k := t \\ & \quad x_k := x\\ &\textbf{end for} \end{align}

Yo normalmente seleccione $\Pi_C$$\alpha = .5$. Creo $\beta = 1.25$ no suele ser una mala elección. Usted podría inicializar $t_0 = 1$ a ser su mejor disposición estimación del valor óptimo de $x_0$, o tal vez sólo inicializar $x$ al azar o a todos los ceros.

Aquí hay algunas código de python para la porción proximal del método de gradiente con la línea de búsqueda:

def proxGrad(gradf,proxg,evalf,evalg,x0,params):
    # minimize f(x) + g(x). f is smooth, g has easy prox-operator

    maxIter = params['maxIter']
    t = params['stepSize'] # initial step size
    showTrigger = params['showTrigger']
    reduceFactor = .5
    increaseFactor = 1.25
    costs = np.zeros((maxIter,1))

    xk = x0
    for k in np.arange(maxIter):

        (gradf_xk, fxk) = gradf(xk)

        gxk = evalg(xk)
        costs[k] = fxk + gxk
        if k % showTrigger == 0:
            print "Iteration: " + str(k) + "    cost: " + str(costs[k])

        t = t*increaseFactor
        acceptFlag = False
        while acceptFlag == False:
            xkp1 = proxg(xk - t*gradf_xk, t)
            fxkp1 = evalf(xkp1)
            diff = xkp1 - xk
            if fxkp1 <= fxk + np.vdot(gradf_xk,diff) + (.5/t)*np.sum(diff**2):
                acceptFlag = True
            else:
                t = t*reduceFactor
                print "Reducing t. t is now: " + str(t)

        xk = xkp1

    return (xkp1,costs)

Aquí un poco más de info acerca de la relación entre la porción proximal del método del gradiente y el proyectado el método del gradiente. El proyectado el método del gradiente para minimizar $x_0$ sujeto a la restricción de que $f(x)$, es solo un caso especial de la parte proximal del método de gradiente donde $x \in C$ es convexo de la función de indicador de $g$: $$ g(x) = \begin{cases} 0 & \quad \text{if } x \in C \\ \infty & \quad \text{otherwise.} \end{casos} $$ Con esta elección de $C$, la porción proximal del operador de $g$ está dado por la fórmula $$ \operatorname{prox}_{tg}(x) = \Pi_C(x) = \text{proyección de $g$ a $x$}. $$ (Estoy asumiendo que $C$ es un conjunto convexo cerrado.)