Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Algoritmos de programación cuadrática

Definición de programación cuadrática

La programación cuadrática es el problema de encontrar un vector x que minimice una función cuadrática, posiblemente sujeto a restricciones lineales:

minx12xTHx+cTx

tal que A·x ≤ b, Aeq·x = beq, l ≤ x ≤ u.

algoritmo de interior-point-convex quadprog

El algoritmo interior-point-convex realiza los siguientes pasos:

Nota

El algoritmo tiene dos rutas de código. Toma uno cuando la matriz Hessiana H es una matriz ordinaria (completa) de dobles, y toma la otra cuando H es una matriz dispersa. Para obtener información detallada sobre el tipo de datos dispersos, consulte Sparse Matrices (MATLAB). Generalmente, el algoritmo es más rápido para problemas grandes que tienen relativamente pocos términos no nulos cuando se especifica H como sparse. Del mismo modo, el algoritmo es más rápido para problemas pequeños o relativamente densos cuando se especifica H como full.

Resolución de problemas

El algoritmo comienza intentando simplificar el problema eliminando redundancias y simplificando las restricciones. Las tareas realizadas durante el paso de presolver incluyen:

  • Compruebe si las variables tienen límites superiores e inferiores iguales. Si es así, compruebe si hay viabilidad y, a continuación, fije y elimine las variables.

  • Compruebe si cualquier restricción de desigualdad lineal implica sólo una variable. Si es así, compruebe si hay viabilidad y cambie la restricción lineal a un límite.

  • Compruebe si cualquier restricción de igualdad lineal implica sólo una variable. Si es así, compruebe si hay viabilidad y, a continuación, fije y elimine la variable.

  • Compruebe si cualquier matriz de restricción lineal tiene cero filas. Si es así, compruebe si hay viabilidad y elimine las filas.

  • Compruebe si los límites y las restricciones lineales son coherentes.

  • Compruebe si las variables aparecen sólo como términos lineales en la función objetiva y no aparecen en ninguna restricción lineal. Si es así, Compruebe la viabilidad y el límite, y fije las variables en sus límites apropiados.

  • Cambie las restricciones de desigualdad lineal a las restricciones de igualdad lineal mediante la adición de variables de holgura.

Si el algoritmo detecta un problema inviable o ilimitado, detiene y emite un mensaje de salida apropiado.

El algoritmo podría llegar a un único punto factible, que representa la solución.

Si el algoritmo no detecta un problema inviable o ilimitado en el paso de presolver, continúa, si es necesario, con los demás pasos. Al final, el algoritmo reconstruye el problema original, deshaciendo todas las transformaciones de Presolución. Este paso final es el paso possolve.

Para obtener más información, consulte Gould y Toint [63].

Generar punto inicial

El punto inicial x0 para el algoritmo es:

  1. Inicialice x0 en ones(n,1), donde n es el número de filas en H.

  2. Para los componentes que tienen una ub de límite superior y un lbde límite inferior, si un componente de x0 no está estrictamente dentro de los límites, el componente se establece en (ub + lb)/2.

  3. Para los componentes que sólo tienen un límite, modifique el componente si es necesario para mentir estrictamente dentro del límite.

  4. Tome un paso predictor (vea Predictor-corrector), con correcciones menores para la viabilidad, no un paso de predictor completo-corrector. Esto coloca el punto inicial más cerca del Ruta central sin implicar la sobrecarga de un paso de predictor-corrector completo. Para detalles de la ruta central, vea Nocedal y Wright [7], Página 397.

Predictor-corrector

Los algoritmos interior-punto-convexos dispersos y llenos diferencian principalmente en la fase del predecir-corrector. Los algoritmos son similares, pero difieren en algunos detalles. Para la descripción básica del algoritmo, vea Mehrotra [47].

Escueto predictor-corrector.  Al igual que el fmincon algoritmo de punto interior, el escaso algoritmo interior-point-convex intenta encontrar un punto en el que las condiciones Karush-Kuhn-Tucker (KKT) se mantienen. Para el problema de programación cuadrática descrito en Definición de programación cuadrática, estas condiciones son:

Hx+cAeqTyA¯Tz=0A¯xb¯s=0Aeqxbeq=0sizi=0, i=1,2,...,ms0z0.

Aquí

  • A¯ es la matriz de desigualdad lineal extendida que incluye límites escritos como desigualdades lineales. b¯ es el vector de desigualdad lineal correspondiente, incluidos los límites.

  • s es el vector de holguras que convierten las restricciones de desigualdad a las igualdades. s tiene mde la longitud, el número de desigualdades y de límites lineares.

  • z es el vector de los multiplicadores de Lagrange correspondientes a s.

  • y es el vector de los multiplicadores de Lagrange asociados a las restricciones de igualdad.

