Mínimos cuadrados no lineales

Los Mínimos cuadrados no lineales es la forma de análisis de mínimos cuadrados que se usa para encajar un conjunto de m observaciones con un modelo que es no lineal en n parámetros desconocidos (m > n). Se utiliza en algunas formas de regresión no lineal. La base del método es para aproximar el modelo por uno lineal y para refinar los parámetros por iteraciones sucesivas. Hay muchas similitudes con mínimos cuadrados lineales, pero también algunas diferencias importantes.

Teoría

Considere un conjunto de ${\displaystyle m}$ observaciones, ${\displaystyle (x_{1},y_{1}),(x_{2},y_{2}),\dots ,(x_{m},y_{m}),}$ y una curva (función del modelo) ${\displaystyle y=f(x,{\boldsymbol {\beta }}),}$ que además de la variable ${\displaystyle x}$ también depende de ${\displaystyle n}$ parámetros, ${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\beta _{2},\dots ,\beta _{n}),}$ con ${\displaystyle m\geq n.}$ Se desea encontrar el vector ${\displaystyle {\boldsymbol {\beta }}}$ de parámetros tales que la curva se ajuste mejor a los datos dados en el sentido de mínimos cuadrados, es decir, la suma de cuadrados

${\displaystyle S=\sum _{i=1}^{m}r_{i}^{2}}$

esta es minimizada cuando los errores r_i están dados por

${\displaystyle r_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }})}$

para ${\displaystyle i=1,2,\dots ,m.}$

El mínimo valor de S se produce cuando el gradiente es cero. Dado que el modelo contiene n parámetros hay n ecuaciones de gradiente:

${\displaystyle {\frac {\partial S}{\partial \beta _{j}}}=2\sum _{i}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}=0\quad (j=1,\ldots ,n).}$

En un sistema no lineal, los derivados ${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}$ son funciones tanto de la variable independiente y los parámetros, por lo que estas ecuaciones gradiente no tienen una solución cerrada. En lugar de ello, los valores iniciales deben ser elegidos para los parámetros. Entonces, los parámetros se refinan iterativamente, es decir, los valores se obtienen por aproximación sucesiva,

${\displaystyle \beta _{j}\approx \beta _{j}^{k+1}=\beta _{j}^{k}+\Delta \beta _{j}.\,}$

Aquí, k es un número de iteración y el vector de incrementos, ${\displaystyle \Delta {\boldsymbol {\beta }}\,}$ que se conoce como el vector de desplazamiento. En cada iteración del modelo se linealiza por aproximación a un primer orden en serie de Taylor de expansión sobre ${\displaystyle {\boldsymbol {\beta }}^{k}\!}$

${\displaystyle f(x_{i},{\boldsymbol {\beta }})\approx f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}{\frac {\partial f(x_{i},{\boldsymbol {\beta }}^{k})}{\partial \beta _{j}}}\left(\beta _{j}-\beta _{j}^{k}\right)\approx f(x_{i},{\boldsymbol {\beta }}^{k})+\sum _{j}J_{ij}\,\Delta \beta _{j}.}$

El jacobiano , J, es una función de las constantes, la variable independiente y los parámetros, por lo que cambia de una iteración a la siguiente. Por lo tanto, en términos del modelo linealizado, ${\displaystyle {\frac {\partial r_{i}}{\partial \beta _{j}}}=-J_{ij}}$ y los residuos se dan por

${\displaystyle r_{i}=\Delta y_{i}-\sum _{s=1}^{n}J_{is}\ \Delta \beta _{s};\ \Delta y_{i}=y_{i}-f(x_{i},{\boldsymbol {\beta }}^{k}).}$

Sustituyendo estas expresiones en las ecuaciones de gradiente, se convierten

${\displaystyle -2\sum _{i=1}^{m}J_{ij}\left(\Delta y_{i}-\sum _{s=1}^{n}J_{is}\ \Delta \beta _{s}\right)=0}$

que, en el reordenamiento, convertido en n ecuaciones lineales simultáneas, las ecuaciones normales

${\displaystyle \sum _{i=1}^{m}\sum _{s=1}^{n}J_{ij}J_{is}\ \Delta \beta _{s}=\sum _{i=1}^{m}J_{ij}\ \Delta y_{i}\qquad (j=1,\dots ,n).\,}$

Las ecuaciones normales se escriben en notación matricial como

${\displaystyle \mathbf {\left(J^{T}J\right)\Delta {\boldsymbol {\beta }}=J^{T}\ \Delta y} .}$

Cuando las observaciones no son igualmente fiable, una suma ponderada de los cuadrados puede ser minimizado,

${\displaystyle S=\sum _{i=1}^{m}W_{ii}r_{i}^{2}.}$

Cada elemento de la matriz de peso diagonal W debería, idealmente, ser igual al recíproco de la varianza de error de la medida.^[1]​ Las ecuaciones normales son entonces:

