mínimos cuadrados no lineales

Cuando las observaciones no son igualmente fiables, se puede minimizar una suma ponderada de cuadrados. $${\displaystyle S=\sum _{i=1}^{m}W_{ii}r_{i}^{2}.}$$

Cada elemento de la matriz de ponderación diagonal W debería, idealmente, ser igual al recíproco de la varianza del error de la medición. ^{[ 1 ]} Las ecuaciones normales son entonces, de forma más general, $${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\mathbf {J} ^{\mathsf {T}}\mathbf {W} \ \Delta \mathbf {y} .}$$

Algunos problemas de mal condicionamiento y divergencia pueden corregirse encontrando estimaciones iniciales de parámetros cercanas a los valores óptimos. Una buena manera de hacerlo es mediante simulación por computadora . Tanto los datos observados como los calculados se muestran en una pantalla. Los parámetros del modelo se ajustan manualmente hasta que la concordancia entre los datos observados y calculados sea razonablemente buena. Aunque se trata de un juicio subjetivo, es suficiente para encontrar un buen punto de partida para el refinamiento no lineal. Las estimaciones iniciales de los parámetros pueden crearse mediante transformaciones o linealizaciones. Mejor aún, los algoritmos evolutivos, como el Algoritmo de Embudo Estocástico, pueden conducir a la cuenca de atracción convexa que rodea las estimaciones óptimas de los parámetros. Se ha demostrado que los algoritmos híbridos que utilizan aleatorización y elitismo, seguidos de métodos de Newton, son útiles y computacionalmente eficientes .

Cualquiera de los métodos que se describen a continuación puede aplicarse para encontrar una solución.

El criterio de sentido común para la convergencia es que la suma de los cuadrados no aumenta de una iteración a la siguiente. Sin embargo, este criterio suele ser difícil de implementar en la práctica, por diversas razones. Un criterio de convergencia útil es $${\displaystyle \left|{\frac {S^{k}-S^{k+1}}{S^{k}}}\right|<0.0001.}$$ El valor 0,0001 es algo arbitrario y puede que sea necesario cambiarlo. En particular, puede que sea necesario aumentarlo cuando los errores experimentales sean grandes. Un criterio alternativo es $${\displaystyle \left|{\frac {\Delta \beta _{j}}{\beta _{j}}}\right|<0.001,\qquad j=1,\dots ,n.}$$

Nuevamente, el valor numérico es algo arbitrario; 0,001 equivale a especificar que cada parámetro debe ajustarse con una precisión del 0,1 %. Esto es razonable cuando es menor que la mayor desviación estándar relativa de los parámetros.

Hay modelos para los que 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 cálculo de${\displaystyle f(x_{i},{\boldsymbol {\beta }})}$para${\displaystyle \beta _{j}}$y${\displaystyle \beta _{j}+\delta \beta _{j}}$. El incremento,${\displaystyle \delta \beta _{j}}$El tamaño debe elegirse de manera que la derivada numérica no esté sujeta a errores de aproximación por ser demasiado grande, ni a errores de redondeo por ser demasiado pequeña.

En la sección correspondiente de la página sobre mínimos cuadrados ponderados se ofrece información al respecto .

Pueden producirse múltiples mínimos en diversas circunstancias, algunas de las cuales son:

• Un parámetro se eleva a una potencia de dos o más. Por ejemplo, al ajustar datos a una curva lorentziana.$${\displaystyle f(x_{i},{\boldsymbol {\beta }})={\frac {\alpha }{1+\left({\frac {\gamma -x_{i}}{\beta }}\right)^{2}}}}$$dónde${\displaystyle \alpha }$es la altura,${\displaystyle \gamma }$es la posición y${\displaystyle \beta }$es la mitad del ancho a la mitad de la altura, hay dos soluciones para la mitad del ancho,${\displaystyle {\hat {\beta }}}$y${\displaystyle -{\hat {\beta }}}$que dan el mismo valor óptimo para la función objetivo.
• Se pueden intercambiar dos parámetros sin cambiar el valor del modelo. Un ejemplo sencillo es cuando el modelo contiene el producto de dos parámetros, ya que${\displaystyle \alpha \beta }$dará el mismo valor que${\displaystyle \beta \alpha }$.
• Un parámetro está en una función trigonométrica, como por ejemplo:${\displaystyle \sin \beta }$, que tiene valores idénticos en${\displaystyle {\hat {\beta }}+2n\pi }$Véase el algoritmo de Levenberg-Marquardt para un ejemplo.