El algoritmo primero predice un paso de la fórmula de Newton-Raphson, luego calcula un paso corrector. El corrector intenta aplicar mejor la restricción no lineal sizi = 0.

Definiciones para el paso predictor:

  • Rd, el residuo dual:

    rd=Hx+cAeqTyA¯Tz.

  • Req, la restricción de igualdad primaria residual:

    req=Aeqxbeq.

  • Rineq, la restricción de la desigualdad primaria residual, que incluye límites y holguras:

    rineq=A¯xb¯s.

  • Rsz, la complementariedad residual:

    Rsz = Sz.

    S es la matriz diagonal de los términos flojos, z es la matriz de la columna de los multiplicadores de Lagrange.

  • Rc, la complementariedad media:

    rc=sTzm.

En un paso de Newton, los cambios en x, s, y, y z, son dados por:

(H0AeqTA¯TAeq000A¯I000Z0S)(ΔxΔsΔyΔz)=(rdreqrineqrsz).(1)

Sin embargo, un paso completo del neutonio pudo ser inviable, debido a las restricciones de la positividad en s y z. Por lo tanto, quadprog acorta el paso, si es necesario, para mantener positividad.

Además, para mantener una posición "centrada" en el interior, en lugar de tratar de resolver sizi = 0, el algoritmo toma un parámetro positivo σ, y trata de resolver

siZi = σrc.

quadprog sustituye Rsz en la ecuación del paso del neutonio con rsz + ΔsΔz – σrc1, donde 1 es el vector de unos. Además, quadprog reordena las ecuaciones de Newton para obtener un sistema simétrico, más estable numéricamente para el cálculo del paso predictor.

Después de calcular el paso corregido del neutonio, el algoritmo realiza más cálculos para conseguir ambos un paso actual más largo, y prepararse para pasos subsecuentes mejores. Estos cálculos de corrección múltiple pueden mejorar tanto el rendimiento como la robustez. Para obtener más información, consulte Gondzio [4].

Predictor completo-corrector.  El algoritmo de predictor-corrector completo no combina límites en restricciones lineales, por lo que tiene otro conjunto de variables de holgura que corresponden a los límites. El algoritmo desplaza los límites inferiores a cero. Y, si sólo hay un límite en una variable, el algoritmo lo convierte en un límite inferior de cero, al negar la desigualdad de un límite superior.

A¯ es la matriz lineal extendida que incluye las desigualdades lineales y las igualdades lineales. b¯ es el vector de igualdad lineal correspondiente. A¯ también incluye los términos para extender el vector x con las variables flojas s que convierten restricciones de la desigualdad a las restricciones de la igualdad:

A¯x=(Aeq0AI)(x0s),

donde x0 significa el vector x original.

Las condiciones KKT son

Hx+cA¯Tyv+w=0A¯x=b¯x+t=uvixi=0, i=1,2,...,mwiti=0, i=1,2,...,nx,v,w,t0.(2)

Para encontrar la solución x, las variables flojas y las variables duales a Ecuación 2, el algoritmo considera básicamente un paso de Newton-Raphson:

(HA¯T0IIA¯0000I0I00V00X000W0T)(ΔxΔyΔtΔvΔw)=(Hx+cA¯Tyv+wA¯xb¯uxtVXWT)=(rdrprubrvxrwt),(3)

donde X, V, Wy T son matrices diagonales correspondientes a los vectores x, v, wy t respectivamente. Los vectores residuales en el extremo derecho de la ecuación son:

  • Rd, el doble residual

  • Rp, el primario residual

  • Rub, el límite superior residual

  • Rvx, la complementariedad residual de menor límite

  • Rwt, el residuo de complementariedad del límite superior

El algoritmo resuelve Ecuación 3 primero convirtiéndolo a la forma simétrica de la matriz

(DA¯TA¯0)(ΔxΔy)=(Rrp),(4)

Donde

D=H+X1V+T1WR=rdX1rvx+T1rwt+T1Wrub.

Todos los inversos de la matriz en las definiciones de D y R son simples de calcular porque las matrices son diagonales.

Para derivar Ecuación 4 de Ecuación 3, observe que la segunda fila de Ecuación 4 es la misma que la segunda fila de matriz de Ecuación 3. La primera fila de Ecuación 4 viene de resolver las dos últimas filas de Ecuación 3 para δv y δw, y entonces resolviendo para δt.