${\displaystyle \mathbf {\left(J^{T}WJ\right)\Delta {\boldsymbol {\beta }}=J^{T}W\ \Delta y} .}$

Estas ecuaciones forman la base para el algoritmo de Gauss-Newton para un problema de mínimos cuadrados no lineal.

Interpretación geométrica

En mínimos cuadrados lineales la función objetivo, S, es una función cuadrática de los parámetros.

${\displaystyle S=\sum _{i}W_{ii}\left(y_{i}-\sum _{j}X_{ij}\beta _{j}\right)^{2}}$

Cuando sólo hay un parámetro, la gráfica de S con respecto a ese parámetro será una parábola. Con dos o más parámetros, los contornos de S con respecto a cualquier par de parámetros serán elipses concéntricas (suponiendo que la matriz de ecuaciones normales ${\displaystyle \mathbf {X^{T}WX} }$ es definida positiva). Los valores de los parámetros mínimos se encuentran en el centro de las elipses. La geometría de la función objetivo general puede describirse como el elíptico paraboloide. En NLLSQ la función objetivo es cuadrática con respecto a los parámetros sólo en una región cercana a su valor mínimo, donde la serie truncada de Taylor es una buena aproximación al modelo.

${\displaystyle S\approx \sum _{i}W_{ii}\left(y_{i}-\sum _{j}J_{ij}\beta _{j}\right)^{2}}$

Cuanto más los valores de los parámetros difieren de sus valores óptimos, más los contornos se desvían de la forma elíptica. Una consecuencia de esto es que las estimaciones de parámetros iniciales deben ser lo más cercanas posible a sus valores óptimos (desconocidos!). También explica cómo la divergencia puede surgir como el algoritmo de Gauss-Newton es convergente sólo cuando la función objetivo es aproximadamente cuadrática en los parámetros.

Computación

Estimaciones de parámetros iniciales

Algunos problemas de mal acondicionamiento y divergencia se pueden corregir mediante la búsqueda de estimaciones de parámetros iniciales que están cerca de los valores óptimos. Una buena manera de hacer esto es por simulación por computadora . Los datos observados y calculados se muestran en una pantalla. Los parámetros del modelo se ajustan a mano hasta que el acuerdo entre los datos observados y calculados es razonablemente bueno. Aunque esto será un juicio subjetivo, basta con encontrar un buen punto de partida para el refinamiento no lineal. Las estimaciones de parámetros iniciales se pueden crear usando transformaciones o linealizaciones. Los mejores algoritmos evolutivos como el Algoritmo del Embudo Estocástico pueden conducir a la cuenca convexa de atracción que rodea las estimaciones de los parámetros óptimos. Se ha demostrado que los algoritmos híbridos que usan la aleatorización y el elitismo, seguidos por los métodos de Newton, son útiles y computacionalmente eficientes.

Solución

Cualquier método entre los descritos a continuación se puede aplicar para encontrar una solución.

Criterios de convergencia

El criterio del sentido común para la convergencia es que la suma de los cuadrados no disminuye de una iteración a la siguiente. Sin embargo, este criterio es a menudo difícil de aplicar en la práctica, por diversas razones. Un criterio útil de convergencia es:

${\displaystyle \left|{\frac {S^{k}-S^{k+1}}{S^{k}}}\right|<0.0001.}$

El valor 0.0001 es algo arbitrario y puede que tenga que cambiarse. En particular, puede ser necesario aumentar cuando los errores experimentales son grandes. Un criterio alternativo es:

${\displaystyle \left|{\frac {\Delta \beta _{j}}{\beta _{j}}}\right|<0.001,\qquad j=1,\dots ,n.}$

Una vez más, el valor numérico es algo arbitrario; 0.001 es equivalente a especificar que cada parámetro debe ser refinado a 0.1% de precisión. Esto es razonable cuando es menor que la desviación estándar relativa más grande en los parámetros.

Cálculo del Jacobiano por aproximación numérica

Existen modelos para los cuales es muy difícil o incluso imposible derivar expresiones analíticas para los elementos del Jacobiano. Entonces, la aproximación numérica.

${\displaystyle {\frac {\partial f(x_{i},{\boldsymbol {\beta }})}{\partial \beta _{j}}}\approx {\frac {\delta f(x_{i},{\boldsymbol {\beta }})}{\delta \beta _{j}}}}$

Se obtiene mediante el cálculo de ${\displaystyle f(x_{i},{\boldsymbol {\beta }})\,}$ para ${\displaystyle \beta _{j}\,}$ y ${\displaystyle \beta _{j}+\delta \beta _{j}\,}$. El tamaño debe ser elegido por lo que la derivada numérica no está sujeta a error de aproximación por ser demasiado grande, o el error de redondeo por ser demasiado pequeño.

Referencias

• C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN 0-89871-433-8. Online copy
• T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). Vieweg+Teubner, ISBN 978-3-8348-1022-9.

Notas