No todos los mínimos múltiples tienen valores iguales de la función objetivo. Los falsos mínimos, también conocidos como mínimos locales, ocurren cuando el valor de la función objetivo es mayor que su valor en el llamado mínimo global. Para asegurar que el mínimo encontrado sea el mínimo global, el proceso de refinamiento debe comenzar con valores iniciales de los parámetros muy diferentes. Si se encuentra el mismo mínimo independientemente del punto de partida, es probable que sea el mínimo global.

Cuando existen múltiples mínimos, se produce una consecuencia importante: la función objetivo tendrá un punto estacionario (por ejemplo, un máximo o un punto de silla ) en algún lugar entre dos mínimos. La matriz de ecuaciones normales no es definida positiva en un punto estacionario de la función objetivo, porque el gradiente se anula y no existe una dirección de descenso única. El refinamiento a partir de un punto (un conjunto de valores de parámetros) cercano a un punto estacionario estará mal condicionado y debe evitarse como punto de partida. Por ejemplo, al ajustar una función lorentziana, la matriz de ecuaciones normales no es definida positiva cuando la semianchura de la lorentziana es cero. ^{[ 2 ]}

En ocasiones, un modelo no lineal puede transformarse en uno lineal. Esta aproximación suele ser aplicable, por ejemplo, en las proximidades del mejor estimador, y constituye uno de los supuestos básicos de la mayoría de los algoritmos de minimización iterativos. Cuando una aproximación lineal es válida, el modelo puede utilizarse directamente para la inferencia mediante mínimos cuadrados generalizados , donde se aplican las ecuaciones del ajuste de plantilla lineal ^{[ 3 ]} .

Otro ejemplo de aproximación lineal sería cuando el modelo es una función exponencial simple, $${\displaystyle f(x_{i},{\boldsymbol {\beta }})=\alpha e^{\beta x_{i}},}$$ que se puede transformar en un modelo lineal tomando logaritmos. $${\displaystyle \log f(x_{i},{\boldsymbol {\beta }})=\log \alpha +\beta x_{i}}$$ Gráficamente esto corresponde a trabajar en un gráfico semilogarítmico . La suma de los cuadrados se convierte en $${\displaystyle S=\sum _{i}(\log y_{i}-\log \alpha -\beta x_{i})^{2}.}$$ Este procedimiento debe evitarse a menos que los errores sean multiplicativos y sigan una distribución log-normal, ya que puede generar resultados engañosos. Esto se debe a que, independientemente de los errores experimentales en y , los errores en log y son diferentes. Por lo tanto, al minimizar la suma de cuadrados transformada, se obtendrán resultados distintos tanto para los valores de los parámetros como para sus desviaciones estándar calculadas. Sin embargo, con errores multiplicativos que siguen una distribución log-normal, este procedimiento proporciona estimaciones de parámetros insesgadas y consistentes.

Otro ejemplo lo proporciona la cinética de Michaelis - Menten , utilizada para determinar dos parámetros.${\displaystyle V_{\max }}$y${\displaystyle K_{m}}$: $${\displaystyle v={\frac {V_{\max }[S]}{K_{m}+[S]}}.}$$ El diagrama de Lineweaver-Burk$${\displaystyle {\frac {1}{v}}={\frac {1}{V_{\max }}}+{\frac {K_{m}}{V_{\max }[S]}}}$$ de${\textstyle {\frac {1}{v}}}$contra${\textstyle {\frac {1}{[S]}}}$es lineal en los parámetros${\textstyle {\frac {1}{V_{\max }}}}$y${\textstyle {\frac {K_{m}}{V_{\max }}}}$pero muy sensible a los errores de datos y con un fuerte sesgo hacia el ajuste de los datos a un rango particular de la variable independiente.${\displaystyle [S]}$.