Para resolver Ecuación 4, el algoritmo sigue los elementos esenciales de Altman y Gondzio [1]. El algoritmo resuelve el sistema simétrico por una descomposición del LDL. Como señalaron autores como Vanderbei y Carpenter [2], esta descomposición es numéricamente estable sin pivotes, por lo que puede ser rápido.

Después de calcular el paso corregido del neutonio, el algoritmo realiza más cálculos para conseguir ambos un paso actual más largo, y prepararse para pasos subsecuentes mejores. Estos cálculos de corrección múltiple pueden mejorar tanto el rendimiento como la robustez. Para obtener más información, consulte Gondzio [4].

El algoritmo de quadprog predictor-corrector completo es en gran medida el mismo que en el algoritmo de 'interior-point' linprog , pero también incluye términos cuadráticos. Véase Predictor-corrector.

Referencias

[1] Altman, Anna and J. Gondzio. Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization. Optimization Methods and Software, 1999. Available for download here.

[2] Vanderbei, R. J. and T. J. Carpenter. Symmetric indefinite systems for interior point methods. Mathematical Programming 58, 1993. pp. 1–32. Available for download here.

Condiciones de parada

El algoritmo predictor-corrector itera hasta que alcanza un punto que es factible (satisface las restricciones dentro de las tolerancias) y donde los tamaños de los pasos relativos son pequeños. Específicamente, defina

ρ=max(1,H,A¯,Aeq,c,b¯,beq)

El algoritmo se detiene cuando se cumplen todas estas condiciones:

rp1+rub1ρTolConrdρTolFunrcTolFun,

Donde

rc=maxi(min(|xivi|,|xi|,|vi|),min(|tiwi|,|ti|,|wi|)).

Rc esencialmente mide el tamaño de los residuales de complementariedad xv y tw, que son cada vectores de ceros en una solución.

Detección de inviabilidad

quadprog calcula un función del mérito φ en cada iteración. La función Merit es una medida de viabilidad. quadprog se detiene si la función de mérito crece demasiado. En este caso, quadprog declara que el problema es inviable.

La función Merit está relacionada con las condiciones KKT para el problema — ver Predictor-corrector. Utilice las siguientes definiciones:

ρ=max(1,H,A¯,Aeq,c,b¯,beq)req=Aeqxbeqrineq=A¯xb¯srd=Hx+c+AeqTλeq+A¯Tλ¯ineqg=xTHx+fTxb¯Tλ¯ineqbeqTλeq.

La notación A¯ Y b¯ significa los coeficientes de desigualdad lineal, aumentados con términos para representar límites para el algoritmo disperso. La notación λ¯ineq también representa los multiplicadores de Lagrange para las restricciones de desigualdad lineal, incluidas las restricciones enlazadas. Esto se llamaba z en Predictor-corrector, y λeq se llamaba y.

La función Merit φ es

1ρ(max(req,rineq,rd)+g).

Si esta función de mérito se vuelve demasiado grande, quadprog declara que el problema es inviable y se detiene con la bandera de salida -2.

algoritmo de trust-region-reflective quadprog

Muchos de los métodos utilizados en Optimization Toolbox™ Solvers se basan en regiones de confianza, un concepto simple pero poderoso en la optimización.

Para entender el enfoque de la región de confianza para la optimización, considere el problema de minimización no restringida, minimice f(x), donde la función toma argumentos vectoriales y devuelve escalares. Suponga que está en un punto x en n-Space y desea mejorar, es decir, pasar a un punto con un valor de función inferior. La idea básica es aproximar f con una función más simple q, que refleja razonablemente el comportamiento de la función f en un barrio N alrededor del punto x. Este vecindario es la región de confianza. Un paso de ensayo s se computa reduciendo al mínimo (o aproximadamente reduciendo al mínimo) el excedente N. Este es el subproblema de la región de confianza,

mins{q(s), sN}.(5)

El punto actual se actualiza para ser x + s si f(x + s) < f(x); de lo contrario, el punto actual permanece inalterado y N, la región de confianza, se encoge y se repite el cálculo del paso de prueba.

Las preguntas clave en la definición de un enfoque específico de la región de confianza para minimizar la f(x) son cómo elegir y calcular la aproximación q (definida en el punto actual x), cómo elegir y modificar la región de confianza N, y Cómo resolver con precisión el subproblema confianza-región. Esta sección se centra en el problema no restringido. Las secciones posteriores discuten complicaciones adicionales debido a la presencia de restricciones en las variables.

