Main Content

Algoritmos de mínimos cuadrados (ajuste de modelos)

Definición de mínimos cuadrados

Mínimos cuadrados, en general, es el problema de encontrar un vector x que es un minimizador local para una función que sea una suma de cuadrados, posiblemente sujeta a algunas restricciones:

minxF(x)22=minxiFi2(x)

de forma que A·xb, Aeq·x = beq, lb ≤ x ≤ ub, c(x) ≤ 0, ceq(x) = 0.

Hay varios solvers de Optimization Toolbox™ disponibles para distintos tipos de F(x) y distintos tipos de restricciones:

SolverF(x)Restricciones
mldivideC·xdNinguna
lsqnonnegC·xdx ≥ 0
lsqlinC·xdLímite, lineal
lsqnonlinF(x) generalLineal y no lineal general
lsqcurvefitF(x, xdata) – ydataLineal y no lineal general

Hay seis algoritmos de mínimos cuadrados en los solvers de Optimization Toolbox, además de los algoritmos usados en mldivide:

  • Interior-point lsqlin

  • Conjunto activo lsqlin

  • Trust-region-reflective (mínimos cuadrados no lineales o lineales, límites de restricción)

  • Levenberg-Marquardt (mínimos cuadrados no lineales, límites de restricción)

  • El algoritmo fmincon 'interior-point' modificado para solvers de mínimos cuadrados no lineales lsqnonlin y lsqcurvefit (restricciones no lineales y lineales generales).

  • El algoritmo utilizado por lsqnonneg

Todos los algoritmos excepto el algoritmo active-set de lsqlin son a gran escala; consulte Algoritmos a gran escala frente a algoritmos a media escala. Para ver una descripción general de métodos de mínimos cuadrados no lineales, consulte Dennis [8]. Puede encontrar detalles específicos sobre el método Levenberg-Marquardt en Moré [28].

Mínimos cuadrados lineales: interior-point o active-set

El algoritmo 'interior-point' de lsqlin utiliza interior-point-convex quadprog Algorithm y el algoritmo 'active-set' de lsqlin utiliza el algoritmo active-set de quadprog. La definición del problema de quadprog es minimizar una función cuadrática

minx12xTHx+cTx

sujeta a restricciones lineales y restricciones de límites. La función lsqlin minimiza la norma euclídea al cuadrado del vector Cx – d sujeta a restricciones lineales y restricciones de límites. En otras palabras, lsqlin minimiza

Cxd22=(Cxd)T(Cxd)=(xTCTdT)(Cxd)=(xTCTCx)(xTCTddTCx)+dTd=12xT(2CTC)x+(2CTd)Tx+dTd.

Esto se ajusta al marco de trabajo de quadprog estableciendo la matriz H en 2CTC y el vector c en (–2CTd). (El término aditivo dTd no tiene efecto en la ubicación del mínimo). Después de esta reformulación del problema de lsqlin, quadprog calcula la solución.

Nota

El algoritmo 'interior-point-convex' de quadprog tiene dos rutas de código. Toma una cuando la matriz hessiana H es una matriz ordinaria (completa) de dobles y toma la otra cuando H es una matriz dispersa. Para obtener detalles del tipo de datos dispersos, consulte Matrices dispersas. Por lo general, el algoritmo es más rápido para problemas grandes que tienen relativamente pocos términos distintos de cero cuando especifica H como sparse. De forma similar, el algoritmo es más rápido para problemas pequeños o relativamente densos cuando especifica H como full.

Mínimos cuadrados trust-region-reflective

Algoritmo de mínimos cuadrados trust-region-reflective

Muchos de los métodos utilizados en los solvers de Optimization Toolbox se basan en regiones de confianza, un concepto sencillo, pero potente de optimización.