Las ecuaciones normales $${\displaystyle \left(\mathbf {J} ^{\mathsf {T}}\mathbf {WJ} \right)\Delta {\boldsymbol {\beta }}=\left(\mathbf {J} ^{\mathsf {T}}\mathbf {W} \right)\Delta \mathbf {y} }$$ puede ser resuelto para${\displaystyle \Delta {\boldsymbol {\beta }}}$mediante descomposición de Cholesky , como se describe en mínimos cuadrados lineales . Los parámetros se actualizan iterativamente. $${\displaystyle {\boldsymbol {\beta }}^{k+1}={\boldsymbol {\beta }}^{k}+\Delta {\boldsymbol {\beta }}}$$ donde k es el número de iteración. Si bien este método puede ser adecuado para modelos simples, fallará si se produce divergencia. Por lo tanto, es fundamental protegerse contra la divergencia.

El mínimo en la suma de cuadrados se puede encontrar mediante un método que no implica la formación de las ecuaciones normales. Los residuos con el modelo linealizado se pueden escribir como $${\displaystyle \mathbf {r} =\Delta \mathbf {y} -\mathbf {J} \,\Delta {\boldsymbol {\beta }}.}$$ El jacobiano se somete a una descomposición ortogonal; la descomposición QR servirá para ilustrar el proceso. $${\displaystyle \mathbf {J} =\mathbf {QR} }$$ donde Q es un ortogonal${\displaystyle m\times m}$matriz y R es una${\displaystyle m\times n}$matriz que se particiona en una${\displaystyle n\times n}$bloquear,${\displaystyle \mathbf {R} _{n}}$y un${\displaystyle (m-n)\times n}$bloque cero.${\displaystyle \mathbf {R} _{n}}$es triangular superior.

$${\displaystyle \mathbf {R} ={\begin{bmatrix}\mathbf {R} _{n}\\\mathbf {0} \end{bmatrix}}}$$

El vector residual se multiplica por la izquierda por${\displaystyle \mathbf {Q} ^{\mathsf {T}}}$.

$${\displaystyle \mathbf {Q} ^{\mathsf {T}}\mathbf {r} =\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} -\mathbf {R} \ \Delta {\boldsymbol {\beta }}={\begin{bmatrix}\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} -\mathbf {R} \ \Delta {\boldsymbol {\beta }}\right)_{n}\\\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{m-n}\end{bmatrix}}}$$