En el método estándar de la región de confianza ([48]), la aproximación cuadrática q se define por los dos primeros términos de la aproximación de Taylor a F en x; la vecindad N es generalmente esférica o elipsoide en forma. Matemáticamente el subproblema de la confianza-región se indica típicamente

min{12sTHs+sTg  such that  DsΔ},(6)

donde g es el gradiente de f en el punto actual x, H es la matriz del hessian (la matriz simétrica de los segundos derivados), D es una matriz de escalamiento diagonal, Δ es un escalar positivo, y ∥. ∥ es el 2-Norm. Existen buenos algoritmos para resolver Ecuación 6 (ver [48]); Estos algoritmos típicamente implican el cómputo de un eigensystem completo y un proceso de Newton aplicado a la ecuación secular

1Δ1s=0.

Tales algoritmos proporcionan una solución exacta a Ecuación 6. Sin embargo, requieren tiempo proporcional a varios factorizaciones de H. Por lo tanto, para problemas a gran escala es necesario un enfoque diferente. Se han propuesto varias estrategias de aproximación y heurística, basadas en Ecuación 6, en la bibliografía ([42] y [50]). El enfoque de aproximación seguido en Optimization Toolbox Solvers es restringir el subproblema de la región de confianza a un S subespacial de dos dimensiones ([39] y [42]). Una vez que el subespacio S se ha calculado, el trabajo para solucionar Ecuación 6 es trivial incluso si se necesita la información completa del valor propio/vector (puesto que en el subespacio, el problema es solamente de dos dimensiones). El trabajo dominante se ha trasladado ahora a la determinación del subespacio.

El subespacio S de dos dimensiones se determina con la ayuda de un proceso de gradiente conjugado precondicionado descrito a continuación. El solucionador define S como el espacio lineal atravesado por s1 y s2, donde s1 está en la dirección del gradiente g, y s2 es una aproximación Newton Dirección, es decir, una solución para

Hs2=g,(7)

o una dirección de curvatura negativa,

s2THs2<0.(8)

La filosofía detrás de esta opción de S es forzar la convergencia global (vía la dirección más escarpada de la pendiente o la dirección negativa de la curvatura) y alcanzar la convergencia local rápida (vía el paso del neutonio, cuando existe).

Un bosquejo de la minimización no restringida usando ideas de la confianza-región es ahora fácil de dar:

  1. Formule el subproblema bidimensional de la confianza-región.

  2. Resuelve Ecuación 6 para determinar el paso de prueba s.

  3. Si f(x + s) < f(x)Entonces x = x + s.

  4. Ajuste Δ.

Estos cuatro pasos se repiten hasta la convergencia. La dimensión de la confianza-región Δ se ajusta según reglas estándares. En particular, se reduce si no se acepta el paso de prueba, es decir, f(x + s) ≥ f(x). Vea [46] y [49] para una discusión de este aspecto.

los solucionadores de Optimization Toolbox tratan algunos casos especiales importantes de f con funciones especializadas: menos cuadrados no lineales, funciones cuadráticas, y mínimos-cuadrados lineares. Sin embargo, las ideas algorítmicas subyacentes son las mismas que para el caso general. Estos casos especiales se discuten en secciones más últimas.

El método de región de confianza subespacial se utiliza para determinar una dirección de búsqueda. Sin embargo, en vez de restringir el paso a (posiblemente) un paso de la reflexión, como en el caso no lineal de la minimización, un tramos la búsqueda de línea reflexiva se realiza en cada iteración. Consulte [45] para obtener más información sobre la búsqueda de líneas.

Método de gradiente de conjugado precondicionado

Una manera popular de resolver los sistemas definidos positivos simétricos grandes de ecuaciones lineares Hp = –g es el método de gradientes de conjugado precondicionado (PCG). Este enfoque iterativo requiere la capacidad de calcular los productos vectoriales matriciales de la forma H·v donde v es un vector arbitrario. La matriz definida positiva simétrica M es un precondicionador para H. Es decir M = C2Donde C–1HC–1 es una matriz bien acondicionada o una matriz con valores propios agrupados.

En un contexto de minimización, puede asumir que la matriz de hessian H es simétrica. Sin embargo, H se garantiza para ser positivo definido solamente en el vecindario de un minimizador fuerte. El algoritmo PCG sale cuando se encuentra una dirección de curvatura negativa (o cero), es decir, dTHd ≤ 0. La dirección de la salida de PCG, p, es una dirección de la curvatura negativa o un aproximado (tol controla cómo aproximado) solución al sistema del neutonio Hp = –g. En cualquiera de los casos, p se utiliza para ayudar a definir el subespacio bidimensional utilizado en el enfoque de la región de confianza discutido en Métodos de la confianza-región para la minimización no lineal.

