Métodos numéricos en optimización

Traballo Fin de Grao Métodos numéricos en optimización Cheyenne Amoroso Sanmiguel 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Métodos numéricos en optimización Cheyenne Amoroso Sanmiguel Julio 2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Matemática Aplicada Título: Métodos numéricos en optimización Breve descrición do contido Implementación de métodos numéricos de optimización con y sin restricciones utilizando Python. Comparación con métodos de librerías de código libre. Recomendacións Buenas capacidades de programación. iii Índice general Resumen viii Introducción xi 1. Métodos de numéricos 1 1.1. MétododeGauss-Newton ............................ 3 1.1.1. Descripción y propiedades . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2. Implementación del método de Gauss-Newton . . . . . . . . . . . . . 6 1.2. Método de Gauss-Newton amortiguado . . . . . . . . . . . . . . . . . . . . . 6 1.2.1. Implementación del método de Gauss-Newton amortiguado . . . . . 7 1.3. Método de Levenberg-Marquardt . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.1. Descripción y propiedades . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2. Implementación.............................. 10 1.4. Métodos Quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4.1. MétodoDogbox.............................. 14 1.4.2. MétodoNL2SOL............................. 17 2. Resultados numéricos 21 2.1. Descripción de los problemas . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.1.1. Test1:Gauss1 .............................. 22 2.1.2. Test2:Hahn1............................... 23 2.1.3. Test3:Bennett5 ............................. 24 2.2. Resultadosobtenidos............................... 26 2.2.1. Resultados Test 1: Gauss1 . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.2. Resultados Test 2: Hahn1 . . . . . . . . . . . . . . . . . . . . . . . . 29 2.2.3. Resultados Test 3: Bennett5 . . . . . . . . . . . . . . . . . . . . . . . 31 3. Conclusiones 35 v vi ÍNDICE GENERAL A. Conceptos básicos 37 B. Quasi-Newton 39 B.1.MétododeBetts ................................. 40 B.2. Método de Bartholomew-Biggs . . . . . . . . . . . . . . . . . . . . . . . . . 40 B.3.MétododeFletcher-Xu.............................. 41 C. Código Python 45 C.1.Gauss-Newton................................... 45 C.2. Gauss-Newton amortiguado . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Bibliografía 49 xiv INTRODUCCIÓN manera: m´ın θ∈Rnf(θ) = m´ın θ∈Rn 1 2r(θ)Tr(θ) = m´ın θ∈Rn 1 2 m X i=1 [ri(θ)]2, m ≥n, siendo ri(θ) = yi−h(xi,θ) , ri:Rn→R . La interpretación geométrica que podemos dar a este tipo de problemas es encontrar un punto en la supercie z=r(ω) en Rm lo más próximo posible al origen. En el caso de la regresión consideramos la supercie z= (h(x1,θ), ..., h(xm,θ)) ∈Rm . El trabajo se estructura en dos capítulos como sigue: en el Capítulo 1 estudiaremos los aspectos más relevantes de algoritmos utilizados en el problema general de estimación de parámetros. Para ello nos guiaremos principalmente por Sun y Yuan (2006), Bjorck (1996) y Gill et al. (1981). En el Capítulo 2 se mostrará el desempeño de los algoritmos sobre una base de datos tomada de Shen et al. (2013), realizándose una comparación de los distintos métodos explicados y las distintas implementaciones de cada uno de ellos. Los cálculos y la programación se llevarán a cabo a través de Python mediante las bibliotecas libres Virtanen et al. (2020), van der Walt et al. (2011) y Hunter (2007). Se incluirán además tres apéndices, en el primero de ellos encontraremos una serie de conceptos básicos para una mejor comprensión de lo expuesto en el primer capítulo. En el segundo se mostrarán algunos métodos más generales para la resolución del problema expuesto. Finalmente, se incluirá en un tercer apéndice el código implementado por la autora de este trabajo para aquellos métodos que no fuera posible encontrar programado en ninguna biblioteca de Python . Capítulo 1 Métodos de numéricos En este capítulo estudiaremos diferentes métodos numéricos en optimización para problemas de suma de cuadrados de funciones no lineales. Los algoritmos para la optimización no lineal sin restricciones con múltiples variables se pueden clasicar en los siguientes tres grupos: 1. Búsqueda sin el uso de derivadas. Ejemplos de esto son el método simplex, el algoritmo de Hooke y Jeeves, el método de Rosenbrock y el método de las direcciones conjugadas. 2. Búsqueda usando información de la derivada primera. En este caso tenemos el método de máximo descenso, el método del gradiente conjugado y los métodos quasi-Newton. 3. Búsqueda usando información de la derivada segunda como, por ejemplo, el método de Newton. A lo largo de este capítulo asumiremos que las componentes ri(ω) son de clase dos. Denotaremos al jacobiano de r(ω) por J(ω)∈Rm×n donde J(ω)ij =∂ri(ω) ∂ωj (1.1) y a las matrices Hessianas de ri(ω) son Gi(ω)∈Rn×n , donde Gi(ω)jk =∂2ri(ω) ∂ωj∂ωk , i = 1, ..., m. (1.2) Por otra parte, el gradiente de f(ω) se denotará con ∇f(ω) = m X i=1 ri(ω)∇ri(ω) = J(ω)Tr(ω), (1.3) 1 2 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS y su Hessiana con ∇2f(ω) = m X i=1 (∇ri(ω)∇ri(ω)T+ri(ω)∇2ri(ω)) = J(ω)TJ(ω) + S(ω), (1.4) donde S(ω) = m X i=1 ri(ω)∇2ri(ω). (1.5) Existen dos formas de ver el problema (4). La primera de ellas, como ya mencionábamos en la introducción, es interpretarlo como la resolución de un sistema sobredeterminado de m ecuaciones no lineales ri(ω)=0, i = 1, ..., m. En este caso podemos aproximar r(ω) por un modelo lineal centrado en un punto dado, ωk , ˜rk(ω) = r(ωk) + J(ωk)(ω−ωk). (1.6) Cuando m > n tenemos el sistema sobredeterminado rk(ω) = 0 , el cual podemos resolver como un problema de mínimos cuadrados lineal para así aproximar una solución. La segunda forma de ver este problema es tratarlo como un caso especial de optimización en Rn . Para ello empleamos un modelo cuadrático en torno a ωk qk(ω) = f(ωk) + ∇f(ωk)T(ω−ωk) + 1 2(ω−ωk)T∇2f(ωk)(ω−ωk). (1.7) Haciendo uso de las ecuaciones (1.3) y (1.4) llegamos a la siguiente expresión qk(ω) = 1 2r(ωk)Tr(ωk)+(J(ωk)Tr(ωk))T(ω−ωk)+1 2(ω−ωk)T(J(ωk)TJ(ωk)+S(ωk))(ω−ωk). (1.8) Empleando el método de Newton, el mínimo de la función qk se puede calcular de la forma ωk+1 =ωk−(J(ωk)TJ(ωk) + S(ωk))−1J(ωk)Tr(ωk), (1.9) el cual converge localmente de forma cuadrática. El principal problema de este método es que al calcular las mn2 derivadas segundas necesarias para obtener el término S(ωk) el coste computacional puede ser relativamente alto. Por lo tanto, puede resultar útil omitir dicho término. De (1.5) deducimos que cuando ri(ω), i = 1, ..., n, es próximo a cero o es levemente no lineal, en cuyo caso ∇2ri(ω) se aproxima a cero, el término S(ω) es pequeño y puede omitirse. A partir de este razonamien- to se siguen los métodos de Gauss-Newton y Levenberg-Marquardt que trataremos más 1.1. MÉTODO DE GAUSS-NEWTON 3 en detalle a continuación. En este capítulo trataremos también los métodos quasi-Newton que sí utilizan información de la derivada segunda, aunque no directamente y sí mediante aproximaciones recursivas. 1.1. Método de Gauss-Newton 1.1.1. Descripción y propiedades El método Gauss-Newton puede pensarse de dos maneras distintas: la primera de ellas es pensar que se obtiene omitiendo el término de segundo grado, S(ω) , en el modelo cuadrático (1.8). De esta forma llegamos a la siguiente aproximación qG k(ω) = 1 2r(ωk)Tr(ωk) + (J(ωk)Tr(ωk))T(ω−ωk) + 1 2(ω−ωk)TJ(ωk)TJ(ωk)(ω−ωk), (1.10) y aplicando el método de Newton a la condición necesaria de mínimo se obtiene la solución mediante ωk+1 =ωk−(J(ωk)TJ(ωk))−1J(ωk)Tr(ωk). (1.11) Cabe destacar que para que esto tenga sentido es necesario que la matriz J(ωk) sea de rango completo. Como hemos explicado anteriormente, este razonamiento será útil cuando ri(ω), i = 1, ..., n , sea próximo a cero o levemente no lineal. Si nos encontramos en este caso es esperable que el comportamiento del método sea similar al método de Newton dado en (1.9). En particular, en el caso de que nos encontremos con un problema consistente tal que r(˜ω)=0 , donde ˜ ω es solución, la velocidad de convergencia local será idéntica con ambos métodos. Sin embargo, en el caso de residuos grandes la velocidad de convergencia local será inferior en el método de Gauss-Newton. La segunda forma de pensar este método es como una serie de aproximaciones de r(ω) de la forma (1.6). Si denotamos por sk a la corrección de la aproximación ωk , dicha corrección se obtendrá al resolver el problema de mínimos cuadrados m´ın s∈Rn 1 2||r(ωk) + J(ωk)s||2 2. (1.12) De esta forma, la nueva aproximación será ωk+1 =ωk+sk . En el caso de que la matriz J(ωk) tenga rango completo la solución al problema lineal de mínimos cuadrados es única y vemos que coincide con la obtenida mediante el razonamiento 4 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS anterior sk=−(J(ωk)TJ(ωk))−1J(ωk)Tr(ωk). (1.13) Tenemos que ∇f(ωk) es distinto de cero, por lo que deducimos que sk se trata de una dirección de descenso ya que sT k∇f(ωk) = sT kJ(ωk)Tr(ωk) = −sT kJ(ωk)TJ(ωk)sk≤0, (1.14) donde la desigualdad será estricta si no se cumple que J(ωk)sk= 0 . Esta última igualdad es equivalente a J(ωk)Tr(ωk) = ∇f(ωk) = 0 . En el caso de que J(ωk) esté mal condicionado o sea singular también podremos calcular sk usando descomposiciones de J(ωk) como, por ejemplo, la factorización QR o la descomposición en valores singulares o SVD ( Singular Value Decomposition ). Si J(ωk) no es de rango completo podemos tomar sk como sk=−J(ωk)Ir(ωk), (1.15) donde J(ωk)I denota la matriz pseudo-inversa de J(ωk) . A continuación analizaremos la convergencia del método. Para ello enunciaremos tres teoremas que nos permitirán concluir lo siguiente: 1. Si S(˜ ω)=0 , el método tiene convergencia cuadrática. 2. Si S(˜ ω) es pequeño comparado con J(˜ ω)TJ(˜ ω) el método es localmente Q-linealmente convergente. 3. Si S(˜ ω) es demasiado grande el método no convergerá. Teorema 1.1.1. Sea f:Rn−→ R, f ∈ C2 . Supongamos que ˜ ω es el mínimo local del problema (4) y que J(˜ ω)TJ(˜ ω) es denida positiva. Supongamos además que la sucesión {ωk} generada por el Algoritmo 1.1.5 converge a ˜ ω . Entonces, si ∇2f(ω) y (J(ω)TJ(ω))−1 son Lipschitzianas en un entorno de ˜ ω , se cumple ||ωk+1 −˜ ω||2≤ ||(J(˜ ω)TJ(˜ ω))−1||2||S(˜ ω)||2||ωk−˜ ω||2+O(||ωk−˜ ω||2 2) (1.16) Demostración. Ver Sun y Yuan (2006) (p. 357-358).  1.1. MÉTODO DE GAUSS-NEWTON 5 Teorema 1.1.2. Sea f: Ω ⊂Rn−→ R, f ∈ C2(Ω) donde Ω es un conjunto convexo y abierto. Sea J(ω) de Lipschitz continua en Ω , es decir, ||J(u)−J(v)||2≤γ||u−v||2,∀u,v∈Ω, (1.17) y ||J(ω)||2≤α, ∀ω∈Ω . Supongamos que existen ˜ ω∈Ω y λ, σ ≥0 tal que, si J(˜ ω)Tr(˜ ω) = 0 , λ es el menor valor propio de J(˜ ω)TJ(˜ ω) y ||(J(ω)−J(˜ ω))Tr(˜ ω)||2≤σ||ω−˜ ω||2,∀ω∈Ω. (1.18) Si σ < λ , entonces, para cualquier c∈(1,λ σ) , existe ε > 0 tal que ∀ω0∈N(˜ ω, ε) la secuencia generada por el método de Gauss-Newton mediante el Algoritmo 1.1.5 está bien denida, converge a ˜ ω y satisface ||ωk+1 −˜ ω||2≤cσ λ||ωk−˜ ω||2+cαγ 2λ||ωk−˜ ω||2 2 ||ωk+1 −˜ ω||2≤cσ +λ 2λ||ωk−˜ ω||2<||ωk−˜ ω||2. (1.19) Demostración. Ver Sun y Yuan (2006) (p. 358-359).  Teorema 1.1.3. Consideremos que las suposiciones del Teorema 1.1.1 o las del Teorema 1.19 se verican. Si r(˜ ω) = 0 , entonces existe ε > 0 tal que para cualquier ω0∈N(˜ ω, ε) la sucesión {ωk} generada por el método de Gauss-Newton converge de forma cuadrática a ˜ ω . Demostración. Ver Sun y Yuan (2006) (p. 360).  El siguiente ejemplo (Fletcher, 1980, pág. 94) muestra que el método de Gauss-Newton funciona bien para problemas con residuos pequeños. Sin embargo, cuando sean grandes o los problemas sean claramente no lineales el método no será localmente convergente. Ejemplo 1.1.4. Consideramos el problema con m =2 y n =1 dado por r1(ω) = ω+ 1, r2(ω) = λω2+ω+ 1 donde λ es un parámetro cualquiera. Tenemos el siguiente problema escrito de forma estándar m´ın ω∈Rnf(ω) = m´ın ω∈Rn 1 2 2 X i=1 ri(ω)2= m´ın ω∈Rn 1 2[(ω+ 1)2+ (λω2+ω+ 1)2], 6 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS de manera que ˜ ω= 0 . Llegamos a que, en este caso particular, ωk+1 =2λ2ω3 k+λω2 k+ 2λωk 1 + (2λωk+ 1)2. Observamos que si λ= 0 entonces ω1=0=˜ ω , lo que muestra que, cuando el problema es lineal, el método de Gauss-Newton lo resuelve en una única iteración. En el caso de que λ6= 0 ωk+1 =λωk+O(||ωk||2 2). Así, cuando λ es sucientemente pequeño la convergencia del método será lineal, mientras que si |λ|>1 el método no será convergente. 1.1.2. Implementación del método de Gauss-Newton El método presentado en la sección anterior no está, bajo mi conocimiento, programado en ninguna biblioteca libre de Python, por lo que ha sido implementado por la autora de este trabajo. El código utilizado podemos encontrarlo en el apéndice C.1. El algoritmo programado, basado en la descripción dada del método de Gauss-Newton, es el siguiente: Algoritmo 1.1.5. (Gauss-Newton) Paso 0. Dados ω0 , ε > 0 , k:= 0 . Paso 1. Si ||∇f(ωk)||2≤ε , stop. Paso 2. Cálculo de sk resolviendo J(ωk)TJ(ωk)sk=−J(ωk)Tr(ωk). Paso 3. ωk+1 =ωk+sk , k:= k+ 1 . Volver al Paso 1.  Como ya hemos comentado, cabe la posibilidad de que J(ωk) no sea una matriz de rango completo. Como solución se propone calcular sk como sk=−J(ωk)Ir(ωk) donde J(ωk)I denota la matriz pseudo-inversa de J(ωk) . 1.2. Método de Gauss-Newton amortiguado Según los problemas observados en el ejemplo mostrado en la sección anterior, en la práctica se utiliza un método de Gauss-newton con búsqueda de línea o también llamado método de Gauss-Newton amortiguado ( damped Gauss-Newton ). Sea αk un escalar y sk la solución de 1.12, que denominaremos dirección de Gauss-Newton, el método es como sigue ωk+1 =ωk+αksk. (1.20) 1.2. MÉTODO DE GAUSS-NEWTON AMORTIGUADO 7 Se verica que sk se mantiene invariante ante transformaciones lineales sobre la variable independiente ω y, como ya vimos anteriormente, es una dirección de descenso. Para que esta modicación del método sea un algoritmo viable tenemos que tener cuidado con la elección de αk . Generalmente existen dos maneras que facilitan la elección del escalar: 1. Tomar αk como el mayor número de la secuencia 1,1 2,1 4, ... tal que se cumpla la desigualdad ||r(ωk)||2 2− ||r(ωk+αksk)||2 2≥1 2αk||J(ωk)sk||2 2 (1.21) 2. Tomar αk como la solución de m´ın α||r(ωk+αk)||2 2 (1.22) Observación 1.2.1 . El primer método explicado para la búsqueda de αk se trata básicamente de la regla de búsqueda de paso Armijo-Goldstein 1 , mientras que el segundo es el cálculo del paso óptimo. Tomando de forma adecuada αk tenemos que este método siempre toma pasos de descenso, por lo tanto es localmente convergente. De hecho, se trata de un método globalmente convergente en la mayoría de los casos pero continúa siendo muy lento cuando el residuo es grande o los problemas son claramente no lineales. 1.2.1. Implementación del método de Gauss-Newton amortiguado Al igual que ocurría con el método de Gauss-Newton este método no está programado en ninguna biblioteca libre de Python, por lo que ha sido implementado por la autora de este trabajo. En el apéndice C.2 podemos encontrar el código utilizado. El algoritmo programado, basado en la descripción dada, es el siguiente: Algoritmo 1.2.2. (Gauss-Newton amortiguado) Paso 0. Dados ω0 , ε > 0 , k:= 0 . Paso 1. Si ||∇f(ωk)||2≤ε , stop. Paso 2. Cálculo de sk resolviendo J(ωk)TJ(ωk)sk=−J(ωk)Tr(ωk). Paso 3. Cálculo de αk . Dos opciones: 1 La regla de Armijo-Goldstein sigue el esquema siguiente: (i) Se elige t > 0, β ∈(0,1), ε ∈(0,1) . (ii) Se toma p= 0 y se comprueba si se verica q(βpt)≤q(0) + βptεq0(0) . (iii) Si se cumple se toma α=βpt , en caso contrario se hace p=p+ 1 y se vuelve a (i). 8 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS 1. Tomar αk= m´ax 1,1 2,1 4, ... tal que se verique ||r(ωk)||2 2− ||r(ωk+αksk)||2 2≥1 2αk||J(ωk)sk||2 2 2. Tomar αk como la solución de m´ın α||r(ωk+αk)||2 2 Paso 4. ωk+1 =ωk+αksk , k:= k+ 1 . Volver al Paso 1.  1.3. Método de Levenberg-Marquardt En la sección 1.2 veíamos que el método Gauss-Newton amortiguado era globalmente convergente en la mayoría de los casos. Sin embargo, cuando J(ω) no sea de rango completo este método tendrá dicultades para trabajar o el algoritmo convergerá a un punto no estacionario. Para solucionar este problema podemos tener en cuenta las segundas derivadas (Sección 1.4), o bien, estabilizar el método Gauss-Newton amortiguado. Para llevar a cabo la segunda estrategia mencionada consideramos la técnica de región de conanza ( trust-region technique ), con la que obtenemos el método de Levenberg- Marquardt. Este algoritmo fue publicado por primera vez por Kenneth Levenberg, Levenberg (1944), y redescubierto en 1963 por Donald Marquardt, Marquardt (1963). 1.3.1. Descripción y propiedades El método de Gauss-Newton podíamos pensarlo de dos maneras distintas: una de ellas era pensarlo como una serie de aproximaciones de r(ω) de la forma (1.6) obteniendo un problema lineal de mínimos cuadrados (1.12). Desafortunadamente, esta linealización no es efectiva para todos los (ω−ωk) . Para solucionar esto empleamos la técnica de región de conanza y añadimos una restricción. Consideramos entonces el siguiente modelo de región de conanza: m´ın ω∈Rn 1 2||r(ωk) + J(ωk)(ω−ωk)||2 2 s.t ||ω−ωk||2≤∆k, (1.23) 1.3. MÉTODO DE LEVENBERG-MARQUARDT 9 que podemos verlo como un problema no lineal de mínimos cuadrados con restricciones. Este modelo puede reescribirse como m´ın ω∈Rn 1 2r(ωk)Tr(ωk) + r(ωk)TJ(ωk)(ω−ωk) + 1 2(ω−ωk)TJ(ωk)TJ(ωk)(ω−ωk) s.t ||ω−ωk||2<∆k. Denotamos s=ω−ωk . La solución al problema (1.23) se obtendrá resolviendo el sistema (J(ωk)TJ(ωk) + µkI)s=−J(ωk)Tr(ωk), (1.24) de donde obtenemos que ωk+1 =ωk−(J(ωk)TJ(ωk) + µkI)−1J(ωk)Tr(ωk). (1.25) El parámetro ∆k≥0 se encarga de controlar las iteraciones y limita el tamaño de sk . Nótese que el parámetro sk está bien denido incluso cuando J(ωk) no es de rango completo. Si se verica que ||(J(ωk)TJ(ωk))−1J(ωk)Tr(ωk)||2≤∆k entonces se cumple que µk= 0 . En otro caso, existe µk>0 tal que la solución sk satisface ||sk||2= ∆k y (J(ωk)TJ(ωk) + µkI)sk=−J(ωk)Tr(ωk). (1.26) Otra forma de ver este método es como una mezcla entre Gauss-Newton y el método de descenso de gradiente ( steepest descent method ), de tal manera que este método permite elegir dos direcciones distintas. En el caso de que µk= 0 la dirección sería la misma que en el Gauss-Newton, mientras que si µk es demasiado grande el sistema (1.24) se reduciría a µkIs=−J(ωk)Tr(ωk). (1.27) y la dirección que obtendríamos estaría muy próxima a la de máximo descenso. A continuación enunciaremos una serie de propiedades del método Levenberg-Marquardt donde denotaremos s=s(µ) como la solución del sistema (J(ω)TJ(ω) + µI)s=−J(ω)Tr(ω). Teorema 1.3.1. Cuando µ aumenta de forma monótona desde cero, ||s(µ)||2 decrecerá de forma estrictamente monótona. Demostración. Ver Sun y Yuan (2006, pág. 363-364).  16 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS donde b=||P C||2 ||C||2∈[0,1] . Trivialmente qk(s(a)) disminuye de forma monótona alcanzando así un valor más bajo para el modelo. Es probable que el algoritmo exhiba una convergencia lenta cuando el rango de J(ω) es menor que el número de variables. Implementación De la descripción dada del método de Dogbox en la subsección anterior se sigue el siguiente algoritmo: Algoritmo 1.4.1. (Dogbox) Paso 0. Elegimos el iterante inicial ω0 , el parámetro inicial de la región de conanza ∆0 y jamos k:= 0 . Paso 1. Construimos el modelo cuadrático en torno a ωk , es decir, m´ın ω∈Rn(J(ωk)Tr(ωk))T(ω−ωk) + 1 2(ω−ωk)TMk(ω−ωk), s.t ||ω−ωk||∞≤∆k, Paso 2. Calculamos sk : if N=−Mk(J(ωk)Tr(ωk)) es factible sk=N else if C=r(ωk)TJ(ωk)J(ωk)Tr(ωk) (r(ωk)TJ(ωk)BkJ(ωk)Tr(ωk)J(ωk)Tr(ωk) es factible encontrar el máximo α tal que C+α(N−C)∈Tk , sk=C+α(N−C) else encontrar el máximo β tal que PC ≡βC ∈Tk y el máximo α tal que PC +α(N−PC)∈Tk , sk=PC +α(N−PC) Paso 3. Calculamos el radio δk=f(ωk)−f(ωk+sk) qk(0)−qk(sk) . Paso 4. Elegimos el nuevo punto ωk+1 : if δk≤0.1 ωk+1 =ωk else ωk+1 =ωk+sk Paso 5. Actualizamos ∆k : if δk<0.25 ∆k+1 =||sk||2/4 1.4. MÉTODOS QUASI-NEWTON 17 else if δk>0.75 y ||sk||2= ∆k ∆k+1 = 2∆k else ∆k+1 = ∆k Paso 6. k:= k+ 1 . Volver al Paso 1.  Al igual que ocurría con el método de Levenberg-Marquardt, en la biblioteca libre Scipy encontramos distintos módulos que implementan este algoritmo. En esta memoria utilizaremos dos: scipy.optimize.least_squares y scipy.optimize.curve_t . Analizaremos su comportamiento y lo compararemos tanto con las implementaciones de los restantes métodos descritos en este trabajo como entre ellos. 1.4.2. Método NL2SOL La aplicación directa de los métodos quasi-Newton al problema de mínimos cuadrados no lineales descrito anteriormente, en ocasiones en la práctica, no resultan tan ecientes. Una de las razones por las que ocurre esto es que estos métodos ignoran la información que aporta J(ωk) , pero a menudo J(ωk)TJ(ωk) es la parte dominante de ∇2f(ω) . Una estrategia posible es, en lugar de aproximar ∇2f(ω) , incluir una aproximación quasi-Newton del término desconocido de la segunda derivada, S(ω) . Sea Bk la aproximación de la matriz S(ωk) , se tiene entonces (Bk+J(ωk)TJ(ωk))sk=−J(ωk)Trk. (1.44) Dado que S(ωk+1) = m X i=1 ri(ωk+1)∇2ri(ωk+1) podemos tomar Bk+1 = m X i=1 ri(ωk+1)(Hi)k+1 como la aproximación de S(ωk+1) , donde (Hi)k+1 es una aproximación de ∇2ri(ωk+1) . Tenemos entonces Bk+1(ωk+1 −ωk) = m X i=1 (ri(ωk+1)(Hi)k+1)(ωk+1 −ωk) = m X i=1 ri(ωk+1)(∇ri(ωk+1)− ∇ri(ωk)) =J(ωk+1)Tr(ωk+1)−J(ωk)Tr(ωk+1) =: yk, (1.45) 18 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS lo cual es una condición quasi-Newton impuesta sobre Bk . En Dennis et al. (1981) se sugiere una regla de aproximación para S(ωk) simple, general y geométrica. La idea es elegir un conjunto de características deseables para la aproximación y posteriormente seleccionar S(ωk+1) como el punto factible más cercano a S(ωk) . Veamos entonces qué propiedades debe vericar S(ωk+1) . Es razonable comenzar con S(ω0)=0 dado que es barato computacionalmente y, además, es lógico pensar que q0=qG 0 . Supongamos entonces que S(ωk) está disponible. Recordamos que queremos aproximar Pm i=1 ri(ω)∇2ri(ω) , por lo tanto es obvio que debe ser simétrico. Dennis et al. (1981) se basa en el siguiente resultado para proporcionar un formula de actualización así como sus propiedades. Notación 1.4.2 . Sea A∈ Mm×n(C) , se dene la norma de Frobenius, llamada también norma de HilbertSchmidt, como ||A||F:= (tr(ATA))1 2=Pm i=1 Pn j=1 |aij|21 2 . Teorema 1.4.3. Sea vT ksk>0 y T∈Rn×n una matriz simétrica y denida positiva tal que TTTsk=vk donde vk=∇f(ωk+1)− ∇f(ωk). Entonces la actualización Bk+1 =Bk+(yk−Bksk)vT k+vk(yk−Bksk)T sT kvk −(yk−Bksk)Tsk (sT kvk)2vkvT k (1.46) es la única solución del problema de minimización m´ın ||T−T(Bk+1 −Bk)T−1||F s.t (Bk+1 −Bk) simétrica , Bk+1sk=yk. Dennis et al. (1981) emplearon la condición quasi-Newton (1.45) y la fórmula de actualización (1.46), y presentaron un nuevo algoritmo quasi-Newton con región de conanza que recibe el nombre de NL2SOL. En cada etapa es necesario resolver el subproblema m´ın ω∈Rn 1 2r(ωk)Tr(ωk)+(J(ωk)Tr(ωk))T(ω−ωk) +1 2(ω−ωk)T(J(ωk)TJ(ωk) + Bk)(ω−ωk) s.t ||ω−ωk||2≤∆k. (1.47) Se cumple que, en el caso en el que el residuo sea cero, la matriz Bk debería desaparecer o, al menos, acercarse a 0 . Sin embargo, en el algoritmo propuesto, la actualización (1.46) no 1.4. MÉTODOS QUASI-NEWTON 19 garantiza que desaparezca la matriz Bk . Cuando esto ocurre, el método de Gauss-Newton es superior al NL2SOL. Además, esto puede interferir con la convergencia superlineal del método. Para solucionar dicho problema utilizaremos una modicación directa del autoescalado. La idea es actualizar γkBk en lugar de Bk para obtener Bk+1 . Podemos tomar como factor de escalado γk=|(ωk+1 −ωk)Tyk| |(ωk+1 −ωk)TBk(ωk+1 −ωk)|. (1.48) Dado que queremos que Bk no sea demasiado grande tomamos γk= m´ın |(ωk+1 −ωk)Tyk| |(ωk+1 −ωk)TBk(ωk+1 −ωk)|,1 (1.49) A partir de los experimentos numéricos reejados en Dennis et al. (1981) observamos que para problemas con residuos grandes, el algoritmo quasi-Newton NL2SOL es signicativamente ventajoso. Sin embargo, para problemas con residuos signicativamente pequeños, el rendimiento de NL2SOL y del algoritmo de Levenberg-Marquardt implementado en Moré (1978) es similar. En el caso de los problemas con residuo nulo seguimos prerendo el método de Gauss-Newton. 20 CAPÍTULO 1. MÉTODOS DE NUMÉRICOS Capítulo 2 Resultados numéricos A lo largo de este capítulo analizaremos el comportamiento de los algoritmos descritos empleando el lenguaje de programación interpretado Python 1 . Como ya hemos comentado, haremos uso tanto de bibliotecas libres como SciPy, donde encontramos implementados los métodos de Levenberg-Marquardt y Dogbox, como de código propio en el caso de los algoritmos restantes. La validación de los algoritmos descritos la realizaremos mediante problemas de distintos grados de dicultad tomados de Statistical Reference Datsets Project (STRD) desarrollado por el Sta Statistical Engineering Division and the Mathematical and Computational Sciences Division del National Institute of Standards and Techonology (NIST), que proporciona bases de datos referenciadas con valores certicados. En primer lugar, haremos una descripción de los distintos problemas para, posteriormente, mostrar y analizar los resultados obtenidos. Encontrar el mejor código es una tarea casi imposible y muy dependiente de la de- nición de mejor. Para cualquiera que sea el criterio utilizado, el procedimiento de prueba debe intentar medir la capacidad del código para encontrar soluciones. Los problemas de mínimos cuadrados no lineales son intrínsecamente difíciles y, generalmente, es posible encontrar un conjunto de datos que haga fallar incluso a los códigos más robustos. Por lo 1 Se trata de un lenguaje de programación multiparadigma, ya que soporta programación orientada a objetos, programación imperativa y, en menor medida, programación funcional. Es un lenguaje interpretado, dinámico y multiplataforma. Se administra por la Python Software Foundation y posee una licencia de código abierto, denominada Python Software Foundation License, que es compatible con la Licencia pública general de GNU a partir de la versión 2.1.1, e incompatible en ciertas versiones anteriores. En este trabajo se ha usado la versión 3.7.1. 21 22 CAPÍTULO 2. RESULTADOS NUMÉRICOS tanto, la mayoría de las evaluaciones de software de mínimos cuadrados no lineales también deben incluir una medida de la conabilidad del código, es decir, el código debe reconocer correctamente cuándo ha encontrado una solución. Los conjuntos de datos empleados son particularmente adecuados para tales pruebas de robustez y conabilidad. Se incluyen problemas de mínimos cuadrados no lineales generados y del mundo real. Los valores certicados son soluciones "mejores disponibles", obtenidas con precisión de 128 bits y conrmadas por, al menos, dos algoritmos y paquetes de software diferentes que utilizan derivadas analíticas. 2.1. Descripción de los problemas En esta sección presentaremos 3 problemas test ordenados por grado de dicultad, según la clasicación del NIST. Emplearemos dichos problemas para analizar tanto el buen comportamiento de los algoritmos descritos en el primer capítulo como su implementación. 2.1.1. Test 1: Gauss1 Los datos empleados en este problema son datos generados que se corresponden con dos normales en una base de decrecimiento exponencial con errores normales de media cero y varianza 6.25. Sea {(xi, yi)}250 i=1 una muestra donde xi, yi∈R son los valores de las variables explicativa y respuesta, respectivamente, correspondientes al i -ésimo individuo. Sea θ∈R8 el vector de parámetros desconocidos. Consideraremos el ajuste no lineal dado por: yi=h(xi,θ) + εi, εi∈N(0,6.25), i = 1, ..., 250, donde h(x, θ) = θ1·exp(−θ2x) + θ3·exp −(x−θ4)2 θ2 5+θ6·exp −(x−θ7)2 θ2 8. El problema a resolver sigue la fórmula general (4): m´ın θ∈R8f(θ) = m´ın θ∈R8 1 2r(θ)Tr(θ) = m´ın θ∈R8 1 2 250 X i=1 [ri(θ)]2, siendo ri(θ) = yi−h(xi,θ) , ri:R8→R , i= 1, ..., 250. En la Figura 2.1 representamos el ajuste benchmark tomando como valores certi- cados θ= (9.8778210871e+01, 1.0497276517e-02, 1.0048990633e+02, 6.7481111276e+01, 2.1. DESCRIPCIÓN DE LOS PROBLEMAS 23 Figura 2.1: Ajuste de Gauss1 con la función del benchmark. 2.3129773360e+01, 7.1994503004e+01, 1.7899805021e+02, 1.8389389025e+01). Empleando dichos valores certicados obtenemos que RSS(θ) = 1.3158222432e+03. Al realizar la implementación de este problema utilizaremos dos iterantes iniciales distintos: θ1 0= (97.0,0.009,100.0,65.0,20.0,70.0,178.0,16.5), θ2 0= (94.0,0.0105,99.0,63.0,25.0,71.0,180.0,20.0). 2.1.2. Test 2: Hahn1 Los datos empleados en este problema son el resultado de un estudio del NIST que involucra la dilatación térmica del cobre. La variable respuesta es el coeciente de dilatación térmica, y la variable independiente es la temperatura medida en grados Kelvin. Sea {(xi, yi)}236 i=1 una muestra donde xi, yi∈R son los valores de las variables explicativa y respuesta, respectivamente, correspondientes al i -ésimo individuo. Sea θ∈R7 el vector de parámetros desconocidos. Consideraremos el ajuste no lineal dado por: yi=h(xi,θ) + εi, i = 1, ..., 236, donde h(x, θ) = θ1+θ2x+θ3x2+θ4x3 1 + θ5x+θ6x2+θ7x3. 24 CAPÍTULO 2. RESULTADOS NUMÉRICOS El problema a resolver sigue la fórmula general (4): m´ın θ∈R7f(θ) = m´ın θ∈R7 1 2r(θ)Tr(θ) = m´ın θ∈R7 1 2 236 X i=1 [ri(θ)]2, siendo ri(θ) = yi−h(xi,θ) , ri:R7→R , i= 1, ..., 236. En la Figura 2.2 representamos el ajuste benchmark tomando como valores certi- cados θ= (1.0776351733e+00, -1.2269296921e-01, 4.0863750610e-03,-1.4262662514e-06, - 5.7609940901e-03, 2.4053735503e-04, -1.2314450199e-07). Empleando dichos valores certi- cados obtenemos que RSS(θ) = 1.5324382854e+00. Figura 2.2: Ajuste de Hahn1 con la función del benchmark. Al realizar la implementación de este problema utilizaremos dos iterantes iniciales distintos: θ1 0= (10,−1,0.05,−0.00001,−0.05,0.001,−0.000001), θ2 0= (1,−0.1,0.005,−0.000001,−0.005,0.0001,−0.0000001). 2.1.3. Test 3: Bennett5 Los datos empleados en este problema son el resultado de un estudio del NIST que involucra el modelado de magnetización de superconductividad, donde la variable respuesta es el magnetismo, y la variable independiente es el tiempo tomando como unidad de medida 2.1. DESCRIPCIÓN DE LOS PROBLEMAS 25 los minutos. Sea {(xi, yi)}154 i=1 una muestra donde xi, yi∈R son los valores de las variables explicativa y respuesta, respectivamente, correspondientes al i -ésimo individuo. Sea θ∈R3 el vector de parámetros desconocidos. Consideraremos el ajuste no lineal dado por: yi=h(xi,θ) + εi, i = 1, ..., 154, donde h(x, θ) = θ1(θ2+x) −1 θ3. El problema a resolver sigue la fórmula general (4): m´ın θ∈R3f(θ) = m´ın θ∈R3 1 2r(θ)Tr(θ) = m´ın θ∈R3 1 2 154 X i=1 [ri(θ)]2, siendo ri(θ) = yi−h(xi,θ) , ri:R3→R , i= 1, ..., 154. En la Figura 2.3 representamos el ajuste benchmark tomando como valores certicados θ= (−2.5235058043e+ 03,4.6736564644e+ 01,9.3218483193e−01) . Empleando dichos valores certicados obtenemos que RSS(θ)=5.2404744073e−04 . Figura 2.3: Ajuste de Bennett5 con la función del benchmark. Al realizar la implementación de este problema utilizaremos dos iterantes iniciales distintos: θ1 0= (−2000,50,0.8), θ2 0= (−1500,45,0.85). 32 CAPÍTULO 2. RESULTADOS NUMÉRICOS mos que sea el más eciente debido a que no se alcanza la solución propuesta. Por lo tanto, a la vista de los resultados expuestos, obtenemos que entre los métodos implementados que sí alcanzan los valores certicados el que tiene un menor coste computacional es el método de Levenberg-Marquardt al emplear el módulo curve_t . Destaca el método de Gauss-Newton amortiguado alacanzando la solución en un tiempo muy competitivo. En la gura 2.8 se representan los ajustes obtenidos. Tiempo Error RSS Gauss-Newton 1.585631e+00 s 2.2892082490e-02 5.2404744072e-04 Gauss-Newton amortiguado 4.452369e-01 s 2.2892082490e-02 5.2404744072e-04 Levenberg-Marquardt (curve_t) 8.676696e-02 s 2.2892082490e-02 5.2404744073e-04 Levenberg-Marquardt (least_squares) 2.393365e-02 s 2.3478732372e-02 5.512508738e-04 Levenberg-Marquardt (leastsq) 1.595163e-02 s 2.4900365520e-02 6.200282030e-04 Dogbox (curve_t) 1.087105e-01 s 2.2892082490e-02 5.2404744073e-04 Dogbox (least_squares) 9.129357e-02 s 2.2892082490e-02 5.2404744073e-04 Cuadro 2.5: Resultados obtenidos para Bennett5 con el iterante inicial θ1 0 . Figura 2.8: Ajustes de Bennett5 partiendo del iterante inicial θ1 0 . 2.2. RESULTADOS OBTENIDOS 33 En el Cuadro 2.6 se presentan los resultados obtenidos empleando como iterante inicial θ2 0 . No se aprecian diferencias de comportamiento entre los métodos utilizados alcanzándose en todos los casos los valores certicados. De nuevo, obtenemos que el método más eciente es Levenberg-Marquardt mediante el módulo curve_t . En la gura 2.9 se representan los ajustes obtenidos. Tiempo Error RSS Gauss-Newton 1.681291e+00 s 2.2892082490e-02 5.2404744073e-04 Gauss-Newton amortiguado 2.098632e-01 s 2.2892082490e-02 5.2404744073e-04 Levenberg-Marquardt (curve_t) 1.097035e-02 s 2.2892082490e-02 5.2404744073e-04 Levenberg-Marquardt (least_squares) 2.493477e-02 s 2.2894309265e-02 5.2414939672e-04 Levenberg-Marquardt (leastsq) 1.595736e-02 s 2.2892082491e-02 5.2404744080e-04 Dogbox (curve_t) 1.994944e-02 s 2.2892082490e-02 5.2404744073e-04 Dogbox (least_squares) 1.499319e-02 s 2.2892082490e-02 5.2404744073e-04 Cuadro 2.6: Resultados obtenidos para Bennett5 con el iterante inicial θ2 0 . Figura 2.9: Ajustes de Bennett5 partiendo del iterante inicial θ2 0 . Según los resultados expuestos en el Cuadro 2.5 y en el Cuadro 2.6 apreciamos que em- 34 CAPÍTULO 2. RESULTADOS NUMÉRICOS pleando el método de Levenberg-Marquardt mediante el módulo curve_t el coste computacional se reduce en una octava parte al utilizar como iterante inicial a θ2 0 respecto a utilizar θ1 0 . Destaca que en ambos casos el método de Gauss-Newton amortiguado mejora al método de Gauss-Newton alcanzando la solución en un tiempo muy competitivo respec- to al resto. Debemos tener siempre en consideración que ambos métodos, Gauss-Newton y Gauss-Newton amortiguado, son de implementación propia. Capítulo 3 Conclusiones El objetivo de este Trabajo de Fin de Grado era analizar el funcionamiento de algunos de los métodos existentes para la resolución del problema de ajuste no lineal y comparar sus distintas implementaciones presentes en Python . Las conclusiones generales tras haber realizado el presente Trabajo de Fin de Grado son las siguientes: Se han complementado los conocimientos adquiridos en el Grado de Matemáticas impartidos en las asignaturas de Métodos Numéricos: 'Cálculo Numérico nunha Variable', 'Análise Numérica Matricial' y 'Métodos Numéricos en Optimización e Ecuacións Diferenciais'. Se han aprendido los fundamentos de una serie de métodos numéricos que permiten la resolución de problemas de optimización de mínimos cuadrados no lineales con el n poder aplicar dichos conocimientos a realizar un análisis de los propios métodos, su implementación en ordenador y aplicarlos a problemas concretos. Se han adquirido nuevos conocimientos en el ámbito de la programación, aprendiendo un nuevo lenguaje que no se incluye en el programa formativo del Grado en Matemáticas como es Python . Se han extraído también una serie de conclusiones particulares tras la realización del trabajo: El método de Levenberg-Marquardt aplicado en los problemas benchmark descritos en el presente Trabajo de Fin de Grado ha demostrado ser signicativamente más eciente en el cálculo de parámetros frente al resto de métodos expuestos. 35 36 CAPÍTULO 3. CONCLUSIONES Aunque en los resultados expuestos no sea posible apreciarlo con claridad, se deduce que en los casos en los que el residuo sea pequeño o el problema no sea claramente no lineal los métodos de Gauss-Newton y Gauss-Newton amortiguado resultan ser muy competitivos. Esto no es posible apreciarlo debido a que ambos han sido implementados por la autora de este trabajo. A juicio de la autora del presente Trabajo de Fin de Grado habría resultado interesante realizar la implementación de los métodos propuestos en el Apéndice B; en particular, el método propuesto por Fletcher-Xu. Apéndice A Conceptos básicos En este apéndice expondremos una serie de conceptos básicos a tener en cuenta para una mejor comprensión del trabajo. Para ello nos guiaremos por Viaño y Burguera (2013) donde podrá encontrarse más información. Sea Ω⊆Rn un conjunto abierto y f: Ω ⊆Rn−→ R . Diremos que un punto u∈Ω es mínimo local o relativo de f si existe ε > 0 tal que B(u, ε)⊂Ω y, además, f(u)≤f(ω),∀ω∈B(u, ε) . Si la desigualdad es estricta diremos que es un mínimo local estricto. Diremos que se trata de un mínimo global o absoluto si f(u)≤f(ω),∀ω∈Ω . De nuevo, si la desigualdad es estricta diremos que se trata de un mínimo global estricto. Por denición, el mínimo global estricto, si existe, es único. Sea ω∈Ω y supongamos que f es derivable en ω . Denotamos por ∇f(ω) al vector gradiente de f en ω . Tenemos la siguiente condición necesaria de mínimo relativo. Teorema A.0.1. Sea Ω⊆Rn un conjunto abierto, no vacío y sea f: Ω ⊆Rn−→ R un funcional con un mínimo relativo en u∈Ω . Si f es derivable en u entonces se verica que ∇f(u) = θ . Supongamos ahora que f es dos veces derivable en ω∈Ω , denotamos por ∇2f(ω) la matriz Hessiana de f en ω . Teorema A.0.2. Sea Ω⊆Rn un conjunto abierto, no vacío y sea f: Ω ⊆Rn−→ R un funcional con un mínimo relativo en u∈Ω . Se supone f derivable en un entorno de u y dos veces derivable en u . Entonces se verica: 1. ∇f(u) = θ . 2. La matriz hessiana de f es semidenida positiva. 37 38 APÉNDICE A. CONCEPTOS BÁSICOS Veamos ahora una serie de condiciones sucientes que nos permiten armar la existencia de un mínimo relativo Teorema A.0.3. Sea Ω⊆Rn un conjunto abierto y no vacío, u∈Ω . Sea f: Ω ⊆Rn−→ R un funcional derivable en un entorno de u y tal que ∇f(u) = θ . Entonces: 1. Si la función f es dos veces derivable en u y su matriz hessiana es denida positiva en u entonces la función f admite un mínimo relativo estricto en u . 2. Si la función f es dos veces derivable en una bola B(u, r)⊆Ω de forma que la matriz hessiana es semidenida positiva en todos los puntos de la bola entonce f admite un mínimo relativo en u . Apéndice B Quasi-Newton En este apéndice se mostrarán algunos métodos numéricos más generales para la resolución de los problemas de minimización de suma de cuadrados de funciones no lineales. Estos métodos podemos incluirlos dentro del marco de los quasi-Newton. Recordemos que la estrategia a seguir era introducir una aproximación quasi-Newton del término S(ωk) que denotábamos por Bk de manera que (Bk+J(ωk)TJ(ωk))sk=−J(ωk)Tr(ωk). Tomando Bk+1 = m X i=1 ri(ωk+1)(Hi)k+1, donde (Hi)k+1 era una aproximación de ∇2ri(ωk+1) , llegábamos a que Bk+1 tenía que satisfacer Bk+1(ωk+1 −ωk) = m X i=1 (ri(ωk+1)(Hi)k+1)(ωk+1 −ωk) = m X i=1 ri(ωk+1)(∇ri(ωk+1)− ∇ri(ωk)) =J(ωk+1)Tr(ωk+1)−J(ωk)Tr(ωk+1) =: yk. De forma análoga, si pedimos que se verique (Bk+1 +J(ωk+1)TJ(ωk+1))sk=J(ωk+1)Trk+1 −J(ωk)Tr(ωk) obtenemos que Bk+1 debe cumplir Bk+1sk=J(ωk+1)Tr(ωk+1)−J(ωk)Tr(ωk)−J(ωk+1)TJ(ωk+1)sk:= ˜ yk. De esta forma llegamos a una nueva condición quasi-Newton. 39 40 APÉNDICE B. QUASI-NEWTON B.1. Método de Betts Siguiendo esta estrategia, en Betts (1976), se describe un nuevo algoritmo donde se propone emplear la siguiente fórmula de rango uno Bk+1 =Bk+(˜ yk−Bksk)(˜ yk−Bksk)T sk(˜ yk−Bksk). Esta fórmula de rango uno ha sido reportada por varios autores, incluyendo entre ellos a Broyden, Murtagh y Sargent. Dado que se utiliza un algoritmo recursivo de rango uno para generar B la estimación completa no estará disponible hasta que no se completen las n iteraciones. La contribución de Bk afecta a la determinación de la dirección sk siempre que las n iteraciones hayan sido realizadas. De lo contrario, la dirección de búsqueda sk se calcula directamente como una dirección de Gauss, lo que equivale a tomar ∇2f(ωk) = J(ωk)TJ(ωk) . Como ocurría con el método NL2SOL es razonable tomar en la primera iteración la aproximación de Gauss. La aproximación de Gauss resulta algo más ventajosa si los residuos son lineales o casi nulos pero este algoritmo continúa siendo muy competitivo respecto al método de Gauss-Newton. La contribución de B se vuelve signicativa para residuos no lineales y puntos donde ri(ω) es grande. El algoritmo resulta efectivo en problemas generales que se caracterizan por residuos distintos de cero, principalmente porque la contribución de B está incluída en la estimación total de la Hessiana. En Levenberg-Marquardt ya se incluía esta matriz B pero se representaba como λI . Aunque no se llegó a probar,los resultados numéricos dados en Betts (1976) tienden a indicar que el algoritmo es cuadráticamente convergente, lo que se debe a la contribución de la matriz B . B.2. Método de Bartholomew-Biggs En Bartholomew-Biggs (1977) se propone un nuevo algoritmo de forma análoga al método NL2SOL expuesto en el Capítulo 1. Con respecto a la cuestión de elegir las fórmulas de actualización necesitamos que estas mantengan a la matriz Bk simétrica, pero no hay razón para suponer que deba ser denida positiva. Sin embargo, es preferible que Bk+1 +JT(ωk+1)J(ωk+1) sí sea denida positiva; de lo contrario (Bk+1 + JT(ωk+1)J(ωk+1))sk+1 =JT(ωk+1)r(ωk+1) puede no dar una dirección de descenso. Las B.3. MÉTODO DE FLETCHER-XU 41 dos modicaciones consideradas fueron la actualización de rango dos dada por Powell Bk+1 =Bk+yk−Bksk)sT k+sk(yk−Bksk)T sT ksk −sksT k(yk−Bksk)Tsk (sT ksk)2 y la fórmula de rango uno propuesta por Wolfe Bk+1 =Bk+(yk−Bksk)(yk−Bksk)T (yk−Bksk)Tsk . Supongamos que la aproximación Bk es exactamente igual a Pm i=1 ri(ωk)∇2ri(ωk) y que todas las subfunciones ri(ωk) son cuadráticas por lo que las matrices ∇2ri(ωk) son constantes. Debido a que, en general, rk+1 6=rk , la verdadera matriz Pm i=1 ri(ωk+1)∇2ri(ωk+1) diere de Bk por una matriz de rango n como resultado de los cambios realizados en ri(ωk) . Por lo tanto, no se puede garantizar, simplemente mediante el uso de una actualización basada en que Bk+1 sea correcto incluso siendo Bk correcto. Se propone entonces una estrategia de actualización basada en el caso especial de que rk+1 =γkrk para algún escalar γk . Si se verican las condiciones expuestas tenemos que Pm i=1 ri(ωk+1)∇2ri(ωk+1) = γkBk . Si denimos γk=rT k+1rk rT krk y reescalamos la matriz Bk tenemos que la matriz Bk+1 debería reejar los cambios realizados en Pm i=1 ri(ωk)∇2ri(ωk) . Se proponen entonces dos modicaciones de las fórmulas de actualización Bk+1 =γkBk+(yk−γkBksk)sT k+sk(yk−γkBksk)T sT ksk −sksT k(yk−γkBksk)Tsk (sT ksk)2 y Bk+1 =γkBk+(yk−γkBksk)(yk−γkBksk)T (yk−γkBksk)Tsk . La cuestión de la mejor actualización para Bk se encuentra todavía abierta. Como mencionábamos al principio, es deseable que la actualización mantenga denida positiva la matriz (Bk+1 +JT(ωk+1)J(ωk+1)) . Sería interesante desarrollar un algoritmo que tuviera la capacidad de cambiar automáticamente de una aproximación Gauss-Newton a una aproximación basada en una de las actualizaciones dadas cuando hay indicios de que el problema a tratar es considerado difícil. B.3. Método de Fletcher-Xu En la Sección 1.1 hemos visto que para un problema de residuos pequeños, el método de Gauss-Newton converge a una velocidad lineal rápida que, con precisión limitada, puede ser preferible a la convergencia superlineal del método BFGS. Fletcher y Al-Baali (1985) proponen un método híbrido, HY1, intentando combinar las mejores características de ambos 48 APÉNDICE C. CÓDIGO PYTHON Bibliografía Bartholomew-Biggs, M. C. (1977). The estimation of the hessian matrix in nonlinear least squares problems with non-zero residual. Mathematical Programming , 12:6780. Betts, J. T. (1976). Solving the nonlinear least-square problem: Application of a general method. Journal of Optimations Theory and Applications , 18(4):469483. Bjorck, A. (1996). Least squares methods. In Ciarlet, P. G. y Lions, J. L., editors, Handbook of Numerical Analysis, Vol. I . North-Holland. Dennis, J. E., Gay, D. M., y Welsch, R. E. (1981). An adaptive nonlinear least-squares algorithm. ACM Transactions on Math. Software , 7(3):348368. Fletcher, R. (1980). Practical Methods of Optimization, Vol. 1, UnconstrainedOptimization . John Wiley and Sons, New York. Fletcher, R. y Al-Baali, M. (1985). Variational methods for nonlinear least squares. Operational Res. Soc. , 36:405421. Fletcher, R. y Xu, C. (1987). Hybrid methods of nonlinear least squares. IMA Journal of Numerical Analysis , 7:371389. Gill, P. E., Murray, W., y Wright, M. H. (1981). Practical Optimization . Academic Press, London. Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science and Engineering , 9:9095. Levenberg, K. (1944). A method for the solution of certain problems in leas squares. Quart. Appl. Math , 2:164168. Marquardt, D. W. (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math , 11(2):431441. 49 50 BIBLIOGRAFÍA Moré, J. J. (1978). The levenberg-marquardt algorithm: implementation and theory. In Watson, G., editor, Lecture Notes in Mathematics 630: Numerical Analysis , Berlín,Heidelberg. Springer-Verlag. Moré, J. J. (1983). Recent developments in algorithms and software for trust region methods. In Bachem, A., Grotschel, M., y Korte, B., editors, Mathematical Programming: The State of the Art , Berlín,Heidelberg. Springer-Verlag. Moré, J. J., Garbow, B. S., y Hilstrom, K. E. (1980). User Guide for Minpack-1 . Argonne National Laboratory, Argonne, Illinois. Powell, M. J. D. (1970). A new algorithm for unconstrained optimization. In Rosen, J., Mangasarian, O., y Ritter, K., editors, Nonlinear Programming , pages 3165, London. Academic Pres. Shen, V. K., Siderius, D. W., Krekelberg, W. P., y Hatch, H. W. (2013). Nist standard reference simulation website. Technical Report 173, National Institute of Standards and Technology, Gaithersburg MD, 20899. http://doi.org/10.18434/T4M88Q . Sun, W. y Yuan, Y.-X. (2006). Optimization Theory and Methods . Springer US. van der Walt, S., Colbert, S. C., y Varoquaux, G. (2011). The NumPy Array: A Structure for Ecient Numerical Computation. Computing in Science and Engineering , 13:2230. Viaño, J. M. y Burguera, M. (2013). Lecciones de Métodos Numérico: 4. Optimización . Andavira Editora, S.L., Santiago de Copostela. Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., y Reddy, T. et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientic Computing in Python. Nature Methods , 17:261272. Voglis, C. y Lagaris, I. E. (2004). A rectangular trust region dogleg approach for unconstrained and bound constrained nonlinear optimization. In Proc. WSEAS International Conference on Applied Mathematics , pages 17.