Introducción a la Programación Convexa Discplinada (DCP). Algunas aplicaciones y su resolución informática

Grado en Matemáticas TRABAJO FIN DE GRADO Introducción a la Programación Convexa Discplinada (DCP). Algunas aplicaciones y su resolución informática. María Cabeza Cobano Sevilla, Junio de 2019 Índice general Resumen....................................... iii Abstract....................................... iv ÍndicedeFiguras.................................. v ÍndicedeCuadros.................................. vii 1. Introducción 1 1.1. Programación Matemática . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1. Resolver el problema . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. Programación convexa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2. Análisis convexo 5 2.1. Conjuntosconvexos.............................. 5 2.1.1. Operaciones que preservan la convexidad . . . . . . . . . . . . . . 5 2.2. Funcionesconvexas.............................. 7 2.2.1. Definiciones.............................. 7 2.2.2. Propiedades básicas . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2.3. Operaciones que preservan la convexidad . . . . . . . . . . . . . . 10 2.3. Problemas de optimización convexa . . . . . . . . . . . . . . . . . . . . . 12 2.3.1. Problemalineal............................ 13 2.3.2. Problema de mínimos cuadrados . . . . . . . . . . . . . . . . . . . 13 2.3.3. Problema cuadrático . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.4. Problema de cono de segundo orden . . . . . . . . . . . . . . . . . 13 2.3.5. Problemas de programación semidefinida . . . . . . . . . . . . . . 14 2.3.6. Problemas de programación geométrica . . . . . . . . . . . . . . . 14 2.4. Algoritmosnuméricos............................. 14 2.4.1. Formasestándar ........................... 14 2.4.2. Código personalizado . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5. Ventajas de los problemas convexos . . . . . . . . . . . . . . . . . . . . . 17 2.6. Convexidad y diferenciabilidad . . . . . . . . . . . . . . . . . . . . . . . . 17 2.7. Verificación de convexidad . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8. Aplicaciones.................................. 19 3. Programación Convexa Disciplinada 21 3.1. Conjunto de reglas de convexidad . . . . . . . . . . . . . . . . . . . . . . 22 3.1.1. Reglas de nivel superior . . . . . . . . . . . . . . . . . . . . . . . 22 3.1.2. Reglas libres de productos . . . . . . . . . . . . . . . . . . . . . . 23 3.1.3. Reglasdesigno ............................ 26 3.1.4. Reglas de composición . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2. Verificación .................................. 29 i 3.3. Crear programas convexos disciplinados . . . . . . . . . . . . . . . . . . . . 31 3.4. Implementación de las funciones básicas . . . . . . . . . . . . . . . . . . . 31 3.4.1. Implementaciones gráficas . . . . . . . . . . . . . . . . . . . . . . 32 3.4.2. Uso de las implementaciones gráficas . . . . . . . . . . . . . . . . 32 4. CVXR 35 4.1. Expresiones .................................. 36 4.2. Funciones ................................... 38 4.2.1. Funciones escalares . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2.2. Funciones element-wise . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.2.3. Funciones vector / matriz . . . . . . . . . . . . . . . . . . . . . . 42 4.3. Construirunproblema............................ 43 4.3.1. Problema vectorial . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.4. Implementación................................ 45 4.5. Mínimoscuadrados.............................. 47 4.6. Aplicaciones.................................. 49 4.6.1. Problema de la Catenaria . . . . . . . . . . . . . . . . . . . . . . 49 4.6.2. Regresión logística . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Conclusiones 59 Bibliografía 61 ii Resumen En este trabajo de fin de grado se explica la Programación Convexa Disciplinada, una metodología para trabajar con problemas de optimización convexa. Su objetivo es simplificar el proceso de verificación de la convexidad de un problema, que en muchos casos es un trabajo intratable. Se basa en una librería de funciones, a partir de las cuales se construyen los problemas y un conjunto de normas impuesto sobre estas funciones, que establecen como combinarlas. No se restringe a un tipo determinado de forma estándar como hacen otras metodologías y no pierde generalidad ya que la librería en la que se basa es extensible. Esta metodología se puede llevar a la práctica a través del paquete CVXR de R. Permite formular los problemas de una forma muy natural e intuitiva. En CVXR es posible formular y resolver problemas para los que no existe un código personalizado. Dado un problema, CVXR emplea la programación convexa disciplinada previamente a resolver, para probar que sea convexo y por tanto garantizar que la solución óptima es global. iii Abstract Convex Disciplined Programming is a methodology for working with convex optimization problems. The objective is to simplify the convexity verification, because this is often an intractable task. It is based on a library of functions, that we use for construct the problems, and a set of rules that constitute a set of sufficient conditions for guarantee convexity. DCP is not restricted to a standard way as other methodologies and it does not lose generality since the library is extensible. This methodology can be put into practice through the package CVXR of R. It allows the user to formulate convex optimization problems in a natural and intuitive mathematical syntax. In CVXR it is possible to formulate and solve problems for which there is no personalized code. Given a problem, CVXR applies Disciplined Convex Programming to verify the problem’s convexity and therefore guarantee that the optimal solution is global. iv Índice de figuras 2.1. Funciónconvexa. ............................... 7 2.2. Funcióncóncava. ............................... 8 2.3. Epigrafo.................................... 9 2.4. Hipografo ................................... 9 4.1. Representación gráfica de la solución al problema de la Catenaria mediante CVXR y solución analítica. . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2. Representación gráfica del problema de la Catenaria con masa no uniforme. 53 4.3. Representación gráfica del problema de la Catenaria con restricciones sobre elsuelo. .................................... 54 v Índice de cuadros 4.1. Funciones escalares disponibles en el paquete CVXR . . . . . . . . . . . . 39 4.2. Funciones element-wise disponibles en el paquete CVXR . . . . . . . . . . 41 4.3. Funciones vector/matriz disponibles en el paquete CVXR . . . . . . . . . 42 4.4. Solución del problema de mínimos cuadrados mediante CVXR y lm. . . . 48 4.5. Solución del probema de mínimos cuadrados con restricciones. . . . . . . 49 4.6. Primero valores de la solución del problema de la Catenaria. . . . . . . . . 51 4.7. Primeros valores de la solución del probema de la catenaria con masa no uniforme. ................................... 52 4.8. Primeros valores de la solución del probema de la catenaria con restricciones sobreelsuelo.................................. 54 4.9. Resultados del problema Logístico mediante CVXR y glm . . . . . . . . 57 vii 2.1. CONJUNTOS CONVEXOS x, y ∈S2⇒ ∀α∈[0,1],αx + (1 −α)y∈S2ya que S2es convexo. Por tanto αx + (1 −α)y∈S1∩S2, luego S1∩S2es convexo. Esta propiedad se extiende a la intersección de un número infinito de conjuntos: si Sα es convexo para cada α∈Aentonces ∩α∈ASαes convexo. Suma La suma de dos conjuntos se define como: S1+S2={x1+x2|x1∈S1, x2∈S2} Si S1yS2son convexos entonces S1+S2es convexo. Demostración Para dos puntos cualesquiera x, y ∈S1+S2se tiene que: x∈S1+S2⇒x=x1+x2con x1∈S1,x2∈S2. y∈S1+S2⇒y=y1+y2con y1∈S1,y2∈S2. Sea α∈[0,1], αx + (1 −α)y=αx1+αx2+ (1 −α)y1+ (1 −α)y2= αx1+ (1 −α)y1+αx2+ (1 −α)y2∈S1+S2 ya que αx1 + (1 −α ) y1∈S1 porque S1 es convexo y αx2 + (1 −α ) y2∈S2 porque S2 es convexo, queda demostrado que S1+S2es convexo. Producto cartesiano El producto cartesiano de dos conjuntos se define como: S1×S2={(x1, x2)|x1∈S1, x2∈S2}(2.2) Si S1yS2son convexos entonces S1×S2es convexo. Demostración Para dos puntos cualesquiera x, y ∈S1×S2se tiene que: x∈S1×S2⇒x= (x1, x2)con x1∈S1,x2∈S2. y∈S1×S2⇒y= (y1, y2)con y1∈S1,y2∈S2. Sea α∈[0,1], αx + (1 −α)y=α(x1, x2) + (1 −α)(y1, y2) = = (αx1+ (1 −α)y1, αx2+ (1 −α)y2)∈S1×S2 ya que αx1 + (1 −α ) y1∈S1 porque S1 es convexo y αx2 + (1 −α ) y2∈S2 porque S2 es convexo, queda demostrado que S1×S2es convexo. 6CAPÍTULO 2. ANÁLISIS CONVEXO 2.2. FUNCIONES CONVEXAS 2.2. Funciones convexas 2.2.1. Definiciones Una función f : Rn→R es convexa si dom f es un conjunto convexo y para todo x, y ∈dom f, y para todo θ∈[0,1] se tiene que: f(θx + (1 −θ)y)≤θf(x) + (1 −θ)f(y)(2.3) 0 2 4 6 8 10 0 5000 10000 20000 x exp Figura 2.1: Función convexa. Una función f : Rn→R es cóncava si dom f es un conjunto convexo y para todo x, y ∈dom f, y para todo θ∈[0,1] se tiene que: f(θx + (1 −θ)y)≥θf(x) + (1 −θ)f(y)(2.4) CAPÍTULO 2. ANÁLISIS CONVEXO 7 2.2. FUNCIONES CONVEXAS 0 2 4 6 8 10 −2 −1 0 1 2 x log Figura 2.2: Función cóncava. Se cumple que si f es cóncava entonces -f es convexa. Las funciones afines son aquellas para las que se da la igualdad, así que son tanto convexas como cóncavas. Es decir, una función f : Rn→R es afín si para todo x, y ∈dom f, y para todo θ∈[0,1] se tiene que: f(θx + (1 −θ)y) = θf(x) + (1 −θ)f(y)(2.5) Grafo El grafo de una función f:Rn→Res un subconjunto de Rn+1 que se define como: {(x, f(x)) |x∈domf}(2.6) Epigrafo El epigrafo de una función f:Rn→Rse define como: Epi(f) = {(x, t)|x∈domf, f(x)≤t}(2.7) Es el conjunto de puntos situado sobre la curva. 8CAPÍTULO 2. ANÁLISIS CONVEXO 2.2. FUNCIONES CONVEXAS Figura 2.3: Epigrafo Hipografo El hipografo de una función f:Rn→Rse define como: Hip(f) = {(x, t)|x∈domf, t ≤f(x)}(2.8) Es el conjunto de puntos situado bajo la curva. Figura 2.4: Hipografo 2.2.2. Propiedades básicas Un principio fundamental del análisis convexo es la relación entre funciones convexas y conjuntos convexos, y se establece a través del epigrafo e hipografo de las funciones. Lema. Una función f : Rn→R es convexa si y solo si su epigrafo es un conjunto convexo. Demostración Supongamos que f es convexa, entonces sean ( x, t ) , ( y, s ) ,∈Epi ( f )se tiene que para todo λ∈[0,1]: f(λx + (1 −λ)y)≤λf(x) + (1 −λ)f(y)≤λt + (1 −λ)s Por tanto, λ(x, t) + (1 −λ)(y, s)∈Epi(f), por tanto se tiene que Epi(f) es convexo. Supongamos ahora que Epi(f) es un conjunto convexo. Dados x, y ∈dom f por definición de Epi(f) se tiene que (x, f(x)),(y, f(y)) ∈Epi(f). CAPÍTULO 2. ANÁLISIS CONVEXO 9 2.2. FUNCIONES CONVEXAS Como Epi(f) es convexo, λ ( x, f ( x )) + (1 −λ )( y, f ( y )) ∈Epi ( f ). Desarrollando esta expresion obtenemos: λ(x, f(x)) + (1 −λ)(y, f(y)) = (λx + (1 −λ)y, λf(x) + (1 −λ)f(y)) ∈Epi(f) Por tanto, f(λx + (1 −λ)y)≤λf(x) + (1 −λ)f(y) es decir, f es convexa. Lema. Una función f : Rn→R es cóncava si y solo si su hipografo es un conjunto convexo. Demostración Supongamos que f es cóncava, entonces sean ( x, t ) , ( y, s ) ∈Hip ( f )se tiene que para todo λ∈[0,1]: −f(λx + (1 −λ)y)≤λ(−1)f(x) + (1 −λ)(−1)f(y)⇒ f(λx + (1 −λ)y)≥λf(x) + (1 −λ)f(y)≥λt + (1 −λ)s Por tanto, λ(x, t) + (1 −λ)(y, s)∈Hip(f), por tanto se tiene que Hip(f) es convexo. Supongamos ahora que Hip(f) es un conjunto convexo. Dados x, y ∈dom f por definición de Hip(f) se tiene que (x, f(x)),(y, f(y)) ∈Hip(f). Como Hip(f) es convexo, λ ( x, f ( x )) + (1 −λ )( y, f ( y )) ∈Hip ( f ). Desarrollando esta expresión obtenemos: λ(x, f(x)) + (1 −λ)(y, f(y)) = (λx + (1 −λ)y, λf(x) + (1 −λ)f(y)) ∈Hip(f) Por tanto, f(λx + (1 −λ)y)≥λf(x) + (1 −λ)f(y) es decir, f es cóncava. Estas relaciones se pueden expresar también de la siguiente forma, suponiendo f una función convexa y g cóncava: f(x) = inf{y|(x, y)∈Epi(f)} g(x) = sup{y|(x, y)∈Hip(g)} 2.2.3. Operaciones que preservan la convexidad En esta sección describiremos algunas operaciones que preservan la convexidad o concavidad de las funciones o que nos permiten construir nuevas funciones convexas o cóncavas. Suma Si f1yf2son funciones convexas, entonces también lo es su suma f1+f2. Demostración 10 CAPÍTULO 2. ANÁLISIS CONVEXO 2.2. FUNCIONES CONVEXAS Como f1 y f2 son convexas, dados x, y ∈ dom ( f1 + f2 ) se tiene que para todo λ∈ [0 , 1]: f1(λx + (1 −λ)y)≤λf1(x) + (1 −λ)f1(y) f2(λx + (1 −λ)y)≤λf2(x) + (1 −λ)f2(y) Entonces, (f1+f2)(λx + (1 −λ)y) = f1(λx + (1 −λ)y) + f2(λx + (1 −λ)y)≤ ≤λf1(x) + (1 −λ)f1(y) + λf2(x) + (1 −λ)f2(y) = =λ(f1(x) + f2(x)) + (1 −λ)(f1(y) + f2(y)) Luego, f1+f2es convexa. Producto por escalar Si f es una función convexa y α≥0, entonces αf es convexa. Demostración Dados x, y ∈dom f se tiene que para todo λ∈[0,1]: f(λx + (1 −λ)y)≤λf(x) + (1 −λ)f(y)⇒ αf(λx + (1 −λ)y)≤α(λf(x) + (1 −λ)f(y)) = λ(αf(x)) + (1 −λ)(αf(y)) Luego, αf es convexa. Transformaciones afines Suponiendo f:Rn→R, A ∈Rn×myb∈R. Se define g:Rm→Rcomo g(x) = f(Ax +b)(2.9) con dom g ={ x|Ax + b∈dom f }. Entonces si f es convexa también lo es g. De la misma forma, cuando f es cóncava también lo es g. Demostración Dados x, y ∈dom g cualesquiera, entonces para todo λ∈[0,1] se tiene: g(λx + (1 −λ)y) = f(A(λx + (1 −λ)y) + b) = (Por definición de g) =f(λ(Ax +b) + (1 −λ)(Ay +b)) ≤(Descomponiendo b=λb + (1 −λ)b) ≤λf(Ax +b) + (1 −λ)f(Ay +b) = (f es convexa) =λg(x) + (1 −λ)g(y) Composición Vamos a estudiar las condiciones de f : R→ ( R∪ + ∞ )y g : Rn→ ( R∪ + ∞ )que garantizan la convexidad o concavidad de h=f◦g. Lema : Si f : R→ ( R∪ + ∞ )es convexa y no decreciente y g : Rn→ ( R∪ + ∞ )es convexa, entonces h=f◦ges convexa. Demostración CAPÍTULO 2. ANÁLISIS CONVEXO 11 2.3. PROBLEMAS DE OPTIMIZACIÓN CONVEXA Dados x, y ∈dom g se tiene que para todo λ∈[0,1]: (f◦g)(λx + (1 −λ)y) = f(g(λx + (1 −λ)y)) ≤ ≤f(λg(x) + (1 −λ)g(y)) (g convexa y f no decreciente) ≤λf(g(x)) + (1 −λ)f(g(y)) = (f convexa) =λ(f◦g)(x) + (1 −λ)(f◦g)(y) Luego, f◦ges convexa. Se obtienen resultados similares para funciones cóncavas y/o no crecientes: Si f : R→ ( R∪ + ∞ )es convexa y no creciente y g : Rn→ ( R∪−∞ )es cóncava, entonces h=f◦ges convexa. Si f : R→ ( R∪−∞ )es cóncava y no decreciente y g : Rn→ ( R∪−∞ )es cóncava, entonces h=f◦ges cóncava. Si f : R→ ( R∪−∞ )es cóncava y no creciente y g : Rn→ ( R∪ + ∞ )es convexa, entonces h=f◦ges cóncava. Además se tienen reglas similares para las funciones con múltiples argumentos. 2.3. Problemas de optimización convexa Un problema de optimización convexa es de la forma: Minimizar f0(x) Sujeto a: fi(x)≤0, i = 1, ..., n aT jx=bj, j = 1, ..., m x∈X (2.10) El vector x son las variables del problema, contenida en un conjunto X⊆Rk . Tanto la función objetivo f0 : Rk→R , como las funciones que describen las restricciones f1, ..., fn:Rk→Rson convexas. Se dice que un punto es factible para el problema si satisface todas las restricciones. El problema es factible si existe al menos un punto factible e infactible en caso contrario. El valor óptimo para el problema que hemos definido viene dado por: p∗=inf{f0(x)|fi(x)≤0, i = 1, ..., m aT ix=bi, i = 1, ..., n}(2.11) Si el problema es infactible consideramos p∗ = ∞ , ya que el ínfimo de un conjunto vacío es ∞. En el caso en que haya puntos factibles xk tal que f0 ( xk ) → −∞ cuando k→ ∞ entonces decimos que el problema es ilimitado. Decimos que x∗es la solución óptima si x∗es factible y f(x∗) = p∗. El conjunto de problemas convexos (CP, convex problem) es un subconjunto estricto de los problemas de programación no lineal (NLP, nonlinear programming). Un problema no 12 CAPÍTULO 2. ANÁLISIS CONVEXO 2.3. PROBLEMAS DE OPTIMIZACIÓN CONVEXA lineal es un problema de optimización donde la función objetivo o alguna de las funciones de las restricciones no son lineales. Dentro de la definición general del problema de optimización convexa se encuentran casos especiales como el problema de mínimos cuadrados (LS, least squares), los problemas de programación lineal (LP, linear programming), los problemas de programación semidefinida (SDP, semidefinite program) [71], los problemas cuadráticos (QP, quadratic program), los problemas de cono de segundo orden (SOCP, second-order cone program) [45], y los problemas de programación geométrica (GP, geometric programming) ([23], [5]). Vamos a definir cada uno de estos tipos de problemas explicitando su formulación. 2.3.1. Problema lineal El problema lineal es el tipo de problema más popular y más estudiado. Se formula de la siguiente forma: Minimizar cTx, Sujeto a: aT ix≤bi, i = 1, ..., n, (2.12) donde c, a1, ..., an∈Rk,b1, ..., bn∈Rson parámetros del problema. 2.3.2. Problema de mínimos cuadrados El problema de mínimos cuadrados es un problema sin restricciones que se formula como: Minimizar kAx +bk2 2(2.13) donde A∈Rn×k,b∈Rnson los parámetros del problema. 2.3.3. Problema cuadrático Un problema de optimización es cuadrático si la función objetivo es cuadrática y las restricciones son afines. Minimizar cTx+1 2xTQx Sujeto a: Ax =b(2.14) donde Q∈Sk + ( Sk + es el conjunto de las matrices cuadradas de dimensión k×k definidas positivas), c∈Rk,A∈Rn×k,b∈Rnson los parámetros del problema. 2.3.4. Problema de cono de segundo orden El problema de cono de segundo orden se formula como: Minimizar fTx Sujeto a: kAix+bik2≤cT ix+di, i = 1, ..., m Fx =g (2.15) CAPÍTULO 2. ANÁLISIS CONVEXO 13 2.4. ALGORITMOS NUMÉRICOS donde A1, ..., An∈Rni×n , F∈Rp×n , b∈Rni , c1, ..., cm, d1, ..., dm∈Rn , g∈Rp son los parámetros del problema. 2.3.5. Problemas de programación semidefinida El problema de programación semidefinida se formula como: Minimizar cTx Sujeto a: x1F1+... +xnFn+G≤0, Ax =b (2.16) donde F1, ..., Fn, G ∈Sk + ( Sk + es el conjunto de las matrices cuadradas de dimensión k×kdefinidas positivas), A∈Rp×k,b∈Rpson los parámetros del problema. 2.3.6. Problemas de programación geométrica Un problema de programación geométrica se formula como: Minimizar f0(x) Sujeto a: fi(x)≤1, i = 1, ..., m hi(x) = 1, i = 1, ..., p (2.17) donde f0, ..., fm son de la forma fi ( x1, ..., xn ) = PL k=1 ckxa1 1...xan n y h1, ..., hp son funciones monomiales. 2.4. Algoritmos numéricos Gran parte del estudio de algoritmos numéricos en las últimas décadas se ha centrado en los métodos de punto interior ([84], [68]). Inicialmente, dicho trabajo se limitó a LPs ([40], [60]), pero pronto se amplió para abarcar otros CPs ([51], [52]). Estos métodos consisten en la búsqueda del óptimo a través de caminos que recorren el interior de la región factible. Existen varios solvers excelentes disponibles, tanto comerciales como de distribución gratuita, vamos a enumerar algunos de ellos. Se dividen en dos clases: los que se basan en las formas estándar y los que se basan en un código personalizado. 2.4.1. Formas estándar La mayoría de los solvers para programación convexa están diseñados para una familia limitada de problemas con una estructura específica, o que obedecen a ciertas convenciones, se conocen como formas estándar. Dos ejemplos comunes pueden ser el problema de mínimos cuadrados o los problemas lineales. La mayoría de los problemas convexos se pueden transformar en una o más de estas formas estándar. Sin embargo, las transformaciones requeridas a menudo están lejos de ser obvias. Programas convexos smooth Se han diseñado varios solvers para resolver problemas convexos en forma estándar de problemas no lineales, bajo la condición de que la función objetivo y las funciones de 14 CAPÍTULO 2. ANÁLISIS CONVEXO 2.4. ALGORITMOS NUMÉRICOS restricciones de desigualdad sean dos veces diferenciables al menos en la región en que el algoritmo desea buscar. A estos problemas los llamaremos CP smooth, mientras que a los CP que no lo cumplan los llamaremos CP nonsmooth. Entre los paquetes de software que resuelven CP smooth encontramos LOQO ([73]) (emplea un método primal-dual), y el paquete comercial MOSEK ([50]) (implementa el algoritmo homogéneo). Una dificultad práctica en el uso de solvers de CP smooth es que el solver debe poder calcular el gradiente y el hessiano del objetivo y las funciones de restricción de desigualdad en los puntos elegidos. En algunos casos, para realizar estos cálculos, puede ser necesario un código personalizado. Muchos marcos de modelado simplifican este proceso en la mayoría de los casos al permitir que las funciones se expresen en forma matemática natural y que calculen las derivadas automáticamente ([24], [13]). Programas cónicos Una familia de formas estándar que se ha vuelto muy común son las forma cónicas primales y duales: Minimizar cTxMinimizar bTy Sujeto a: Ax =b o Sujeto a: ATy+z=c x∈K z ∈K∗ (2.18) donde los conjuntos K y K∗ son conos convexos cerrados. La forma cónica más común es la LP donde K=K∗es: K=K∗=Rn +={x∈Rn|xi≥0, i = 1, ..., n}(2.19) Se puede demostrar que prácticamente cualquier CP puede representarse en forma cónica, con la elección adecuada de K o K∗ ([52]). En la práctica, dos formas cónicas (además del LP) dominan los recientes estudios e implementaciones. Uno es el programa semidefinido (SDP), para el cual K = K∗ es un isomorfismo del cono de matrices semidefinidas positivas. Sn +={X=XT∈Rn×n|λmin(X)≥0}(2.20) Otro es el programa de cono de segundo orden, para el cual K= K∗ es el producto cartesiano de uno o más conos de segundo orden (o de Lorentz). K=Qn1×···×Qnk, Qn={(x, y)∈Rn×R| k xk2≤y}(2.21) Para SDP y SOCP se han descubierto muchas aplicaciones y su geometría admite ciertas optimizaciones algorítmicas útiles ([53], [31]). Los solvers disponibles públicamente para SDP y SOCP incluyen SeDuMi [65], CDSP [11], SDPA [27], SDPT3 [69] y DSDP [8]. Son generalmente bastante eficientes, confiables y están completamente basados en datos: es decir, no requieren código externo para realizar cálculos de funciones. Programación geométrica Un problema de programación geométrica es en realidad un caso un tanto único, ya que en realidad no es convexo, pero una transformación simple produce un problema CAPÍTULO 2. ANÁLISIS CONVEXO 15 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD 3.1. Conjunto de reglas de convexidad El conjunto de reglas de convexidad establece como las variables, parámetros, funciones y conjuntos pueden ser combinados para formar DCPs. Los DCPs son un subconjunto estricto de los CPs generales, por lo que otra manera de decir esto es que el conjunto de reglas impone un conjunto de convenciones o restricciones a los CPs. El conjunto de reglas se puede dividir en cuatro categorías: reglas de nivel superior, reglas sin productos, reglas de signo y reglas de composición. Antes de enunciar las reglas, vamos a definir dos conceptos que utilizaremos, expresión constante y aserción. Una expresión constante es una expresión que involucra solo valores numéricos y / o parámetros mientras que una expresión no constante depende del valor de al menos una variable del problema. Obviamente una expresión constante es trivialmente afín, convexa y cóncava. Estas expresiones deben estar bien definidas, es decir, que se verifiquen para cualquier conjunto de parámetros que satisfagan las aserciones del problema. Una aserción se asemeja a una restricción, pero contiene solo expresiones constantes, está determinada completamente por los valores numéricos proporcionados para los parámetros de un modelo antes del comienzo de cualquier algoritmo de optimización numérico. Las aserciones sirven como condiciones previas, garantizando que un problema es numéricamente válido o físicamente significativo. Hay varias razones por las cuales se puede necesitar una aserción; por ejemplo: Representar los límites físicos dictados por el modelo. Por ejemplo, si un parámetro w representa el peso físico de un objeto, una afirmación w > 0impone el hecho de que el peso debe ser positivo. Para asegurar que está bien definido. Por ejemplo, si x, y, z son variables y a, b, c son parámetros, entonces la restricción de desigualdad ax + by + z/c ≤ 1está bien planteada solo si c es distinta de cero; esto puede ser asegurado por una aserción como c6= 0 oc > 0. Para garantizar el cumplimiento de las condiciones previas adjuntas a una función o conjunto en la librería. Por ejemplo, una función f ( x ) = kxkp , está parametrizada por un valor p≥ 1. Si p se proporciona como un parámetro, entonces se requeriría una aserción como p≥ 1para garantizar que la función se utiliza correctamente. Para asegurar el cumplimiento con las reglas de signo o las reglas de composición que veremos a continuación. 3.1.1. Reglas de nivel superior Como su nombre indica, las reglas de nivel superior establecen la estructura superior de los DCPs. Estas reglas son más descriptivas que restrictivas, en el sentido de que casi todos los problemas convexos siguen estas convenciones, pero se deben explicitar. Tipos de problemas. Un DCP válido puede ser: T1 un problema de minimización: un objetivo convexo y cero o más restricciones convexas. 22 CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD T2 un problema de maximización: un objetivo cóncavo y cero o más restricciones convexas. T3 un problema de factibilidad: no hay objetivo, y una o más restricciones convexas. Un DCP válido también puede incluir cualquier número de aserciones. Restricciones. Las restricciones válidas incluyen: T4 una restricción de igualdad con expresiones afines a izquierda y derecha. T5 una desigualdad menor o menor o igual que ( < , ≤ ), con una expresión convexa a la izquierda y una expresión cóncava a la derecha. T6 una desigualdad mayor 0 mayor o igual que ( > , ≥ ), con una expresión cóncava a la izquierda y una expresión convexa a la derecha. T7 una restricción de pertenencia establecida ( lexp1, ..., lexpm ) ∈cset , donde m≥ 1, lexp1, ..., lexpmson expresiones afines, y cset es un conjunto convexo. Las restricciones de no igualdad ( 6 =) y las restricciones de no pertenencia no están permitidas, porque son convexas solo en casos excepcionales, y la filosofía de la programación convexa disciplinada no considera casos especiales. Sin embargo, en las aserciones estas dos operaciones pueden ser empleadas libremente ya que las aserciones no están restringidas de la manera en que lo están las restricciones reales. Expresiones y afirmaciones constantes. T8 Cualquier expresión bien definida que conste únicamente de valores numéricos y parámetros es una expresión constante válida. T9 Cualquier expresión booleana que realice pruebas o comparaciones en expresiones constantes válidas es una aserción válida. T10 Si una función o conjunto está parametrizado, entonces esos parámetros deben ser expresiones constantes válidas. Si consideramos por ejemplo esta función y conjunto: fp:Rn→R, fp(x) =kxkp Bp={x∈Rn| k xkp≤1}(3.1) Esta última regla simplemente impone que los parámetros como p deben ser constantes. Por supuesto esto es generalmente asumido, pero debemos explicitarlo. 3.1.2. Reglas libres de productos La convención más destacada impuesta por la programación convexa disciplinada es la prohibición de productos entre expresiones no constantes, lo que incluye además operaciones relacionadas, como la exponencial. La necesidad de esta imposición se debe a que no existe un principio en el análisis convexo que pueda determinar la curvatura en tal caso. Podemos observarlo a partir de ejemplos. Supongamos que x es una variable escalar, entonces: La expresión x·x, un producto de dos expresiones afines, es convexa. CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 23 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD La expresión x·logx , un producto entre una expresión afín y una expresión cóncava, es convexa. La expresión x·ex , un producto entre una expresión afín y una expresión convexa, no es convexa ni cóncava. Por otro lado, a partir de los principios más básicos del análisis convexo tenemos las siguientes afirmaciones, probadas en el capítulo anterior: La suma de dos o más expresiones convexas (cóncavas, afines) es convexa (cóncava, afín). El producto de una expresión convexa (cóncava) y de una expresión constante no negativa es convexa (cóncava). El producto de una expresión convexa (cóncava) y una expresión constante no positiva, o la simple negación de la primera, es cóncava (convexa). El producto de una expresión afín y cualquier constante es afín. Transformamos estos principios en un conjunto de reglas y como resultado tenemos las reglas, apropiadamente llamadas, libres de productos: Reglas libre de productos: todas las expresiones numéricas válidas deben estar libres de productos; tales expresiones incluyen: PN1 Una variable simple. PN2 Una expresión constante. PN3 Una llamada a una función en la librería. Cada argumento de la función debe ser una expresión libre de producto. PN4 La suma de dos o más expresiones libres de productos. PN5 La diferencia de expresiones libres de producto. PN6 La negación de una expresión libre de producto. PN7 El producto de una expresión libre de producto y una expresión constante. PN8 La división de una expresión libre de producto por una expresión constante. Asumimos en cada una de estas reglas que los resultados están bien definidos; por ejemplo, que las dimensiones son compatibles. En el caso escalar, un modo de reformular estas reglas es decir que una expresión numérica válida puede reducirse a la forma: a+ n X i=1 bixi+ L X j=1 cjfj(argj,1, argj,2, ..., argj,mj)(3.2) donde a, bi , cj son constantes; xi son las variables problema; y fj : Rmj→R son funciones de la librería, y sus argumentos argj,k son expresiones libres de productos. Podemos observar ciertos casos especiales: 1. Si L= 0 es una expresión afín simple. 2. Si, además, b1=b2=... =bn= 0, entonces es una expresión constante. 24 CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD Para una ilustración del uso de estas reglas libres de productos, supongamos que a, b y c son parámetros; x, y, z son variables; y f ( · ) , g ( ·,· ) y h ( ·,·,· )son funciones de la librería. Entonces la expresión: af(x) + y+ (h(x, bg(y, z), c)−z+b)/a satisface la regla libre de producto, puede ser reescrita siguiendo forma 3.2: (b/a) + y+ (−1/a)z+af(x) + (1/a)h(x, bg(y, z), c) Por otro lado, la expresión siguiente: axy/2−f(x)g(y, z) + h(x, yb, z, c) no se puede escribir como 3.2, hay un producto entre dos variables (xy) y un producto entre dos funciones (f(x)g(y,z)). Para expresiones de conjuntos, las reglas libres de producto son similares. Utilizamos también las operaciones vistas en el capítulo anterior. Reglas libres de producto para expresiones de conjunto: todas las expresiones de conjunto válidas usadas en restricciones deben estar libres de producto; tales expresiones incluyen: PSI Una llamada a un conjunto convexo en la librería. PS2 Una llamada a una función en la librería. Cada argumento de la función debe ser una expresión de conjunto libre de producto o una expresión constante (numérica). PS3 La suma de conjuntos libres de producto, o de un conjunto libre de producto y una expresión constante, o viceversa. PS4 La diferencia de un conjunto libre de productos y una expresión constante, o viceversa. PS5 La negación de un conjunto libre de productos. PS6 El producto de un conjunto libre de productos y una expresión constante. PS7 La división de un conjunto libre de producto por una expresión constante. PS8 La intersección de dos o más conjuntos libres de productos. PS9 El producto cartesiano de dos o más conjuntos libres de productos. También asumimos en cada una de estas reglas que los resultados están bien definidos; por ejemplo, que las dimensiones son compatibles. En otras palabras llamamos conjuntos válidos a expresiones que se pueden reducir a la forma: a+ n X i=1 biSi+ L X j=1 cjfj(argj,1, argj,2, ..., argj,mj)(3.3) donde a, bi , cj son constantes, y Sk son conjuntos de la librería, o intersecciones y / o productos cartesianos de conjuntos válidos. Las funciones fj : Rm j→R son funciones de la librería, y sus argumentos argj,k son expresiones de conjuntos libres de productos. CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 25 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD Las expresiones de conjuntos están más restringidas que las expresiones numéricas, ya que las funciones fjdeben ser afines, lo veremos a continuación. Hemos visto en el capítulo anterior que la suma, intersección y producto cartesiano de conjuntos convexos es convexa; mientras que las uniones y las diferencias de conjuntos generalmente no son convexas. Lo que puede no estar claro es por qué las PS8-PS9 se consideran reglas “sin productos”. Al examinar estas reglas en términos de funciones de indicadores, el enlace se vuelve claro. Consideremos, por ejemplo, el problema. Minimizar ax +by, Sujeto a: (x, y)∈(S1×S2)∪S2 si φ1 y φ2 son funciones indicador convexas para los conjuntos S1 y S2 respectivamente, es decir son de la forma: φi(x) =      0si x ∈Si ∞si x /∈Si entonces el problema se puede reducir a: Minimizar ax +by + (φ1(x) + φ2(y))φ2(x, y) El producto escalar se ha convertido en la suma de las funciones indicador y la unión se ha convertido en un producto. Ahora la función objetivo viola la regla libre de producto debido a la operación generada por la unión, luego no se permite. Si consideramos ahora el problema: Minimizar ax Sujeto a: x∈S1∩S2 se puede reducir a: Minimizar ax +φ1(x) + φ2(x) La intersección se ha convertido también en la suma de las funciones indicador y por tanto es válida. 3.1.3. Reglas de signo Las reglas de signo se denominan de esta manera porque su objetivo es determinar el signo de c1...cL en 3.2 . Son un conjunto de reglas que establecen las condiciones suficientes para determinar si las expresiones son convexas, cóncavas o afines. Podemos establecer de manera concisa las reglas de signo para las expresiones numéricas de la siguiente manera. Reglas de signo para expresiones numéricas. Dada una expresión libre de producto, lo siguiente debe ser cierto en su forma reducida 3.2: 26 CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD SN1 Para que la expresión 3.2 sea convexa, cada término cjfj ( ... )debe ser convexo; por lo tanto, uno de los siguientes puntos debe ser cierto: •fj(argj,1, ..., argj,mj)es afín. •fj(argj,1, ..., argj,mj)es convexo y cj≥0. •fj(argj,1, ..., argj,mj)es cóncavo y cj≤0. SN2 Para que la expresión 3.2 sea cóncava, cada término cjfj ( ... )debe ser cóncavo; por lo tanto, uno de los siguientes puntos debe ser cierto: •fj(argj,1, ..., argj,mj)es afín. •fj(argj,1, ..., argj,mj)es cóncavo y cj≥0. •fj(argj,1, ..., argj,mj)es convexo y cj≤0. SN3 Para que la expresión 3.2 sea afín, cada función fj debe ser afín, al igual que cada uno de sus argumentos argj,1, ..., argj,mj. SN4 Para que la expresión 3.2 sea constante, entonces se debe cumplir que L = 0 y b1=b2=... =b= 0. Todos los argumentos de las funciones deben obedecer estas reglas también, según su curvatura esperada dictada por las reglas de composición. Por ejemplo, supongamos que se espera que la expresión af(x) + y+ (h(x, bg(y, z), c)−z+b)/a sea convexa y que la librería indica que la función f ( · )es convexa, g ( ·,· )es cóncava y h(·,·,·)es convexa. Entonces la regla del signo dicta que: Para que af(x)sea convexa se debe tener a≥0 Para que (1/a)h(x, bg(y, z), c)sea convexa se debe tener 1/a ≥0 Además se requiere que el segundo argumento de h sea convexo, es decir, bg ( y, z ) debe ser convexo, esto implica que b≤0 Es responsabilidad del modelador asegurarse de que los valores de los coeficientes c1, ..., cL obedecen las reglas de signo. Esto se puede lograr agregando las aserciones apropiadas al modelo. Solo hay una regla de signo para las expresiones de conjuntos: Regla de signos para conjuntos. Dada una expresión de conjunto libre de producto, lo siguiente debe ser cierto para su forma reducida 3.3: SS1 Cada función fj, y cualquier función utilizada en sus argumentos, debe ser afín. Por ejemplo, para una expresión en la forma reducible con L= 0, x∈a+Pn i=0 biSi↔ ∃(t1, t2, ..., tn)∈S1×S2×... ×Sntal que x=a+Pn i=0 biti CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 27 3.1. CONJUNTO DE REGLAS DE CONVEXIDAD Cuando L > 0, también se realizan sustituciones similares de forma recursiva en los argumentos de la función, produciendo un resultado similar: una serie de restricciones de pertenencia a conjuntos simples de la forma tk∈Sk , y una restricción de igualdad única. Por lo tanto, para asegurar que esta restricción de igualdad implícita sea convexa, las expresiones de conjuntos (específicamente, aquellas usadas en restricciones) deben reducirse a combinaciones afines de conjuntos. Por supuesto, las expresiones de conjuntos utilizadas en las aserciones no están restringidas de esta manera. 3.1.4. Reglas de composición Un principio básico del análisis convexo es que la composición de una función convexa con una transformación afín sigue siendo convexa. De hecho, bajo ciertas condiciones, también es cierto para composiciones con transformaciones no lineales. Se han planteado estos resultados en el capítulo anterior. Supongamos que tenemos dos funciones convexas f, g y queremos realizar una composición. Por el resultado visto en el capítulo 2, sabemos que para asegurar que f◦g es convexa necesitamos que f sea no decreciente. Se puede generalizar este razonamiento a funciones con múltiples argumentos, creando un conjunto de reglas denominadas reglas de composición. Las reglas de composición. Considera una expresión numérica de la forma f(argj,1, ..., argj,mj)(3.4) donde f es una función de la librería. La expresión debe satisfacer exactamente una de las siguientes reglas: C1-C3 Para que la expresión sea convexa, entonces f debe ser afín o convexa, y una de las siguientes afirmaciones debe ser cierta para cada k= 1, ..., m: •C1 fes no decreciente en el argumento k y argkes convexo. •C2 fes no creciente en el argumento k y argkes cóncavo. •C3 argkes afín. C4-C6 Para que la expresión sea cóncava, entonces f debe ser afín o cóncavo, y una de las siguientes afirmaciones debe ser cierta para cada k= 1, ..., m: •C4 fes no decreciente en el argumento k y argkes cóncavo. •C5 fes no creciente en el argumento k y argkes convexo. •C6 argkes afín. C7 Para que la expresión sea afín, entonces f debe ser afín, y cada argumento es afín para todo k= 1, ..., m. Estas reglas pueden ser demasiado conservativos y no dar como válida una composición que si lo sea. Es un pequeño precio que hay que pagar. 28 CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 3.2. VERIFICACIÓN 3.2. Verificación Para resolver un problema mediante programación convexa disciplinada (DCP) el primer paso será la verificación, es decir, comprobar que es un problema convexo disciplinado válido, probar que incluye solo funciones y conjuntos de la librería y los combina de acuerdo al conjunto de reglas. La prueba de validez es de naturaleza jerárquica, refleja la estructura del problema y sus expresiones. Para ilustrar el proceso, consideremos el siguiente problema de optimización: Minimizar ex+y Sujeto a: ax2−3≤√y log(x+y)≥3 (3.5) donde x, y son variables. Vamos a dividir la prueba en dos etapas: En la primera etapa se verificará si cada una de las expresiones del problema está libre de productos. Vamos a representar la jerarquía de la prueba, además de indicar la regla empleada en cada una de ellas. ex+yes libre de productos, por (PN4) •exes libre de productos (PN3) ◦xes libre de productos (PN1) •yes libre de productos (PN1) ax2−3es libre de producto, por (PN5) •ax2es libre de producto (PN7) ◦a es libre de producto, por (PN2) ◦x2es libre de producto, por (PN3) xes libre de producto, por (PN1) •3 es libre de producto, por (PN2) √yes libre de producto, por (PN3) •y es libre de producto, por (PN1) log(x+y)es libre de producto (PN3) •x+yes libre de producto (PN4) ◦xes libre de producto, por (PN1) ◦yes libre de producto, por (PN1) 3 es libre de producto, por (PN2) En la segunda etapa se verificará si el problema cumple las reglas de nivel superior, signo y composición. Lo representaremos de la misma forma. Como veremos a continuación, el problema de minimización es válido si a≥ 0por (T1) La función objetivo es válida, por (T1) CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 29 3.2. VERIFICACIÓN ex+yes convexa (SN1) •exes convexa, por (c1) ◦e.es convexa, no decreciente (librería) xes convexa, (SN1) •yes convexa, (SN1) La primera restricción es válida si a≥0, por (T1) ax2−3≤√yes valida si a≥0, por (T5) •ax2−3es convexa si a≥0, por (SN1) ◦ax2es convexa si a≥0, por (SN1) x2es convexa, por(C1) square(.) es convexa, creciente (librería) xes convexa, (SN1) ◦ −3es convexa, (SN1) •√yes cóncava, por (C4) ◦√.es cóncava y no decreciente (librería) ◦yes cóncava (SN2) La segunda restricción es válida, por (T1) log(x+y)≥3es válida por(T6) •log(x+y)es cóncavo, (C4) ◦x+yes cóncava (SN2) xes cóncava, (SN2) yes cóncava, (SN2) •3 es convexa, (SN1) También es posible examinar gráficamente la estructura del problema y su validez mediante un árbol de expresiones del problema. Cada hoja representaría una operación o variable y se incluyen también las reglas usadas en cada paso. El proceso de verificación garantiza una de estas tres conclusiones: Válido: las reglas están totalmente satisfechas. Condicionalmente válido: las reglas se cumplirán completamente si se cumplen una o más condiciones previas adicionales en los parámetros. No válido: una o más de las reglas de convexidad han sido violadas. En este ejemplo, llegamos a la conclusión de que es condicionalmente válido: a través del análisis hemos visto que debe cumplirse una condición adicional, a≥ 0. Si esta condición previa estuviera asegurada de alguna manera, entonces la prueba habría determinado de manera concluyente que el problema es DCP válido. Esto se puede lograr agregando la aserción a≥ 0a la lista de restricciones. Si, por otro lado, tuviéramos que agregar una aserción a < 0, se violaría la regla de signos SN2. 30 CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 3.3. CREAR PROGRAMAS CONVEXOS DISCIPLINADOS 3.3. Crear programas convexos disciplinados Como hemos visto, verificar el conjunto de reglas de convexidad es condición suficiente, no necesaria, para asegurar que el problema sea convexo. Un ejemplo muy común e importante son las formas cuadráticas, sabemos que son convexas y sin embargo no verifican las reglas libres de productos. Para poder resolver problemas con formas cuadráticas se incluye en la librería la función square(x) que significa x2. Es bastante sencillo encontrar problemas que son convexos pero no DCPs. Si un problema no verifica las reglas de convexidad no significa que no se pueda resolver en el marco de programación convexa disciplinada, a veces simplemente hay que reescribirlo. Por ejemplo, vamos a considerar el problema de maximizar la entropía: Maximizar −Pn i=1 xilogxi Sujeto a: Ax =b 1Tx= 1 x≥0 (3.6) donde x∈Rn es la variable del problema y A∈Rm×n , b∈Rm son parámetros; y log ( . ) está definida en la librería. La expresión xilogxi no verifica la regla libre de productos PN7, ya que contiene un producto entre dos expresiones no constantes. Luego este problema no es DCP válido, sin embargo sabemos que es un problema convexo. La forma más sencilla de resolverlo en el marco de la programación convexa disciplinada es agregando una nueva función a la librería. En este caso necesitaríamos añadir la siguente función: entr ( x ) = −xlogx para x > 0. Es una función cóncava, con signo y monotonía desconocidos. A partir de esta función, el problema se puede reescribir de la siguiente manera como un problema dDCP válido: Maximizar −Pn i=1 entr(xi) Sujeto a: Ax =b 1Tx= 1 x≥0 (3.7) Una vez que se define e implementa una función, se puede reutilizar libremente para otros problemas y se pueden compartir con otros usuarios. El esfuerzo invertido en agregar una nueva función a la librería se amortiza significativamente. 3.4. Implementación de las funciones básicas La implementación es una descripción fácil de usar de una función o conjunto que permite su uso en algoritmos de búsqueda de solución numérica óptima. Es importante tener en cuenta que no se utiliza para determinar si un problema en particular es o no DCP válido, entra en juego solo después de la verificación, al calcular una solución numérica. Existen una gran variedad de métodos o solvers que se pueden emplear para resolver problemas convexos: métodos primal/dual, métodos de barrera, métodos de plano de corte, etc. No los podemos tratar todos en detalle en este trabajo, pero para más información se puede consultar ([33]). Todos estos métodos pueden adaptarse a la programación CAPÍTULO 3. PROGRAMACIÓN CONVEXA DISCIPLINADA 31 4.2. FUNCIONES 4.2. Funciones Vamos a decribir las funciones que pueden ser aplicadas en CVXR. Operadores Los operadores + ,−,∗, % ∗ %, / se tratan como funciones. + y - son funciones afines. ∗ y / son afines en CVXR porque expr1 ∗ expr2 y expr1% ∗ % expr2 solo se permiten cuando una de las expresiones es constante y expr1 / expr2 solo se permite cuando expr2 es una constante escalar. Indexación y corte Todas las expresiones no escalares se pueden indexar utilizando la sintaxis “expr[i, j]”. La indexación es una función afín. La sintaxis “expr[i]” se puede usar como una abreviatura para “expr[i, 1]” cuando expr es un vector columna. De manera similar, “expr[i]” es una abreviatura de “expr[1, i]” cuando expr es un vector fila. Las expresiones no escalares también se pueden dividir utilizando la sintaxis de indexado de R estándar. Por ejemplo, “expr[i: j, r]” selecciona las filas i a j de la columna r y devuelve un vector. CVXR admite la indexación avanzada utilizando listas de índices o matrices booleanas. La semántica es la misma que en R. Siempre que R devuelva un vector numérico, CVXR devuelve un vector columna. Transponer La transposición de cualquier expresión también es una función afín. Se puede obtener utilizando la sintaxis t(expr). 4.2.1. Funciones escalares Una función escalar toma uno o más escalares, vectores o matrices como argumentos y devuelve un escalar. 38 CAPÍTULO 4. CVXR 4.2. FUNCIONES Cuadro 4.1: Funciones escalares disponibles en el paquete CVXR Función Significado Dominio Signo Curvatura Monotonía geo_mean(x,p) p∈Rn +, p 6= 0 x1/n 1···x1/n n(xp1 1··· xpn n) 1 1Tp x∈Rn +Positivo Cóncava ↑ harmonic_mean(x) n 1 x1+...+1 xn x∈Rn +Positivo Cóncava ↑ lambda_max(X) λmin(X)X∈SnDesconocido Convexa Ninguna lambda_min(X) λmin(X)X∈SnDesconocido Cóncava Ninguna lambda_sum_largest(X,k) k=1,...,n suma de los k mayores autovalores de X x∈SnDesconocido Convexa Ninguna lambda_sum_smallest(X,k) k=1,...,n suma de los k menores autovalores de X X∈SnDesconocido Cóncava Ninguna log_det(X) log(det(X)) X∈Sn +Desconocido Cóncava Ninguna log_sum_exp(X) C∈Rp×q log(Pij exij )X∈Rm×nDesconocido Convexa ↑ matrix_frac(x,P) xTP−1x X ∈Rn P∈Sn ++ Positivo Convexa Ninguna max_entries(X) maxij{Xij}X∈Rm×nSigno de X Convexa ↑ min_entries(X) minij{Xij}X∈Rm×nSigno de X Cóncava ↑ mixed_norm(X,p,q) (Pk(Pl|xk,l|p)q/p)1/q X∈Rm×nPositivo Convexa Ninguna cvxr_norm(x) cvxr_norm(x,2) pPix2 ix∈RnPositivo Convexa ↑para xi≥0, ↓para xi≤0 cvxr_norm(x,"fro") qPij x2 ij X∈Rm×nPositivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 cvxr_norm(x,1) Pij |xij|X∈Rm×nPositivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 cvxr_norm(x,"inf") maxij{|xij |} X∈Rm×nPositivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 cvxr_norm(x,"nuc") tr((XTX)1/2)X∈Rm×nPositivo Convexa Ninguna cvxr_norm(X) cvxr_norm(X,2) pλmax(XTX)X∈Rm×nPositivo Convexa Ninguna p_norm(X,p) p≥1op=∞kXkp= (Pij |xij|p)1/p X∈Rm×nPositivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 p_norm(X,p) p < 1op6= 0 kXkp= (Pij xp ij)1/p X∈Rm×n +Positivo Cóncava ↑ quad_form(x,P) p∈Sn +constante xTPx X ∈RnPositivo Convexa ↑para xi≥0, ↓para xi≤0 quad_form(x,P) p∈Sn −constante xTPx X ∈RnNegativo Cóncava ↑para xi≥0, ↓para xi≤0 quad_form(c,X) c∈Rnconstante cTXc X ∈Rn×nDepende de c,X Afín Depende de c quad_over_lin(X,y) (Pij x2 ij)/y x ∈Rn +Positivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 ↓en y sum_entries(X) Pij xij x∈Rn +Signo de X Afín ↑ sum_largest(X,k) k= 1,2, ... Suma de los k mayores Xij X∈SnSigno de X Convexa ↑ sum_smallest(X,k) k= 1,2, ... Suma de los k menores Xij X∈SnSigno de X Cóncava ↑ sum_squares(X) Pij x2 ij X∈Rm×nPositivo Convexa ↑para Xij ≥0, ↓para Xij ≤0 matrix_trace(X) tr(X) X∈SnSigno de X Afín ↑ tv(x) Pi|xi+1 −xi|X∈RnPositivo Convexa Ninguna tv(X) Pij ||hXi+1,j −Xij Xi,j+1 −Xij i||2X∈Rm×nPositivo Convexa Ninguna tv(X1, ..., Xk)Pij ||        X(1) i+1,j −X(1) ij X(1) i,j+1 −X(1) ij . . . X(k) i+1,j −X(k) ij X(k) i,j+1 −X(k) ij        ||2X(i)∈ Rm×n Positivo Convexa Ninguna CAPÍTULO 4. CVXR 39 4.2. FUNCIONES Aclaraciones El dominio Sn se refiere al conjunto de matrices simétricas. Los dominios Sn + y Sn − se refieren al conjunto de matrices semi-definidas positivas y semi-definidas negativas, respectivamente. De manera similar, Sn ++ y Sn −− se refieren al conjunto de matrices definidas positivas y negativas, respectivamente. Para una expresión vectorial x, “cvxr_norm(x)” y “cvxr_norm(x, 2)” dan la norma euclídea ( kxk2 = qx2 1+x2 2+... +x2 n ), mientras que para una expresión matricial X “cvxr_norm(X)” y “cvxr_norm(X, 2)” dan la norma espectral ( kAk2 = qλmax(A∗A) , donde A∗es la matriz traspuesta conjugada de A). La función “cvxr_norm(X,”fro“)” se denomina norma de Frobenius ( kAkF = tr ( A∗A ) 1/2 = ( Pm i=1 Pn j=1 |Ai,j|2 ) 1/2 ) y cvxr_norm (X, “nuc”) la norma nuclear( kAkN = tr (( ATA ) 1/2 )). La norma nuclear también se puede definir como la suma de los valores singulares de X. Las funciones max_entries y min_entries dan la entrada más grande y más pequeña, respectivamente, en una sola expresión. Estas funciones no deben confundirse con max_elemwise y min_elemwise que dan el máximo o el mínimo de una lista de expresiones escalares. La función sum_entries suma todas las entradas en una sola expresión, mientras que “sum” suma una lista de expresiones. Funciones a lo largo de un eje Las funciones sum_entries, cvxr_norm, max_entries y min_entries se pueden aplicar a lo largo de un eje. Dada una expresión expr de dimensión (m,n) la sintaxis func(expr, axis = 1) aplica func a cada fila, devolviendo una expresión de dimensión (m,1). La sintaxis func(expr, axis = 2) aplica func a cada columna, devolviendo una expresión de dimensión (1,n). En el siguiente ejemplo se construyen dos expresiones, en la primera se suman las variables de la matriz por filas y en la segunda por columnas. library(CVXR, warn.conflicts=FALSE) X <- Variable(5,4) row_sums <- sum_entries(X, axis=1)# Tiene dimensión (5, 1) col_sums <- sum_entries(X, axis=2)# Tiene dimensión (1, 4) 40 CAPÍTULO 4. CVXR 4.2. FUNCIONES 4.2.2. Funciones element-wise Estas funciones operan en cada elemento de sus argumentos. Por ejemplo, si X es una matriz de variables con dimensión (5,4), entonces “abs(X)” es una matriz de dimensión (5 ,4), es decir, “abs(X)[1, 2]” es equivalente a “abs(X [1, 2])”. Las funciones elementwise que toman múltiples argumentos, como “max_elemwise” y “mul_elemwise”, operan en los elementos correspondientes de cada argumento. Por ejemplo, si X e Y son matrices de variables de dimensión (3,3), entonces “max_elemwise (X, Y)” es una matriz 3 por 3. “max_elemwise (X, Y) [2, 1]” es equivalente a “max_elemwise (X [2, 1], Y [2, 1])”. Esto significa que todos los argumentos deben tener las mismas dimensiones o ser escalares. Cuadro 4.2: Funciones element-wise disponibles en el paquete CVXR Función Significado Dominio Signo Curvatura Monotonía abs(x) |x|x∈RPositivo Convexa ↑para x≤0, ↓para x≥0 entr(x) −xlog(x)x > 0Desconocido Cóncava Ninguno exp(x) exx∈RPositivo Convexa ↑ huber(x,M), M≥0exx∈RPositivo Convexa ↑para x≤0, ↓para x≥0 inv_pos(x) 1/x x > 0Positivo Convexa ↓ kl_div(x,y) xlog(x/y)−x+y x > 0,y > 0Positivo Covexa Ninguno log(x) log(x)x > 0Desconocido Cóncava ↑ log1p(x) log(x+ 1) x > −1Signo de x Cóncava ↑ logistic(x) log(x+ex)x∈RPositivo Convexa ↑ max_elewise( x1 , ..., xk) max{x1, ..., xk}xi∈Rmax(sign(xi)) Convexa ↑ min_elewise( x1 , ..., xk) min{x1, ..., xk}xi∈Rmin(sign(xi)) Cóncava ↑ mul_elewise( x1 , ..., xk)c∈R c∗x x ∈Rsign(cx)Afín Depende de c neg(x) max{−x, 0}x∈RPositivo Convexa ↓ pos(x) max{x, 0}x∈RPositivo Convexa ↑ power(x,0) 1 x∈RPositivo Constante constante power(x,1) x x∈RLa misma de x Afín ↑ power(x,p) p=2,4,8,... xpx∈RPositivo Convexa ↑para x≤0 ↓para x≥0 power(x,p) p<0 xpx > 0Positivo Convexa ↓ power(x,p) 0<p<1 xp≥0Positivo Cóncava ↑ power(x,p) p>1 p6= 1,4,8, ... xpx≥0Positivo Convexa ↑ scalene(x, α, β)α≥0β≥ 0 αpos(x) + βneg(x)x∈RPositivo Convexa ↑para x≤0 ↓para x≥0 sqrt(x) √x x ≥0Positivo Cóncava ↑ square(x) x2x∈RPositivo Convexa ↑para x≤0 ↓para x≥0 CAPÍTULO 4. CVXR 41 4.2. FUNCIONES 4.2.3. Funciones vector / matriz Una función de vector / matriz toma uno o más escalares, vectores o matrices como argumentos y devuelve un vector o matriz. Cuadro 4.3: Funciones vector/matriz disponibles en el paquete CVXR Función Significado Dominio Signo Curvatura Monotonía bmat([[X11(, ..., X1q] , ..., [Xq1, ..., Xqq]])   X(1,1) ... X(1,q) . . .. . . X(p,1) ... X(p,q)  x(i,j)∈ Rmi×nj signo(Pij X(i,j) 11 )Afín ↑ conv(c,x) , c∈R c*x X∈RnSigno(c1x1)Afín Depende de C cumsum_axis(X,axis=1) Suma acumulada a lo largo del eje dado X∈Rm×nSigno de X Afín ↑ diag(x), M≥0      X1 ... Xn      X∈RnSigno de X Afín ↑ diag(X)   X1 . . . Xn  X∈Rn×nSigno de X Afín ↑ diff(X,k=1,axis=1) k∈0,1,2, ... Diferencias de orden k (el argumento k se denomina diferencias y también se puede usar lag a lo largo del eje dado X∈Rm×nSigno de X Afín ↑ hstack(X1, ..., Xk)X(1) ... X(k)X(i)∈Rm×nisigno(PiX(i) 11 )Afín ↑ kronecker(C,X) C∈Rp×q   C11X ... C1qX . . .. . . Cp1X ... CpqX  X∈Rm×nSigno(c11x11)Afín Depende de C reshape_expr(X, n0, m0)X’∈Rm0×n0X∈Rm×n m0n0=mn Signo de X Afín ↑ vec(X) x’∈Rmn X∈Rm×nSigno de X Afín ↑ vstack(X1, ..., Xk)  X(1) . . . X(n)  X(i)∈Rmi×nsigno(PiX(i) 11 )Afín ↑ Aclaraciones La entrada de bmat es una lista de listas de expresiones CVXR y construye una matriz de bloques. Los elementos de cada lista interna se apilan horizontalmente, y luego las matrices de bloques resultantes se apilan verticalmente. La salida de conv(c,x) es de tamaño n+m-1 y se define como y[k] = Pk j=0 c[j]x[k−j]. La salida de vec (X) es la matriz X aplanada en orden de columna principal en un vector. Formalmente, x0 i=Xi mod m,[i/m] La salida X de reshape_expr (X, m’, n’) es la matriz X convertida en una matriz m’xn’. Las entradas se toman de X por columnas y se almacenan en X’ por columnas también. Formalmente, X0 ij =vec(X)m0j+i. 42 CAPÍTULO 4. CVXR 4.3. CONSTRUIR UN PROBLEMA 4.3. Construir un problema Empecemos con un problema simple. Sean x, y, z variables escalares, definimos: Minimizar x+ 2y−0.1z Sujeto a: 2x−y≥0 −x+ 3y≥0 x+y≤5 z≥1.1 x, y, z ≥0 z∈Z (4.1) Es un problema convexo, concretamente lineal, con f ( x, y ) = x + 2 y− 0 . 1 z función objetivo y con las restricciones f1 ( x ) = − 2 x + y , f2 ( x, y ) = x− 3 y , f3 ( x, y ) = x + y− 5, f4(x, y)=1.1−z,f5(x, y) = −x,f6(x, y) = −y,f7(x, y) = −z,z∈Z. Vamos a describir detalladamente el código CVXR necesario para resolverlo. Se carga la librería CVXR y empezamos por declarar las variables. Podemos imponer que alguna de las variables sea entera mediante Int(): library(CVXR, warn.conflicts=FALSE) x <- Variable() y <- Variable() z <- Int() Ahora se define el objetivo, se utiliza Maximize() para funciones cóncavas y Minimize() para funciones convexas. En este caso empleamos Maximize(): objetivo <- Maximize(x +2*y-0.1*z) A continuación vamos a definir las restricciones de igualdad y desigualdad: restricciones <- list(x>=0, y>=0, z>=0, x+y<=5,2*x-y>=0, -x+3*y>=0,z>=1.1) Entonces, utilizando el objetivo y las restricciones como objetos de entrada, ya podemos definir el problema mediante la función “Problem()”: problema <- Problem(objetivo,restricciones) Los objetos problema son bastante flexibles, pueden tener restricciones o no. Una vez definido el problema, comprobamos si es dcp. is_dcp(problema) ## [1] TRUE Por último, para resolver el problema se utiliza la función solve: solucion <- solve(problema) Esta llamada traduce el problema a un formato que un solver convexo puede entender, pasa el problema al solver y luego devuelve los resultados en una lista. Veamos como recuperar estos resultados: CAPÍTULO 4. CVXR 43 4.3. CONSTRUIR UN PROBLEMA solucion$status ## [1] "optimal" Devuelve el estado del problema, que en este caso es óptimo por lo que devuelve “optimal”. Es decir, se ha encontrado solución óptima para el problema. Se encontraron valores para las variables de decisión que satisfacen todas las restricciones y maximizan (en otros casos minimizan) el objetivo. En general, al aplicar el metodo solve() puede ocurrir que el problema sea óptimo como este, pero también nos podemos encontrar otros dos casos: Si solucion$status == “infeasible”, el problema no se ha resuelto porque no existe una combinación de variables de decisión que pueda satisfacer todas las restricciones. Si solucion$estado == “unbounded”: el problema no se resolvió porque el objetivo se puede hacer arbitrariamente pequeño para un problema de minimización o arbitrariamente grande para un problema de maximización. Por lo tanto, no hay una solución óptima porque, para cualquier solución dada, siempre es posible encontrar otra aún mejor. solucion$value ## [1] 8.133333 Devuelve el valor óptimo de la función objetivo. Si el problema es infactible devuelve +Inf en los problemas de mínimos y -Inf en los problemas de máximos. Si el problema es ilimitado, devuelve -Inf para un problema de minimización y + Inf para un problema de maximización, solucion$getValue(x) ## [1] 1.666667 solucion$getValue(y) ## [1] 3.333333 solucion$getValue(z) ## [1] 2 Devuelve los valores de las variables para los que se alcanza el óptimo. Si el problema es infactible o ilimitado devuelve NA. 4.3.1. Problema vectorial Vamos a ver ahora como construir un problema con variables vectoriales. Consideremos el siguiente problema. Minimizar P3 i=1 xi Sujeto a: Ax =b x≥0 con A=   4 2 4 1 2 3 1 1 2   , b=   4 4 4   . 44 CAPÍTULO 4. CVXR 4.4. IMPLEMENTACIÓN El procedimiento será el mismo que en el caso anterior, empezamos declarando la variable. library(CVXR, warn.conflicts=FALSE) X <- Variable(3,1) Ahora se define el objetivo, en este caso vamos a emplear “Minimize()” y la función de la librería “sum_entries()”: objetivoVec <- Minimize(sum_entries(X)) Se definen las restricciones: A = matrix(c(4,2,4,1,2,3,1,1,2), nrow=3) b = matrix(c(4,4,4), nrow=3) restriccionesVec <- list(A%*%X==b, X>=0) A partir del objetivo y las restricciones podemos definir el problema. problemaVec <- Problem(objetivoVec,restriccionesVec) Comprobamos que es dcp. is_dcp(problemaVec) ## [1] TRUE Por último, resolvemos mediante solve. solucionVec <- solve(problemaVec) solucionVec$status ## [1] "infeasible" El estado del problema es infactible. solucionVec$value ## [1] Inf Como es un problema de mínimo y es infactible devuelve +Inf. solucionVec$getValue(X) ## [1] NA Vemos que al ser el problema infactible devuelve NA. 4.4. Implementación CVXR representa las funciones, variables, restricciones y otras partes de un problema de optimización utilizando objetos R de clase S4. S4 nos permite sobrecargar operaciones matemáticas estándar para que CVXR se combine perfectamente con el script R original y otros paquetes. Cuando se invoca una operación en una variable, se crea un nuevo objeto que representa el árbol de expresión correspondiente con el operador como nodo raíz y CAPÍTULO 4. CVXR 45 4.4. IMPLEMENTACIÓN los argumentos como hojas. Este árbol crece automáticamente a medida que se agregan más elementos, lo que nos permite encapsular la estructura de una función objetivo o restricción. Una vez que el usuario llama a la orden solve, ocurre la verificación DCP. CVXR recorre el árbol de la expresión de forma recursiva, determinando el signo y curvatura de cada subexpresión basándose en las propiedades de las funciones básicas que las componen. Si se considera que el problema se ajusta a las reglas, es transformado en un programa de cono equivalente utilizando implementaciones gráficas de funciones convexas. Luego, CVXR pasa la descripción del problema a la librería CVXcanon C ++ [49], que genera los datos para el programa de cono y envía estos datos a la interfaz de R específica del solver. Los resultados del solver se devuelven al usuario en una lista. Este diseño e infraestructura orientados a objetos de CVXR se tomaron en gran parte de CVXPY (versión escrita en lenguaje Phyton). CVXR interactúa con los solvers de conos de código abierto ECOS [22] y SCS [55] a través de sus respectivos paquetes R. ECOS es un solver de punto interior, que logra una alta precisión para problemas pequeños y medianos, mientras que SCS es un solver de primer orden que es capaz de manejar problemas más grandes y restricciones semidefinidas. Ambos solvers ejecutan un solo hilo a la vez, el soporte para SCS multihilo se agregará en el futuro. Además, la versión 0.99 de CVXR proporciona soporte para los solvers comerciales MOSEK [4] y GUROBI [37] a través de paquetes binarios R publicados por los proveedores respectivos. No es difícil conectar solvers adicionales siempre que el solver tenga una API que pueda comunicarse con R. Los usuarios que deseen emplear un solver personalizado pueden obtener la combinación de los datos canonizados para un problema y solver directamente con get_problem_data (problema, solver). Cuando un problema puede ser resuelto por más de un solver, el argumento del solver de la función de resolución puede usarse para indicar una preferencia. Los solvers disponibles, según los paquetes instalados en una sesión, se devuelven a través de installed_solvers (). Si ejecutamos dicha función, installed_solvers () ## [1] "ECOS" "ECOS_BB" "SCS" obtenemos los solvers mencionados, ECOS y SCS, y además aparece ECOS_BB, que es una extensión reciente de ECOS para resolver programas enteros mixtos o booleanos mixtos. También podemos obtener qué solver ha utilizado para resolver un cierto problema. Podemos verlo utilizando el ejemplo aterior, solucion$solver ## [1] "ECOS_BB" Veamos ahora como podemos elegir el solver usado para resolver el problema, Si queremos indicar ECOS_BB (en este caso es el que ha usado). solucionECOS <- solve(problema,"ECOS_BB") Si le indicamos que use, por ejemplo ECOS, 46 CAPÍTULO 4. CVXR 4.5. MÍNIMOS CUADRADOS # solucionECOS <- solve(problema,"ECOS") Nos devuelve error porque el solver ECOS no es válido para problemas mixtos enteros. Además podemos pedir que nos devuelva el tiempo empleado en cargar y resolver el problema. solucion$setup_time ## [1] 2.559999e-05 solucion$solve_time ## [1] 1.024e-05 Tenemos una gran librería de funciones que deberían ser suficiente para modelar la mayoría de los problemas de optimización convexa. Sin embargo, es posible añadir nuevas funciones a la librería. El proceso conlleva crear en R una clase S4 para la función, la sobrecarga de métodos que caracterizan sus propiedades DCP y representar su implementación gráfica como una lista de operadores lineales que especifican la factibilidad del problema. 4.5. Mínimos cuadrados Para mostrar algunas de las ventajas de CVXR sobre otros paquetes de R vamos a tomar como ejemplo el problema de mínimos cuadrados, el cual puede ser resuelto en R utilizando lm. Minimizar ky−Xβ k2 2(4.2) Se traduce a CVXR directamente como se muestra a continuación: # Generar los datos: n <- 100 p <- 10 beta <- -4:5 x <- matrix(rnorm(n*p), nrow = n) colnames(x) <- paste("beta_",beta) y <- x%*%beta+rnorm(n) # Resolver el problema usando CVXR library(CVXR, warn.conflicts=FALSE) # Definir la variable y el objetivo betah <-Variable(p) objetivo <- Minimize(sum((y-x%*%betah)^2)) CAPÍTULO 4. CVXR 47 4.6. APLICACIONES Cuadro 4.8: Primeros valores de la solución del probema de la catenaria con restricciones sobre el suelo. x y 0.0000000 1.0000000 0.0025235 0.9801594 0.0051234 0.9603286 0.0078044 0.9405084 0.0105718 0.9207003 0.0134311 0.9009051 restricciones3[[4]] <- (y[m] == 0.5) prob3 <- Problem(objetivo, restricciones3) result3 <- solve(prob3) xs2 <- result3$getValue(x) ys2 <- result3$getValue(y) 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 x y Figura 4.3: Representación gráfica del problema de la Catenaria con restricciones sobre el suelo. En el cuadro 4.8 tenemos los primeros valores de la solución obtenida y en la gráfica 4.3 podemos observar su representación. 54 CAPÍTULO 4. CVXR 4.6. APLICACIONES 4.6.2. Regresión logística Supongamos que queremos resolver un modelo de predicción para una variable binaria yi∈ { 0 , 1 } , utilizamos la regresión logística (Cox 1958 [16]; Freedman 2009 [25]) . Modelamos la respuesta condicional como y|x∼Bernoulli ( gβ ( x )), donde gβ ( x ) = 1 1+e−xTβ es la función logística, y maximizamos el logaritmo de la función de probabilidad, generando el problema de optimización: Minimizar m X i=1 (yilog(gβ(xi)) + (1 −yi)log(1 −gβ(xi))) (4.4) Podemos simplificar esta expresión: yilog(gβ(xi)) = yilog(1 1+e−xT iβ) = −yilog(1 + e−xT iβ) (1 −yi)log(1 −gβ(xi)) = (1 −yi)log(1 −1 1+e−xT iβ) = (1 −yi)log(1+e−xT iβ−1 1+e−xT iβ) = = (1 −yi ) xT iβ− (1 −yi ) log (1+ e−xT iβ ) = (1 −yi ) xT iβ−log (1+ e−xT iβ )+ yilog (1+ e−xT iβ ) Finalmente nos queda: Minimizar m X i=1 (1 −yi)xT iβ−log(1 + e−xT iβ)(4.5) Lo natural sería escribir directamente en la sintaxis convencional de R la función log (1+ exp (X%*%beta)). Pero esta representación no es DCP válida, viola la regla de composición C4. La función log(x) es cóncava no decreciente, para verificar la regla el argumento debería ser cóncavo y, sin embargo, exp(x) es convexa. Por tanto, aunque el objetivo sea convexo, CVXR rechazará el problema. Se necesita encontrar otra escritura del problema. CVXR proporciona la función logistic como camino directo para expresar f(z) = log(1 + ez). Entonces nos quedaría Pm i=1(1 −yi)xT iβ−logistic(−xT iβ). Entonces nuestro problema completo es expresado en CVXR como: library(CVXR) # Crear una muestra con 20 variables y 1000 observaciones. n <- 20 m <- 1000 sigma <- 45 DENSITY <- 0.2 set.seed(183991) # Se genera una muestra normal de tamaño n beta_true <- stats::rnorm(n) # Se toma una muestra de tamaño (1-DENSITY)*n de numeros entre 1 y n # sin repetición idxs <- sample(n, size = floor((1-DENSITY)*n), replace = FALSE) CAPÍTULO 4. CVXR 55 4.6. APLICACIONES # Hace 0 el valor las componentes dados por idxs beta_true[idxs] <- 0 # Se crea una matriz m*n cuyas entradas vienen dadas por una muestra # normal X <- matrix(stats::rnorm(m*n, 0,5), nrow = m, ncol = n) # Genera los datos y: y=x*beta+eps con eps~N(0,sigma) y <- sign(X %*% beta_true +stats::rnorm(m, 0, sigma)) # Definir el problema beta <- Variable(n) obj <- -sum(X[y <= 0, ] %*% beta) -sum(logistic(-X%*% beta)) prob <- Problem(Maximize(obj)) Podemos observar que para la función objetivo se ha empleado la función logistic, como se ha explicado anteriormente, pero además no se ha escrito de forma directa (1 −yi ) xT iβ . Hemos definido y a partir de la función signo, luego toma valores en {-1,0,1}. Necesitamos que y sea binaria y que tome únicamente valores en {0,1}, para ello cuando yi = − 1vamos a considerar que vale cero. La expresión (1 −yi ) xT iβ es 0 si yi = 1 y es xT iβ para yi = 0, luego lo que hacemos es tomar las xital que yi≤0, es decir, tal que yi= 0 oyi=−1. # Comprobar que el problema es convexo is_dcp(prob) ## [1] TRUE #Resolver result <- solve(prob) beta_log <- result$getValue(beta) #Podemos ver que solver ha empleado en este caso result$solver ## [1] "ECOS" Ahora vamos a resolver el problema usando gml para comparar los resultados: d <- data.frame(y = as.numeric(y >0), X=X) glm <- stats::glm(formula = y~0+X, family = "binomial",data = d) En el cuadro 4.9 tenemos los valores de la solución obtenida por CVXR y glm. 56 CAPÍTULO 4. CVXR 4.6. APLICACIONES Cuadro 4.9: Resultados del problema Logístico mediante CVXR y glm CVXR.est GLM.est beta_1 -0.0003345 -0.0003345 beta_2 -0.0149849 -0.0149849 beta_3 0.0144564 0.0144564 beta_4 0.0172645 0.0172645 beta_5 0.0206029 0.0206029 beta_6 -0.0200824 -0.0200824 beta_7 -0.0137066 -0.0137066 beta_8 -0.0009806 -0.0009806 beta_9 -0.0246409 -0.0246409 beta_10 0.0101345 0.0101345 beta_11 -0.0113936 -0.0113936 beta_12 0.0187613 0.0187613 beta_13 0.0051938 0.0051938 beta_14 0.0015573 0.0015573 beta_15 0.0202801 0.0202801 beta_16 -0.0176528 -0.0176528 beta_17 0.0191530 0.0191530 beta_18 0.0131130 0.0131130 beta_19 -0.0265186 -0.0265186 beta_20 0.0039912 0.0039912 CAPÍTULO 4. CVXR 57 Conclusiones En este trabajo se ha expuesto una metodología para la optimización convexa denominada Programación Convexa Disciplinada. En el capítulo 2 tenemos una larga lista de aplicaciones de la programación convexa, pero existen razones para creer que la cantidad de aplicaciones aumentará drásticamente si se puede hacer más fácil el hecho de crear, analizar y resolver CPs. Una fuente prometedora de nuevas aplicaciones para la programación convexa es la extensión y mejora de las aplicaciones existentes para la programación lineal, y si nos fijamos en su historia vemos que no experimentó un gran desarrollo hasta que se hicieron grandes avances computacionales que facilitaban su resolución. Por todo esto, la importancia de este método se debe a que en muchos casos verificar si un problema dado es convexo puede ser un trabajo intratable, pero utilizando DCP se puede llevar a cabo de forma automática. En el capítulo 3 hemos descrito las normas que esta metodología impone sobre la construcción de un problema para así, simplificar el proceso de verificar si es convexo. De forma que la verificación de la convexidad consistirá simplemente en probar que el problema está construido a partir de funciones de la librería de acuerdo a las normas establecidas. Es necesario destacar que el cumplimiento de estas reglas es una condición suficiente de convexidad, no necesaria, es decir, la programación Convexa Disciplinada no se puede utilizar para rechazar que un problema sea convexo. Además, ofrece otras ventajas como el hecho de que, mediante las implementaciones gráficas, permite trabajar con funciones no diferenciables sin la pérdida de rendimiento que normalmente va asociada a ellas, como hemos podido comprobar en el capítulo 3 mediante un ejemplo con la función mínimo. Otra de sus principales ventajas es su generalidad, puede trabajar con cualquier tipo de problema, no se restringe a un tipo determinado de formas estándar como hacen otras metodologías. Esta generalidad no se pierde por el hecho de que existan un conjunto de reglas determinado, ya que hemos visto que la librería en la que se basa es extensible. En el último capítulo hemos llevado esta metodología a la práctica a través del paquete CVXR de R. Como se ha visto en los ejemplos y aplicaciones planteados, dado un problema, CVXR emplea la programación convexa disciplinada previamente a resolver, para probar que sea convexo. Una vez probada la convexidad, utiliza las implementaciones gráficas para transformar el problema y lo pasa a un solver, como ECOS o SCS, para resolverlo. Permite formular los problemas de una forma muy natural e intuitiva, utilizando la amplia lista de funciones que tiene implementada. Además, es posible formular y resolver problemas para los que no exista un código personalizado. Como se ha mostrado en la 59 4.6. APLICACIONES resolución del problema de mínimos cuadrados, puede no ser el paquete óptimo para resolver ciertos problemas comunes para los cuales ya exista un paquete exclusivo, ya que CVXR probablemente sea más lento que un algoritmo personalizado y serían necesarias más líneas de código, pero estos paquetes están limitados. Si necesitaramos añadir nuevas restricciones CVXR sería el paquete ideal, debido a su flexibilidad. Como hemos visto en las aplicaciones, es muy sencillo incluir nuevas restricciones o transformar un problema añadiendo pocas líneas de código. 60 CAPÍTULO 4. CVXR Bibliografía [1] Achtziger, W., Bendsoe, M., Ben-Tal, A. and J Zowe 1992. Equivalent displacement based formulations for maximum strength truss topology design. Impact of Computing in Science and Engineering, [2] Alizadeh, F. 1995. Interior point methods in semidefinite programming with applications to combinatorial optimization. SIAM Journal on Optimization. [3] Allaire, J., Xie, Y., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., Wickham, H., Cheng, J., Chang, W. and Iannone, R. 2018. Rmarkdown: Dynamic documents for r. [4] Andersen, E. and Andersen, K. 2000. The mosek interior point optimizer for linear programming: An implementation of the homogeneous algorithm. In High Performance Optimization. [5] Avriel, M., Dembo, R. and Passy, U. 1975. Solution of generalized geometric programs. International Journal for Numerical Methods in Engineering. [6] Ben-Tal, A. and Bendsoe, M. 1993. A new method for optimal truss topology design. SIAM J. Optim. [7] Bendsoe, M., A Ben-Tal and Zowe, J. 1994. Optimization methods for truss geometry and topology design. Structural Optimization. [8] Benson, S. 2002. DSDP 4.5: A daul scaling algorithm for semidefinite programming. Web site: http://www-mix.mcs.an1.gov/benson/dsdp/. [9] Bertsimas, D. and Nino-Mora, J. 1999. Optimization of multiclass queuing networks with changeover times via the achievable region approach: Part ii, the multistation case. Mathematics of Operations Research. [10] Biswas, P. and Ye, Y. 2004. Semidefinite programming for ad hoc wireless sensor network localization. Technical report, Stanford University. [11] Borchers, B. 1999. CDSP, a c library for semidefinite programming. Optimization Methods and Software. [12] Boyd, S. and Vandenberghe, L. 2004. Convex optimization. Cambridge University Press. [13] Brooke, A., Kendrick, D., Meeraus, A. and Raman, R. 1998. GAMS: A user’s guide. The Scientific Press. Web site: http: //www . gams . com/docs/gams/GAMSUsersGuide.pdf. [14] Calafiore, G. and Indri, M. 2000. Robust calibration and control of robotic manipulator. In American Control Conference. [15] Chinneck, J. 2003. MProbe 5.0 (software package). Web site: http://www.sce.carleton.ca/ 61 4.6. APLICACIONES Bibliografía faculty/chinneck/mprobe.html. [16] Cox, D.R. 1958. The regression analysis of binary sequences. Journal of the Royal Statistical Society. [17] Crusius, C. 2002. A parser/solver for convex optimization problems. PhD thesis, Stanford University. [18] Dantzig, G.B. 1963. Linear programming and extensions. Princeton University Press. [19] Davidson, T., Luo, Z. and Wong, K. 2000. Design of orthogonal pulse shapes for communications via semidefinite programming. IEEE Transactions on Communications, [20] Diamond, S. and Boyd, S. 2016. CVXPY: A python-embedded modeling language for convex optimization. Journal of Machine Learning Research. [21] Doherty, A., Parrilo, P. and Spedalieri, F. 2002. Distinguishing separable and entangled states. Physical Review Letters. [22] Domahidi, A., Chu, E. and Boyd, S. 2013. ECOS: An socp solver for embedded systems. In Proceedings of the European Control Conference. [23] Dufin, R. 1970. Linearizing geometric programs. SIAM Review. [24] Fourer, R., Gay, D. and Kernighan, B. 1999. AMPL: A modeling language for mathematical programming. Duxbury Press. [25] Freedman, D.A. 2009. Statistical models: Theory and practice. Cambridge University Press. [26] Fu, Anqi, Narasimhan, B. and Boyd, S. 2017. CVXR: An r package for disciplined convex optimization. [27] Fujisawa, K., Kojima, M., Nakata, K. and Yamashita, M. 2002. SDPA (semi-definite programming algorithm) user’s manual-version 6.00. Technical report, Tokyo Insitute of Technology. [28] Gelfand, I.M. and Fomin, S.V. 1963. Calculus of variations. Prentice-Hall. [29] GenBeta https://www.genbeta.com/desarrollo/r-un-lenguaje-y-entorno-de- programacion-para-analisis-estadistico. [30] Ghaoui, L.E., Commeau, J.-L., Delebecque, F. and Nikoukhah, R. 1999. LMITOOL 2.1 (software package).". Web site: http://robotics.eecs.berkeley.edu/ elghaoui/lmitool/lmitool.html. [31] Giiler, 0. and Hauser, R. 2002. Self-scaled barrier functions on symmetric cones and their classification. Foundations of Computational Mathematics. [32] Goldfarb, D. and Iyengar, G. 2002. Robust portfolio selection problems. Technical report, Computational Optimization Research Center, Columbia University. [33] Grant, M. 2004. Disciplined convex programming. PhD thesis, Department of Electrical Engineering, Stanford University. [34] Grant, M. and Boyd, S. 2014. CVX: MATLAB software for disciplined convex 62 CAPÍTULO 4. CVXR Bibliografía 4.6. APLICACIONES programming, version 2.1. [35] Grant, M., Boyd, S. and Ye, Y. 2006. Disciplined convex programming. Global Optimization: From Theory to Implementation, Nonconvex Optimization and its Applications. [36] Griva, I.A. and Vanderbei, R.J. 2005. Case studies in optimization: Catenary problem. Optimization and Engineering. [37] Gurobi Optimization, I. 2016. Gurobi optimizer reference manual. Disponible en http://www.gurobi.com. [38] Hershenson, M., Mohan, S., Boyd, S. and Lee, T. 1999. Optimization of inductor circuits via geometric programming. In Proceedings 36th IEEE/ACM Integrated Circuit Design Automation Conference. [39] Huaizhong, L. and Fu, M. 1997. A linear matrix inequality approach to robust h, filtering. IEEE Transactions on Signal Processing. [40] Karmarkar, N. 1984. A new polynomial-time algorithm for linear programming. Combinatorica. [41] Kortanek, K., Xu, X. and Ye, Y. 1997. An infeasible interior-point algorithm for solving primal and dual geometric progams. Mathematical Programming. [42] Lasserre, J. 2002. Bounds on measures satisfying moment conditions. Annals of Applied Probability. [43] Lasserre, J. 2001. Global optimization with polynomials and the problem of moments. SIAM Journal of Optimization. [44] Lenstra, J.K., Rinnooy-Kan, A.H.G. and Schrijver, A. 1991. History of mathematical programming: A collection of personal reminiscenses. Elsevier Science Ltd. [45] Lobo, M., Vandenberghe, L., Boyd, S. and Lebret, H. 1998. Applications of secondorder cone programming. Linear Algebra and its Applications. [46] Lofberg, J. 2004. YALMIP: A toolbox for modeling and optimization in matlab. In Proceedings of the IEEE International Symposium on Computed Aided Control Systems Design. [47] Lovasz, L. 1986. An algorithmic theory of numbers, graphs and convexity. Regional Conference Series in Applied Mathematics. SIAM, Philadelphia. [48] Luque-Calvo, P.L. 2017. Escribir un trabajo fin de estudios con r markdown. Disponible en http://destio.us.es/calvo. [49] Miller, J., Quigley, P. and Zhu, J. 2015. CVXcanon: Automatic canonicalization of disciplined convex programs. https://stanford.edu/class/ee364b/projects/2015projects/ reports/miller_quigley_zhu_report.pdf. [50] MOSEK-ApS 2001. Mosek (software package). Web site: http://www.mosek.com. [51] Nesterov, Y. and Nemirovsky, A. 1988. A general approach to polynomial-time algorithms design for convex programming. Technical report, Centr. Econ. & Math. Inst., USSR Acad. Sci., Moscow, USSR. [52] Nesterov, Y. and Nemirovsky, A. 1993. Interior-point polynomial algorithms in convex programming: Theory and algorithms. Society of Industrial and Applied Mathematics CAPÍTULO 4. CVXR 63