Restricciones de igualdad lineal

Las restricciones lineales complican la situación descrita para la minimización no restringida. Sin embargo, las ideas subyacentes descritas anteriormente se pueden llevar a cabo de una manera limpia y eficiente. Los métodos de la región de confianza en los solucionadores de Optimization Toolbox generan iteraciones estrictamente factibles.

La igualdad lineal general limitada problema de minimización se puede escribir

min{f(x)  such that  Ax=b},(9)

donde A es una matriz m-por-n (m ≤ n). Algunos Optimization Toolbox solucionadores preprocesan A para eliminar dependencias lineales estrictas utilizando una técnica basada en la factorización LU de UnT [46]. Aquí A se asume para ser de mespeso.

El método utilizado para resolver Ecuación 9 difiere del enfoque sin restricciones en dos formas significativas. Primero, un punto factible inicial x0 se computa, utilizando un paso de mínimos cuadrados dispersos, de modo que Ax0 = b. En segundo lugar, el algoritmo PCG es substituido por gradientes conjugados pre-condicionados reducidos (RPCG), vea [46], para calcular un paso reducido aproximado del neutonio (o una dirección de la curvatura negativa en el espacio nulo de A). El paso de álgebra lineal clave consiste en resolver sistemas de la forma

[CA˜TA˜0][st]=[r0],(10)

Donde A˜ aproxima A (los pequeños no ceros de A se fijan a la fila proporcionada cero no se pierde) y C es una aproximación positiva-definida simétrica escasa a H, es decir, C = H. Vea [46] para más detalles.

Restricciones de caja

El problema de la caja restringida es de la forma

min{f(x)  such that  lxu},(11)

donde l es un vector de límites inferiores, y u es un vector de límites superiores. Algunos (o todos) de los componentes de l pueden ser iguales a – ∞ y algunos (o todos) de los componentes de u pueden ser iguales a ∞. El método genera una secuencia de puntos estrictamente factibles. Se utilizan dos técnicas para mantener la viabilidad y lograr un comportamiento de convergencia robusto. En primer lugar, un paso de Newton modificado a escala sustituye el paso de Newton sin restricciones (para definir el subespacio Sde dos dimensiones). Segundo las reflexiones se utilizan para aumentar el tamaño del escalón.

El paso modificado escalado del neutonio se presenta de examinar las condiciones necesarias de Kuhn-Tucker para Ecuación 11,

(D(x))2g=0,(12)

Donde

D(x)=diag(|vk|1/2),

y el vector v(x) se define a continuación, para cada 1 ≤ i ≤ n:

  • Si gi < 0 Y ui < ∞ entonces vi = xi – ui

  • Si gi ≥ 0 Y li > –∞ entonces vi = xi – li

  • Si gi < 0 Y ui = ∞ entonces vi = –1

  • Si gi ≥ 0 Y li = –∞ entonces vi = 1

El sistema no lineal Ecuación 12 no es diferenciable por todas partes. La no diferenciabilidad ocurre cuando vi = 0. Puede evitar estos puntos manteniendo una viabilidad estricta, es decir, restringiendo l < x < u.

El escalón modificado de Newton sk para el sistema no lineal de ecuaciones dadas por Ecuación 12 se define como la solución al sistema lineal

M^DsN=g^(13)

en la iteración kTH, donde

g^=D1g=diag(|v|1/2)g,(14)

Y

M^=D1HD1+diag(g)Jv.(15)

Aquí Jv desempeña el papel del jacobiano de | v|. Cada componente diagonal de la matriz diagonal Jv es igual a 0, – 1 o 1. Si todos los componentes de l y u son finitos, Jv = diag(sign(g)). En un punto donde gi = 0, Vi puede que no sea diferenciable. Jiiv=0 se define en un punto. La no diferenciabilidad de este tipo no es motivo de preocupación porque, para dicho componente, no es significativo el valor Vi Toma. Más, |Vi| seguirá siendo discontinua en este punto, pero la función |vigi es continuo.

Segundo las reflexiones se utilizan para aumentar el tamaño del escalón. Un paso de reflexión (Single) se define como sigue. Dado un paso p que intersecta una restricción enlazada, considere la primera restricción enlazada cruzada por p; Supongamos que es la restricción iTH Bound (ya sea el iTH Upper o iTH inferior Bound). A continuación, el paso de reflexión pR = p excepto en el componente iTH, donde pRi = –pi.