Para comprender el enfoque trust-region de la optimización, considere el problema de minimización no restringida, minimice f(x), donde la función toma argumentos de vector y devuelve escalares. Suponga que está en un punto x en el espacio n y desea mejorar, es decir, moverse a un punto con un valor de función inferior. La idea básica es aproximar f con una función más sencilla q, que refleja razonablemente el comportamiento de la función f en un entorno N alrededor del punto x. Este entorno es la región de confianza. Se calcula un paso de prueba s minimizando (o minimizando aproximadamente) en N. Este es el subproblema trust-region,

mins{q(s), sN}.(1)

El punto actual se actualiza para ser x + s si f(x + s) < f(x); de lo contrario, el punto actual permanece sin cambiar y N, la región de confianza, disminuye y se repite el cálculo del paso de prueba.

Las principales cuestiones al definir un enfoque específico trust-region para minimizar f(x) son cómo elegir y calcular la aproximación q (definida en el punto actual x), cómo elegir y modificar la región de confianza N y cómo resolver de forma precisa el subproblema trust-region. Esta sección se centra en el problema no restringido. En secciones posteriores, se abordan otras complicaciones causadas por la presencia de restricciones en las variables.

En el método estándar trust-region ([48]), la aproximación cuadrática q se define por los dos primeros términos de la aproximación de Taylor a F en x; el entorno N suele tener forma de esfera o de elipsoide. En términos matemáticos, el subproblema trust-region se indica, por lo general, como

min{12sTHs+sTg  such that  DsΔ},(2)

donde g es el gradiente de f en el punto actual x, H es la matriz hessiana (la matriz simétrica de segundas derivadas), D es una matriz de escalado diagonal, Δ es un escalar positivo y ‖ . ‖ es la norma euclídea. Existen buenos algoritmos para resolver Ecuación 2 (consulte [48]); estos algoritmos incluyen, por lo general, el cálculo de todos los valores propios de H y un proceso de Newton aplicado a la ecuación secular

1Δ1s=0.

Estos algoritmos proporcionan una solución precisa para Ecuación 2. Sin embargo, requieren tiempo proporcional a varias factorizaciones de H. Por lo tanto, para los problemas trust-region es necesario un enfoque diferente. Varias estrategias heurísticas y de aproximación, basadas en Ecuación 2, se han propuesto en la literatura ([42] y [50]). El enfoque de aproximación seguido en los solvers de Optimization Toolbox es restringir el subproblema trust-region a un subespacio bidimensional S ([39] y [42]). Una vez calculado el subespacio S, el trabajo para resolver Ecuación 2 es trivial incluso si se necesita información completa del valor propio/vector propio (ya que, en el subespacio, el problema solo es bidimensional). El trabajo dominante se ha desplazado ahora a la determinación del subespacio.

El subespacio bidimensional S se determina con la ayuda de un proceso de gradiente conjugado precondicionado descrito a continuación. El solver define S como el espacio lineal que abarcan s1 y s2, donde s1 está en la dirección del gradiente g y s2 es una dirección de Newton aproximada, es decir, una solución para

Hs2=g,(3)

o una dirección de curvatura negativa,

s2THs2<0.(4)

La filosofía sobre la que se basa esta elección de S es forzar la convergencia global (mediante la dirección de descenso más pronunciado o la dirección de curvatura negativa) y lograr una convergencia local rápida (mediante el paso de Newton, cuando exista).

Ahora es fácil dar una aproximación de minimización no restringida utilizando ideas trust-region:

  1. Formule el subproblema trust-region bidimensional.

  2. Resuelva Ecuación 2 para determinar el paso de prueba s.

  3. Si f(x + s) < f(x), entonces x = x + s.

  4. Ajuste Δ.

Estos cuatro pasos se repiten hasta la convergencia. La dimensión Δ de la región de confianza se ajusta de acuerdo con las reglas estándar. En particular, disminuye si el paso de prueba no se acepta, es decir, f(x + s) ≥ f(x). Consulte [46] y [49] para ver una exposición sobre este aspecto.