Esto no tiene efecto sobre la suma de cuadrados ya que${\displaystyle S=\mathbf {r} ^{\mathsf {T}}\mathbf {Q} \mathbf {Q} ^{\mathsf {T}}\mathbf {r} =\mathbf {r} ^{\mathsf {T}}\mathbf {r} }$porque Q es ortogonal . El valor mínimo de S se alcanza cuando el bloque superior es cero. Por lo tanto, el vector de desplazamiento se encuentra resolviendo $${\displaystyle \mathbf {R} _{n}\ \Delta {\boldsymbol {\beta }}=\left(\mathbf {Q} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{n}.}$$

Estas ecuaciones se resuelven fácilmente ya que R es triangular superior.

Una variante del método de descomposición ortogonal implica la descomposición en valores singulares , en la que R se diagonaliza mediante transformaciones ortogonales adicionales.

$${\displaystyle \mathbf {J} =\mathbf {U} {\boldsymbol {\Sigma }}\mathbf {V} ^{\mathsf {T}}}$$ dónde${\displaystyle \mathbf {U} }$es ortogonal,${\displaystyle {\boldsymbol {\Sigma }}}$es una matriz diagonal de valores singulares y${\displaystyle \mathbf {V} }$es la matriz ortogonal de los autovectores de${\displaystyle \mathbf {J} ^{\mathsf {T}}\mathbf {J} }$o equivalentemente los vectores singulares derechos de${\displaystyle \mathbf {J} }$En este caso, el vector de desplazamiento viene dado por $${\displaystyle \Delta {\boldsymbol {\beta }}=\mathbf {V} {\boldsymbol {\Sigma }}^{-1}\left(\mathbf {U} ^{\mathsf {T}}\ \Delta \mathbf {y} \right)_{n}.}$$

La relativa simplicidad de esta expresión resulta muy útil en el análisis teórico de mínimos cuadrados no lineales. La aplicación de la descomposición en valores singulares se analiza en detalle en Lawson y Hanson. ^{[ 6 ]}

En la literatura científica existen numerosos ejemplos donde se han utilizado diferentes métodos para problemas de ajuste de datos no lineales.

• Inclusión de segundas derivadas en el desarrollo en serie de Taylor de la función modelo. Este es el método de Newton en optimización .$${\displaystyle f(x_{i},{\boldsymbol {\beta }})=f^{k}(x_{i},{\boldsymbol {\beta }})+\sum _{j}J_{ij}\,\Delta \beta _{j}+{\frac {1}{2}}\sum _{j}\sum _{k}\Delta \beta _{j}\,\Delta \beta _{k}\,H_{jk_{(i)}},\ H_{jk_{(i)}}={\frac {\partial ^{2}f(x_{i},{\boldsymbol {\beta }})}{\partial \beta _{j}\,\partial \beta _{k}}}.}$$La matriz H se conoce como matriz hessiana . Si bien este modelo presenta mejores propiedades de convergencia cerca del mínimo, su rendimiento es mucho peor cuando los parámetros se alejan de sus valores óptimos. El cálculo de la matriz hessiana incrementa la complejidad del algoritmo. Este método no es de uso generalizado.
• Método de Davidon-Fletcher-Powell . Este método, una variante del método pseudo-Newton, es similar al anterior, pero calcula la matriz hessiana mediante aproximaciones sucesivas, para evitar tener que utilizar expresiones analíticas para las segundas derivadas.
• Descenso más pronunciado . Si bien se garantiza una reducción en la suma de cuadrados cuando el vector de desplazamiento apunta en la dirección del descenso más pronunciado, este método suele tener un rendimiento deficiente. Cuando los valores de los parámetros están lejos de ser óptimos, la dirección del vector de descenso más pronunciado, que es normal (perpendicular) a los contornos de la función objetivo, es muy diferente de la dirección del vector de Gauss-Newton. Esto hace que la divergencia sea mucho más probable, especialmente porque el mínimo a lo largo de la dirección del descenso más pronunciado puede corresponder a una pequeña fracción de la longitud del vector de descenso más pronunciado. Cuando los contornos de la función objetivo son muy excéntricos, debido a que existe una alta correlación entre los parámetros, las iteraciones del descenso más pronunciado, con corte de desplazamiento, siguen una trayectoria lenta y en zigzag hacia el mínimo.
• Búsqueda de gradiente conjugado . Este es un método mejorado basado en el descenso más pronunciado con buenas propiedades de convergencia teórica, aunque puede fallar en computadoras digitales de precisión finita incluso cuando se usa en problemas cuadráticos. ^{[ 8 ]}

Los métodos de búsqueda directa se basan en la evaluación de la función objetivo para diversos valores de parámetros y no utilizan derivadas. Ofrecen alternativas al uso de derivadas numéricas en el método de Gauss-Newton y los métodos de gradiente.

• Búsqueda de variables alternas. ^{[ 4 ]} Cada parámetro se varía sucesivamente mediante la adición de un incremento fijo o variable, y se conserva el valor que produce una reducción en la suma de cuadrados. El método es simple y efectivo cuando los parámetros no están altamente correlacionados. Presenta propiedades de convergencia muy deficientes, pero puede ser útil para obtener estimaciones iniciales de los parámetros.
• Búsqueda de Nelder-Mead (símplex) . En este contexto, un símplex es un politopo de n  +  1 vértices en n dimensiones; un triángulo en un plano, un tetraedro en el espacio tridimensional, etc. Cada vértice corresponde a un valor de la función objetivo para un conjunto particular de parámetros. La forma y el tamaño del símplex se ajustan variando los parámetros de tal manera que el valor de la función objetivo en el vértice más alto siempre disminuya. Si bien la suma de cuadrados puede disminuir rápidamente al principio, puede converger a un punto no estacionario en problemas cuasiconvexos, como demuestra el ejemplo de MJD Powell.

En la obra Numerical Recipes se ofrecen descripciones más detalladas de estos y otros métodos , junto con código informático en varios lenguajes.