PONTIFICIA UNIVERSIDAD CATÓLICA DEL PERÚ Escuela de Posgrado Análisis e Implementación del Método de Levenberg-Marquardt para la Estimación de Parámetros Mecánicos de un Altavoz Tesis para obtener el grado académico de Magíster en Física Aplicada que presenta: Víctor Raúl Medina Chávez Asesor: Jorge Nestor Moreno Ruiz Lima, 2017 Resumen Si bien el modelo elástico tradicional de la parte mecánica de un altavoz de radiación directa, consistente en un sistema masa-resorte-amortiguador, es ampliamente utilizado en aplicaciones electroacústicas debido a su simplici- dad y aceptable concordancia (para fines prácticos) con datos experimentales, existen características en el comportamiento del altavoz que dicho modelo no reproduce. Una de esas características, relacionada al fenómeno de creep que presenta la suspensión mecánica viscoelástica de los altavoces, es el incre- mento de los valores de la resistencia mecánica en el rango de las frecuencias bajas. El presente trabajo presenta un modelo viscoelástico de la parte me- cánica de un altavoz de radiación directa y muestra como este resulta ser un modelo más completo respecto al modelo elástico tradicional al reproducir el mencionado comportamiento del altavoz asociado a la presencia del creep. Este resultado comparativo se valida mediante el ajuste del modelo a las cur- vas de datos experimentales obtenidos de la medición de distintos altavoces. El método de Levenberg–Marquardt, empleado para el ajuste del modelo por mínimos cuadrados no lineales, es estudiado con detalle. Introducción El presente trabajo de tesis se desarrolla en dos planos. Uno general que aborda el tratamiento matemático del criterio de los mínimos cuadrados y del método de Levenberg–Marquardt, y otro específico que muestra la apli- cación de dicho método al caso particular de la estimación experimental de los parámetros de un modelo viscoelástico de la suspensión mecánica de un altavoz de radiación directa. Por lo general, la aplicación de los mínimos cuadrados presentada a los estudiantes de ingeniería se limita al problema trivial de ajustar una recta a un conjunto de datos sujetos a errores de medición, lo que en última instancia se reduce a resolver un sistema de ecuaciones lineales. No es poco común, además, que los mínimos cuadrados se introduzcan como axioma, o que sus fundamentos matemáticos se discutan de forma superficial. En este contexto, raramente se tratan problemas reales que involucren modelos no lineales, y cuando esto ocurre la solución recae en el uso de software comercial. Este trabajo intenta cubrir esta brecha exponiendo algunos aspectos relevantes del tema y ofreciendo referencias a literatura complementaria. La justificación matemática de los mínimos cuadrados como criterio para el ajuste óptimo de un modelo a un conjunto de datos experimentales se analiza en el capítulo 2. Los primeros intentos por darle un marco lógico a los mínimos cuadrados se basaron en el análisis probabilístico del error aleatorio de medición. Ejemplo de ello son las deducciones de Gauss y Laplace tratadas en este capítulo, las cuales se basan en la ley normal del error y el teorema del límite central respectivamente. En la última sección de este capítulo se desarrolla una expresión general para la discrepancia entre el modelo y los datos medidos y se muestra que los minimos cuadrados son un caso particular. Este desarrollo se basa en ciertos requerimientos de naturaleza elemental y no considera las propiedades del iv error de medición. Una vez establecidos los mínimos cuadrados como criterio de ajuste el si- guiente paso consiste en diseñar un procedimiento que implemente tal ajuste. Dos métodos precursores que cumplen con tal fin, el del máximo descenso y el de Gauss–Newton, son tratados con regular detalle en el capítulo 3. Sin embargo, no es poco frecuente que dichos métodos y sus posteriores variantes presenten problemas de convergencia cuando los modelos son complejos. Un método que hereda las ventajas de los dos anteriores y supera el problema de convergencia es el de Levenberg–Marquardt analizado en el capítulo 4. Por su eficiencia y confiabilidad, este método se ha convertido en la herramienta estándar para el ajuste de modelos no lineales. Yendo al plano específico, los capítulos 1 y 5 muestran una aplicación del método de Levenberg–Marquardt al campo de la electroacústica. Este mate- rial puede ser de utilidad tanto para quien busque un complemento práctico de la teoría, como para el acústico interesado en los detalles específicos que involucra la medición de los elementos presentes en un modelo viscoelástico de la parte mecánica de un altavoz de radiación directa. Dada la naturaleza de los materiales usados en la construcción de los alta- voces, la resistencia mecánica asociada a las pérdidas de energía varía con la frecuencia. Este hecho, comprobado experimentalmente, sugiere la necesidad de ampliar el modelo tradicional de resistencia mecánica constante. El mode- lo viscoelástico propuesto en el capítulo 1 cumple este propósito y reproduce con mayor fidelidad el comportamiento observado experimentalmente. El ejemplo de aplicación se completa en el capítulo 5 donde se muestran los resultados obtenidos al emplear el método de Levenberg–Marquardt en la estimación de los parámetros mecánicos de un altavoz de muestra según el modelo viscoelástico presentado en el primer capítulo. En términos generales, la importancia del método de Levenberg–Marquardt radica en el proceso de ajuste de modelos, el cual permite al investigador teó- rico contrastar la exactitud de sus modelos contra los datos experimentales para así saber qué partes de estos se comprueban correctas y que otras ne- cesitan ser refinadas. Al investigador práctico, por otro lado, le ofrece la posibilidad de extraer de forma eficiente información específica respecto a los parámetros o elementos de un sistema a partir de un conjunto de datos v medidos. En cuanto al rango de aplicación, el método es potencialmente aplicable a cualquier sistema susceptible de ser modelado matemáticamente. Son di- versas las áreas en las que el método es comúnmente empleado. La industría química, por ejemplo, extrae información respecto a las cantidades óptimas de los componentes que intervienen en una reacción química basándose en modelos ajustados experimentalmente. Otros ejemplos se dan en campos co- mo las finanzas, en las que se ajustan los modelos de las series de tiempo correspondientes a los precios de las acciones. El método también se aplica en los modelos de riesgos de las empresas aseguradoras, en la optimización de técnicas agropecuarias, en sistemas de diagnóstico clínico de ciertas en- fermedades crónicas, y en general en casi todos los ámbitos de las ciencias e ingeniería. El método de Levenberg–Marquardt es sin duda uno de los implementos más importantes dentro de la caja de herramientas matemáticas con que cuenta la física aplicada. Índice general Resumen II Introducción III 1. Parámetros Mecánicos de un Altavoz según el Modelo TGM 1 1.1. Modelo Tradicional del Altavoz . . . . . . . . . . . . . . . . . 2 1.2. Modelo TGM de la suspensión . . . . . . . . . . . . . . . . . . 5 1.3. Arreglo Experimental . . . . . . . . . . . . . . . . . . . . . . . 8 2. Mínimos Cuadrados como Criterio para el Ajuste Óptimo de Modelos 10 2.1. Formulación del Problema . . . . . . . . . . . . . . . . . . . . 10 2.2. Deducción basada en la Teoría de la Probabilidad . . . . . . . 13 2.2.1. Modelo Lineal y Reducción de Modelos No Lineales . . 14 2.2.2. Enunciado de Gauss de la Media Aritmética . . . . . . 17 2.2.3. Demostración de Gauss de la Ley Normal del Error . . 18 2.2.4. Deducción de Gauss de los Mínimos Cuadrados . . . . 20 2.2.5. Deducción Alterna de Laplace . . . . . . . . . . . . . . 22 2.3. Interpretación Geométrico Vectorial del Método . . . . . . . . 25 2.3.1. Obtención de las Ecuaciones Normales Mediante Mé- todos Vectoriales . . . . . . . . . . . . . . . . . . . . . 26 2.3.2. Evaluación de la Confiabilidad de los Parámetros Es- timados . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4. Mínimos Cuadrados y Modelos Aproximados . . . . . . . . . . 33 2.4.1. Expresión General para la Función de Discrepancia . . 34 2.4.2. Ejemplos Particulares de la Función de Discrepancia . 37 vii 3. Algoritmos Precursores: Métodos del Máximo Descenso y de Gauss-Newton 39 3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2. Método del Máximo Descenso . . . . . . . . . . . . . . . . . . 42 3.2.1. Elección de la Dirección de Minimización . . . . . . . . 42 3.2.2. Reescalamiento de la Función a Minimizar . . . . . . . 45 3.2.3. Cálculo de la Longitud del Paso . . . . . . . . . . . . . 50 3.2.4. Demostración de la Convergencia del Método . . . . . 52 3.3. Método de Gauss-Newton . . . . . . . . . . . . . . . . . . . . 61 3.3.1. Elección de la Dirección de Minimización . . . . . . . . 61 3.3.2. Cálculo de la Longitud del Paso . . . . . . . . . . . . . 63 3.3.3. Demostración de la Convergencia del Método . . . . . 64 4. Método de Levenberg-Marquardt 68 4.1. Enfoque de Levenberg . . . . . . . . . . . . . . . . . . . . . . 69 4.1.1. Construcción del Método . . . . . . . . . . . . . . . . . 69 4.1.2. Demostración de la Utilidad del Método . . . . . . . . 74 4.1.3. Cálculo del Vector de Paso . . . . . . . . . . . . . . . . 79 4.2. Enfoque de Marquardt . . . . . . . . . . . . . . . . . . . . . . 80 4.2.1. Base Teórica del Método . . . . . . . . . . . . . . . . . 80 4.2.2. Escalamiento del Espacio de Parámetros . . . . . . . . 84 4.2.3. Construcción del Algoritmo . . . . . . . . . . . . . . . 85 4.3. Diferencia entre los Dos Enfoques . . . . . . . . . . . . . . . . 88 5. Estimación de Parámetros y Conclusiones 90 5.1. Procedimiento de Estimación de Parámetros . . . . . . . . . . 90 5.2. Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Bibliografía 95 Capítulo 1 Parámetros Mecánicos de un Altavoz según el Modelo TGM El altavoz, visto como sistema electro-mécano-acústico formado por blo- ques que agrupan elementos de un mismo comportamiento y función, puede ser caracterizado por un modelo discreto de parámetros concentrados cuya complejidad dependerá de la exactitud que se busque obtener. La mayoría de estos modelos son válidos únicamente bajo dos condiciones de operación: 1. Régimen lineal del altavoz, obtenido a pequeños niveles de señal de entrada, con lo cual la presencia de las características no lineales tanto del factor de fuerza como de la compliancia mecánica son mínimas y despreciables. 2. Rango de frecuencias de pistón rígido, ubicado por debajo de las fre- cuencias donde empiezan a presentarse los modos del cono, lo que per- mite aproximar al sistema móvil del altavoz como un cuerpo rígido. En el presente capítulo se introduce y analiza el modelo TGM (Truncated Generalized Maxwell) de la suspensión mecánica de un altavoz, el cual amplía el modelo tradicional incorporando el carácter dependiente de la frecuencia que presenta la parte real de la impedancia mecánica. Asímismo, se describe el procedimiento para estimar el valor de los parámetros que mejor ajusten el modelo a los datos experimentales de un altavoz en particular. Finalmen- te, se exponen detalles técnicos del arreglo de medición, y se comparan los resultados experimentales con aquellos simulados por el modelo. 1 2 1.1. Modelo Tradicional del Altavoz Antes de tratar el modelo TGM mencionado, conviene hacer una breve descripción del modelo tradicional. Esto permite introducir al modelo TGM como una ampliación del modelo tradicional. El modelo tradicional lineal del altavoz ha demostrado ser, a pesar de su simplicidad, una buena aproximación en la mayoría de los casos. El pequeño número de parámetros que contiene, la evidente interpretación física de los mismos, así como la disponibilidad de numerosos métodos para medirlos, tanto en el dominio del tiempo como de la frecuencia [28], [29], [18] y [24], hacen que el uso de este modelo sea común en aplicaciones prácticas. En la figura 1.1 se muestra el esquema de un altavoz dinámico con el fin de facilitar la identificación de los elementos que conforman el modelo. El funcionamiento de este tipo de altavoces se puede simplificar de la siguiente manera: la corriente aplicada a través de los terminales de la bobina circula perpendicularmente a las líneas de campo magnético presentes en el entre- hierro que completa el circuito magnético conformado por el polo, el imán, y las placas frontal y posterior. Las fuerzas de Lorentz resultantes actúan sobre las cargas que fluyen por la bobina. Esto genera el movimiento axial del diafragma formado por el cono y la tapa, ambos adheridos a la bobina y centrados por el sistema de suspensión compuesto por la araña y el anillo. Este movimiento se transmite a las partículas de aire circundantes y de ahí se propaga al resto del medio en forma de onda acústica. Cabe mencionar que el aire ejerce sobre el diafragma una fuerza de resistencia al movimiento. araña cono tapa antipolvo anillo canasta bobina placa frontal placa posterior imán poloducto de ventilación entrehierro Figura 1.1: Esquema de la sección transversal de un altavoz dinámico. 3 cuerpo rígido suspensión sistema móvil resorte ideal (k) amortiguador ideal (r) xfuerza restitutiva = - k x fuerza disipativa = - r x 0 fuerza motriz = (Bl) i . Figura 1.2: Modelo tradicional de la parte mecánica del altavoz. Este efecto se toma en cuenta introduciendo en el modelo lo que se conoce como impedancia de radiación, que por simplificación muchas veces es des- cartada sin que esto afecte apreciablemente el desempeño del modelo, ya que su influencia es pequeña en comparación con los otros parámetros. En el modelo tradicional [1] se asume que la fuerza motriz desarrollada en la bobina tiene un valor igual a la corriente que la atraviesa multipli- cada por una constante denominada factor de fuerza (Bl), donde ‘B’ es la densidad de flujo magnético en el entrehierro y ‘l’ la longitud de la bobina. Asimismo, se asume que las fuerzas restitutivas y disipativas ejercidas por la suspensión sobre el diafragma son respectivamente proporcionales al despla- zamiento y la velocidad del mismo. Para la parte eléctrica el modelo asume un comportamiento ideal de la bobina y de su resistencia intrínseca. Bajo estos supuestos, la parte mecánica del modelo tradicional está con- formada por un resorte más un amortiguador lineales y sin masa que carac- terizan la suspensión, y por un cuerpo rígido que caracteriza al sistema móvil que incluye al diafragma, a la bobina y a la parte de la suspensión sujeta a movimiento. La conexión entre estos elementos se indica en la figura 1.2. Considerando algunas leyes físicas básicas (ley de voltajes de Kirchhoff y segunda ley de Newton) es posible derivar la expresión matemática que R L m r c (Bl ) x + - (Bl ) i + - v + - i Bl x Figura 1.3: Circuito eléctrico equivalente. 4 1 R+(jω)L V( )ω Bl (j )m + r + 1/ ω (jω)c 1 Bl X( )ωI( )ω F ( )ω V ( )ω b b Figura 1.4: Diagrama de bloques del modelo tradicional. describe, de acuerdo al modelo, los procesos eléctricos y mecánicos presentes en el altavoz. Esta expresión está compuesta por dos ecuaciones diferenciales, donde las funciones se pueden conocer experimentalmente y los coeficientes son las incógnitas a determinar: v(t) = R i(t) + L di(t) dt + vb(t) vb(t) = Bl dx(t) dt , (1.1) fb(t) = m d2x(t) dt2 + r dx(t) dt + 1 c x(t) fb(t) = Bl i(t). (1.2) La definición de los símbolos se puede ver en el cuadro 1.1. Estas ecuaciones son idénticas a las que se obtienen del circuito eléctrico equivalente mostrado en la figura 1.3. En este circuito se utiliza un girador, el cual es el componente ideal que separa la parte eléctrica de la mecánica para así enfatizar la interdependencia entre ambas. Dicho en términos generales, un girador es un elemento discreto eléctrico ideal de dos puertos cuya ley de funcionamiento es la siguiente: el voltaje presente en los terminales del pri- mer puerto es proporcional a la corriente que fluye por el segundo puerto y, similarmente, el voltaje presente en los terminales del segundo puerto es pro- porcional a la corriente que fluye por el primer puerto. Este comportamiento puede observarse en la figura 1.3. Vistas en el dominio de la frecuencia, dichas ecuaciones se pueden repre- sentar mediante el diagrama de bloques de la figura 1.4, donde las mayúsculas denotan la transformada de Fourier; e.g., X(ω) = F{x(t)}. La ventaja del diagrama de bloques es que muestra en forma sencilla la relación entre las señales. Esto permite que tanto las funciones de transferen- cia como las impedancias puedan ser derivadas directamente. Así por ejemplo 5 v(t) Voltaje aplicado [V] i(t) Corriente en la bobina [A] R Resistencia de la bobina [Ω] L Inductancia de la bobina [H] Bl Factor de fuerza [N/A] vb(t) Voltaje inducido en la bobina [V] fb(t) Fuerza motriz desarrollada en la bobina [N] x(t) Desplazamiento de la bobina [m] m Masa del sistema móvil más la carga de aire [kg] r Coeficiente de amortiguamiento viscoso de la suspensión [N s/m] c Compliancia de la suspensión [m/N] k Coeficiente de elasticidad de la suspensión (= 1/c) [N/m] t Tiempo [s] ω Frecuencia angular [rad/s] Cuadro 1.1: Nomenclatura (modelo tradicional). se obtiene la función de transferencia del altavoz I(ω) V (ω) = 1 R + (jω)L+ (Bl)2 (jω)m+ r + 1/(jω)c , (1.3) y la impedancia mecánica Fb(ω) Ẋ(ω) = ZM(ω) = (jω)m+ r + 1/(jω)c. (1.4) 1.2. Modelo TGM de la suspensión En la sección anterior se vio cómo el modelo tradicional posee un amor- tiguador viscoso como único elemento disipador de energía (elemento que ejerce una fuerza que se opone al movimiento). En consecuencia, la expre- sión de la resistencia mecánica, definida como la parte real de (1.4)) y por ende responsable de las pérdidas de energía en el sistema, es una constante r cuyo valor es el coeficiente de amortiguación viscoso. Sin embargo, esta aproximación básica adoptada por el modelo tradicio- nal respecto a la resistencia mecánica no reproduce satisfactoriamente los resultados obtenidos experimentalmente, como se menciona en [13], donde se 6 Muestra Modelo del Altavoz Diámetro 1 Peerless 831727 10 in 2 JBL 2241H 18 in 3 JBL 2245H 18 in 4 JBL 2206H 12 in 5 JBL 2226H 15 in Figura 1.5: Resistencia Mecánica (dividida entre Bl) de los altavoces en el cuadro con- tiguo. midió la resistencia mecánica dividida entre el factor de fuerza de distintos altavoces. El resultado de esas mediciones se muestra en la figura 1.5. En ella se observa que en realidad la resistencia mecánica es dependiente de la frecuencia y se aleja notoriamente del valor constante asumido. Esta limitación motiva la elaboración de un modelo más complejo que incluya propiedades físicas adicionales a la de restitución y amortiguamiento que presenta la suspensión. Una de tales propiedades es la que se conoce como creep, comúnmente presente en materiales viscoelásticos como los que m r2 c1 F - r1 c2 + Figura 1.6: Circuito eléctrico equivalente. F Fuerza motriz [V] m Masa del sistema móvil más la carga de aire [kg] r1 Resistencia mecánica de la rama Voigt del modelo TGM [N·s/m] c1 Compliancia mecánica de la rama Voigt del modelo TGM [m/N] r2 Resistencia mecánica de la rama Maxwell del modelo TGM [N·s/m] c2 Compliancia mecánica de la rama Maxwell del modelo TGM [m/N] Cuadro 1.2: Nomenclatura (modelo TGM). 7 0 50 100 150 200 250 5 10 15 20 25 30 Mechanical Impedance, Real Part Hz K g /s 0 50 100 150 200 250 -300 -200 -100 0 100 200 Mechanical Impedance, Imaginary Part Hz K g /s standard elastic model visco-elastic model standard elastic model visco-elastic model Figura 1.7: Comparación de los resultados obtenidos para la impedancia mecánica cal- culada por el método tradicional y el método TGM. componen la suspensión. El fenómeno de creep se manifiesta como un conti- nuo y lento desplazamiento del diafragma bajo una fuerza constante aplicada, hecho que impide alcanzar el estado estacionario esperado. La presencia de creep en los altavoces fue reportada por primera vez en [6] y modelada según ciertas estructuras en [14]. Un modelo particular, entre los diversos modelos existentes para caracteri- zar materiales que exhiben creep, que ha probado ser apropiado para modelar la suspensión (resistencia y compliancia mecánica) es el modelo generalizado de Maxwell truncado, o TGM por sus siglas en inglés. Reemplazando por tanto el modelo tradicional de la suspensión por el modelo TGM, se consigue un modelo mejorado de la impedancia mecánica del altavoz. Para este caso, el circuito eléctrico equivalente de la parte mecánica, ya sin incluir el girador, se muestra en la figura 1.6. En la expresión correspondiente a la impedancia mecánica del modelo TGM ZM = r1 + r2 1 + (ωc2r2)2 + j [ (ωm− 1 ωc1 )− ωc2r 2 2 1 + (ωc2r2)2 ] , (1.5) puede observarse que la parte real ya no es constante sino que depende de la frecuencia. Nótese en la figura 1.7 que el efecto de los parámetros r2 y c2, introducidos por el modelo TGM (a diferencia del modelo tradicional), se manifiesta en el rango de las frecuencias bajas de la parte real de la impe- 8 dancia mecánica, reproduciendo el comportamiento asociado al creep que se desea modelar. El efecto de r2 y c2 es en cambio irrelevante en todo el rango de frecuencias en la parte imaginaria de la impedancia mecánica, la cual es dominada por los parámetros del modelo tradicional. Para ilustrar las diferencias entre el modelo tradicional y el TGM, en la figura 1.7 se compara tanto la parte real como la imaginaria de las impe- dancias mecánicas simuladas según cada modelo. La impedancia medida del altavoz no se muestra en la figura. De las gráficas puede observarse que a diferencia del modelo tradicional el modelo TGM presenta un incremento de la resistencia mecánica en el rango de las frecuencias bajas similar al que presentan los datos experimentales (figura 1.5). Esta característica, como se verá más adelante, permite una mejor coincidencia entre los datos simulados con el modelo TGM y los datos obtenidos experimentalmente. 1.3. Arreglo Experimental En la figura 1.8 se muestra el arreglo experimental empleado en este trabajo para la toma de datos. El altavoz se sujeta verticalmente mediante una abrazadera de bronce que evita fugas de flujo magnético. Como señal de entrada se utiliza un barrido sinusoidal generado por el analizador y transmitido al altavoz a través de un amplificador de potencia. El nivel de esta señal debe ser suficientemente pe- queño como para no excitar las no linealidades del altavoz y suficientemente grande como para alcanzar una buena relación señal a ruido. Tres señales, conectadas a las entradas del analizador multicanal, son medidas simultá- neamente: el voltaje en los terminales del altavoz, el voltaje a través de una resistencia de 1Ω en serie con el altavoz y por ende la corriente que circu- la por la bobina, y la velocidad del cono mediante un transductor láser de velocidad. Las componentes de frecuencia de estas tres señales son medidas por el analizador mediante la técnica SSR (Steady State Response) que determina la amplitud y fase de cada señal medida tomando como referencia la señal armónica de entrada y descartando los primeros periodos para conseguir que los transitorios se extingan. 9 Figura 1.8: Arreglo experimental Las impedancias eléctrica y mecánica del altavoz se calculan a partir de las componentes de frecuencia ya medidas de la velocidad, corriente, y voltaje. Asimismo, es posible calcular el factor de fuerza siguiendo el procedimiento descrito en [24]. La definición teórica de la impedancia mecánica del altavoz es la razón entre la fuerza ejercida en la bobina y la velocidad del cono. Sin embargo, en este caso la fuerza no se mide directamente sino que se calcula asumiendo idealmente que la misma es igual a la corriente multiplicada por el factor de fuerza. La validez de este supuesto ya se ha verificado en anteriores investi- gaciones. Para ello se mide la impedancia mecánica sin hacer uso del motor del altavoz, i.e., aplicando y midiendo directamente una fuerza externa al sistema móvil del altavoz. En [25] se dan detalles de un arreglo experimental usado con ese propósito. Capítulo 2 Mínimos Cuadrados como Criterio para el Ajuste Óptimo de Modelos En las diversas ramas de la ciencia se presenta con frecuencia el pro- blema consistente en ajustar un modelo matemático parametrado a datos observados experimentalmente de cierto sistema. Para llevar a cabo el ajus- te es necesario establecer un criterio que evalúe el grado de concordancia, o discrepancia, entre el modelo y los datos experimentales. El presente capítulo analiza uno de esos criterios en particular, denomi- nado método de los mínimos cuadrados, y justifica el uso del mismo para conseguir el ajuste óptimo de modelos. En la primera sección se introduce una descripción del método; en la se- gunda, se dan detalles de la derivación del método hecha por Gauss [9] con base en la teoría matemática de la probabilidad; en la tercera, se presenta una interpretación del método propuesta por Kolmogorov [15] dentro del campo del álgebra lineal vectorial; y finalmente en la cuarta, se propone un enfo- que distinto del problema, apartándolo de su interpretación probabilística convencional y asociándolo con algunos conceptos de espacios normados. 2.1. Formulación del Problema En términos generales se puede definir un sistema como una caja que admite un conjunto de señales de entrada, las procesa, y produce como res- 10 11 puesta un grupo de señales de salida. Cuando este proceso está libre de elementos aleatorios, produce siempre las mismas salidas en respuesta a unas mismas entradas y se dice que el sistema es determinista. En este caso, las entradas y salidas del sistema están relacionadas de manera precisa mediante una expresión matemática que representa el modelo exacto del sistema. Contar con una expresión matemática que describa el comportamiento de cierto sistema resulta muy útil en el proceso de diseño, análisis y optimización del mismo. Esta expresión o modelo matemático, al ser una forma equivalente del sistema, permite mediante operaciones simples simular la respuesta del sistema ante situaciones hipotéticas. Sin embargo, conocer el modelo exacto de un sistema es solo un caso ideal, pues no es posible conocer la totalidad de elementos que intervienen en un sistema real. En la práctica se requiere un modelo que aproxime al sistema con un grado de exactitud suficiente. El proceso de construcción de un modelo implica identificar los elementos representativos que conforman un sistema. Un modelo es ajustable cuando posee parámetros que pueden asumir va- lores arbitrarios que permitan adaptar el modelo a distintos sistemas particu- lares. Así, en términos más precisos, el proceso de ajuste consiste en estimar los valores particulares de los parámetros que produzcan un mínimo grado de discrepancia entre las respuestas observadas del sistema y las reproducidas por el modelo. La figura 2.1 ilustra el proceso de generación de las respuestas reales (y∗), medidas (y) y predichas (ŷ). Dado que todo proceso de medición introdu- ce inevitablemente una componente de error aleatorio (∆), las respuestas medidas poseen la forma y = y∗ + ∆. 12 MODELO EXACTO Parámetros de ajuste b = (b ,b ,…,b )1 2 k SISTEMA entrada x = (x ,x ,…,x )1 2 m MODELO APROXIMADO Parámetros de ajuste b = (b ,b ,…,b )1 2 k A B MEDICIÓN Error de medición Δ respuesta real *y CÁLCULO CÁLCULO DATOS MEDIDOS (y ) i = 1,2,...,ni DATOS PREDICHOS (ŷ ) i = 1,2,...,ni DATOS PREDICHOS (ŷ ) i = 1,2,...,ni y = *y + Δ ŷ ŷ discrepancia por error aleatorio discrepancia por error sistemático Figura 2.1: Proceso de generación de datos medidos (y) y predichos (ŷ). Dos casos se muestran en esta figura: Caso A: se tiene el modelo exacto del sistema, por lo que es posible conse- guir que las respuestas predichas coincidan con las reales al asignarles valores apropiados a los parámetros de ajuste. En este caso, la discre- pancia entre las respuestas medidas y predichas se debe únicamente al error aleatorio de medición y, por ende, el criterio para evaluarla debe fundarse en la naturaleza aleatoria del error. Caso B: se tiene un modelo aproximado del sistema, por lo que aun cuan- do se estimen los valores óptimos de los parámetros de ajuste, habrá siempre una diferencia entre las respuestas reales y predichas. En este caso, la discrepancia se debe al error sistemático introducido por el uso de un modelo aproximado y, por ende, el criterio para evaluarla debe considerar la naturaleza sistemática del error. Para expresar matemáticamente lo anteriormente dicho, se considera el mo- delo ŷi = f(x1i, x2i, . . . , xmi; b1, b2, . . . , bk) (2.1) = f(xi, b) = fi(b), donde xi = (x1i, x2i, . . . , xmi) es el i-ésimo vector de entradas o variables independientes, b= (b1, b2, . . . , bk) es el vector de parámetros de ajuste, e ŷi es la i-ésima respuesta escalar predicha por el modelo. Asimismo, se define 13 el residuo (ri) para cada uno de los n datos medidos yi = y∗i + ∆i, i = 1, 2, . . . , n, como la diferencia entre las correspondientes respuestas medida y predicha ri = yi − ŷi, i = 1, 2, . . . , n; (2.2) y se asume que todo criterio de ajuste legítimo define la discrepancia entre los datos medidos y predichos en función de tales residuos. En particular, el método de los mínimos cuadrados le asigna a la discre- pancia Φ un valor igual a la sumatoria de los cuadrados de los residuos Φ = n∑ i=1 r2 i , (2.3) de modo que el vector de parámetros de ajuste bmin que minimice Φ debe ser solución del sistema de ecuaciones ∂Φ ∂bj = 2 n∑ i=1 ri ∂ri ∂bj = 0, j = 1, 2, . . . , k. (2.4) En este punto es importante anotar dos hechos. Primero, si algún b0 es solución del sistema de ecuaciones ri(b1, b2, . . . , bk) = 0, i = 1, 2, . . . , n, (2.5) dicho b0 es también solución de (2.4). Segundo, para el caso en que los ri son una combinación lineal de las componentes de b, se tiene que (2.5) posee infinitas soluciones si n < k. En dicho caso, n ≥ k es condición necesaria para que (2.4) posea una solución única bmin. 2.2. Deducción de Gauss basada en la Teoría de la Probabilidad La primera mención del método de los mínimos cuadrados aparece en una publicación de Legendre de 1806 [19], en la que el método es propuesto sin 14 fundamento matemático y a manera de principio. Tres años después, Gauss publica la primera deducción del método [9], fundada en principios de natu- raleza más elemental, en la que conecta el método con la teoría matemática de la probabilidad. En esta sección se analizan los detalles de la deducción de Gauss y se expone al final una deducción alterna sugerida por Laplace. 2.2.1. Modelo Lineal y Reducción de Modelos No Lineales a la Forma Lineal Como punto de partida se presenta un ejemplo de aplicación del método. Se considera para ello el caso más sencillo en el que el modelo, expresado en (2.1), es lineal en los parámetros que se busca estimar. Se tiene entonces fi(b) = c1ib1 + c2ib2 + · · ·+ ckibk, (2.6) de donde se desprende la forma de los residuos ri = yi − k∑ j=1 cjibj. (2.7) El vector de parámetros bmin es en este caso la solución única del sistema de ecuaciones conocido como ecuaciones normales [c1c1]b1 + [c1c2]b2 + . . . + [c1ck]bk = [c1y] [c1c2]b1 + [c2c2]b2 + . . . + [c2ck]bk = [c2y] ... ... . . . ... ... [c1ck]b1 + [c2ck]b2 + . . . + [ckck]bk = [cky]  (2.8) obtenido al introducir (2.7) en (2.4). Cabe señalar que los corchetes [ ] en la notación de Gauss representan el producto interno de los vectores contenidos en ellos, i.e., [cjch] = n∑ i=1 cjichi, j, h = 1, 2, . . . , k. En la sección 2.3.1 se analizan las condiciones que deben cumplir los vectores c1, c2, . . . , ck para garantizar que (2.8) posea una única solución. Aplicando técnicas de álgebra elemental es posible obtener la solución 15 de (2.8). Una de estas técnicas, basada en el uso de determinantes, permite expresar de forma explícita cada componente (bmin)j del vector solución en función de los coeficientes presentes en el sistema de ecuaciones. En particular para obtener (bmin)1 se define A =  [c1c1] [c1c2] . . . [c1ck] [c1c2] [c2c2] . . . [c2ck] ... ... . . . ... [c1ck] [c2ck] . . . [ckck]  y g =  [c1y] [c2y] ... [cky]  , (2.9) con lo cual (2.8) toma la forma vectorial A b = g. (2.8a) Si b = bmin es solución de esta ecuación, entonces también lo será de la ecuación A b = b1 (bmin)1 g, (2.10) puesto que para b=bmin ambas ecuaciones son idénticas al ser b1/(bmin)1 =1. Definiendo A1 =  [c1c1]− [c1y] (bmin)1 [c1c2] . . . [c1ck] [c1c2]− [c2y] (bmin)1 [c2c2] . . . [c2ck] ... ... . . . ... [c1ck]− [cky] (bmin)1 [c2ck] . . . [ckck]  , se expresa (2.10) en forma matricial A1b = 0, (2.10a) la cual hace evidente que aparte de la solución trivial b= 0 existen infinitas soluciones de la forma b = α bmin, α ∈ R, por ende A1 es singular y tanto su determinante como la determinante de su 16 transpuesta AT 1 son cero. Se tiene entonces la ecuación |AT 1 | = ∣∣∣∣∣∣∣∣∣∣∣∣∣ [c1c1]− [c1y] (bmin)1 [c1c2]− [c2y] (bmin)1 . . . [ckck]− [cky] (bmin)1 [c1c2] [c2c2] . . . [c2ck] ... ... . . . ... [c1ck] [c2ck] . . . [ckck] ∣∣∣∣∣∣∣∣∣∣∣∣∣ = 0, (2.11) de donde se despeja el valor de (bmin)1 (bmin)1 = ∣∣∣∣∣∣∣∣∣∣∣∣ [c1y] [c2y] . . . [cky] [c1c2] [c2c2] . . . [c2ck] ... ... . . . ... [c1ck] [c2ck] . . . [ckck] ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ [c1c1] [c1c2] . . . [c1ck] [c1c2] [c2c2] . . . [c2ck] ... ... . . . ... [c1ck] [c2ck] . . . [ckck] ∣∣∣∣∣∣∣∣∣∣∣∣ . (2.12) Las demás componentes del vector solución bmin se pueden obtener siguiendo el mismo procedimiento. Para el caso general en el que el modelo es no lineal en los parámetros, el tratamiento estándar consiste en reducir el modelo fi(b) a la forma lineal 〈fi(b)〉 mediante su aproximación de Taylor de primer grado alrededor de un punto inicial b0. Introduciendo δ = (b− b0) se tiene 〈yi〉 = 〈fi(b)〉 = fi(b0) + k∑ j=1 ∂fi(b0) ∂bj δj, (2.13) donde 〈yi〉 es el i-ésimo dato predicho por el modelo linealizado. El residuo correspondiente a esta versión aproximada del modelo se define como 〈ri〉 = [yi − fi(b0)]− k∑ j=1 ∂fi(b0) ∂bj δj, (2.14) el cual, al igual que (2.7), posee forma lineal. A partir de este punto es posible seguir el procedimiento descrito al inicio 17 para obtener el vector solución bmin; para ello basta hacer la identificación cji = ∂fi(b0) ∂bj . 2.2.2. Enunciado de Gauss de la Media Aritmética Para deducir la ley de probabilidad del error aleatorio de medición, Gauss utilizó el siguiente enunciado: dado cualquier número de medidas y1, y2, y3, . . ., de una cantidad desconocida η, el valor η∗ con mayor probabilidad de ser di- cho η desconocido, es la media aritmética de las medidas. Cabe señalar que el valor real de η, denotado por y∗, no es necesariamente igual a η∗. Este enunciado se puede derivar a partir de propiedades y axiomas de naturaleza más evidente, formulados de la siguiente manera: Propiedad I - Las diferencias entre el valor más probable y las medidas individuales no dependen de la posición del punto cero desde el cual son calculadas. Propiedad II - La razón entre el valor más probable y cualquier medida individual no depende de la unidad de medida utilizada. Axioma I - El valor más probable es independiente del orden en el que se realizan las medidas, por lo tanto es una función simétrica de las medidas. Axioma II - El valor más probable, considerado como una función de las medidas individuales, posee primeras derivadas continuas con respecto a ellas. Si se expresa el valor más probable η∗ en términos de las smedidas y1, y2, . . . , ys mediante la función g(y1, y2, . . . , ys), entonces por el teorema del valor medio, aplicable gracias al axioma II, se tiene g(y) = g(y1, y2, . . . , ys) = g(0, 0, . . . , 0) + y1 [ ∂g ∂y1 ] + y2 [ ∂g ∂y2 ] + . . .+ ys [ ∂g ∂ys ] , (2.15) donde los corchetes denotan que las derivadas parciales están evaluadas en el punto ty, con t ∈ [0, 1]. 18 Según la propiedad II se tiene que la ecuación g(ky1, ky2, . . . , kys) = k g(y1, y2, . . . , ys) debe ser válida para todo k 6= 0. Evaluando en particular en el límite k → 0, y considerando la continuidad de g, se obtiene g(0, 0, . . . , 0) = 0, y así (2.15) se transforma en g(y1, y2, . . . , ys) = α1(y)y1 + α2(y) y2 + . . .+ αs(y) ys. De acuerdo al axioma I, todos los α′s deben ser iguales, por lo tanto g(y1, y2, . . . , ys) = α(y)[y1 + y2 + . . .+ ys]. Finalmente, la propiedad I se traduce en g(y1 + h, y2 + h, . . . , ys + h) = g(y1, y2, . . . , ys) + h, de donde se desprende α(y) = 1 s . De estos resultados se concluye que la forma de g es η∗ = g(y1, y2, . . . , ys) = 1 s (y1 + y2 + . . .+ ys), (2.16) que corresponde a la media aritmética de las medidas. 2.2.3. Demostración de Gauss de la Ley Normal del Error A partir del enunciado de la media aritmética es posible deducir que la densidad de probabilidad del error de medición ∆ sigue una ley normal. Si p(∆) es la función de densidad de probabilidad de ∆, y ε es la me- nor cantidad a la cual el instrumento de medición es sensible, entonces la probabilidad de obtener un error de valor ∆0 es p(∆0)ε. 19 Si para una cierta cantidad η, cuyo valor real es y∗, se tienen s medidas y1, y2, . . . , ys en las que se introducen los errores ∆1 =y1−y∗, ∆2 =y2−y∗, . . . , ∆s=ys−y∗, entonces la probabilidad de que ocurran dichas medidas es εsp(y1 − y∗)p(y2 − y∗) . . . p(ys − y∗). Si se asume, antes del proceso de medición, que todos los valores de η tienen la misma probabilidad de ser el valor real a medir, entonces por el teorema de Bayes de la probabilidad inductiva [16] se concluye que, una vez realizadas las mediciones, la probabilidad de que el valor real de η esté en el intervalo [y∗, y∗+ dy∗] es εsp(y1 − y∗)p(y2 − y∗) . . . p(ys − y∗) dy∗ ∞∫ −∞ εsp(y1 − η)p(y2 − η) . . . p(ys − η) dη , y por consiguiente la hipótesis más probable, con respecto al verdadero valor de η, es aquel valor de y∗ que maximice la expresión p(y1 − y∗)p(y2 − y∗) . . . p(ys − y∗). Dicho y∗, simbolizado por η∗, satisface por ende la ecuación s∑ q=1 d dη ln [p(yq − η∗)] = 0. (2.17) Por otro lado, según el postulado de la media aritmética, se tiene que s∑ q=1 (yq − η∗) = 0. (2.16a) Condición suficiente para que exista compatibilidad entre (2.17) y (2.16a) es que la función de densidad p satisfaga la ecuación diferencial d dη ln [p(yq − η)] = c(yq − η), (2.18) donde c es una constante arbitraria. Efectuando el cambio de variable ∆ = (yq−η), se obtiene la solución de 20 (2.18) p(∆) = Ae− 1 2 c∆ 2 . Para determinar el valor particular de la constante de integración A, se con- sidera que p(∆) posee, por ser función de densidad de probabilidad, la pro- piedad ∞∫ −∞ p(∆) d∆ = 1, por lo tanto A ∞∫ −∞ e− 1 2 c∆ 2 d∆ = 1; y puesto que ∞∫ −∞ e−x 2 dx = √ π, entonces A = √ c 2π . Finalmente, introduciendo σ= √ 2/c se obtiene p(∆) = 1√ πσ e− ∆2 σ2 , (2.19) la cual es la función de densidad normal. 2.2.4. Deducción de Gauss del Criterio de los Mínimos Cuadrados Como paso previo a la deducción del criterio de los mínimos cuadrados, es conveniente introducir el concepto de ‘peso’ asociado a cada medida. En términos generales, el peso de una medida indica el grado de precisión in- volucrado en el proceso de medición del que resulta dicha medida. De esta manera, si una medida está menos propensa a error, mayor es su peso. En conformidad con lo arriba expresado, algunas posibles definiciones del peso w asociado a una cierta medición, que introduce un error ∆ con densidad de probabilidad p(∆), son: 21 - la densidad de probabilidad de que la medida esté libre de errores w = p(0), - la probabilidad de que el error caiga dentro de un cierto rango de tole- rancia |∆| ≤ ∆0 w = ∆0∫ −∆0 p(∆) d∆, - la inversa del tamaño del intervalo centrado en cero dentro del cual los errores caen con cierta probabilidad P0 P0 = 1/2w∫ −1/2w p(∆) d∆. Sin embargo, el peso se define formalmente como la inversa del cuadrado de la desviación estándar del error de medición w =  ∞∫ −∞ ∆2p(∆) d∆ −1 . (2.20) Esta definición es adecuada puesto que la desviación estándar, al ser la raiz cuadrada del error cuadrático medio, constituye una medida de la presencia del error, el cual está en relación inversa con la precisión y por tanto con el peso. Para el caso tratado, de densidad de error normal, se tiene a partir de (2.19) y (2.20) w = σ−2. (2.21) Dadas n respuestas con valores reales y∗1, y∗2, . . . , y∗n, la probabilidad de ob- tener, a partir de estas, sendas medidas y1, y2, . . . , yn es P = 1√ πσ1 e − (y1−y ∗ 1)2 σ2 1 1√ πσ2 e − (y2−y ∗ 2)2 σ2 2 . . . 1√ πσn e − (yn−y∗n)2 σ2 n , donde se considera que los valores de la desviación estándar del error intro- ducido en el proceso de medición de cada respuesta son, en general, distintos entre sí. 22 El criterio que define el método de los mínimos cuadrados establece que el ajuste óptimo de un modelo a un conjunto de medidas y1, y2, . . . , yn se consigue cuando se maximiza la probabilidad de obtener dichas medidas, adoptando los correspondientes valores predichos por el modelo ŷ1, ŷ2, . . . , ŷn como los valores reales de los que se obtienen las medidas, i.e., cuando el vector de parámetros de ajuste es tal que maximice la expresión P = 1√ πσ1 e − (y1−ŷ1)2 σ2 1 1√ πσ2 e − (y2−ŷ2)2 σ2 2 . . . 1√ πσn e − (yn−ŷn)2 σ2 n , o, en forma equivalente al tomar logaritmos y considerar (2.21), que minimice Φ = w1(y1 − ŷ1)2 + w2(y2 − ŷ2)2 + · · ·+ wn(yn − ŷn)2 = n∑ i=1 wi(yi − ŷi)2. (2.22) En este caso general, al método se le denomina mínimos cuadrados pondera- dos, sin embargo la mayor cantidad de aplicaciones del método se dan para el caso especial en que todos los procesos de medición poseen la misma precisión o peso wi=w (i=1, 2, . . . , n), de manera que (2.22) se convierte en Φ = n∑ i=1 (yi − ŷi)2, que no es otra cosa que (2.3). 2.2.5. Deducción Alterna de Laplace Posteriormente a la primera deducción de Gauss, tratada en las subsec- ciones precedentes, el método de los mínimos cuadrados fue derivado de una manera completamente distinta por Laplace [17], cuyo enfoque influenció subsecuentes investigaciones en el tema; entre ellas, un trabajo de Gauss [10] en el que expone su segunda deducción del método. El principio en el que se basa esta deducción se puede enunciar en los siguientes términos: Para las n respuestas predichas por el modelo mediante las expresiones 23 lineales ŷ1 = c11b1 + c21b2 + . . . + ck1bk ŷ2 = c12b1 + c22b2 + . . . + ck2bk ... ... ... . . . ... ŷn = c1nb1 + c2nb2 + . . . + cknbk,  se tienen respectivamente las estimaciones y1, y2, . . . , yn obtenidas mediante medición, donde n se supone mayor que el número de parámetros k. Para la cantidad µ1(c11b1 + c21b2 + · · ·+ ck1bk) + µ2(c12b1 + c22b2 + · · ·+ ck2bk) + . . . + µn(c1nb1 + c2nb2 + · · ·+ cknbk) (A) se tiene entonces la estimación µ1y1 + µ2y2 + · · ·+ µnyn. (B) Si los µi se escogen de forma tal que en la expresión (A) el coeficiente re- sultante de b1 es uno y los coeficientes resultantes de b2, b3, . . . , bk son cero, i.e., c11µ1 + c12µ2 + . . . + c1nµn = 1 c21µ1 + c22µ2 + . . . + c2nµn = 0 ... ... . . . ... ... ck1µ1 + ck2µ2 + . . . + cknµn = 0,  (2.23) entonces la expresión (B) es una estimación de b1. Puesto que existen infinitos juegos de µi que son solución de (2.23), el problema consiste en establecer condiciones adicionales para determinar el juego de valores µi que asegure que (B) sea la mejor estimación posible de b1. Esto se consigue considerando la noción de peso de una estimación y definiendo como mejor estimación a aquella que posea el mayor peso. Si se tiene que los pesos asociados a las estimaciones y1, y2, . . . , yn son respectivamente w1, w2, . . . , wn, entonces es posible demostrar [7] que el peso W asociado a la expresión lineal (B) satisface 1 W = µ 2 1 w1 + µ 2 2 w2 + · · ·+ µ 2 n wn . (2.24) 24 Para determinar el juego de valores µi que minimicen la inversa del peso (1/W ) y, simultáneamente, cumplan con el sistema de restricciones (2.23), se emplea el método de los multiplicadores de Lagrange (véase sección 4.2.1), de cuya aplicación se obtiene el sistema de ecuaciones µ1 + w1c11λ1 + w1c21λ2 + . . . + w1ck1λk = 0 µ2 + w2c12λ1 + w2c22λ2 + . . . + w2ck2λk = 0 ... ... ... . . . ... ... µn + wnc1nλ1 + wnc2nλ2 + . . . + wncknλk = 0,  (2.25) que junto al sistema de ecuaciones (2.23) determinan unívocamente el juego de valores µi con los que se calcula la mejor estimación (bmin)1 del parámetro b1 mediante la expresión (bmin)1 = µ1y1 + µ2y2 + · · ·+ µnyn. (2.26) Despejando los µi en (2.25) e introduciéndolos en (2.23) y (2.26) se obtienen las ecuaciones (bmin)1 + [wc1y]λ1 + [wc2y]λ2 + . . . + [wcky]λk = 0 1 + [wc1c1]λ1 + [wc1c2]λ2 + . . . + [wc1ck]λk = 0 [wc1c2]λ1 + [wc2c2]λ2 + . . . + [wc2ck]λk = 0 ... ... . . . ... ... [wc1ck]λ1 + [wc2ck]λ2 + . . . + [wckck]λk = 0,  a partir de las cuales, siguiendo la técnica descrita en la sección 2.2.1, se deduce ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ 1 [wc1y] [wc2y] . . . [wcky] 1/(bmin)1 [wc1c1] [wc1c2] . . . [wc1ck] 0 [wc1c2] [wc2c2] . . . [wc2ck] ... ... ... . . . ... 0 [wc1ck] [wc2ck] . . . [wckck] ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ = 0, 25 de donde finalmente se despeja el valor buscado (bmin)1 = ∣∣∣∣∣∣∣∣∣∣∣∣ [wc1y] [wc2y] . . . [wcky] [wc1c2] [wc2c2] . . . [wc2ck] ... ... . . . ... [wc1ck] [wc2ck] . . . [wckck] ∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣ [wc1c1] [wc1c2] . . . [wc1ck] [wc1c2] [wc2c2] . . . [wc2ck] ... ... . . . ... [wc1ck] [wc2ck] . . . [wckck] ∣∣∣∣∣∣∣∣∣∣∣∣ . (2.27) Para el caso especial en que w1 = w2 = · · · = wn, se puede apreciar que el valor resultante de (bmin)1 es idéntico al obtenido en (2.12), el cual es la solución de las ecuaciones normales ordinarias dictadas por el método de los mínimos cuadrados; por lo tanto, el enfoque alterno de Laplace conduce al establecimiento de dicho método. 2.3. Interpretación Geométrico Vectorial del Método según Kolmogorov La originalidad e importancia del trabajo de Kolmogorov [15], concernien- te al método de los mínimos cuadrados, radica no en la propuesta de algún nuevo procedimiento de deducción del método, sino en su interpretación en- marcada en el campo de la geometría vectorial n-dimensional. En esta sección se presentan los aspectos más relevantes del mencionado trabajo en el que mediante un ejemplo del caso más simple, correspondien- te a un modelo lineal, se muestra cómo la aplicación de métodos generales del álgebra lineal vectorial, tales como el concepto de ortogonalidad, permi- te obtener de manera mucho más transparente todos los resultados básicos del método, así como los instrumentos para evaluar la confiabilidad de las estimaciones involucradas. 26 2.3.1. Obtención de las Ecuaciones Normales Mediante Métodos Vectoriales Considérese un sistema cuyo modelo exacto f(x, b) es lineal en los pará- metros bj, j = 1, 2, . . . , k. Por definición, para dicho modelo existe un vector de parámetros exactos β, a priori desconocidos, con los cuales se calculan las respuestas del sistema (y∗i ) mediante la expresión lineal y∗i = k∑ j=1 cjiβj, i = 1, 2, . . . , n, (2.28) donde los coeficientes cji, al estar en función de las entradas xi, son conocidos. Asimismo, asúmase que dados los valores de y∗i y xi, los βj están deter- minados de manera única por (2.28). Esto implica que el rango de la matriz [cji]n×k no es menor que k, por lo tanto se debe cumplir que n ≥ k. Sin embargo, la única forma de conocer las respuestas del sistema (y∗i ) es mediante un proceso de medición que inevitablemente introduce error. Así en lugar de obtener los verdaderos valores (y∗i ), se obtienen experimentalmente los valores yi = y∗i + ∆i, i = 1, 2, . . . , n, (2.29) lo que elimina la posibilidad de obtener los βj a partir de (2.28), ya que los errores ∆i se introducen en (2.28) como incógnitas adicionales. Ante la imposibilidad de conocer los parámetros exactos βj, se utilizan parámetros con valores alternos arbitrarios bj, por lo que el modelo ya no entrega las respuestas reales (y∗i ), sino las respuestas predichas ŷi = k∑ j=1 cjibj, i = 1, 2, . . . , n, (2.30) las cuales, en general, son distintas a las primeras. Como ya se indicó, el método de los mínimos cuadrados define los residuos ri = yi − ŷi, i = 1, 2, . . . , n, (2.31) y establece que los parámetros bj que mejor substituyen a los parámetros exactos βj en lograr que el modelo caracterice al sistema, son aquellos que 27 minimizan la norma del residuo [rr] = n∑ i=1 riri, (2.32) donde los corchetes denotan el producto escalar de los vectores contenidos. Es posible demostrar, aplicando cálculo elemental, que los bj que mini- mizan (2.32) están determinados unívocamente por el sistema de ecuaciones normales k∑ j=1 [chcj ]bj = [chy], h = 1, 2, . . . , k. (2.33) Sin embargo, el uso del cálculo no es la única manera de establecer la equi- valencia entre cumplir la condición de minimizar (2.32) y ser solución de las ecuaciones (2.33); también es posible obtener dicho resultado mediante el uso de métodos vectoriales. Con ese fin, considérense y∗i , cji, yi, ∆i, ŷi, ri, (i = 1, 2, . . . , n), como las componentes de los vectores n-dimensionales y∗, cj , y, ∆, ŷ, r, con lo cual (2.28) - (2.31) se reescriben en forma vectorial y∗ = k∑ j=1 βjcj , (2.28a) y = y∗ + ∆, (2.29a) ŷ = k∑ j=1 bjcj , (2.30a) r = y − k∑ j=1 bjcj . (2.31a) Si se denota con V al subespacio lineal generado dentro del espacio vectorial Rn por el conjunto de vectores {cj}, se desprende entonces de (2.30a) que el vector ŷ está contenido en V (ŷ ∈ V ). En lenguaje vectorial, cumplir la condición de minimizar [rr] equivale a que ŷ sea la proyección ortogonal de y sobre V y r sea el vector complemen- tario ortogonal a V , con lo cual [rcj ] = 0, j = 1, 2, . . . , k. (2.34) La figura 2.2 ilustra el trazado, en conformidad con el método, de los 28 ξ η ζ c1 c2 *y y ŷ Δ r V Figura 2.2: Representación vectorial del método para n = 3 y k = 2. vectores mencionados para un caso arbitrario en que el número de muestras es n = 3 y el número de parámetros es k = 2. Nótese que en este caso el subespacio lineal V corresponde a un plano en R3 que pasa por el origen. Reemplazando en (2.34) la expresión para r dada por (2.31a) se obtienen las ecuaciones k∑ j=1 [chcj ]bj = [chy], h = 1, 2, . . . , k, (2.33) correspondientes al sistema de ecuaciones normales formulado en (2.33). Este sistema se puede reescribir en forma matricial (véase (2.9), en la página 15), A b = g, (2.35) donde el determinante del sistema, |A| = |[chcj ]|, es el determinante de Gram de los vectores c1, c2, . . . , ck. Bajo el supuesto, enunciado líneas arriba, que la matriz [cji]n×k posee ran- go k, se deduce que los vectores c1, c2, . . . , ck son linealmente independientes y, en consecuencia, su determinante de Gram es distinto de cero (|A| 6= 0). Por lo tanto, al ser A invertible, (2.33) determina unívocamente los paráme- tros bj. 29 2.3.2. Evaluación de la Confiabilidad de los Parámetros Estimados Hasta este punto se ha visto cómo el método estima los parámetros bj a partir de las respuestas medidas (sujetas a error) ante la inaccesibilidad de las respuestas reales necesarias para calcular los parámetros verdaderos βj. Es importante entonces establecer un criterio para evaluar el grado con que se desvían los parámetros estimados respecto a los parámetros reales; lo cual será un indicativo de la confiabilidad de los resultados. Esto se consigue analizando la esperanza y la varianza de los parámetros bj. La esperanza E{bj} indica el valor al que converge la media de los parámetros, por lo que un valor distinto de βj (de conocerse βj) señalaría la presencia de error sistemático en los resultados. La varianza V {bj} por otro lado, al ser el valor cuadrático medio de la desviación de bj respecto a su media, determina la dispersión o, visto a la inversa, la exactitud de los resultados. Antes de proceder al cálculo de E{bj} y V {bj} es necesario asumir las siguientes condiciones respecto a la distribución de probabilidad de los errores de medición ∆i: C1) ∆i y ∆l, con i 6= l, son variables aleatorias independientes. C2) E{∆i} = 0. C3) V {∆i} = E{∆ 2 i } = s2, donde s2 es finito e independiente de i. C4) E{∆i∆l} = 0 para i 6= l. Como punto de partida, defínase el sistema de vectores u1, u2, . . . , uk, en el subespacio lineal V , de forma que sea δk-ortogonal al sistema c1, c2, . . . , ck, i.e., que cumpla [cjuh] = I (k) jh , (2.36) donde I(k) es la matriz identidad de k × k. Puesto que cada uno de los sistemas {c1, c2, . . . , ck} y {u1, u2, . . . , uk} es una base generadora de V , se tiene que uj = k∑ λ=1 Qjλcλ, cj = k∑ λ=1 Ajλuλ. 30 Tomando el producto escalar de la primera igualdad con uh y el producto escalar de la segunda con ch, y considerando (2.36), se tiene Qjh = [ujuh], Ajh = [cjch], (2.37) con lo cual, uj = k∑ h=1 [ujuh]ch, (2.38) cj = k∑ h=1 [cjch]uh. (2.39) De estas dos últimas ecuaciones y de (2.37), se desprende que la matriz Q es la inversa de la matriz A. Tomando ahora el producto escalar de (2.28a) y uj , se obtiene βj = [y∗uj ], (2.40) y en forma similar, (2.30a) implica que bj = [ŷ uj ]. (2.41) Por otro lado, puesto que r es un vector ortogonal a V , se cumple que [r uj ] = 0, con lo cual, de (2.41) y (2.31a), se determina que bj = [y uj ]. (2.42) De esta última ecuación, junto con (2.29a) y (2.40), se deduce que bj = βj + [∆uj ], (2.43) y por ende, considerando C2), se llega al resultado E{bj} = βj, (2.44) 31 que establece que la esperanza de cada uno de los parámetros estimados (aproximados), es igual al correspondiente parámetro verdadero. Se dedu- ce entonces que los valores de los parámetros estimados no contienen error sistemático sí y solo sí se cumple la condición C2). Obtenido el valor de la esperanza de los parámetros bj, resta determinar la varianza de los mismos, que indicará la precisión con la que cada parámetro bj aproxima en promedio a su respectivo βj. Para ello se calcula primero la expresión más general correspondiente a E{(bj − βj)(bh − βh)}, la cual, empleando (2.43), C3) y C4), da como resultado E{(bj − βj)(bh − βh)} = E{[∆uj ][∆uh]} = n∑ i=1 n∑ l=1 E{∆i∆l}ujiuhl = s2 n∑ i=1 ujiuhi = [ujuh]s2, que, en vista de (2.37), se puede expresar como E{(bj − βj)(bh − βh)} = Qjhs 2. (2.45) Finalmente, para el caso particular h = j, se obtiene el valor buscado de la varianza de bj V {bj} = E{(bj − βj)2} = Qjjs 2. (2.46) Nótese que la varianza de los bj es proporcional a la varianza de los errores de medición (V {∆i} = s2). Por esta razón es necesario, en caso de no conocerse a priori el valor de s, establecer un valor σ que sirva como aproximación del mismo. Con ese fin, se define el vector ∆∗ = ŷ − y∗, (2.47) el cual, por (2.29a) y (2.31a), se puede expresar como ∆ = ∆∗ + r. (2.48) Puesto que ∆∗ está contenido en V junto con ŷ e y∗, y r es ortogonal a V , se concluye que esta última ecuación representa la descomposición de ∆ en su proyección ortogonal sobre V y el complemento ortogonal a dicha proyección. 32 Si se define la base ortonormal del espacio vectorial n-dimensional eτ = (eτ1, eτ2, . . . , eτn), τ = 1, 2, . . . , n, [eτeτ ′ ] = I (n) ττ ′ , de forma tal que los k primeros vectores estén contenidos en V y los n − k restantes sean ortogonales a V ; entonces, considerando ∆̃τ = [∆eτ ], (2.49) se obtiene ∆ = n∑ τ=1 ∆̃τeτ , (2.50) ∆∗ = k∑ τ=1 ∆̃τeτ , (2.51) r = n∑ τ=k+1 ∆̃τeτ . (2.52) A partir de la propiedad de ortonormalidad de los vectores eτ de la base, se deduce que [rr] = n∑ τ=k+1 ∆̃τ∆̃τ , (2.53) y puesto que por (2.49) y C3) E{∆̃τ∆̃τ} = E{ n∑ i=1 ∆ieτi n∑ i′=1 ∆′ieτi′} = n∑ i=1 n∑ i′=1 E{∆i∆′i}eτieτi′ = s2 n∑ i=1 eτieτi = s2[eτ ] = s2, (2.54) se concluye finalmente que E{[rr]} = n∑ τ=k+1 E{∆̃τ∆̃τ} = n∑ τ=k+1 s2 = (n− k) s2. (2.55) Definiendo σ = √ [rr]/(n− k), (2.56) 33 se cumple que E{σ2} = s2, (2.57) y además es posible demostrar, siguiendo un procedimiento similar al em- pleado, que V {σ2} = E{(σ2 − s2)} = 2s2/(n− k), (2.58) lo que significa que σ2 tiene un valor muy cercano a s2, con probabilidad muy cercana a uno, para valores de (n− k) grandes. Se concluye así que conforme (n−k) crece, la aproximación de s2 mediante σ2 mejora, y para estos casos se hace válido aproximar la varianza de los parámetros bj dada en (2.46) mediante la expresión V {bj} ≈ Qjjσ 2. (2.59) 2.4. Aplicabilidad del Criterio de los Mínimos Cuadrados cuando el Modelo es Aproximado y no Existe Error de Medición Dado que hasta este punto se asumió conocido el modelo exacto de un sistema, el estudio del criterio de los mínimos cuadrados ha estado siempre confinado al caso en que el error aleatorio de medición es la única fuente de discrepancia entre los datos observados de un sistema y aquellos predichos por su modelo. Sin embargo, es virtualmente imposible conocer todos los detalles internos de un sistema, lo que impide elaborar un modelo exacto del mismo. En la práctica solo se cuenta con un modelo aproximado incapaz de reproducir la verdadera respuesta del sistema que se obtiene a partir de un proceso de medición libre de error. Para un modelo aproximado no existe un juego de parámetros de ajuste exactos o verdaderos, ya que la discrepancia entre el sistema y el modelo es ahora de naturaleza sistemática. En vista que los argumentos que justifican el uso del criterio de los mínimos cuadrados tienen su fundamento en el error aleatorio de medición, cabe cuestionarse qué tan correcto es aplicar dicho 34 criterio en el ajuste de modelos aproximados donde el error es sistemático y se asume la ausencia de error de medición. La presente sección busca discutir tal cuestionamiento. Para ello se pro- pone el diseño de una función Ψ a partir de ciertas propiedades elementales adecuadas para cuantificar la discrepancia entre el modelo aproximado y el sistema. La aplicación del criterio de los mínimos cuadrados en este nuevo contexto será entonces más o menos correcto según el grado de concordancia con dicha función Ψ. 2.4.1. Expresión General para la Función de Discrepancia Como punto de partida hay que establecer que el error o residuo entre el valor real y∗i (coincidente en este caso con el valor medido yi) y el valor ŷi predicho por el modelo, debe evaluarse según una escala que esté en función a su diferencia. Este residuo puede expresarse matemáticamente en forma general como ri = r(yi, ŷi), (2.60) donde valores positivos o negativos significan error por exceso o defecto res- pectivamente. Un caso particular de (2.60), conocido como error relativo, se presenta cuando la medida del error es proporcional a la desviación pero normalizada respecto al valor medido. La forma del residuo en este caso es ri = yi − ŷi yi . (2.61) Otro caso particular se da cuando el error se define igual a la desviación. En este caso el residuo toma su forma más sencilla y se le conoce como error absoluto, ri = yi − ŷi. (2.62) En este trabajo se adopta una definición del residuo independiente del signo de la desviación: ri = |yi − ŷi|. (2.63) Una vez definido el error o residuo individual al representar cierto valor yi 35 por medio de otro ŷi, resta definir el error o residuo global entre las n-tuplas (y1, y2, . . . , yn) y (ŷ1, ŷ2, . . . , ŷn). La función de discrepancia es justamente la que cuantifica dicho residuo global mediante una expresión matemática (a determinar) en función de los residuos individuales. Por consiguiente, cualquier intento por definir una estrategia para calcular el residuo global debe expresarse matemáticamente en función de los residuos individuales como una función de discrepancia Ψ = Ψ(r1, r2, . . . , rn). (2.64) El diseño de Ψ se basa en las siguientes propiedades o requerimientos elementales: P1) El residuo global no depende del orden de los residuos individuales en Ψ(r1, r2, . . . , rn) y por convención se adoptará un orden ascendente de los mismos, i.e., ri+1 ≥ ri. P2) La expresión general del residuo global contiene un número infinito de argumentos. El caso particular para n argumentos toma la forma Ψ(r1, r2, . . . , rn) = Ψ(r1, r2, . . . , rn, 0, 0, . . . ). Nótese que los residuos individuales con valor cero no contribuyen al residuo global. P3) El residuo global para el caso particular de n = 1 es igual al residuo individual, i.e., Ψ(r1) = r1. P4) El residuo global es creciente respecto a cada uno de sus argumentos, i.e., ∂Ψ/∂ri > 0 en todos los casos. P5) El residuo global posee la misma dimensión física que los residuos in- dividuales. La propiedad P5) significa que Ψ y los ri se expresan en la misma unidad de medida. Por lo tanto, si se utiliza una unidad de medida contenida k veces en la unidad original, se cumple Ψ(kr1, kr2, . . . , krn) = kΨ(r1, r2, . . . , rn), k > 0. (2.65) Esto revela que Ψ es una función homogénea de primer grado. Introduciendo k = r−1 n en (2.65) se obtiene la forma general que poseen todas las funciones 36 homogéneas de primer grado Ψ = rn f( r1 rn , r2 rn , . . . , rn−1 rn ), (2.66) donde f es una función arbitraria cualquiera. Derivando Ψ respecto a cada ri y eliminando f en (2.66) resulta la ecua- ción diferencial denominada relación de homogeneidad de Euler r1 ∂Ψ ∂r1 + r2 ∂Ψ ∂r2 + · · ·+ rn ∂Ψ ∂rn = Ψ, (2.67) cuyo conjunto solución es la totalidad de funciones homogéneas de primer grado [4]. La solución general de esta ecuación se puede construir mediante la com- binación lineal de todas las soluciones particulares obtenidas por el método de la separación de variables [26]. Una función con variables separadas Ψ = R1(r1)R2(r2) · · ·Rn(rn), (2.68) es solución de (2.67) si y solo si r1 R1 ∂R1 ∂r1 + r2 R2 ∂R2 ∂r2 + · · ·+ rn Rn ∂Rn ∂rn = 1. (2.69) Como cada sumando en el miembro izquierdo de esta ecuación depende de una variable distinta, estos deberán ser constantes αi de valor arbitrario que cumplan con la restricción α1 + α2 + · · ·+ αn = 1. Por lo tanto, (2.69) es equivalente a n ecuaciones diferenciales ordinarias de primer grado ri ∂Ri ∂ri − αiRi = 0, i = 1, 2, . . . , n, cuyas soluciones, donde los ci representan constantes de integración, son Ri(ri) = cir αi i , i = 1, 2, . . . , n. 37 Considerando (2.66) y (2.68) se obtiene la solución particular Ψ = c rn ( r1 rn )α1 ( r2 rn )α2 · · · ( rn−1 rn )αn−1 , (2.70) donde c es una constante arbitraria. Finalmente, la solución general de (2.67) se construye combinando lineal- mente las soluciones particulares (2.70) generadas al considerar todas las posibles (n− 1)-tuplas (α1, α2, . . . , αn−1). La solución general es entonces la integral Ψ = rn ∞∫ −∞ ∞∫ −∞ . . . ∞∫ −∞ c(α1, α2, . . . , αn−1) ( r1 rn )α1 ( r2 rn )α2 · · · · · · ( rn−1 rn )αn−1 dα1 dα2 . . . dαn−1 (2.71) donde c(α1, α2, . . . , αn−1) es una función arbitraria compatible con las pro- piedades P1) – P5). Nótese además que (2.71) es una forma distinta pero equivalente de expresar (2.66). 2.4.2. Ejemplos Particulares de la Función de Discrepancia Para identificar dentro de la expresión general (2.71) qué función de dis- crepancia es la apropiada para cierto problema de ajuste de modelos, es necesario derivar la forma particular de c(α1, α2, . . . , αn−1) a partir a los requerimientos específicos de dicho problema. Sin embargo, con el fin de mostrar algunos ejemplos de funciones de dis- crepancia particulares se seleccionarán arbitrariamente algunos c especiales. Uno de los casos más sencillos se da cuando c(α1, α2, . . . , αn−1) = δ(α1)δ(α2) . . . δ(αn−1) + δ(α1 − 1)δ(α2) . . . δ(αn−1) + δ(α1)δ(α2 − 1) . . . δ(αn−1) + · · · · · · + δ(α1)δ(α2) . . . δ(αn−1 − 1), (2.72) 38 dando como resultado la función Ψ = r1 + r2 + . . . + rn, (2.73) correspondiente al criterio denominado de los mínimos errores absolutos. Un caso más general y complejo se presenta cuando c(α1, α2, . . . , αn−1) = ∞∑ m=0 1/β m ) m∑ m1=0 m−m1∑ m2=0 · · · m−m1−...−mn−3∑ mn−2=0 m! m1!m2! · · ·mn−1! δ(α1−βm1)δ(α2−βm2) . . . δ(αn−1−βmn−1), (2.74) donde mn−1 = m− (m1 +m2 + · · ·+mn−2). Considerando la serie binomial y el teorema del multinomio es posible demostrar que cuando se cumple la condición de convergencia rβ1 + rβ2 + . . . + rβn < 2 rβn, (2.75) se obtiene la siguiente forma para la función de discrepancia Ψ = (rβ1 + rβ2 + . . . + rβn)1/β, (2.76) la cual coincide con la β-norma de un espacio vectorial de dimensión finita n. Finalmente, nótese que el criterio de los mínimos cuadrados corresponde al caso especial β = 2. Capítulo 3 Algoritmos Precursores: Métodos del Máximo Descenso y de Gauss-Newton El objetivo del campo de la optimización consiste en el desarrollo de algo- ritmos para evaluar el punto bmin que minimice (o maximice) cierta función suave Φ : Rk → R en un subconjunto dado A de su dominio, i.e., Φ(bmin) = mı́n b∈A⊂Rk Φ(b). Un grupo especial de estos algoritmos lo constituyen aquellos diseñados para minimizar la función particular Φ(b) = n∑ i=1 r 2 i (b) = ‖r‖2, estudiada en el capítulo anterior. En el presente capítulo se analizan dos de estos algoritmos, conocidos co- mo: ‘del máximo descenso’ y ‘Gauss-Newton’ respectivamente. Si bien ambos son de naturaleza elemental y fueron concebidos en una etapa inicial, su im- portancia radica en que han sido precursores de algoritmos más sofisticados y de mejor desempeño, tales como el de ‘Levenberg-Marquardt’ tratado en el siguiente capítulo. 39 40 3.1. Introducción Todos los algoritmos de minimización requieren como punto de partida un tanteo inicial del mínimo, denotado por b(0), a partir del cual calculan, mediante un proceso iterativo, una secuencia (b(1), b(2), . . . , b(s), . . . ) de pun- tos tentativos que mejoran progresivamente el tanteo inicial y convergen en el mínimo bmin. Para calcular un nuevo punto tentativo b(s+1), estos algoritmos utilizan información sobre Φ, y sus derivadas, en el punto tentativo actual b(s) y, eventualmente, en los anteriores (. . . , b(s−2), b(s−1)). La condición fundamental a cumplir en este caso es que Φ(b(s+1)) < Φ(b(s)), aunque existen excepciones como en el caso de los llamados algoritmos no monótonos, en los que la condición es que la función decrezca luego de un cierto número p de iteraciones, Φ(b(s+p)) < Φ(b(s)). El proceso de iteración finaliza cuando se obtiene un punto cuya distan- cia al punto anterior es menor a cierto umbral preestablecido, o cuando la diferencia entre los valores de Φ en los puntos anterior y actual es igualmente menor a cierto umbral. La rapidez con la que el algoritmo arriba a dicho punto final depende no solo de la eficiencia del mismo, reflejada en su razón de convergencia como medida de su velocidad en alcanzar el mínimo, sino también de qué tan cerca o lejos del mínimo esté el tanteo inicial. Existen dos estrategias fundamentales para construir procedimientos ite- rativos que busquen el mínimo: de búsqueda en línea (line search) y de región de confianza (trust region). Cronológicamente, los algoritmos de búsqueda en línea, como los del máximo descenso y de Gauss-Newton, se desarrollaron an- tes que los algoritmos de región de confianza, los cuales tienen su génesis en el algoritmo de Levenberg-Marquardt. En la estrategia de búsqueda en línea se elige en cada iteración un vector dirección unitario δ̂(s), con el que se construye un vector de paso que conduce al nuevo punto tentativo b(s+1) a partir del actual b(s). La longitud αs del paso se obtiene resolviendo el problema de minimización unidimensional Φ(b(s) + αs δ̂ (s)) = mı́n α> 0 Φ(b(s) + α δ̂ (s)). (3.1) En la práctica, sin embargo, los algoritmos de búsqueda en línea no optan por 41 resolver este problema, por ser este un proceso de alto costo computacional, sino que generan un número limitado de longitudes de paso de prueba hasta conseguir una que se aproxime suficientemente a la longitud de paso óptimo que resulta de (3.1). En la estrategia de región de confianza se extrae en cada iteración in- formación de Φ en el punto tentativo actual b(s) para crear, alrededor del mismo, un modelo ms aproximado que reemplace a Φ en la búsqueda del mínimo. La vecindad Vs en la cual ms constituye una buena aproximación de Φ se conoce como región de confianza. El paso δ(s) que conduce al nuevo punto tentativo b(s+1) se calcula resolviendo aproximadamente el problema de minimización ms(b(s) + δ(s)) = mı́n δ ∈Vs ms(b(s) + δ). (3.2) Si la iteración no produce una disminución aceptable de Φ, es necesario en- coger la región de confianza Vs y volver a resolver (3.2). Por lo general, Vs es una bola definida por ‖δ‖ ≤ ∆(s), donde al escalar ∆(s) se le denomina radio de la región de confianza. En la mayoría de los casos, el modelo ms lo constituye la función cuadrá- tica correspondiente a la serie de Taylor truncada de Φ alrededor de b(s) ms(b(s) + δ) = Φ(s) + (∇TΦ(s)) δ + 1 2 δ T B(s) δ, donde para abreviar: Φ(s) = Φ(b(s)), ∇Φ(s) =∇Φ(b(s)), y B(s) es el Hessiano en b(s) B(s) =∇(∇TΦ(s)) =  ∂2Φ(b(s)) ∂b 2 1 ∂2Φ(b(s)) ∂b1∂b2 · · · ∂2Φ(b(s)) ∂b1∂bk ∂2Φ(b(s)) ∂b1∂b2 ∂2Φ(b(s)) ∂b 2 2 · · · ∂2Φ(b(s)) ∂b2∂bk ... ... . . . ... ∂2Φ(b(s)) ∂b1∂bk ∂2Φ(b(s)) ∂b2∂bk · · · ∂2Φ(b(s)) ∂b 2 k  , o alguna aproximación del mismo. Antes de pasar a ver en detalle los dos algoritmos de minimización men- cionados, es conveniente resaltar la importancia de la expansión en series de 42 Taylor como herramienta matemática fundamental empleada en la minimi- zación de funciones continuamente diferenciables. Para el caso particular de interés, correspondiente a funciones de más de una variable independiente a ser expandidas hasta el término de segundo grado, el teorema de Taylor toma la forma enunciada a continuación. Teorema de Taylor. Sea Φ(b) : Rk → R una función continua y doble- mente diferenciable en todos los puntos b de la bola ‖ b − b0 ‖ < δ0. Si se tiene además que δ ∈ Rk con ‖δ‖ = δ < δ0, entonces existe un t ∈ [0, 1] para el cual se cumple Φ(b0 + δ) = Φ(b0) + (∇TΦ(b0)) δ + 1 2 δ T ∇(∇TΦ(b0 + t δ)) δ. (3.3) La demostración de esta versión multivariable del teorema se puede lle- var a cabo generalizando el procedimiento seguido en la demostración del caso univariable que aparece en diversos textos de cálculo. El lector intere- sado puede consultar en particular el artículo de Pringsheim [27], que ofrece también un recuento histórico del desarrollo del tema. 3.2. Método del Máximo Descenso Una manera sencilla de visualizar el método del máximo descenso, cono- cido también como ‘steepest descent’, es mediante un símil con la trayectoria que sigue una esfera que rueda sobre una superficie irregular por acción única de la gravedad hasta que alcanza un valle, o mínimo relativo, en el que queda encajonada y detenida. En la presente sección se analizan en detalle las tres etapas que componen este método: la elección de la dirección de minimización, el reescalamiento de la función a minimizar, y el cálculo de la longitud del paso en cada iteración. La sección concluye con un análisis de la convergencia del método. 3.2.1. Elección de la Dirección de Minimización En cada nueva iteración se busca elegir la dirección δ̂ en la cual la función Φ a minimizar decrezca de forma más abrupta a partir del punto tentativo actual, i.e., la dirección de máximo descenso. Esta dirección, evidentemente, 43 está dada por el negativo de la gradiente en dicho punto, que apunta hacia donde la derivada o razón de cambio de Φ adquiere su menor valor. En particular, para el caso de la (s+1)-ava iteración, la búsqueda del nuevo punto tentativo b(s+1) se lleva a cabo a partir del actual b(s), siguiendo la dirección δ̂ (s) = − ∇Φ(s) ‖∇Φ(s)‖ . Si bien lo anteriormente expresado fundamenta la elección del negativo de la gradiente como dirección de minimización a ser usada por el método, es conveniente analizar, con cierto detalle, las limitaciones de tal elección. Para tal fin, y por simplicidad, se considera el caso de Φ definida en R2. Definiendo una familia de circunferencias concéntricas a cierto punto b en el dominio de Φ, e introduciendo luego coordenadas polares con origen en b, se tiene que el valor que toma Φ en los puntos rêθ = r(cos θ, sen θ) que conforman la circunferencia de radio r, es función únicamente de θ y, de acuerdo a lo establecido por el teorema de Taylor en (3.3), está dado por f(θ) = Φ(b+rêθ) = Φ(b)+(∇TΦ(b)) êθ r+ [1 2 ê T θ ∇(∇TΦ(b+ t rêθ)) êθ ] r2. (3.4) El valor θr correspondiente al punto de la circunferencia en el que Φ toma el menor valor, se obtiene como solución de f ′(θ) = 0, siempre y cuando se cumpla la condición f ′′(θr) > 0. Por lo tanto, igualando a cero la derivada de (3.4) se llega a la expresión que define a θr sin(θr − θ0) = g(r, θr) ‖∇Φ(b)‖ r, (3.5) donde se ha considerado que ∇Φ(b) = ‖∇Φ(b)‖(cos θ0, sin θ0), con la con- dición adicional que ‖∇Φ(b)‖ 6= 0, y a g(r, θr) como la derivada parcial respecto a θ de la expresión entre corchetes en (3.4). Además, por la condición de doble diferenciabilidad impuesta a Φ, se puede afirmar que tanto ‖∇Φ(b)‖ como g(r, θr) poseen valores finitos. Esto trae como consecuencia que al aplicar en (3.5) el límite cuando r → 0, el lado derecho se haga cero, con lo cual se tiene ĺım r→0 θr = θ0 + π. (3.6) 44 b θ = 0θo Φ = Φ(b) C Δ Φ(b) r = ½ a r = a r = 2a - (a,θa) mínimos de en cada Φ circunferencia Figura 3.1: Curva C que une los puntos mínimos de Φ en cada circunferencia. Este resultado indica que el mínimo sobre la circunferencia límite, cuyo radio tiende a cero, se encuentra en la dirección del negativo de la gradiente. En la figura 3.1 se muestra la curva C formada por los puntos que mini- mizan Φ en cada circunferencia, la cual es tangente a −∇Φ(b) en el punto b. En ella se puede también observar que conforme las circunferencias se ale- jan de b, los mínimos sobre estas se apartan en general de la dirección dada por −∇Φ(b). Sin embargo, es posible advertir que en una pequeña vecindad de b, delimitada por una circunferencia de radio menor a cierto umbral, es posible emplear la dirección del negativo de la gradiente como dirección de máximo descenso con un error de aproximación despreciable. En el caso lími- te r → 0 el negativo de la gradiente es efectivamente, como ya se demostró, la dirección de máximo descenso. Como comentario final es oportuno señalar que al recorrer la curva C partiendo de b, el valor de Φ decrece hasta alcanzar un punto, de existir uno, en el que la tendencia decreciente se revierte y Φ empieza a crecer. Es posible demostrar que este punto de cambio es precisamente el punto mínimo bmin requerido. Otra observación, no evidente a simple vista, es que la curva C es en general distinta a la curva trazada desde b siguiendo progresivamente la dirección del negativo de la gradiente en cada uno de sus puntos, i.e., a la curva que parte de b y es ortogonal a las curvas o superficies de valor constante de Φ, denominadas isolíneas o isosuperficies dependiendo de si el dominio de Φ está en R2 o R3 respectivamente. Cabe mencionar que la 45 gradiente en cada punto de C está en la dirección del radio que va del centro b (de la familia de circunferencias) a dicho punto. 3.2.2. Reescalamiento de la Función a Minimizar Uno de los factores que influye en el desempeño global del método del máximo descenso lo constituye la topología de la función Φ a minimizar, que puede en casos extremos degradar la eficiencia del método hasta el punto de hacerlo improductivo. Por este motivo es conveniente previamente ‘mejorar’ la topología de Φ mediante un proceso de reescalamiento. Para ilustrar esta influencia, considérese una función Φ, definida por sim- plicidad en R2, con isolíneas que formen una familia de circunferencias con- céntricas a su punto mínimo, e.g., Φ = x2 + y2. Esta topología de Φ permite alcanzar el mínimo desde cualquier punto siguiendo la recta en la dirección del negativo de su gradiente. Desde el punto de vista del método, se dice que esta topología es perfecta, pues se obtiene la máxima eficiencia al alcanzar el mínimo en un solo paso. Generalizando, es hasta cierto punto válido afirmar que mientras más difiera una topología del caso perfecto, más pasos serán necesarios para alcanzar el mínimo de su respectiva función. En particular, esta afirmación es correcta en la vecindad de un punto mínimo, donde la topología de cualquier función es aproximadamente la de su serie de Taylor truncada al segundo orden, i.e., una familia de elipses concéntricas al mínimo, las cuales se alejan más del caso perfecto, compuesto por circunferencias, conforme la razón entre sus ejes mayor y menor aumente. La figura 3.2 muestra las isolíneas elípticas de cierta función Φ en una región cercana al mínimo bmin. Para esta topología se cumple que desde cualquier punto, excepto aquellos sobre los ejes de las elipses, la dirección del negativo de la gradiente no conduce directamente a bmin. Para analizar la relación entre la forma de las elipses y el número de pasos requeridos para alcanzar bmin, es necesario tener en cuenta la forma en que el método construye cada paso. En la versión convencional del método, cada paso parte de un punto y recorre la dirección del negativo de su gradiente hasta alcanzar un mínimo. En este mínimo, que es el punto de partida del siguiente paso, la isolínea y la recta recorrida son tangentes. Por consiguiente, los sucesivos pasos que lleva a cabo el método son ortogonales entre sí. 46 eje menor eje mayor bmin x y b(0) b(1) isol neas de , aproximadamente í Φ e ípticas alrededor de bminl Figura 3.2: Pasos que lleva a cabo el método en una vecindad cercana del mínimo bmin. Volviendo al caso mostrado en la figura, donde las elipses son alargadas y el primer paso, que va de b(0) a b(1), forma un ángulo cercano a 45◦ con el eje mayor, se observa que el método alcanza el mínimo bmin luego de varios pequeños pasos decrecientes en zigzag. Nótese que si estas isolíneas fuesen aun más alargadas y angostas, i.e., más alejadas de la topología perfecta, los pasos en general serían más cortos y se necesitaría un mayor número de ellos para alcanzar el mínimo. Tres de las técnicas de reescalamiento más útiles para superar las dificul- tades descritas se reseñan a continuación. Una primera técnica, bastante efectiva cuando el punto tentativo actual b(s) ha alcanzado suficiente cercanía al mínimo bmin, consiste en reemplazar Φ por su serie de Taylor truncada de segundo grado alrededor de b(s), denotada por Φ∗, para luego elegir, convenientemente, como dirección de minimización a aquella que conduce al centro de las isolíneas elípticas de Φ∗. La figura 3.3 muestra precisamente el caso en que bmin se encuentra en la región, alrededor de b(s), donde Φ∗ constituye una buena aproximación de Φ. Por consiguiente, la separación entre bmin y el mínimo de Φ∗, que coincide con el centro b0 = (x0, y0) de sus isolíneas elípticas, es razonablemente pequeña y tiende, además, a disminuir conforme b(s) se acerca más a bmin. En consecuencia, cuando b(s) está suficientemente cerca a bmin, dar un paso en la dirección δ̂(s), que conecta b(s) con b0, representa una mejor op- ción para el objetivo de alcanzar bmin, que darlo en la dirección del máximo descenso −∇Φ(s), tal como se aprecia en la figura. 47 b (s) Δ Φ(s)- Λ δ (s) bo=(xo,yo)bmin ≈ isolíneas de Φ* región de validez Φ ≈ Φ* x y Φ = Φ (s)* Figura 3.3: Dirección de paso modificada δ̂ (s) , que apunta a b0 ≈ bmin y reemplaza la dirección original −∇Φ(s). Como se mencionó, el punto b0 es el centro de las isolíneas elípticas de la función Φ∗(b) = Φ(s) +∇Φ(s)(b− b(s)) + 1 2(b− b(s))∇(∇TΦ(s))(b− b(s)), (3.7) que es la versión aproximada de Φ alrededor del punto tentativo actual b(s). De esta familia de elipses, aquella que pasa por b(s) está dada por la ecuación ∇Φ(s)(b− b(s)) + 1 2(b− b(s))∇(∇TΦ(s))(b− b(s)) = 0, (3.8a) que en coordenadas cartesianas, b = (x, y), toma la forma 1 2Ax 2 + 1 2By 2 + Cxy +Dx+ Ey + F = 0, (3.8b) donde A = ∂2Φ(s) ∂x2 , B = ∂2Φ(s) ∂y2 , C = ∂2Φ(s) ∂x∂y , D = ∂Φ(s) ∂x − (A,C) · b(s) y E = ∂Φ(s) ∂y − (C,B) · b(s). El centro de esta elipse, que lo es igualmente de toda la familia de elipses, es el único punto respecto al cual (3.8) es simétrica. Resolviendo la ecuación 48 correspondiente a esta condición se obtienen las coordenadas del centro x0 = CE −BD AB (1− C2 AB ) , y0 = CD − AE AB (1− C2 AB ) , (3.9) cuyos valores están bien definidos, ya que (1 − C2 AB ) = 0 se da únicamente cuando (3.8) representa la ecuación de una recta o parábola y no, como en este caso, la de una elipse, para la cual se cumple siempre que (1− C2 AB ) > 0. Conocidas las coordenadas de b0, es posible definir un vector δ(s) arbitra- rio en la dirección que va de b(s) a b0, el cual por simplicidad se elige δ(s) = (1− C2 AB )(b0 − b(s)), (3.10) con respectivas componentes δ(s) x = ( 1 A ,− C AB ) · (−∇Φ(s)) = − 1 A ∂Φ(s) ∂x + C AB ∂Φ(s) ∂y δ(s) y = (− C AB , 1 B ) · (−∇Φ(s)) = − 1 B ∂Φ(s) ∂y + C AB ∂Φ(s) ∂x . (3.11) Según este resultado, δ(s) se puede interpretar como una modificación de la dirección original −∇Φ(s), en un grado directamente relacionado a C. Este factor C está asociado a la medida del alargamiento que posee la elipse, así como al ángulo de rotación que presenta la misma con respecto a la orientación de referencia en la que el eje mayor es horizontal. Explícitamente, C tiene la forma C = K (1− r) sin(2α), donde r es la razón entre los ejes menor y mayor de la elipse, α el ángulo de rotación de la elipse, y K una constante de escalamiento positiva. De esta expresión se desprende que para valores de r cercanos a uno, correspondientes a elipses con formas cuasi circulares, la modificación de −∇Φ(s) es mínima y δ(s) apunta prácticamente en la dirección original; lo cual significa que la aplicación de este procedimiento tuvo un efecto despreciable. Sin embargo, esta técnica de reescalamiento sí es relevante cuando las elipses son muy alargadas y, por consiguiente, r es pequeño. Asumiendo este caso y con el propósito de ahorrar operaciones, es posible omitir el cálculo de 49 C (consistente en aproximar numéricamente ∂ 2Φ(s) ∂x∂y ) y reemplazarlo por un valor aproximado que solo dependa de los coeficientes A y B. Dicha aproxi- mación dependerá de si el ángulo de rotación de la elipse es tal que origine que sus ejes se desvíen un ángulo pequeño respecto a los ejes coordenados, o de si este hace que los ejes de la elipse formen un ángulo cercano a los 45◦ con los ejes coordenados. De un análisis más completo de la ecuación general de una elipse rotada se puede deducir que cuando A B + B A > 4, el ángulo de rotación corresponde al primer tipo descrito y C ≈ 0, con lo cual δ(s) x = − 1 A ∂Φ(s) ∂x , δ(s) y = − 1 B ∂Φ(s) ∂y , (3.12a) mientras que cuando A B + B A ≤ 4, el efecto de la rotación es apreciable y C ≈ ∓A+B 3 , de modo que δ(s) x = − 1 A ∂Φ(s) ∂x ∓ 1 3( 1 A + 1 B )∂Φ(s) ∂y , δ(s) y = − 1 B ∂Φ(s) ∂y ∓ 1 3( 1 A + 1 B )∂Φ(s) ∂x , (3.12b) donde los signos − y + corresponden a ángulos de rotación agudos y obtusos, respectivamente. Cabe indicar que estos resultados, obtenidos para el caso bidimensional, pueden ser generalizados a n dimensiones. Un segundo procedimiento [22] consiste en mapear el espacio k-dimensional de los vectores b = (b1, b2, . . . , bk) que conforman el dominio de Φ hacia un nuevo espacio b∗ mediante la transformación b ∗1 = arctan(b1) b ∗2 = arctan(b2) ... b ∗k = arctan(bk)  , (3.13) que convierte las isosuperficies de Φ en el espacio original, en superficies apro- ximadamente esféricas en el nuevo espacio, mediante un proceso de acerca- miento al origen, o compresión, de las partes distantes de las isosuperficies, responsables de alargarlas y apartarlas del caso esférico. Al tomar Φ(b) la forma Φ∗(b∗) en el nuevo espacio, se tiene que las deri- 50 vadas parciales en ambos espacios cumplen la relación ∂Φ ∂bj = ∂Φ∗ ∂b ∗j ∂b ∗j ∂bj , (3.14) y puesto que de (3.13) resulta ∂b ∗j ∂bj = 1 1 + b2 j , (3.15) se obtiene finalmente la expresión de la j-ésima derivada parcial de Φ en el espacio transformado a partir de su similar en el espacio original, ∂Φ∗ ∂b ∗j = (1 + b2 j) ∂Φ ∂bj . (3.16) Dado que la gradiente en este nuevo espacio está, en general, mejor di- rigida hacia el mínimo que la gradiente en el espacio original, la presente técnica propone una dirección de minimización modificada que explota esta ventaja. Como resultado, las componentes del vector de paso δ̂(s) para la s-ésima iteración toman la forma δ̂ (s) j = ( 1 + (b(s) j )2 ) ∂Φ(s) ∂bj k∑ h=1 [( 1 + (b(s) h )2 ) ∂Φ(s) ∂bh ]21/2 . (3.17) Una tercera técnica, que toma en cuenta el comportamiento estadístico de los datos medidos que el modelo busca reproducir, es tratada en el siguiente capítulo. 3.2.3. Cálculo de la Longitud del Paso La última etapa del método consiste en calcular, para cada iteración, la longitud αs del paso en la dirección de minimización δ̂(s) que conduce del punto tentativo actual b(s) al siguiente b(s+1) según la fórmula b(s+1) = b(s) + αs δ̂ (s) . (3.18) 51 Como ya se mencionó al inicio del capítulo, el beneficio de resolver de manera exacta (3.1) para obtener el valor ideal de αs no justifica su alto costo compu- tacional sino que degrada la eficiencia del método, por lo que es preferible la opción de calcular αs mediante un procedimiento de tanteo y verificación [2] que ha probado tener un buen desempeño. Este procedimiento toma en cuenta el hecho de que cuando la longitud del paso es mucho menor que la ideal, el vector de paso de la siguiente iteración es aproximadamente colineal al de la iteración actual, mientras que cuando la longitud del paso es mayor que la ideal, los ángulos entre ambos vectores de paso son mayores a 90◦. La corrección de la longitud del paso se lleva a cabo entonces de acuerdo a este esquema, aumentándola o disminuyéndola según sea el caso. Evidentemente, en el caso que la longitud del paso sea la ideal, el ángulo entre dichos vectores de paso es exactamente 90◦. Para la s+1 -ésima iteración, el procedimiento de cálculo de la longitud del paso αs es el siguiente: 1. Calcular Φ(s) y δ̂(s). 2. Verificar que Φ(s) < Φ(s−1). Si esto no se cumple, significa que el valor de αs−1 usado en la iteración anterior fue muy grande. En ese caso dividir αs−1 entre (por ejemplo) 4 y regresar al punto 1. 3. Calcular el coseno del ángulo θ que forma el vector de paso actual con el de la iteración anterior mediante: cos θ = δ̂ (s) · δ̂ (s−1). 4. Verificar que cos θ ≥ 0 para evitar oscilaciones indeseables en el camino hacia el mínimo. Si cos θ < 0, entonces el valor de αs−1 fue muy grande. En ese caso dividir αs−1 entre (por ejemplo) 4 y regresar al punto 1. 5. Calcular por tanteo el valor de αs en función al valor de cos θ. La estra- tegia consiste en asignarle a αs un valor mayor al de αs−1 si cos θ ≈ 1 y, correspondientemente, un valor menor si cos θ ≈ 0. La siguiente fórmula a demostrado ser útil para este fin: αs = (d1 + d2 cos4 θ)αs−1, donde d1 y d2 son parámetros arbitrarios que deben cumplir con 0 < d1 < 1 52 y (1 − d1) < d2 ≤ 1. En particular, los valores d1 = 0,5 y d2 = 1 pro- ducen resultados satisfactorios. 6. Calcular a partir de (3.18) el nuevo punto tentativo b(s+1) y ejecutar la siguiente iteración partiendo del punto 1. Una de las principales ventajas de este procedimiento radica en que la longitud del paso en cada iteración se ajusta automáticamente a la topología local de Φ sin que sea necesario resolver un sistema de ecuaciones u otro tipo de operaciones complejas que no sean el simple cálculo de Φ(s) y δ̂(s). Otra ventaja importante de este procedimiento es que garantiza, tal como se muestra a continuación, la convergencia al mínimo del proceso iterativo que compone el método del máximo descenso. 3.2.4. Demostración de la Convergencia del Método La demostración de la convergencia del método del máximo descenso que aquí se presenta está basada en uno de los primeros trabajos que abordaron este tema [5]. Dado que la demostración expuesta en dicho trabajo es bastante breve y se fundamenta en una serie de proposiciones presentadas sin mayor detalle ni prueba, ha sido necesario presentar aquí la deducción de cada una de dichas proposiciones para subsanar esa carencia y proporcionar un análisis más completo. Si bien esta demostración de convergencia es válida cuando en cada itera- ción el vector de paso se toma en la dirección del negativo de la gradiente o máximo descenso, es posible (con ligeras modificaciones) extender su validez a los casos en que el vector de paso es una versión modificada de la dirección del máximo descenso producto de alguna técnica de reescalamiento. Considérese una región S en cuyo interior y contorno se asume Φ con pri- meras derivadas parciales continuas. Asimismo, sea C el camino segmentado que parte de algún punto tentativo inicial b(0) al interior de S y sigue por el segmento de recta que lo une con el siguiente punto tentativo b(1) hasta alcanzar primero ya sea el borde de S o el punto b(1), y así sucesivamente con los puntos tentativos siguientes. Se asume además que a lo largo del camino C, producido por la secuencia de los b(s), Φ es monótona decreciente. 53 Bajo este escenario existen tres posibilidades: (1) el camino C termina en algún punto del borde de S, (2) el camino C termina en un punto con gradiente cero, i.e., en un punto estacionario de Φ, (3) el proceso continúa indefinidamente. La primera posibilidad queda excluída si se elije b(0) de manera que el valor de Φ en dicho punto, denotado por Φ(0), sea menor que el valor de Φ en cualquier punto del borde de S. Para la segunda posibilidad, la prueba de convergencia es trivial. El caso relevante es por lo tanto el tercero, en que el proceso continúa indefinidamente. De la condición previa sobre la diferenciabilidad de Φ y del supuesto que esta en el punto b(0) interior a S toma un valor menor a los del contorno de S, se desprende la existencia de un punto estacionario, bmin ∈ S, en el cual Φ adquiere su valor mínimo Φmin en S. Establecido esto, la demostración de convergencia consiste en probar que el método converge a un punto estacio- nario de Φ, que en general puede corresponder a un mínimo relativo o a un punto de ensilladura (saddle point). El paso inicial de esta demostración consiste en determinar la conver- gencia de la secuencia de puntos tentativos b(0), b(1), b(2), . . ., que al mismo tiempo genera la serie Φ(0),Φ(1),Φ(2), . . ., la cual por ser monótona decreciente y acotada inferiormente por Φmin posee un valor límite Φ0 = ĺıms→∞Φ(s). Es- to significa que los puntos b(s) ‘convergen’ hacia la isosuperficie Φ(b) = Φ0, lo que implica la existencia de un umbral Nε a partir del cual los puntos b(Nε), b(Nε+1), b(Nε+2), . . ., están contenidos en la región Ωε limitada por las isosuperficies Sε : Φ(b) = Φ0 + ε y S0 : Φ(b) = Φ0, donde ε > 0 es un pará- metro que puede hacerse arbitrariamente pequeño para así conseguir que las isosuperficies Sε y S0 tiendan a ser paralelas y separadas una distancia dε infinitamente pequeña. Afirmar que la secuencia de los b(s) no converge en un punto equivale a afirmar que la serie {Ds = ‖b(s+1) − b(s)‖} de las distancias entre un punto tentativo y el siguiente no convergen en cero. En este caso siempre será posible encontrar, para un ε suficientemente pequeño, un N ′ε ≥ Nε tal que DN ′ε > dε, 54 lo que significa que b(N ′ε) y b(N ′ε+1) están al interior de Ωε y separados una distancia mayor a su grosor dε. Sin embargo, el paso que conduce de b(N ′ε) a b(N ′ε+1), según el método del máximo descenso, es ortogonal a las isosuperficies aproximadamente paralelas Sε y S0 que delimitan Ω. Esto implica necesariamente que b(N ′ε+1) quede fuera de Ω, ya que la longitud del paso DN ′ε es mayor que la distancia dε entre Sε y S0. Se observa entonces que asumir que la secuencia de los b(s) no es convergente conlleva a conclusiones que se contradicen entre sí. Queda así finalmente establecido que la secuencia b(0), b(1), b(2), . . ., con- verge en el punto b(∞), para el cual se cumple Φ(∞) < Φ(s), s = 0, 1, 2, . . . (3.19) Como paso final de la demostración, resta probar que b(∞) es un punto esta- cionario de Φ, i.e., que se cumple la condición ∇Φ(∞) = 0. Adóptese la siguiente notación para el módulo y la dirección del negativo de la gradiente respectivamente, h(b) = ∥∥∥∇Φ(b) ∥∥∥, δ̂(b) = − ∇Φ(b) ‖∇Φ(b)‖ , (3.20) de modo que ∇Φ(b) = −h(b)δ̂(b) = −δ(b). (3.21) Del supuesto h(b(∞)) = h(∞) 6= 0, contrario a lo que se pretende probar, y de la continuidad de ∇Φ(b) en b(∞), se desprende la existencia de una región U consistente en una bola de radio R(ε) centrada en b(∞) en la cual se cumple ∥∥∥∇Φ(b)−∇Φ(∞) ∥∥∥ < εh(∞) (3.22) para todos los b en U . Del teorema de la desigualdad triangular aplicado al triángulo formado por los vectores ∇Φ(b) y ∇Φ(∞) se obtiene ∣∣∣ ∥∥∥∇Φ(b) ∥∥∥−∥∥∥∇Φ(∞) ∥∥∥ ∣∣∣ < ∥∥∥∇Φ(b)−∇Φ(∞) ∥∥∥ < ∥∥∥∇Φ(b) ∥∥∥+ ∥∥∥∇Φ(∞) ∥∥∥, (3.23) 55 de cuya parte izquierda, y utilizando (3.22), se concluye ∣∣∣h(b)− h(∞) ∣∣∣ < εh(∞). (3.24) Para derivar una desigualdad análoga a (3.24) que acote la desviación de los δ̂(b) en U respecto a δ̂(∞), se reescribe (3.24) en la forma ∣∣∣∣h(b) h(∞) − 1 ∣∣∣∣ < ε, (3.25) y se divide (3.22) entre h(∞) para obtener ∥∥∥∥h(b) h(∞) δ̂(b)− δ̂(∞) ∥∥∥∥ < ε. (3.26) Esta desigualdad se expresa convenientemente en la forma equivalente ∥∥∥∥(h(b) h(∞) − 1 ) δ̂(b) + ( δ̂(b)− δ̂(∞))∥∥∥∥ < ε, (3.27) lo que implica ∣∣∣∣ ∥∥∥δ̂(b)− δ̂(∞)∥∥∥− ∣∣∣h(b) h(∞) − 1 ∣∣∣ ∣∣∣∣ < ε, (3.28) y a la vez ∣∣∣h(b) h(∞) − 1 ∣∣∣− ε < ∥∥∥δ̂(b)− δ̂(∞)∥∥∥ < ∣∣∣h(b) h(∞) − 1 ∣∣∣+ ε, (3.29) de cuya parte derecha y de (3.25) finalmente se concluye ∥∥∥δ̂(b)− δ̂(∞)∥∥∥ < 2ε. (3.30) En general, para dos b1 y b2 cualesquiera en U se tienen entonces las expresiones ∥∥∥δ̂(b1)− δ̂(∞)∥∥∥ < 2ε y ∥∥∥δ̂(b2)− δ̂(∞)∥∥∥ < 2ε, (3.31) que sumadas miembro a miembro resultan ∥∥∥δ̂(b1)− δ̂(∞)∥∥∥+ ∥∥∥δ̂(b2)− δ̂(∞)∥∥∥ < 4ε, (3.32) 56 á θ ngulo max delimita K, se cumple: cos max = θ ε θmaxθmax U K b ( )∞ bo δ ∞ ( )ê(bo) ˆ C Figura 3.4: Región cónica K ⊂ U de todos los puntos b que cumplen ê(b) · δ̂ (∞) ≥ ε. Se demuestra que Φ(b) < Φ(∞), ∀ b ∈ K. lo que implica, por el teorema de desigualdad triangular, que ∥∥∥δ̂(b1)− δ̂(b2) ∥∥∥ < 4ε. (3.33) Para dos puntos tentativos consecutivos cualesquiera, b(s) y b(s+1), se cumple, debido a la ortogonalidad entre δ̂(s) y δ̂(s+1), que ‖δ̂ (s+1) − δ̂ (s) ‖ = √ 2. (3.34) Suponiendo que dichos b(s) y b(s+1) estén contenidos en U , debe cumplirse por (3.33) que ‖δ̂(s+1) − δ̂ (s) ‖ < 4ε, lo cual contradice (3.34) para valores de ε menores a √ 2/4. De esto se concluye que b(s) y b(s+1) no pueden estar simultáneamente contenidos en U . Definiendo asimismo para cada punto b en U el vector ê(b) = b− b(∞) ‖b− b(∞)‖ , (3.35) correspondiente a la dirección de su posición respecto a b(∞), es posible deli- mitar una región cónica K ⊂ U , como se muestra en la figura 3.4, compuesta por todos los puntos b ∈ U en los que el ángulo θ(b) entre su dirección ê(b) y la dirección fija δ̂(∞) es menor o igual a un ángulo máximo θmax cuyo coseno 57 es ε, i.e., cos θ(b) = ê(b) · δ̂(∞) ≥ ε = cos θmax, siempre que 0 < ε ≤ 1. (3.36) El valor de Φ para un punto b0 cualquiera enK se puede calcular mediante la integral Φ(b0) = Φ(∞) − ∫ C (−∇Φ(b)) · ê(b) dc, (3.37) donde el camino de integración C es una línea recta que une b(∞) con b0. De (3.21) y (3.22) se deriva ∇Φ(b) = −h(∞)δ̂ (∞) + γ(b), con ∥∥∥γ(b) ∥∥∥ < εh(∞), (3.38) con lo cual el integrando I(b) de la integral en (3.37) adquiere conveniente- mente la forma I(b) = [ h(∞)δ̂ (∞) − γ(b) ] · ê(b), (3.39) y puesto que γ(b) · ê(b) ≤ ∥∥∥γ(b) ∥∥∥ < εh(∞), (3.40) se obtiene para I(b) la cota inferior I(b) > h(∞) [ δ̂ (∞) · ê(b)− ε ] . (3.41) De esta cota inferior y de la condición en (3.36), que satisfacen los b en K, resulta I(b) > 0. (3.42) A partir de (3.37), y considerando que el integrando es positivo para todos los puntos en K (incluídos los que forman el camino de integración C), se concluye para todos los puntos b0 contenidos en K que Φ(b0) < Φ(∞). (3.43) La última parte de la demostración de convergencia requiere la existencia de una región V construída a partir de una bola, de cierto radio R(ε′), con- céntrica e interior a U a la que se le extrae la región cónica correspondiente a su intersección con K. Dicha región V debe cumplir la condición de que el 58 θmaxθmax b ∞( ) ∞( ) δ̂ ( ) ( ) h∞ ∞ε' h b δ̂ ( )∞ δ̂(b) δ(b) δ∞( ) φmax (b)γ superficie esf rica que confina el punto terminal de todos los é δ(b) en la bola de radio ε'R( ) se cumple: 1) cos θmax = ε 2) sen φmax = ε' Figura 3.5: Construcción geométrica de δ(b), cuyo punto terminal pertenece a la bola de radio ε′h(∞). El ángulo máximo ϕmax se obtiene cuando δ(b) es tangente a la superficie exterior de la bola. rayo en la dirección δ̂(b) desde cualquier punto b en V atraviese K. Antes de probar la existencia de dicha región V , es conveniente determinar el valor máximo ϕmax que puede asumir el ángulo ϕ(b) entre δ̂(b) y δ̂(∞) en cualquier punto b de la bola de radio R(ε′). De la definición de continuidad utilizada para deducir (3.22), se deduce igualmente para los puntos de dicha bola que δ(b) = h(∞)δ̂ (∞) + γ(b), con ∥∥∥γ(b) ∥∥∥ < ε′h(∞), (3.44) La interpretación geométrica de esta expresión, ilustrada en la figura 3.5, indica que si se coloca en b el punto inicial de δ(∞) y δ(b), entonces el punto terminal de este último estará en la bola con centro en el punto terminal de δ(∞) y radio ε′h(∞). Esto significa que el caso de ángulo máximo, ϕ(b) = ϕmax, se presenta cuando algún δ(b) particular es tangente a la bola. Por lo tanto senϕmax = ε′. (3.45) Cabe indicar que si el valor de ε′ es tal que ocasiona que ϕmax > θmax en la bola de radio R(ε′) correspondiente, puede entonces existir en dicha 59 bola algún punto b desde el cual el rayo en la dirección δ̂(b) posea un ángulo ϕ(b) > θmax respecto a δ̂(∞) que ocasione que este no atraviese a K. Tal bola, en consecuencia, no puede dar origen a la región V , pues no presenta la propiedad que por definición debe poseer V . Para evitar esta situación se introduce la condición ε′ < √ 1− ε2, con lo cual : ϕmax < θmax. (3.46) Considérese a continuación el conoQ construído tomando como generatriz la familia de rayos que forman un ángulo ϕmax con δ̂ (∞) y como vértice algún punto sobre el rayo que parte de b(∞) siguiendo la dirección −δ̂(∞), de manera tal que su base coincida con la base del cono que delimita K. Para el valor umbral ε′ = √ 1− ε2, definido en (3.46), se tiene que el vér- tice de Q está en b(∞) y por ende Q coincide con el cono de K. Si a partir de dicho umbral se disminuye progresivamente el valor de ε′, ocurre simul- táneamente que la generatriz de Q se torna cada vez más vertical, causando que el vértice del mismo se aleje progresivamente de b(∞), mientras por otro lado la región V se contrae al disminuir el radio R(ε′) que la limita. En el transcurso de este proceso se llega a un estado, mostrado en la figura 3.6a, en el que V resulta tangente e interior a Q. El valor de ε′ que produce este estado se obtiene resolviendo la ecuación √ 1− ε2 √ 1− ε′2 − ε ε′ = R(ε′)/R(ε). (3.47) Es posible probar gráficamente que la región V construída a partir de la bola de radio ε′ obtenido de (3.47), efectivamente cumple el requisito de que desde cualquier punto interior b partan rayos con respectivas direcciones δ̂(b) que atraviesen K. Considérese para ello el caso crítico en el que los rayos parten de los puntos ubicados en la circunferencia de contacto C1 entre V y Q. Probar que dichos rayos atraviesan K es suficiente para garantizar que ocurra lo mismo para el resto de puntos en V . En la figura 3.6b se muestra el segmento de una línea generatriz de Q que va del punto b1 en C1 al punto b2 contenido en el plano de la circunferencia de intersección C2 entre Q y K. Dado que el ángulo que forma dicho segmento con b(∞) es ϕmax, todos los rayos que partan de b1 y posean ese mismo 60 b ∞( ) ∞( ) δ̂ θmax φmax R(ε') R(ε) θmax - φmax U K V Q C1 C2 se cumple: (a) √(1-ε²) (1-ε'²) - = R( )/√ εε' ε' R(ε) b ( )∞ b1 b2 C1 C2 b1 b2 (b) Figura 3.6: (a) Construcción de la región V a partir de un ε′ que la hace tangente e interior a Q, garantizando así que todos los rayos que partan de algún b en V con dirección δ̂(b) atraviesen K. (b) Vista de V , Q y K en la dirección de δ̂ (∞) , donde se observa el caso crítico en el que los rayos que parten de cualquier punto en la circunferencia de contacto C1 entre V y Q caen necesariamente dentro de la región circular sombreada ubicada en el plano que contiene la circunferencia de intersección C2 entre Q y K. ángulo máximo intersectarán el plano de C2 formando una circunferencia circundante a la región circular sombreada que aparece en la figura. Todos los demás rayos con ángulos menores a ϕmax atravesarán necesariamente dicha región sombreada. Esto significa que los rayos que partan de b1 en todas las posibles direcciones permitidas atravesarán K puesto que la región sombreada es interior a K. 61 Aunque esta prueba gráfica asume un espacio geométrico tridimensional, sus resultados pueden extenderse al caso general en k dimensiones. De esta manera queda probada la existencia de la región V . Finalmente, dado que V y K rodean b(∞), y este a su vez es el punto de convergencia de la secuencia b(0), b(1), b(2), . . ., existe entonces algún punto b(s) de la secuencia contenido en V , desde el cual el rayo con dirección δ̂(s) atraviesa necesariamente K. Asimismo, sea ζ un punto sobre la parte de dicho rayo interior a K. Este punto ζ se encuentra antes del punto b(s+1) en el recorrido a lo largo del camino de minimización C seguido por el método, pues b(s+1) está impedido por (3.33) y (3.34) de pertenecer a U . Del caracter monótono decreciente de Φ en C y de (3.37) se tiene por lo tanto Φ(s+1) < Φ(ζ) < Φ(∞), resultado que contradice el supuesto inicial dado en (3.19). Esta contradicción, al ser consecuencia de asumir h(∞) 6= 0, demuestra que efectivamente b(∞) es un punto estacionario de Φ en el cual se cumple ∇Φ(∞) = 0. 3.3. Método de Gauss-Newton A diferencia del método de la sección anterior en el que las característi- cas del modelo f a ajustar intervienen solo indirectamente, en el método de Gauss-Newton se utiliza una versión lineal de f correspondiente a su expan- sión en series de Taylor de primer grado alrededor del punto tentativo b(s) de la iteración anterior para determinar la dirección de minimización sobre la cual se buscará el punto tentativo b(s+1) de la iteración actual. 3.3.1. Elección de la Dirección de Minimización Como se mencionó, la esencia del método se basa en sustituir el modelo no lineal f de k parámetros de ajuste b = (b1, b2, . . . , bk) y m variables inde- pendientes xi = (x1,i, x2i, . . . , xmi) que representan cada uno de los n juegos de datos de entrada (i = 1, 2, . . . , n), ŷi = f(xi, b) = fi(b), 62 donde ŷi es la respuesta predicha por f correspondiente a la i-ésima respuesta medida yi, por su expansión en serie de Taylor truncada en los términos lineales 〈yi〉 = 〈fi(b(s) + δ)〉 = fi(b(s)) + k∑ j=1 ∂fi(b(s)) ∂bj δj, (3.48) que en notación vectorial toma la forma 〈y〉 = f (s) 0 + P (s)δ, (3.48a) donde evidentemente f (s) 0 [n×1] = ( fi(b(s)) ) , i = 1, 2, . . . , n, y P (s) [n×k] = ∂fi(b(s)) ∂bj ) , i = 1, 2, . . . , n; j = 1, 2, . . . , k. Los corchetes 〈 〉 se emplean para distinguir las respuestas predichas por la aproximación lineal del modelo de aquellas predichas por el modelo no lineal original. La función de discrepancia original Φ dada en (2.3) es aproximada utili- zando esta versión lineal del modelo mediante la expresión 〈Φ〉 = n∑ i=1 (yi − 〈yi〉)2, (3.49) con lo cual el valor buscado es ahora aquel δ que minimice 〈Φ〉. Dado que δ aparece linealmente en (3.48), el cálculo del mismo se puede llevar a cabo directamente resolviendo el sistema de ecuaciones ∂〈Φ〉 ∂δj = 0, j = 1, 2, . . . , k, (3.50) cuya forma explícita es equivalente a la ecuación matricial A(s)δ = g(s), (3.51) donde A (s) [k×k] = [P (s)]TP (s) y g (s) [k×1] = [P (s)]T (y − f (s) 0 ). 63 El valor así obtenido, δ = [A(s)]−1g(s), (3.52) define la dirección de minimización δ̂ (s) = δ̂ a emplearse en la iteración actual. En la práctica ha probado ser conveniente considerar únicamente una fracción de δ como vector de paso δ(s); caso contrario la extrapolación podría estar fuera de la región donde (3.48) constituye una buena aproximación de f , lo que eventualmente originaría la divergencia del proceso iterativo. Como consecuencia, el método establece, para la (s+ 1)-ava iteración, el vector de paso δ(s) = κ0 δ, 0 < κ0 ≤ 1. (3.53) 3.3.2. Cálculo de la Longitud del Paso De (3.53) se desprende la longitud del vector de paso αs = κ0 ‖δ‖, (3.54) cuyo cálculo se reduce a determinar el valor κ0 que minimice, dentro del intervalo ]0, 1], a la función g(κ) = Φ(b(s) + κδ) κ ∈ ]0, 1], (3.55) obtenida al evaluar Φ en el punto tentativo resultante para cada valor de κ. Dado que llevar a cabo la tarea de determinar de manera exacta dicho valor κ0 es impráctico por el alto costo computacional involucrado, se pre- fiere resolver el problema de manera aproximada. Un procedimiento simple y confiable para tal fin consiste reemplazar g(κ) por la parábola g∗(κ) que pasa por los puntos (0, g(0)), (1/2, g(1/2)) y (1, g(1)), i.e., g∗(κ) = 2 [g(0)− 2g(1/2) + g(1)]κ2 − [3g(0)− 4g(1/2) + g(1)]κ+ g(0), (3.56) 64 y aproximar el valor κ0 buscado mediante el valor κ0 = 1 2 + 1 4 g(0)− g(1) g(0)− 2g(1/2) + g(1) (3.57) correspondiente al mínimo de la parábola. Dado que este mínimo debe estar en el intervalo mencionado, se le asigna el valor uno en caso este sobrepase dicho límite. Finalmente, se debe comprobar que Φ(s+1) < Φ(s), caso contrario habría que reducir la longitud del paso a la mitad. 3.3.3. Demostración de la Convergencia del Método La demostración de la convergencia del método de Gauss-Newton que aquí se presenta es básicamente una adaptación del procedimiento empleado en la demostración descrita en [12]. Como paso previo se asume que el modelo f(xi, b) = fi(b) satisface las siguientes condiciones: a) Para todos los vectores de entrada xi (i = 1, 2, . . . , n), tanto las prime- ras derivadas de fi(b) respecto a los bj ∂fi(b) ∂bj , j = 1, 2, . . . , k, como las segundas derivadas ∂2fi(b) ∂bj∂bh , j = 1, 2, . . . , k, h = 1, 2, . . . , k, son funciones continuas de b. Este supuesto permite no solo la existencia de las primeras derivadas de la función de discrepancia original Φ ∂Φ ∂bj (b) = −2 n∑ i=1 (yi − fi(b)) ∂fi(b) ∂bj , sino también la aproximación lineal del modelo según (3.48) y, con ello, la definición del sistema de ecuaciones correspondiente a los mínimos cuadrados lineales dado en (3.51). 65 b) Para todo vector u = (u1, u2, . . . , uk), con ‖u‖ 6= 0, y todo b dentro de cierto subconjunto convexo S del espacio de parámetros, se cumple n∑ i=1  k∑ j=1 uj ∂fi(b) ∂bj 2 > 0, (3.58) lo que significa que el conjunto de n vectores vi = ∂fi(b) ∂b1 , ∂fi(b) ∂b2 , . . . , ∂fi(b) ∂bk ) , i = 1, 2, . . . , n, contiene una base del espacio de parámetros Rk. Este supuesto garantiza la existencia de un único vector solución δ del sistema de ecuaciones correspondiente a los mínimos cuadrados lineales en (3.51), ya que de él se desprende que el rango de la matriz A(s) [k×k] de dicho sistema de ecuaciones es justamente k, lo que equivale a que A(s) sea invertible. Una manera sencilla de derivar (3.58) consiste en considerar (3.48) y (3.49) para reescribir las ecuaciones (3.51) en la forma (y1−〈y1〉) ∂f1(b) ∂b1 + (y2−〈y2〉) ∂f2(b) ∂b1 + . . . + (yn−〈yn〉) ∂fn(b) ∂b1 = 0 (y1−〈y1〉) ∂f1(b) ∂b2 + (y2−〈y2〉) ∂f2(b) ∂b2 + . . . + (yn−〈yn〉) ∂fn(b) ∂b2 = 0 ... ... . . . ... ... (y1−〈y1〉) ∂f1(b) ∂bk + (y2−〈y2〉) ∂f2(b) ∂bk + . . . + (yn−〈yn〉) ∂fn(b) ∂bk = 0  , (3.59) donde se evidencia que estas serán linealmente independientes si y solo si para todo u, con ‖u‖ 6= 0, se cumple (y1 − 〈y1〉) ( u1 ∂f1(b) ∂b1 + u2 ∂f1(b) ∂b2 + . . . + uk ∂f1(b) ∂bk ) + (y2 − 〈y2〉) ( u1 ∂f2(b) ∂b1 + u2 ∂f2(b) ∂b2 + . . . + uk ∂f2(b) ∂bk ) + ... (yn − 〈yn〉) ( u1 ∂fn(b) ∂b1 + u2 ∂fn(b) ∂b2 + . . . + uk ∂fn(b) ∂bk ) 6= 0, (3.60) expresión de la cual es posible deducir el requerimiento (3.58). 66 c) Dado el valor Φ0 = ı́nf S̄ Φ(b), donde S̄ es el complemento de S, es posible encontrar un vector b(0) al interior de S tal que Φ(b(0)) < Φ0. Para los fines de esta demostración es conveniente expresar (3.51) como un sistema de ecuaciones en forma explícita 2 k∑ j=1 n∑ i=1 ∂fi(b(s)) ∂bh ∂fi(b(s)) ∂bj ) δj = −∂Φ(b(s)) ∂bh , h = 1, 2, . . . , k, (3.61) el cual por la condición b) siempre posee un único vector solución δ para cada b(s) en S. Como se indicó en la subsección anterior, para todo b(s) se tiene un vector de paso δ(s) expresado en función de la correspondiente solución δ como δ(s) = κ0 δ, 0 < κ0 ≤ 1, donde el factor κ0, introducido en (3.53), se calcula de manera que en el si- guiente punto tentativo b(s+1) se tenga Φ(s+1) < Φ(s). Por lo tanto, la secuen- cia Φ(0),Φ(1),Φ(2), . . ., es decreciente. Asimismo, si b(0) se escoge al interior de S se concluye entonces por la condición c) que todos los vectores de la secuencia b(0), b(1), b(2), . . ., están igualmente contenidos en S. La secuencia de los Φ(s) (s = 0, 1, 2, . . .) converge a un valor finito por ser decreciente y acotada inferiormente. Bajo esta premisa es posible demostrar, siguiendo el mismo procedimiento utilizado para el caso análogo en el método del máximo descenso (ver pág. 53), que la secuencia de los b(s) (s = 0, 1, 2, . . .) converge a un vector denotado como b(∞), i.e., se cumple ĺım s→∞ b(s) = b(∞) (3.62) e igualmente, por la continuidad de Φ respecto a b, ĺım s→∞ Φ(s) = Φ(b(∞)) = Φ(∞). (3.63) 67 Se requiere demostrar ahora que en el punto b(∞) todas las derivadas parciales de Φ son cero. Si se introduce δ∗ como la solución del sistema de ecuaciones 2 k∑ j=1 n∑ i=1 ∂fi(b(∞)) ∂bh ∂fi(b(∞)) ∂bj ) δ ∗j = −∂Φ(b(∞)) ∂bh , h = 1, 2, . . . , k, (3.64) y se asume, contrario a lo que se quiere demostrar, que no todas las derivadas parciales de Φ en b(∞) son cero (∇Φ(b(∞)) 6= 0), entonces se concluye k∑ j=1 ( δ ∗j )2 > 0, puesto que la matriz de dicho sistema de ecuaciones es invertible. Considerando tanto este resultado como la condición c), y tomando el producto escalar de cada miembro de (3.64) con δ∗, se obtiene k∑ j=1 ∂Φ(b(∞)) ∂bj δ ∗j = −2 n∑ i=1  k∑ j=1 δ ∗j ∂fi(b(∞)) ∂bj 2 < 0, (3.65) lo que significa que en el punto b(∞) la derivada de Φ en la dirección δ̂∗ = δ̂ (∞) es negativa. Asimismo, la continuidad de∇Φ(b) y δ(b) en S implica que para cada b(s) en una pequeña vecindad alrededor de b(∞) la derivada de Φ en la respectiva dirección de minimización δ̂(s) es también negativa. A partir de este resultado se puede demostrar que para todos los b(s) en la vecindad mencionada se cumple que el valor de Φ(s+1) está por debajo del valor de Φ(s) en al menos una cierta cantidad fija ε. Por lo tanto, la secuencia de los Φ(s) tiende a un valor infinito negativo, lo cual contradice (3.63) que establece que tal secuencia converge en el valor finito Φ(∞). Dado que el supuesto inicial (∇Φ(b(∞)) 6= 0) conduce a una contradicción, queda entonces demostrado que todas las derivadas parciales de Φ en el punto límite b(∞) son cero. Capítulo 4 Método de Levenberg-Marquardt Cuando el problema de mínimos cuadrados aumenta en complejidad, los métodos de Gauss-Newton (de la serie de Taylor) y del máximo descenso (de la gradiente) se vuelven inaplicables o muy ineficientes debido a sus limita- ciones. En el caso del método de Gauss-Newton la principal deficiencia es la potencial divergencia de los sucesivos parámetros tentativos. Cuando la ite- ración produce un nuevo punto tentativo lejano al anterior alrededor del cual se linealizó el modelo, el supuesto de linealidad deja de ser válido y el valor de Φ en el nuevo punto tentativo podría aumentar respecto al anterior. En ese caso, la secuencia de puntos tentativos varían de forma errática, produciendo una convergencia muy lenta o incluso llegando a divergir. Por otro lado, la principal deficiencia del método del máximo descenso es la convergencia considerablemente lenta que presenta luego de las primeras iteraciones. Como es bien sabido, cuando el modelo es lineal las isosuperficies de Φ son elipsoides concéntricos, mientras que cuando el modelo es no lineal, las isosuperficies se distorsionan según el grado de la no linealidad. En este último caso, la dirección de minimización (i.e., el negativo de la gradiente) en puntos lejanos al mínimo puede no ser adecuada. Independientemente de que el modelo sea o no lineal, las isosuperficies en puntos cercanos al mínimo son prácticamente elipsoides, que en la mayoría de los casos poseen formas alargadas. Por consiguiente, si algún punto tentativo se encontráse en los extremos de dichas elipsoides alargadas, los pasos subsiguientes serían 68 69 progresivamente más pequeños, ralentizándose así la convergencia (véase la fig. 3.2 en la pág. 46). Varias modificaciones han sido sugeridas para limitar la longitud del paso en el método de Gauss-Newton y para acondicionar mediante reescalamien- to la topología de Φ en el método del máximo descenso. Sin embargo, los problemas de divergencia por un lado, y de convergencia lenta por el otro, no han podido ser satisfactoriamente resueltos. Esto genera la necesidad de contar con un método que supere las limitaciones de ambos métodos y al mismo tiempo conserve sus ventajas. Uno de los métodos que cumple tales requisitos es el llamado método de Levenberg-Marquardt. En el presente capítulo se exponen los razonamientos seguidos tanto por Levenberg [20] como por Marquardt [21] al desarrollar, cada uno por su cuenta, el método que lleva sus nombres. 4.1. Enfoque de Levenberg El principal objetivo de Levenberg consistía en limitar o atenuar la longi- tud del paso dado por el método de Gauss-Newton estándar en cada iteración. Un paso que conduzca a un punto tentativo dentro de la región de validez de la aproximación lineal, evita que el valor de Φ se incremente con cada iteración, y por lo tanto que ocurra divergencia. 4.1.1. Construcción del Método Como se vió en el capítulo anterior, el método de Gauss-Newton están- dar resuelve iterativamente el problema de mínimos cuadrados no lineales consistente en estimar el vector de parámetros bmin que minimice Φ = n∑ i=1 r2 i , donde ri = yi − fi(b), es el residuo entre el i-ésimo dato medido yi y su correspondiente valor pre- dicho por el modelo no lineal fi evaluado en algún vector arbitrario de pará- metros b. 70 El proceso iterativo mencionado consiste en una secuencia de problemas de mínimos cuadrados lineales cuyas soluciones se van acercando progresiva- mente a la solución buscada bmin. En estos problemas se reemplaza el modelo no lineal fi por una versión lineal 〈fi〉 correspondiente a su serie de Taylor de primer grado alrededor de el punto tentativo actual b(s). En este punto, que representa la aproximación de bmin obtenida tras la s-ésima iteración, se asume que Φ no posee un valor estacionario. Tras introducir el vector de incrementos δ = b− b(s) se obtiene 〈fi(b(s) + δ)〉 = fi(b(s)) + ∂fi(b(s)) ∂b1 δ1 + ∂fi(b(s)) ∂b2 δ2 + · · ·+ ∂fi(b(s)) ∂bk δk, (4.1) de donde se desprende la versión lineal de los residuos 〈ri〉 = yi − 〈fi(b(s) + δ)〉. (4.2) Dicho esto, el problema lineal consiste entonces en estimar el vector δ∗ que minimice 〈Φ(δ)〉 = n∑ i=1 〈ri〉2, (4.3) el cual se obtiene resolviendo el sistema de ecuaciones lineales resultante de igualar a cero las derivadas parciales de 〈Φ〉 respecto a cada una de las componentes de δ, sistema al que se conoce como ecuaciones normales lineales 1 2 ∂〈Φ〉 ∂δ1 = [c1c1]δ1 + [c1c2]δ2 + . . . + [c1ck]δk + [c1c0] = 0 1 2 ∂〈Φ〉 ∂δ2 = [c2c1]δ1 + [c2c2]δ2 + . . . + [c2ck]δk + [c2c0] = 0 ... ... ... . . . ... ... ... 1 2 ∂〈Φ〉 ∂δk = [ckc1]δ1 + [ckc2]δ2 + . . . + [ckck]δk + [ckc0] = 0  (4.4) donde los corchetes [ ], introducidos en la pág. 14, representan el producto interno de los vectores contenidos, de forma tal que [c1c0] = n∑ i=1 ∂fi(b(s)) ∂b1 [ fi(b(s))− yi ]) , [c1c1] = n∑ i=1 ∂fi(b(s)) ∂b1 ∂fi(b(s)) ∂b1 ) , [c1c2] = n∑ i=1 ∂fi(b(s)) ∂b1 ∂fi(b(s)) ∂b2 ) , . . . (4.5) 71 En las versiones modificadas del método de Gauss-Newton el vector de paso se define como δ(s) = κ0δ ∗, donde 0 < κ0 < 1, mientras en la versión estándar del mismo se toma κ0 = 1. De esta manera, si el módulo de δ∗ es lo suficientemente grande, la aproximación lineal (4.1) no sería válida alrededor del siguiente punto tentativo, con lo cual una disminución en 〈Φ〉 podría no corresponder a una disminución en Φ. En tales casos, Levenberg propone limitar o atenuar el módulo del vector de incrementos δ con el fin de mejorar la aproximación de Taylor de primer grado en (4.1) y simultáneamente minimizar el valor de 〈Φ〉 en (4.3) bajo estas condiciones de atenuación. Asimismo, propone que para conseguir que tanto el módulo de δ como el valor de 〈Φ〉 sean pequeños, se emplee la idea fundamental de los mínimos cuadrados. Como resultado, se origina el llamado método de Levenberg-Marquardt en el que la expresión a minimizar es Φ̄w = w〈Φ(δ)〉+Q(δ), con : Q(δ) = ‖δ‖2, (4.6) donde w es un factor de ponderación positivo que expresa la importancia relativa entre 〈Φ〉 y Q en el proceso de minimización. Por comodidad se opta por trasladar el origen de coordenadas del espacio de parámetros al punto b(s). Por consiguiente, el punto en el cual Φ̄w alcanza su valor mínimo es el vector de incrementos δw que minimiza (4.6). Este vector de incrementos es el vector de paso δ(s) calculado por método. Para obtener δw, las derivadas parciales de Φ̄w respecto a las distintas componentes del vector de incrementos δ se igualan a cero, con lo cual ∂Φ̄w ∂δj = w ∂〈Φ〉 ∂δj + 2δj = 0, j = 1, 2, . . . , k. (4.7) Al dividir cada una de estas ecuaciones entre 2w, y substituir las expresiones de las derivadas parciales de 〈Φ〉 por (4.4), resultan las ecuaciones normales 72 atenuadas ([c1c1] + w−1)δ1 + [c1c2]δ2 + . . . + [c1ck]δk = −[c1c0] [c2c1]δ1 + ([c2c2] + w−1)δ2 + . . . + [c2ck]δk = −[c2c0] ... ... . . . ... ... [ckc1]δ1 + [ckc2]δ2 + . . . + ([ckck] + w−1)δk = −[ckc0]  (4.8) Se puede apreciar que estas ecuaciones son iguales a las ecuaciones normales ordinarias (4.4), excepto por los coeficientes de la diagonal principal, los cua- les están incrementados en w−1. Al igual que para (4.4), la solución de (4.8) se puede obtener mediante métodos simplificados que explotan la simetría de la matriz de coeficientes del sistema. Asimismo, puede verse que un valor de w que tiende al infinito, da co- mo resultado el método de Gauss-Newton estándar, que por lo tanto viene a ser un caso especial del método de Levenberg-Marquardt, el cual puede interpretarse como un método de mínimos cuadrados con atenuación. Para ilustrar la dependencia de este nuevo método respecto al parámetro w, la representación geométrica del problema de minimización de Φ̄w resulta ser de gran utilidad. Con ese fin, considérense simultáneamente la familia de isosuperficies esféricas concéntricas en b(s) y la de isosuperficies elipsoidales concéntricas en b(s) +δ∗, correspondientes a Q y 〈Φ〉 respectivamente. Nótese además que los valores de Q y 〈Φ〉 aumentan conforme las isosuperficies se van alejando de sus centros. Asimismo, las isosuperficies de w〈Φ〉 son las mismas que las de 〈Φ〉 para todos los valores de w. Sea Q0 una isosuperficie cualquiera de Q, y G el grupo de isosuperficies de 〈Φ〉 que la intersectan. Para determinar el punto de Q0 en el que Φ̄w toma el menor valor, basta con identificar en G la isosuperficie 〈Φ〉0 en la que 〈Φ〉 es menor a las demás, ya que el valor deQ es el mismo en todos los puntos deQ0. Tal 〈Φ〉0 es, de acuerdo a lo expresado en el párrafo anterior, la isosuperficie en G más cercana al centro común. En este punto cabe observar que si uno se desplaza secuencialmente por las isosuperficies de G de manera que estas se vayan acercando progresivamente al mínimo, se tiene que el límite de dicha secuencia es 〈Φ〉0. De esto se concluye que 〈Φ〉0 hace contacto con Q0 en un único punto, en el cual además ambas son tangentes. El anterior análisis cualitativo puede corroborarse analíticamente de ma- 73 C LM C GN C MD wa w b w = 0 w δo= 0 δ∞= δ* (b(s)) C LM : Levenberg-Marquardt C GN : Gauss-Newton C MD : M ximo Descensoá se cumple: wa < w b isosuperficies de Q isosuperficies de ‹Φ› =∞ Figura 4.1: Camino de minimización CLM del método de Levenberg-Marquardt, formado por los puntos que minimizan Φ̄w, para todos los valores de w. El camino CLM parte de δ0 y termina en δ∞ pasando por los puntos donde las isosuperficies de Q y 〈Φ〉 son tangentes. Se puede afirmar que CLM es en cierto sentido una interpolación entre los caminos CGN y CMD correspondientes a los métodos de Gauss-Newton y del máximo descenso respectivamente. nera simple. Sea δw el punto donde Φ̄w = w〈Φ(δ)〉+Q(δ) es mínimo. Con- secuentemente, w∇〈Φ(δw)〉 = −∇Q(δw). (4.9) Las gradientes de 〈Φ〉 y Q son entonces paralelas en δw para todos los valores de w, lo que implica que las isosuperficies de 〈Φ〉 y Q que pasan por el mínimo de Φ̄w, para cualquier w, son tangentes. La línea que forman los puntos mínimos δw correspondientes a cada valor de w, es el camino de minimización de el cual el método de Levenberg- Marquardt selecciona el vector de paso δ(s). En este caso particular, dicho camino es una curva que une los puntos δ0 = 0 (b(s) en el sistema de coordena- das original) y δ∞ = δ∗ correspondientes a w = 0 y w =∞ respectivamente. Por otro lado, los caminos de minimización de los métodos de Gauss-Newton y del máximo descenso son líneas rectas. La figura 4.1 muestra la construcción del camino de minimización CLM , correspondiente al método de Levenberg-Marquard, a partir de los puntos de tangencia entre las isosuperficies de 〈Φ〉 y Q, tal como se detalló ante- 74 riormente. Asimismo se muestran, a efectos de comparación, los caminos de minimización del método de Gauss-Newton (CGN) y del método del máximo descenso (CMD). Nótese que en el método de Levenberg-Marquardt, el pará- metro w que controla la longitud del paso, define al mismo tiempo su direc- ción. Dicha dirección coincide con la dirección de minimización del método del máximo descenso cuando w tiende a cero, y varía continuamente has- ta coincidir con la dirección de minimización del método de Gauss-Newton cuando w tiende al infinito. En la figura se puede también observar que para cualquier valor de w fini- to, el método de Levenberg-Marquardt satisface los dos objetivos propuestos. Por un lado mejora la solución δ∞ del método de Gauss-Newton estándar por medio del vector δw con menor módulo, y por otro consigue disminuir el valor inicial de 〈Φ〉, esto es, 〈Φ(δw)〉 < 〈Φ(δ0)〉. Otra cualidad que se desprende de la figura es que el módulo de δw disminuye conforme w disminuye. 4.1.2. Demostración de la Utilidad del Método Las conclusiones recién obtenidas por inspección visual carecen evidente- mente de rigor matemático. Por lo tanto, es necesario presentar las demos- traciones analíticas que garanticen que el método cumple los objetivos de minimización propuestos. Como punto de partida, nótese que en el caso trivial δw = δ0 = δ∞ = 0 (con w arbitrario), el vector solución del método se mantiene en el punto inicial, por lo que no se produce minimización alguna. Este caso trivial se presenta si y solo si los términos [c1c0], [c2c0], . . . , [ckc0] en (4.4) son todos cero. Sin embargo, dichos términos nunca son todos cero puesto que vienen a ser las componentes de la gradiente de Φ en el punto b(s), la cual se asumió no estacionaria en dicho punto. De esta manera queda excluído el único caso en el que no aplican las demostraciones dadas a continuación. Por definición se tiene w〈Φ(δw)〉+Q(δw) = Φ̄w(δw) < Φ̄w(δ∞) = w〈Φ(δ∞)〉+Q(δ∞), y además, dado que δ∞ es el mínimo de 〈Φ〉, se tiene w〈Φ(δ∞)〉+Q(δ∞) < w〈Φ(δw)〉+Q(δ∞), 75 de forma que por transitividad resulta Q(δw) < Q(δ∞). (4.10) Esta desigualdad muestra que para cualquier w, el método de Levenberg- Marquardt produce un vector de incrementos δw cuyo módulo es menor a su correspondiente obtenido por el método de Gauss-Newton estándar. Asimismo, por definición se cumple w〈Φ(δw)〉 < w〈Φ(δw)〉+Q(δw) = Φ̄w(δw) < Φ̄w(δ0) = w〈Φ(δ0)〉+Q(δ0), y puesto que δ0 = 0, se cumple además w〈Φ(δ0)〉+Q(δ0) = w〈Φ(δ0)〉, con lo cual se obtiene 〈Φ(δw)〉 < 〈Φ(δ0)〉. (4.11) Esta expresión indica que para cualquier w, el vector de incrementos δw obtenido como resultado del método de Levenberg-Marquardt disminuye el valor inicial de 〈Φ〉. Sin embargo, para garantizar su utilidad, es necesario demostrar que el método es capaz de reducir el valor que toma Φ en el punto inicial b(s) repre- sentado por el vector de incrementos δ = δ0 = 0. Con ese fin, es conveniente primero resolver las ecuaciones (4.8). Así, utilizando la técnica de los deter- minantes vista en la sección 2.2.1, resulta para la primera componente (δw)1 = −[c1c0]w1−k + · · · w−k + · · · = −[c1c0]w + · · · , (4.12) con lo cual, en la vecindad w = 0, se tiene d(δw)1 dw ∣∣∣∣∣ w=0 = −[c1c0]. (4.13) Para el resto de componentes se sigue este mismo procedimiento. Por otra parte, dada la relación definida entre un vector de parámetros b 76 y su respectivo vector de incrementos δ, es posible demostrar que dΦ dw = ∂Φ ∂b1 · d(δw)1 dw + ∂Φ ∂b2 · d(δw)2 dw + · · ·+ ∂Φ ∂bk · d(δw)k dw ; (4.14) y de la definición del símbolo de sumatoria [ ], es posible igualmente demos- trar que ∂Φ ∂b1 ∣∣∣∣∣ b=b(s) = 2[c1c0], ∂Φ ∂b2 ∣∣∣∣∣ b=b(s) = 2[c2c0], . . . ∂Φ ∂bk ∣∣∣∣∣ b=b(s) = 2[ckc0]. (4.15) Por consiguiente, sustituyendo (4.13) y (4.15) en (4.14) resulta dΦ dw ∣∣∣∣∣ w=0 = −2 ( [c1c0]2 + [c2c0]2 + . . .+ [ckc0]2 ) . (4.16) Esta derivada es negativa puesto que las derivadas parciales en (4.15) no son todas cero de acuerdo al supuesto de que Φ no es estacionaria en b(s). De este modo, Φ en función de w es decreciente en w = 0, asegurando así que existan valores de w para los que el valor de Φ se reduce respecto al valor Φ(b(s)). Existe otra forma de llegar a esta última conclusión. Por definición, la gradiente de Φ̄∆w es cero en el punto δ∆w, esto es ∆w∇〈Φ(δ∆w)〉+∇Q(δ∆w) = 0, (4.17) y dado que δ0 = 0, dicho punto puede expresarse como δ∆w = [ dδw dw (0) + (∆w) ] ∆w, (4.18) donde (∆w) es un vector cuyo límite es cero cuando ∆w tiende a cero. Por otro lado, de (4.6) se desprende ∇Q(δ∆w) = 2 δ∆w, (4.19) de modo que introduciendo (4.18) en (4.19), y este a su vez en (4.17), resulta dδw dw (0) = −1 2∇〈Φ(δ∆w)〉 − (∆w). (4.20) 77 Puesto que esta expresión es válida para cualquier valor positivo de ∆w, entonces también lo será en el límite cuando tiende a cero, con lo cual se obtiene dδw dw (0) = −1 2∇〈Φ(0)〉, (4.21) que no es otra cosa que (4.13). Si bien ya se demostró que 0 = ‖δ0‖ < ‖δw‖ < ‖δ∞‖ = ‖δ∗‖, es conveniente analizar cómo varía ‖δw‖ según se incrementa w. Con ese fin, tómese como punto de partida la siguiente expresión obtenida por definición: (w + ∆w)∇〈Φ(δw+∆w)〉+∇Q(δw+∆w) = 0, (4.22) la cual para valores infinitamente pequeños de ∆w toma la forma (w + ∆w)∇ 〈 Φ ( δw + dδw dw ∆w )〉 +∇Q ( δw + dδw dw ∆w ) = 0, que es a la vez equivalente a (w + ∆w) { ∇〈Φ(δw)〉+ [ ∇∇T 〈Φ(δw)〉 ]dδw dw ∆w } + { ∇Q(δw) + [ ∇∇TQ(δw) ]dδw dw ∆w } = 0. Introduciendo (4.6) y (4.9) en esta última ecuación, y llevándola al límite ∆w = 0, se obtiene [ w∇∇T Φ̄w(δw) ] dδw dw −∇Q(δw) = 0, que, tras aplicarle el producto interno con dδw/dw a ambos miembros, se transforma en dδw dw · ∇Q(δw) = dδw dw · [ Ā dδw dw ] , donde : Ā = w∇∇T Φ̄w(δw). (4.23) La matriz Ā es el producto entre el factor positivo w y el hessiano de Φ̄w 78 evaluado en un mínimo. Por otro lado, la segunda derivada de una función multivariable arbitraria G en un punto P y una dirección v̂ es d2G dv2 (P ) = v̂TH(P )v̂, (4.24) donde H es el hessiano de G. Por lo tanto, si en P la función G tiene un mínimo, necesariamente (4.24) es positiva para cualquier vector dirección v̂. Cabe aquí indicar que a toda matriz M que cumple la condición v̂TM v̂ > 0, para cualquier v̂, se le denomina matriz positiva definida(positive definite) [11]. Luego, como el producto de una matriz positiva definida y una cons- tante positiva sigue siendo una matriz positiva definida, se concluye que Ā es positiva definida, de manera que (4.23) implica dδw dw · ∇Q(δw) > 0. (4.25) Para determinar la variación de ‖δw‖ respecto a w, basta con analizar el comportamiento de su derivada a lo largo de su dominio. Una manera de simplicar el cálculo de la misma, consiste en expresarla en función de la derivada de ‖δw‖2, así se tiene d‖δw‖2 dw = 2 ‖δw‖ d‖δw‖ dw , (4.26) y puesto que ‖δw‖2 = δw · δw, (4.27) resulta d‖δw‖2 dw = 2 dδw dw · δw. (4.28) Introduciendo (4.28) en (4.26) se obtiene la expresión buscada d‖δw‖ dw = dδw dw · δ̂w, (4.29) 79 la cual, tras considerar que ∇Q(δw) es paralela a δ̂w, toma la forma d‖δw‖ dw = ‖∇Q(δw)‖−1 dδw dw · ∇Q(δw). (4.29a) De este resultado y de (4.25) se deduce d‖δw‖ dw > 0, (4.30) para todo w positivo. Con esto queda finalmente demostrado que el módulo del vector δw es estrictamente creciente en w. 4.1.3. Cálculo del Vector de Paso Que el método logre reducir el valor de Φ en cada iteración depende directamente del valor de w que se elija. Por esta razón es importante analizar qué valores de w son los adecuados. Teóricamente, el mejor valor de w a elegir es aquel que satisface dΦ dw (bw) = 0 donde : bw = b(s) + δw; (4.31) sin embargo, resolver esta ecuación en la práctica es un problema complejo. Una estrategia, utilizada por Cauchy [3] para abordar este tipo de pro- blemas, consiste en plantear la aproximación Φ(bw) ≈ Φ(b(s)) + w ( dΦ dw )∣∣∣∣∣ w=0 . (4.32) Bajo el supuesto de que b(s) es tal que origina que el valor de Φ(bw) sea pequeño, el miembro de la izquierda en (4.32) puede considerarse igual a cero, con lo cual resulta la fórmula w ∼= − Φ(b(s)) (dΦ/dw) w=0 = 1 2 Φ(b(s)) [c1c0]2 + [c2c0]2 + . . .+ [ckc0]2 , (4.33) que constituye una buena opción para w bajo las condiciones impuestas. 80 4.2. Enfoque de Marquardt Posterior a Levenberg y de manera independiente, Marquardt concibió prácticamente el mismo método propuesto por Levenberg, aunque siguiendo una línea de razonamiento distinta. El desarrollo de esta sección es hasta cierto punto una reproducción del trabajo original de Marquardt [21]. Cabe señalar que la notación adoptada en este y en los dos capítulos anteriores fue planeada con el fin de lograr la mayor compatibilidad posible con la notación utilizada por Marquardt. 4.2.1. Base Teórica del Método Los siguientes tres teoremas constituyen el fundamento teórico en el que se basa el método. En ellos se busca hacer énfasis en las relaciones geométricas involucradas. Teorema 1. Si λ es un real positivo arbitrario y δ0 satisface la ecuación (A+ λI)δ = g, (4.34) entonces δ0 minimiza a 〈Φ〉 en la esfera cuyo radio ‖δ‖ satisface ‖δ‖2 = ‖δ0‖2. (4.35) Demostración. Antes de abordar la demostración de este caso específico es conveniente analizar el problema más general. Este consiste en determinar el punto que minimiza la función G(δ) sobre la isosuperficie H(δ) = 0, donde G(δ) : Rk → R y H(δ) : Rk → R son doblemente diferenciables. Condición necesaria para el punto mínimo buscado es que la gradiente de G en dicho punto sea ortogonal a la isosuperficie de H en él o, expresado de otra manera, que sea paralela a la gradiente de H en tal punto. Por lo tanto, el punto mínimo debe ser solución de ∇G(δ) + λ∇H(δ) = 0 y H(δ) = 0, donde el valor particular del parámetro λ debe seleccionarse de manera que implique el cumplimiento de la segunda ecuación de restricción. Para calcu- lar tal valor se obtiene δ en función de λ a partir de la primera ecuación. 81 Este resultado se introduce en la ecuación de restricción que toma la forma H(δ(λ)) = 0, de donde finalmente puede despejarse λ. A este procedimiento de minimización con restricción se le denomina mé- todo de los multiplicadores de Lagrange y al parámetro λ, multiplicador de Lagrange. Aplicando este método al caso específico en el que la función a minimizar es 〈Φ(δ)〉 = ‖y − f 0 − Pδ‖2, (4.36) bajo la restricción ‖δ‖2 = ‖δ0‖2, (4.37) se obtiene para el punto mínimo la condición necesaria ∂Φ̃ ∂δ1 = ∂Φ̃ ∂δ2 = · · · = ∂Φ̃ ∂δk = 0, ∂Φ̃ ∂λ = 0, (4.38) donde Φ̃(δ, λ) = ‖y − f 0 − Pδ‖2 + λ(‖δ‖2 − ‖δ0‖2). (4.39) Tomando las derivadas indicadas se tiene −[P T (y − f 0)− P TPδ] + λδ = 0 (4.40) y ‖δ‖2 − ‖δ0‖2 = 0. (4.41) Si se denota la solución de (4.40) como δ0, entonces se cumple (P TP + λI)δ0 = P T (y − f 0), (4.42) que es equivalente a (A+ λI)δ0 = g. (4.43) Demostrar que δ0 cumple (4.41) es trivial. Finalmente, el que el punto estacionario δ0 sea efectivamente un mínimo se concluye del hecho que A es positiva definida, por ser proporcional al hessiano de 〈Φ〉, y λ≥0 (véase pág. 78). Teorema 2. Si δ(λ) es la solución de (4.34) para un valor de λ dado, en- 82 tonces ‖δ(λ)‖2 es una función monótona decreciente continua de λ, tal que cuando λ→∞, ‖δ(λ)‖2 → 0. Demostración. Siempre que los eigenvectores de una matriz A sean lineal- mente independientes, existe una transformación de similaridad (similarity transformation) que reduce A a una matriz diagonal D. Esto es, existe una matriz no singular M cuyas columnas son precisamente los eigenvectores de A, de manera que M−1AM = D, donde los elementos de D son los eigenvalores de A [30]. Puesto que en el presente caso A es simétrica, sus eigenvalores son distintos entre sí, originando que la matriz M sea ortonormal y represente una rotación de coordenadas, con lo cual MTM = I. Más aún, al ser A positiva definida, se cumple que todos los elementos de D son positivos. Premultiplicando (4.34) por M−1 resulta M−1(A+ λI)MM−1δ = M−1g, de modo que δ = M(D + λI)−1M−1g. (4.44) Luego, definiendo v = M−1g, se tiene ‖δ(λ)‖2 = vT [(D + λI)−1]TMTM(D + λI)−1v = vT [(D + λI)2]−1v = k∑ j=1 v2 j (Dj + λ)2 , (4.45) la cual claramente es una función decreciente en λ, (con λ ≥ 0), tal que cuando λ→∞, ‖δ(λ)‖2 → 0. La transformación ortonormal hacia una matriz diagonal ha sido explíci- tamente mostrada con el fin de facilitar la prueba del siguiente teorema. Teorema 3. Si γ es el ángulo entre δ(λ) y −∇Φ, entonces γ es una función monótona decreciente continua de λ, tal que cuando λ→∞, γ → 0. Puesto que −∇Φ es independiente de λ, significa que δ(λ) rota hacia −∇Φ conforme λ→∞. 83 Demostración. Obsérvese que −∇Φ = 2g, de manera que γ es igualmente el ángulo que forman δ(λ) y g. Por definición se tiene cos γ = δ(λ)Tg (‖δ‖)(‖g‖) = vT (D + λI)−1v (vT [(D + λI)2]−1v)1/2(gTg)1/2 = k∑ j=1 v2 j (Dj + λ) k∑ j=1 v2 j (Dj + λ)2 1/2 (gTg)1/2 . (4.46) Derivando y simplificando, resulta d dλ cos γ =  k∑ j=1 v2 j (Dj + λ)  k∑ j=1 v2 j (Dj + λ)3 −  k∑ j=1 v2 j (Dj + λ)2 2  k∑ j=1 v2 j (Dj + λ)2 3/2 (gTg)1/2 (4.47) =  k∑ j=1 v2 j ∏ 1j  k∑ j=1 v2 j ∏ 3j −  k∑ j=1 v2 j ∏ 2j 2  k∑ j=1 v2 j (Dj + λ)2 3/2  k∏ j=1 (Dj + λ)2 2 (gTg)1/2 , (4.48) donde ∏ 1j = k∏ h=1 h6=j (Dh + λ), ∏ 2j = k∏ h=1 h6=j (Dh + λ)2 ∏ 3j = k∏ h=1 h6=j (Dh + λ)3. Puesto que el denominador en (4.48) está compuesto por factores positivos, el signo de d(cos γ)/dλ es entonces el signo del numerador. Notando que ∏ 1j ∏ 3j = ( ∏ 2j)2, 84 el numerador se puede escribir como  k∑ j=1 (vj ∏ 1/2 1j )2  k∑ j=1 (vj ∏ 1/2 3j )2 −  k∑ j=1 (vj ∏ 1/2 1j )(vj ∏ 1/2 3j ) 2 . (4.49) De la desigualdad de Schwarz se desprende que (4.49) es mayor que cero, por lo tanto d(cos γ)/dλ es siempre positiva para cualquier λ > 0. De este resultado se concluye que γ es una función monótona decreciente de λ. Por otro lado, para valores muy grandes de λ, la matriz (A + λI) es dominada por la diagonal λI. Por lo tanto de (4.34) se observa que cuando λ → ∞, δ(λ) → g/λ, de manera que en el límite el ángulo entre δ y g tiende a cero. Para el caso en que λ = 0 en (4.34), los vectores δ y g forman un ángulo 0 < γ < π/2, excepto en el caso trivial cuando A es la matriz identidad. 4.2.2. Escalamiento del Espacio de Parámetros Un aspecto numérico importante a considerar en el procedimiento de solución de sistemas de ecuaciones lineales del tipo A δ = g es el grado de sensibilidad que poseen las matrices A y g respecto a perturbaciones en sus elementos. En ese sentido se dice que A y g son muy sensibles, o están pobrememnte acondicionadas, cuando una pequeña perturbación en ellas produce una gran desviación en el vector solución del sistema δ. Con el fin de lograr que el vector de parámetros δ obtenido en cada iteración posea un buen nivel de inmunidad a las posibles perturbaciones producidas por errores presentes en las matrices del sistema A y g, y puesto que el método es invariante respecto a transformaciones lineales de coorde- nadas, es conveniente llevar a cabo un escalamiento adecuado del espacio de parámetros b. En particular, escalando cada componente bj de b en unidades de la des- viación estándar de las derivadas ∂fi/∂bj tomadas sobre todas las muestras (i = 1, 2, . . . , n), se obtiene una transformación que convierte a la matriz A en la matriz de coeficientes de correlación normalizados de los ∂fi/∂bj. La ventaja de este escalamiento es que minimiza la sensibilidad de A y g, consiguiendo así un acondicionamiento óptimo del sistema de ecuaciones. En el nuevo espacio de parámetros b? se tiene que la matriz escalada A? 85 y el vector escalado g? poseen la forma A? = (a?jj′) = ajj′√ ajj √ aj′j′ ) , (4.50) g? = (g?j ) = gj√ ajj ) , (4.51) de manera que el vector de paso δ? según el método de Gauss-Newton en el nuevo espacio de parámetros se obtiene resolviendo A? δ? = g?. (4.52) El correspondiente vector de paso δ en el espacio de parámetros original se obtiene a partir de su versión δ? mediante la relación δj = δ?j√ ajj . (4.53) 4.2.3. Construcción del Algoritmo El bosquejo general del algoritmo para implementar el método queda finalmente delineado. Específicamente, para la s-ésima iteración se construye la ecuación (A?(s) + λ(s)I) δ?(s) = g?(s), (4.54) la cual equivale a minimizar 〈Φ〉 en el espacio de parámetros original, pero esta vez no sobre la esfera dada en (4.35), sino sobre una elipse de la forma k∑ j=1 ajjδ 2 j = cte. (4.55) A partir del δ(s) obtenido al introducir en (4.53) la solución δ?(s) de la ecua- ción (4.54), se genera el nuevo punto tentativo b(s+1) = b(s) + δ(s), (4.56) 86 el cual conduce a un nuevo valor de Φ(s+1). Cabe mencionar que es esencial seleccionar un valor de λ(s) que conduzca al resultado Φ(s+1) < Φ(s). (4.57) De la teoría precedente se puede afirmar que siempre existirá un λ(s) suficientemente grande como para satisfacer (4.57), a menos que b(s) sea ya un mínimo de Φ. Para determinar el valor de λ(s) que satisfaga (4.57), y que a la vez produzca una rápida convergencia del algoritmo, se requiere un tanteo del tipo ensayo error. En cada iteración se desea minimizar Φ dentro de la máxima vecindad posible sobre la cual la versión linealizada del modelo ofrezca una adecuada representación de la función no lineal. Por consiguiente, la estrategia pa- ra elegir λ(s) debe buscar utilizar un valor pequeño de λ(s) siempre que se den la condiciones que permitan una buena convergencia del método están- dar de Gauss-Newton. Esto es especialmente relevante en los últimos pasos del procedimiento de convergencia, cuando los puntos tentativos están en la vecindad inmediata del mínimo bmin, donde los contornos de Φ son asintó- ticamente elípticos, y la expansión lineal del modelo necesita ser una buena aproximación solo sobre una pequeña región. Por otro lado, valores grandes de λ(s) solo deben ser utilizados cuando sean necesarios para satisfacer (4.57). Si bien es cierto que Φ(s+1) como función de λ posee un mínimo, y que la elección de dicho λ en la (s + 1)− ésima iteración maximiza (Φ(s)−Φ(s+1)), tal elección óptima localmente resulta ser una estrategia pobre globalmente, ya que esta por lo general requiere un valor substancialmente mayor de λ del que es necesario para satisfacer (4.57). Esta estrategia hereda por ende la propiedad del método del máximo descenso de tener un progreso inicial rápido seguido de uno progresivamente más lento. En consecuencia, la estrategia se define de la siguiente manera: 1. Definir el factor de incremento/decremento ν > 1. 2. Asignar un valor inicial a λ; por ejemplo: λ(0) = 10−2. 3. Para la (s+ 1)− ésima iteración: a) Calcular Φ(λ(s)) y Φ(λ(s)/ν). 87 b) Si Φ(λ(s)/ν) ≤ Φ(s), tomar λ(s+1) = λ(s)/ν. c) Si Φ(λ(s)/ν) > Φ(s) y Φ(λ(s)) ≤ Φ(s), tomar λ(s+1) = λ(s). d) Si Φ(λ(s)/ν) > Φ(s) y Φ(λ(s)) > Φ(s), incrementar λ multiplicán- dolo sucesivamente por ν hasta que para el menor u se cumpla: Φ(λ(s)νu) ≤ Φ(s). Tomar λ(s+1) = λ(s)νu. En relación al punto d), es necesario hacer algunos comentarios complemen- tarios. En ocasiones, en problemas en los que la correlación entre los sucesivos puntos tentativos b(s) es extremadamente alta (> 0,99), puede ocurrir que el valor de λ se incremente hasta valores inaceptablemente altos. Para tales casos es preferible modificar la prueba en d) como se detalla a continuación. Sea b(s+1) = b(s) +K(s)δ(s), K(s) ≤ 1. (4.58) Notando que el ángulo γ(s) (definido en la pág 82) es función decreciente de λ(s), se selecciona un ángulo umbral γ0 < π/2 y se fija K(s) = 1, si γ(s) ≥ γ0. (4.59) Si no es posible superar la prueba d) a pesar de haber incrementado λ(s) hasta haber conseguido que se cumpla γ(s) < γ0, entonces se deja de incrementar λ(s) y se toma un K(s) lo suficientemente pequeño como para que se cumpla Φ(s+1) < Φ(s). Esto es siempre posible puesto que γ(s) < γ0 < π/2. Una buena elección para el ángulo umbral es γ0 = π/4. Nótese que el hecho que cos γ sea positivo cuando λ = 0 solo se puede garantizar cuando la matriz A? es positiva definida. En presencia de niveles de correlación muy grandes, el grado en que A? es positiva definida es tenue y puede degradarse notoriamente debido a errores de redondeo presentes en las operaciones de punto flotante realizadas en los cálculos por computadora. En tales casos, el método puede llegar a divergir indistintamente del valor de K(s) utilizado. Un beneficio asociado al hecho de que el método le sume λ a la diagonal de A?, es que se garantiza que cos γ sea positivo aun en el caso que A? sea débilmente positiva definida. A pesar de que no siempre es posible prever la ocurrencia de altas co- rrelaciones entre los parámetros tentativos b(s) relacionados a modelos no 88 lineales, con frecuencia sí es posible reducir substancialmente dichos niveles de correlación mediante un mejor diseño experimental en la producción de los datos medidos. Numéricamente, se considera que el proceso de iteración ha convergido cuando se cumple la condición |δ(s) j | τ + |b(s) j | < ε, j = 1, 2, . . . , k (4.60) para algún ε > 0 suficientemente pequeño (como 10−5) y algún τ adecuado (como 10−3). Si bien la elección de ν es arbitraria, se ha encontrado en la práctica que ν = 10 es una buena elección. 4.3. Similitud y Diferencia entre los Enfoques de Levenberg y Marquardt Si bien ambos enfoques fueron concebidos mediante razonamientos dis- tintos, los dos son versiones de un mismo método, lo cual se evidencia al considerar la equivalencia w−1 = λ entre los parámetros de magnificación- atenuación ‘w’ y ‘λ’ utlizados por Levenberg y Marquardt respectivamente. Asimismo, la diferencia entre ambos enfoques radica en las estrategias distintas utilizadas en la determinación de los parámetros w y λ, con el fin de minimizar Φ en cada iteración. Esto produce que la razón de convergencia no sea la misma para cada caso. En el enfoque de Levenberg, la estrategia consiste en minimizar Φ local- mente como una función del parámetro w−1 que se le suma a la diagonal, con lo cual se consigue una razón de convergencia lineal ‖b(s+1) − bmin‖ ≤ C‖b(s) − bmin‖, C < 1, (4.61) como la que presenta el método del máximo descenso. Por otro lado, en el enfoque de Marquardt, la estrategia consiste en mini- mizar Φ en la mayor vecindad posible sobre la cual la aproximaxión lineal sea una representación adecuada. Esta estrategia busca por ende el menor valor de λ que satisfaga (4.57), con lo cual se obtiene una razón de convergencia 89 cuadrática ‖b(s+1) − bmin‖ ≤ C‖b(s) − bmin‖2, C < 1, (4.62) similar a la que posee el método de Gauss-Newton. Capítulo 5 Estimación de Parámetros y Conclusiones El presente capítulo muestra los resultados y conclusiones de la aplicación del método de Levenberg–Marquardt al caso específico del modelo TGM introducido en el primer capítulo para la estimación de sus parámetros de mejor ajuste según el criterio de los mínimos cuadrados. 5.1. Procedimiento de Estimación de Parámetros El procedimiento de estimación de parámetros se lleva a cabo en dos etapas. En la primera etapa se ajusta la impedancia mecánica del modelo a los datos experimentales. El algoritmo utilizado para ello es una versión optimi- zada del método de Levenberg-Marquardt [23], [8], analizado en el capítulo anterior. De (1.5) se puede reconocer que los parámetros r2 y c2 introducidos por el modelo TGM tienen una mayor influencia en la parte real de la impedancia mecánica a bajas frecuencias. Tomando en cuenta este hecho, los parámetros r1, r2, y c2 son estimados ajustando el modelo de la parte real de la impe- dancia mecánica a los datos medidos en el rango que va de los 5 Hz hasta dos veces la frecuencia de resonancia. Los otros dos parámetrosm y c1 se obtienen ajustando el modelo de la parte imaginaria de la impedancia mecánica a los 90 91 datos medidos en el intervalo de frecuencias cuyos límites resultan de dividir y multiplicar la frecuencia de resonancia por un factor de 1.3. Este angosto rango de frecuencias alrededor de la frecuencia de resonancia se encontró que era el óptimo para determinar estos parámetros, de forma tal que el punto de cruce por cero del modelo de la parte imaginaria de la impedancia mecánica coinciden con la frecuencia de resonancia real. En la segunda etapa se refinan los parámetros obtenidos en la etapa an- terior con el fin de conseguir que el modelo coincida con los datos medidos en ciertos puntos claves de la curva de magnitud de la mobilidad mecánica (definida como la inversa de la magnitud de la impedancia mecánica), tales como el punto de resonancia o de máxima amplitud y los puntos que están 3 dB por debajo de este. Para este propósito se define la función Fi = 1√√√√ r1 + r2 1 + (ωic2r2)2 )2 + (ωim− 1 ωic1 )− ωic2r 2 2 1 + (ωic2r2)2 )2 − xi, (5.1) donde el par (ωi, xi) compuesto por la magnitud de la mobilidad xi medida a la frecuencia ωi representa un punto de coincidencia. Los parámetros para lograr la coincidencia son entonces ceros de la función Fi. Ya que se tienen tres puntos de coincidencia y cinco parámetros, calcular estos últimos consiste en resolver un sistema de ecuaciones subdeterminado y en general no lineal. 90 95 100 105 110 0.8 0.85 0.9 0.95 1 1.05 1.1 Hz m /s /N Mechanical Mobility, Magnitude measured standard elastic model visco-elastic model 90 95 100 105 110 0.8 0.85 0.9 0.95 1 1.05 1.1 Hz m /s /N Mechanical Mobility, Magnitude measured standard elastic model visco-elastic model Figura 5.1: Mobilidad mecánica medida y simulada alrededor de la frecuencia de reso- nancia. Izquierda: sin proceso de refinamiento de parámetros. Derecha: con parámetros refinados 92 Para reducir el número de variables es conveniente excluir aquellas cuyas variaciones producen un menor efecto en (5.1). Tras una serie de pruebas, se determinó que la mayor sensibilidad se presentaba ante cambios en los parámetros m, c1 y c2. Al final se determinan los ceros mediante un algoritmo basado en el mé- todo de Newton–Raphson. La mejora alcanzada en esta etapa se muestra en la figura 5.1. La implementación práctica del método descrito en este trabajo se llevó a cabo mediante un programa especialmente desarrollado para tal fin. El pro- grama controla el analizador de señales, automatiza el proceso de medición, procesa los datos medidos, realiza el ajuste según el Método de Levenberg– Marquardt, y muestra los parámetros resultantes. En la figura 5.2 se pueden ver los resultados medidos y simulados de un altavoz JBL 2206, en los que se observa una muy buena concordancia desde las bajas frecuencias hasta una frecuencia de aproximadamente 1,3 veces la frecuencia de resonancia. Esto demuestra que el modelo TGM es adecuado. Los resultados, tal y como los muestra el programa, pueden verse en la figura 5.3. Nótese que el recuadro superior corresponde a los parámetros Thiele–Small definidos por el modelo tradicional (también calculados por el programa), mientras que el inferior corresponde a los parámetros definidos por el modelo TGM. 0 50 100 150 200 250 0 0.05 0.1 0.15 0.2 Mechanical Mobility, Magnitude Hz m /s /N 0 50 100 150 200 250 -100 -50 0 50 100 Mechanical Mobility, Phase Hz D e g r measured modeled measured modeled 0 50 100 150 200 250 0 10 20 30 40 Hz K g /s Mechanical Impedance, Real Part 0 50 100 150 200 250 -300 -200 -100 0 100 200 Hz K g /s Mechanical Impedance , Imaginary Part measured modeled measured modeled Figura 5.2: Mobilidad mecánica medida y simulada de un altavoz JBL 2206. Izquierda: magnitud y fase. Derecha: parte real y parte imaginaria. 93 Figura 5.3: Parámetros estimados de un altavoz JBL 2206. 5.2. Conclusiones A manera de resumen y epílogo, se pueden mencionar las siguientes con- clusiones: - El modelo TGM de la suspensión mecánica de un altavoz es una ex- tensión del modelo tradicional que describe el incremento que exhibe la resistencia mecánica de un altavoz conforme disminuye la frecuencia en el rango de las frecuencias bajas, logrando con ello un mejor ajuste a los datos medidos que el que se obtiene con la resistencia mecánica constante del modelo tradicional. - El criterio de los mínimos cuadrados cuantifica el grado de discrepancia entre el modelo y los n datos medidos mediante el cuadrado de la norma del vector de residuos (diferencia entre los datos modelados y medidos) tal como esta se define en el espacio euclideano Rn. Bajo ciertos supuestos es posible demostrar que este criterio es el óptimo cuando los datos medidos están contaminados con errores de medición. Existen varios métodos para la estimación de parámetros según este criterio. - El método de Gauss–Newton realiza una aproximación de primer orden de los residuos y estima iterativamente los parámetros que minimizan la sumatoria de los cuadrados de los residuos linealizados. La ventaja de este método es su rápida convergencia una vez alcanzada la vecindad 94 de los parámetros de ajuste óptimo y su desventaja es la divergencia que presenta cuando los parámetros iniciales están fuera de la región de convergencia. - El método del Máximo Descenso utiliza la dirección de la gradiente en el espacio de parámetros para aproximar iterativamente los parámetros de ajuste óptimo. La ventaja de este método es su convergencia segura aun cuando los parámetros iniciales estén fuera de la región de convergencia y su desventaja es la lenta convergencia que presenta una vez alcanzada la vecindad de los parámetros de ajuste óptimo. - El método de Levenberg–Marquardt realiza una interpolación óptima entre el método de Gauss–Newton y el método del Máximo Descenso. Esta interpolación está basada en la máxima vecindad en la cual la aproximación lineal introducida en el método de Gauss–Newton cons- tituye una buena representación del problema no lineal original. - La aplicación del método de Levenberg–Marquardt al problema parti- cular de la estimación de parámetros mecánicos según el modelo TGM demostró ser adecuada por el buen ajuste obtenido (figura 5.2) y efi- ciente por el tiempo, en el orden de los milisegundos, necesario para llevar a cabo la estimación. - Tal como se puede observar en la figura 1.5, el incremento que exhibe la resistencia mecánica del altavoz conforme aumenta la frecuencia, en el rango de las frecuencias altas, se debe a la contribución de la impedancia acústica generada por las cavidades en el circuito magnético del altavoz. El desarrollo de una extensión del modelo que reproduzca dicho efecto es materia de posteriores investigaciones. Bibliografía [1] L. L. Beranek, “Acoustics,” American Institute of Physics, New York, 1986, p. 185. [2] R. R. Brown, J. B. Dennis, y C. Kingsley, “Design of Systems Using Digital Computers,” WADC Tech. Note 56-3383, Servo. Lab., M.I.T., 1956. [3] A. L. Cauchy, “Méthode Générale pour la Résolution des Systèmes D’Équations Simultanées,” Comptes Rendus Acad. Sci. Paris, vol. 25, pp. 536–538, 1847. [4] R. Courant y D. Hilbert, “Methoden der mathematischen Physik II,” 2da ed., Springer-Verlag, Berlin, 1968, p. 13. [5] H. B. Curry, “The Method of Steepest Descent for Non-Linear Minimiza- tion Problems,” Quarterly of Applied Mathematics, vol. 2, pp. 258–261, 1944. [6] B. J. Elliot, “On the Measurement of Low-Frequency Parameters of Moving-Coil Piston Transducers,” presentado en la AES 58th conven- tion, 1977, preprint 1299. [7] R. L. Ellis, “On the Method of Least Squares,” Trans. Camb. Phil. Soc., vol. 8, pp. 204–219, 1894. [8] R. Fletcher, “A Modified Marquardt Subroutine for Non-Linear Least Squares,” Report AERE - R.6799, U.K. Atomic Energy Authority, 1971. [9] C. F. Gauss, “Theorie der Bewegung der Himmelskörper welche in Ke- gelschnitten die Sonne Umlaufen,” Carl Meyer, Hannover, 1865. 95 96 [10] C. F. Gauss, “Theorie der den kleinsten Fehlern unterworfenen Combi- nation der Beobachtungen,” en A. Börsch y P. Simon (eds.), “Abhand- lungen zur Methode der kleinsten Quadrate”, Stankiewicz, Berlin, 1887, pp. 1–53. [11] G. H. Golub y C. F. Van Loan, “Matrix Computations,” 3ra ed., The Johns Hopkins University Press, Baltimore, 1996, p. 140. [12] H. O. Hartley, “The Modified Gauss-Newton Method for the Fitting of Non-Linear Regression Functions by Least Squares,” Technometrics, vol. 3, pp. 269–280, 1961. [13] S. Jønsson, J. N. Moreno, y H. Bøg “Measurement of Closed Box Louds- peaker System Parameters Using a Laser Velocity Transducer and an Audio Analyzer,” presentado en la AES 92nd convention, 1992, preprint 3325. [14] M. H. Knudsen y J. G. Jensen, “Low-Frequency Loudspeaker Models that Include Suspension Creep,” J. Audio Eng. Soc., vol. 41, pp.3–18, 1993. [15] A. N. Kolmogorov, “Justification of the Method of Least Squares,” en A. N. Shiryayev (ed.), “Selected Works of A. N. Kolmogorov”, Kluwer Academis Publishers, Dordrecht, 1992, vol. II, pp. 285–302. [16] A. N. Kolmogorov, “Foundations of the Theory of Probability,” 2da ed., Chelsea Publishing Company, New York, 1956, p. 6. [17] P. S. Laplace, “Théorie Analytique des Probabilités,” Courcier, Paris, 1812, livre II, chap. IV. [18] W. M. Leach, R. W. Schafer, y T. P. Barnwell, “Time-Domain Measure- ments of Loudspeaker Driver Parameter,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-1 27, pp. 734–739, 1979. [19] A. M. Legendre, “Sur la Méthode des Moindres Quarrés,” en A. M. Legendre, “Nouvelles Méthodes pour la Détermination des Orbites des Comètes”, Didot, Paris, 1805, pp. 72–80. 97 [20] K. Levenberg, “A Method for the Solution of Certain Non-Linear Pro- blems in Least Squares,” Quarterly of Applied Mathematics, vol. 2, pp. 164–168, 1944. [21] D. W. Marquardt, “An Algorithm for Least-Squares Estimation of Non- linear Parameters,” J. Soc. Indust. Appl. Math., vol. 11, pp. 431–441, 1963. [22] D. W. Marquardt, “Solution of Nonlinear Chemical Engineering Mo- dels,” Chemical Engineering Progress, vol. 55, pp. 65–70, 1959. [23] J. Moré, “Levenberg-Marquardt Algorithm: Implementation and Theory,” en G. A. Watson (ed.), “Numerical Analysis”, Springer-Verlag, Berlin, 1978. [24] J. N. Moreno, “Measurement of Loudspeaker Parameters Using a Laser Velocity Transducer and Two-Channel FFT Analysis,” J. Audio Eng. Soc., vol. 39, pp. 243–249, 1991. [25] J. N. Moreno y V. R. Medina, “Measurement of Loudspeaker Parameters Considering a Better Fitting for the Mechanical Impedance,” presentado en la AES 51st international conference, Helsinki, 2013. [26] P. M. Morse y H. Feshbach, “Methods of Theoretical Physics,” McGraw- Hill, New York, 1953, p. 494. [27] A. Pringsheim, “Zur Geschichte des Taylorschen Lehrsatzes,” Bibliothe- ca Mathematica, ser. 3, vol. 1, pp. 433–479, 1900. [28] R. H. Small, “Simplified Loudspeaker Measurements at Low Frequen- cies,” J. Audio Eng. Soc., vol. 20, pp. 28–33, 1972. [29] A. N. Thiele, “Measurement of the Thiele/Small Parameters of Twee- ters,” J. Elec. and Electron. Eng., Australia, vol. 9, pp. 186–199, 1989. [30] J. H. Wilkinson, “Algebraic Eigenvalue Problem,” Oxford University Press, Oxford, 1965, p. 2.