Los solvers de Optimization Toolbox tratan unos pocos casos importantes especiales de f con funciones especializadas: mínimos cuadrados no lineales, funciones cuadráticas y mínimos cuadrados lineales. Sin embargo, las ideas algorítmicas subyacentes son las mismas que para el caso general. Estos casos especiales se abordarán en secciones posteriores.

Mínimos cuadrados no lineales a gran escala

Un caso especial importante para f(x) es el problema de mínimos cuadrados no lineales

minxifi2(x)=minxF(x)22,(5)

donde F(x) es una función de valor vectorial con el componente i de F(x) igual a fi(x). El método básico usado para resolver este problema es el mismo que en el caso general descrito en Métodos trust-region para la minimización no lineal. Sin embargo, se aprovecha la estructura del problema de mínimos cuadrados no lineales para mejorar la eficiencia. En particular, una dirección de Gauss-Newton aproximada, es decir, una solución s para

minJs+F22,(6)

(donde J es la matriz jacobiana de F(x)) se utiliza para ayudar a definir el subespacio bidimensional S. No se utilizan las segundas derivadas de la función del componente fi(x).

En cada iteración, se utiliza el método de gradientes conjugados precondicionados para resolver aproximadamente las ecuaciones normales, es decir,

JTJs=JTF,

aunque las ecuaciones normales no están formadas explícitamente.

Mínimos cuadrados lineales a gran escala

En este caso, la función f(x) que se desea resolver es

f(x)=Cx+d22,

sujeta posiblemente a restricciones lineales. El algoritmo genera iteraciones factibles estrictamente que convergen, en el límite, en una solución local. Cada iteración implica la solución aproximada de un sistema lineal amplio (de orden n, donde n es la longitud de x). Las matrices de iteración tienen la estructura de la matriz C. En particular, se utiliza el método de gradientes conjugados precondicionados para resolver aproximadamente las ecuaciones normales, es decir,

CTCx=CTd,

aunque las ecuaciones normales no están formadas explícitamente.

Se utiliza el método de región de confianza de subespacio para determinar una dirección de búsqueda. Sin embargo, en lugar de restringir el paso a (posiblemente) un paso de reflejo, como en el caso de la minimización no lineal, se realiza una búsqueda de recta de reflejo por tramos en cada iteración, como en el caso cuadrático. Consulte [45] para obtener detalles de la búsqueda de recta. Por último, los sistemas lineales representan un enfoque de Newton capturando las condiciones de optimalidad de primer orden en la solución, lo que implica tasas de convergencia local fuertes.

Función de multiplicación de matriz jacobiana.  lsqlin puede resolver el problema de mínimos cuadrados restringidos linealmente sin utilizar la matriz C explícitamente. En su lugar, utiliza una función de multiplicación de matriz jacobiana jmfun,

W = jmfun(Jinfo,Y,flag)

que se proporciona. La función debe calcular los siguientes productos para una matriz Y:

  • Si flag == 0, entonces W = C'*(C*Y).

  • Si flag > 0, entonces W = C*Y.

  • Si flag < 0, entonces W = C'*Y.

Esto puede ser útil si C es grande, pero contiene suficiente estructura para que pueda escribir jmfun sin formar C explícitamente. Para ver un ejemplo, consulte Jacobian Multiply Function with Linear Least Squares.

Método de Levenberg-Marquardt

El problema de mínimos cuadrados minimiza una función f(x) que es una suma de cuadrados.

minxf(x)=F(x)22=iFi2(x).(7)

Los problemas de este tipo se producen en un gran número de aplicaciones prácticas, especialmente las que incluyen ajustar funciones modelo a datos, como una estimación de parámetros no lineales. Este tipo de problema también aparece en sistemas de control, donde el objetivo es que la salida y(x,t) siga una trayectoria modelo continua φ(t) para el vector x y el escalar t. Este problema se puede expresar como

minxnt1t2(y(x,t)φ(t))2dt,(8)

donde y(x,t) y φ(t) son funciones escalares.

Discretice la integral para obtener una aproximación

minxnf(x)=i=1m(y(x,ti)φ(ti))2,(9)

donde ti están espaciadas uniformemente. En este problema, el vector F(x) es

F(x)=[y(x,t1)φ(t1)y(x,t2)φ(t2)...y(x,tm)φ(tm)].

En este tipo de problema, el valor residual F(x)‖ puede ser pequeño en el valor óptimo, dado que la práctica general es establecer trayectorias de destino alcanzables de manera realista. Aunque puede minimizar la función en Ecuación 7 utilizando una técnica de minimización general no restringida, tal y como se describe en Conceptos básicos de la optimización sin restricciones, a menudo se pueden aprovechar determinadas características del problema para mejorar la eficiencia de iteración del procedimiento de resolución. El gradiente y la matriz hessiana de Ecuación 7 tienen una estructura especial.

Denotando la matriz jacobiana de m por n de F(x) como J(x), el vector gradiente de f(x) como G(x), la matriz hessiana de f(x) como H(x) y la matriz hessiana de cada Fi(x) como Di(x),

G(x)=2J(x)TF(x)H(x)=2J(x)TJ(x)+2Q(x),(10)

donde

Q(x)=i=1mFi(x)Di(x).

Una propiedad de la matriz Q(x) es que cuando el valor residual F(x)‖ tiende a cero a medida que xk se acerca a la solución, Q(x) también tiende a cero. Así, cuando F(x)‖ es pequeño en la solución, un método eficaz es utilizar la dirección de Gauss-Newton como base para un procedimiento de optimización.

En cada iteración principal k, el método de Gauss-Newton obtiene una dirección de búsqueda dk que es una solución del problema lineal de mínimos cuadrados

mindknJ(xk)dk+F(xk)22.(11)

La dirección derivada de este método es equivalente a la dirección de Newton cuando los términos de Q(x) = 0. El algoritmo puede usar la dirección de búsqueda dk como parte de una estrategia de búsqueda de recta para garantizar que la función f(x) disminuye en cada iteración.

El método Gauss-Newton a menudo se encuentra con problemas cuando el término de segundo orden Q(x) es significante. El método de Levenberg-Marquardt resuelve este problema.

El método de Levenberg-Marquardt (consulte [25] y [27]) utiliza una dirección de búsqueda que es una solución del conjunto lineal de ecuaciones

(J(xk)TJ(xk)+λkI)dk=J(xk)TF(xk),(12)

o, de manera opcional, de las ecuaciones

(J(xk)TJ(xk)+λkdiag(J(xk)TJ(xk)))dk=J(xk)TF(xk),(13)

donde el escalar λk controla la magnitud y la dirección de dk y diag(A) significa la matriz de términos diagonales en A. Establezca la opción ScaleProblem en 'none' para elegir Ecuación 12 o establezca ScaleProblem en 'Jacobian' para elegir Ecuación 13.

Establezca el valor inicial del parámetro λ0 usando la opción InitDamping. Ocasionalmente, el valor predeterminado 0.01 de esta opción puede no ser adecuado. Si encuentra que el algoritmo Levenberg-Marquardt hace poco progreso inicial, intente establecer InitDamping en un valor diferente al predeterminado, como 1e2.

Cuando λk es cero, la dirección dk es idéntica a la del método de Gauss-Newton. A medida que λk tiende al infinito, dk tiende a la dirección de descenso más pronunciado, con una magnitud que tiende a cero. Por consiguiente, para algunos λk suficientemente grandes, el término F(xk + dk)‖ < ‖ F(xk)‖ se mantiene verdadero, donde ‖ representa la norma euclidiana. Por tanto, puede controlar el término λk para garantizar el descenso incluso cuando el algoritmo encuentra términos de segundo orden, que restringen la eficiencia del método de Gauss-Newton. Cuando el paso se realiza correctamente (da un valor de función inferior), el algoritmo establece λk+1 = λk/10. Cuando el paso no se realiza correctamente, el algoritmo establece λk+1 = λk*10.

Internamente, el algoritmo Levenberg-Marquardt utiliza una tolerancia de optimalidad (criterio de detención) de 1e-4 veces la tolerancia de la función.

El método de Levenberg-Marquardt, por tanto, utiliza una dirección de búsqueda que es un cruce entre la dirección de Gauss-Newton y la dirección de descenso más pronunciado.

Otra ventaja del método de Levenberg-Marquardt es cuando la matriz jacobiana J es de rango deficiente. En este caso, el método de Gauss-Newton puede tener problemas numéricos, dado que el problema de minimización en Ecuación 11 está mal planteado. En cambio, el método de Levenberg-Marquardt tiene un rango completo en cada iteración y, por tanto, evita estos problemas.

La siguiente figura muestra las iteraciones del método de Levenberg-Marquardt cuando minimiza una función de Rosenbrock, un problema de minimización de conocida dificultad que está en formato de mínimos cuadrados.

El método de Levenberg-Marquardt en la función de Rosenbrock

Para una descripción completa de esta figura, incluidos los scripts que generan los puntos iterativos, consulte Minimización de la función banana.

Restricciones de límites en el método de Levenberg-Marquardt

Cuando el problema contiene restricciones de límites, lsqcurvefit y lsqnonlin modifican las iteraciones de Levenberg-Marquardt. Si un punto iterativo propuesto x se encuentra fuera de los límites, el algoritmo proyecta el paso hasta el siguiente punto factible más cercano. En otras palabras, con P definida como el operador de proyección que proyecta puntos no factibles hasta la región factible, el algoritmo modifica el punto propuesto x a P(x). Por definición, el operador de proyección P opera en cada componente xi independientemente de acuerdo con

P(xi)={lbiif xi<lbiubiif xi>ubixiotherwise

o, de manera equivalente,

P(xi)=min(max(xi,lbi),ubi).

El algoritmo modifica el criterio de detención para la medida de optimalidad de primer orden. Deje que x sea un punto iterativo propuesto. En el caso no restringido, el criterio de detención es

f(x)tol,(14)

donde tol es el valor de la tolerancia de optimalidad, que es 1e-4*FunctionTolerance. En el caso acotado, el criterio de detención es

xP(xf(x))2tolf(x).(15)

Para entender este criterio, observe que si x está en el interior de la región factible, el operador P no tiene efecto. Así, el criterio de detención se convierte en

xP(xf(x))2=f(x)2tolf(x),

que es igual que el criterio de detención no restringido original f(x)tol. Si la restricción de límites está activa, lo que significa que x – ∇f(x) no es factible, en un punto donde el algoritmo debería detenerse, el gradiente en un punto del límite es perpendicular al límite. Por tanto, el punto x es igual a P(x – ∇f(x)), la proyección del paso de descenso más pronunciado, como muestra la siguiente figura.

Condición de detención del algoritmo Levenberg-Marquardt

Sketch of x minus the projection of x minus gradient of f(x)

Algoritmo fmincon modificado para mínimos cuadrados restringidos

Para gestionar restricciones lineales y no lineales, el algoritmo 'interior-point' de los solvers lsqnonlin y lsqcurvefit internamente llama a una versión modificada del algoritmo fmincon 'interior-point'. La versión modificada utiliza una aproximación de matriz hessiana diferente a la aproximación utilizada cuando se ejecuta fmincon. De forma predeterminada, fmincon utiliza la aproximación de matriz hessiana BFGS. Más adelante en esta sección puede encontrar la aproximación de matriz hessiana diferente que utilizan los solvers de mínimos cuadrados. Para obtener detalles del algoritmo fmincon'interior-point' y su aproximación de matriz hessiana, consulte Algoritmo interior-point de fmincon.

En el contexto de mínimos cuadrados, la función objetivo f(x) es un vector, no un escalar. La suma de los cuadrados que se desea minimizar es

fT(x)f(x).

El problema que se desea resolver es

minxf(x)22

de forma que

c(x) ≤ 0

ceq(x) = 0.

El lagrangiano de este problema es

L(x,λE,λI)=f(x)22+λETceq(x)+λITc(x).

El gradiente con respecto a x del lagrangiano es

xL(x,λE,λI)=2FT(x)f(x)+AET(x)λE+AIT(x)λI.

En este caso, F es la jacobiana de la función objetivo

Fi,j(x)=fi(x)xj,

y AE y AI son las jacobianas de la función de restricción de igualdad ceq(x) y la función de restricción de desigualdad c(x), respectivamente.

Las condiciones de Karush-Kuhn-Tucker (KKT) del problema son

2FT(x)f(x)+AET(x)λE+AIT(x)λI=0SλIμe=0ceq(x)=0c(x)+s=0.

En este caso, las variables s y λI son no negativas, S = diag(s), μ es el parámetro de barrera que tiende a cero a medida que aumentan las iteraciones y e es un vector de unos.

Para resolver estas ecuaciones, el software realiza iteraciones. Llame al número de iteración k. El algoritmo intenta resolver las ecuaciones dando un paso Δx, Δs, ΔλE, ΔλI usando esta ecuación de actualización:

[H0AETAIT0Sdiag(λi)0SAE000AIS00][ΔxS1ΔsΔλEΔλI]=[2FT(x)f(x)+AET(x)λE+AIT(x)λISλIμeceq(x)c(x)+s].

En este caso, H es la matriz hessiana de la función objetivo. Usando la aproximación de matriz hessiana BFGS, la matriz hessiana aproximada Bk+1 del paso k + 1 es

Bk+1=BkBkskskTBkskTBksk+ykykTykTsk,(16)

donde

sk=xk+1xkyk=xL(xk+1,λk+1)xL(xk,λk+1).

El valor inicial de esta aproximación de matriz hessiana es B0 = I.

Para ver una exposición sobre este enfoque en fmincon, consulte Actualización de la matriz hessiana.

La actualización de matriz hessiana BFGS de mínimos cuadrados tiene la siguiente diferencia en comparación con la actualización de matriz hessiana de fmincon. La matriz hessiana de la función objetivo de mínimos cuadrados es

xx2f(x)2=2FT(x)F(x)+2ifi(x)xx2fi(x).

La actualización de matriz hessiana de mínimos cuadrados utiliza el primer término, 2FTF, como una cantidad separada, y actualiza la aproximación de matriz hessiana utilizando la fórmula BFGS sin el término 2FTF. En otras palabras, la actualización de matriz hessiana BFGS utiliza la fórmula de actualización de matriz hessiana Ecuación 16 para Bk, donde la aproximación de matriz hessiana para el lagrangiano (función objetivo más multiplicadores de Lagrange por las restricciones) es 2FTF + Bk y las cantidades sk y yk son

sk=xk+1xkyk=xL(xk+1,λk+1)xL(xk,λk+1)2FT(xk+1)F(xk+1)(xk+1xk)=(AET(xk+1)AET(xk))λEk+1+(AIT(xk+1)AIT(xk))λIk+1+2FT(xk+1)f(xk+1)2FT(xk)f(xk)2FT(xk+1)F(xk+1)(xk+1xk).(17)

Esta aproximación de matriz hessiana modificada tiene en cuenta la diferencia entre las iteraciones de fmincon y lsqnonlin o las iteraciones de lsqcurvefit para el mismo problema. Internamente, los solvers de mínimos cuadrados utilizan esta aproximación hessiana Ecuación 16 y Ecuación 17 como una función hessiana personalizada en el solver fmincon.

Consulte también

| | | |

Temas relacionados