Indice

1 El problema de cuadrados mínimos 3 1.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Cuadrados mínimos lineales . . . . . . . . . . . . . . . . . . . 4

1.2.1 Caracterización de la solución del problema de cuadra-

dos mínimos lineales . . . . . . . . . . . . . . . . . . . 4

1.3 Métodos directos para problemas sparse . . . . . . . . . . . . 8 1.3.1 Problemas de cuadrados mínimos de matrices banda. . 9 1.3.2 Almacenamiento de matrices banda. . . . . . . . . . . 11 1.3.3 Ecuaciones normales para matrices banda. . . . . . . . 11 1.3.4 Descomposición QR para problemas con matrices banda. 11

1.4 Cuadrados mínimos lineales con restricciones . . . . . . . . . 13

1.4.1 Problemas de cuadrados mínimos con restricciones de

igualdad (CMI). . . . . . . . . . . . . . . . . . . . . . 13 1.4.2 Método de eliminación directa. . . . . . . . . . . . . . 14 1.4.3 Método del espacio nulo. . . . . . . . . . . . . . . . . . 15 1.4.4 El problema de cuadrados mínimos con restricciones

de desigualdad. . . . . . . . . . . . . . . . . . . . . . . 15

2 El problema de cuadrados mínimos no lineales 18 2.1 Métodos de tipo Gauss-Newton . . . . . . . . . . . . . . . . . 20 2.2 Métodos de tipo Newton . . . . . . . . . . . . . . . . . . . . . 24

3 Cuadrados mínimos no lineales con restricciones 27

4 El problema de cuadrados mínimos no lineales con restricciones de caja 30 4.1 Descripción del problema . . . . . . . . . . . . . . . . . . . . 30

4.1.1 Definiciones y notación. . . . . . . . . . . . . . . . . . 30 4.2 Algoritmo básico . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.2.1 Construcción del modelo cuadrático. . . . . . . . . . . 32 4.2.2 Resolución del subproblema. . . . . . . . . . . . . . . 32 4.2.3 Aceptación del paso. . . . . . . . . . . . . . . . . . . . 32 4.2.4 El algoritmo. . . . . . . . . . . . . . . . . . . . . . . . 33 4.3 Convergencia del algoritmo . . . . . . . . . . . . . . . . . . . 34 4.4 Criterios de convergencia . . . . . . . . . . . . . . . . . . . . 38 4.5 Detalles de implementación . . . . . . . . . . . . . . . . . . . 39

4.5.1 Tabla 1. . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5.2 Tabla 2. . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.6 Conclusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

Capítulo 1

El problema de cuadrados mínimos

1.1 Introducción

Un problema básico que se presenta en muchas aplicaciones, consiste en aproximar observaciones sujetas a errores por un modelo matemático. Es lógico pensar que a mayor cantidad de observaciones, será posible obtener mayor precisión en el modelo. Esto lleva a la necesidad resolver un sistema sobredeterminado, lineal o no lineal de ecuaciones. Puede probarse que la solución que minimiza la suma de los cuadrados de los errores es óptima en cierto sentido y este es básicamente el problema de cuadrados mínimos.

Fue Gauss quien se adjudicó el planteo del método de Cuadrados Mínimos en 1795, a los 18 años de edad y cabe señalar el comentario de Adrien Legendre en “Novelles Methodes pour la détermination des orbites de comètes”: ‘De tous les principes qu ’on peut proposer pour cet objet, je pense qu’il n’en est pas de plus general, de plus exact, ni d’une application plus facile que celui qui consiste àrendre minimum la somme carrée des erreurs’.^{1}

En primer lugar, hay cálculos de gran interés en el análisis de datos estadísticos que requieren de un algoritmo capaz de calcular eficientemente el ajuste de los mismos. Si bien esto en sí mismo es una justificación para estudiar el problema, otra razón de peso es su contribución al desarrollo de la optimización numérica.

En muchas ramas de aplicación se hace referencia al problema desde diferentes ópticas: los matemáticos tienen en mente hallar el punto más cercano a una superficie, los estadísticos lo utilizan en análisis de regresión y los ingenieros llegan a él cuando trabajan en estimación de parámetros, filtros o procesos de identificación. Dada la facilidad con que en la actualidad se dispone de datos para determinados procesos es rutinaria la resolución de problemas de gran tamaño; y la aplicación a procesamiento de señales ha creado gran demanda de métodos estables y eficientes para modificar soluciones de cuadrados mínimos, cuando se agregan o suprimen datos.

1De todos los principios que se pueden proponer para este objetivo, pienso que no hay uno más general, más exacto ni de una aplicación más fácil que el de minimizar al suma de los cuadrados de los errores.

En el este capítulo se hace una revisión del problema de cuadrados mínimos lineales (CML), una descripción de los métodos de resolución y se trata el problema de CML con restricciones. En el segundo capítulo se presenta el problema no lineal (CMNL) y se introduce la notación. Se describen dos métodos iterativos diferentes basados en sendos modelos y se efectúan comparaciones entre ambos.

En el capítulo 3 se comenta el problema general de cuadrados mínimos con restricciones y en el cuarto se describe el problema con restricciones de caja, se presenta un nuevo algoritmo para resolverlo, resultados numéricos logrados con el mismo y comparaciones con los obtenidos por otros autores con otras técnicas.

1.2 Cuadrados mínimos lineales

Este problema proviene de la necesidad de aproximar por un modelo matemático lineal un conjunto de datos u observaciones y así surge la necesidad de resolver un sistema sobredeterminado de ecuaciones lineales.

En términos matriciales podemos plantearlo de la siguiente manera: Dado un vector b ∈ ℝ^{n} y una matriz A ∈ ℝ^{m×n}, m > n, se desea hallar

x ∈ ℝ^{n} tal que Ax es la ‘mejor aproximación’ a b.

Hay muchas maneras de definir la expresión ‘la mejor aproximación’. Una elección que se basa en consideraciones estadísticas [2] y conduce a un problema diferenciable, es buscar la solución de

min_{x} ∥Ax − b∥_{2}^{2}  A ∈ ℝ^{m×n}, b ∈ ℝ^{n}  (1.1)

donde ∥∥_{2} denota la norma euclídea.

El vector r = b − Ax se llama vector residual. Por lo que hemos men-

cionado, es evidente que la solución minimiza ∥r∥_{2}^{2} = ∑_{i=1}^{m} r_{i}^{2} es decir, la suma

de los cuadrados de los residuos. Si rank(A) < n, entonces la solución x de (1.1) no es única, sin embargo, de todas las soluciones existe una de norma mínima, como veremos más adelante.

1.2.1 Caracterización de la solución del problema de cuadrados mínimos lineales

Sea S = {x/ ∥Ax − b∥ = min} el conjunto solución del problema (1.1) .

Puede probarse que x ∈ S si y sólo si e verifica la siguiente condición de ortogonalidad:

A^{T} (b − Ax) = 0.  (1.2)

De aquí resulta que la solución del problema de cuadrados mínimos lineales satisface las ecuaciones normales:

A^{T} Ax = A^{T} b.  (1.3)

Estas ecuaciones son consistentes, ya que si definimos,

R(A) = {z = Ax/x ∈ ℝ^{n}}

resulta

A^{T} b ∈ R(A^{T}) = R(A^{T} A).

Más aún, es válido el siguiente teorema:

Teorema 1 La matriz A^{T} A es definida positiva sí y sólo sí A es de rango completo.

En tal caso, es evidente que la solución de (1.1) es única.

En caso contrario, el subespacio nulo de A, N (A) es no trivial y por lo tanto la solución no es única . Si x^{∗} es una solución del problema, entonces el conjunto del todas las soluciones es

S = {x=x^{∗} + z/z ∈ N (A)} .

Si x^{∗} ⊥ N (A) entonces, es la única solución del problema de norma mínima.

El problema de computar la solución de norma mínima de un sistema sobredeterminado de ecuaciones, es decir,

min_{A^{T} y=c} ∥y∥_{2}

con A ∈ ℝ^{m×n}, aparece como un subproblema en algoritmos de optimización. En el caso de matrices de rango completo el sistema A^{T} y = c es consistente y la única solución satisface las ecuaciones normales de segunda clase [2]:

A^{T} Az = c, y = Az,

es decir, y = A(A^{T} A)^{−1}c.

Otra caracterización de la solución de (1.1) está dada por el llamado sistema aumentado.

En el caso que rank(A) = n, el sistema aumentado

[ \begin{matrix}I & A \\ A^{T} & 0 \end{matrix} ] [ \begin{matrix}y \\ x \end{matrix} ] = [ \begin{matrix}b \\ c \end{matrix} ]

es no singular y da la condición para las solución de los problemas primal y dual de cuadrados mínimos:

min_{x} {∥Ax − b∥_{2}^{2} + 2c^{T} x} ,  (1.4)
min_{y} ∥y − b∥_{2}^{2} ,  A^{T} y = c,  (1.5)

pues haciendo c = 0 en (1.4) se obtiene la solución de (1.1). Haciendo b = 0 en (1.5) se obtiene la solución de norma mínima del sistema lineal A^{T} y = c.

Para resolver (1.3) teóricamente puede utilizarse cualquier método para resolución de sistemas lineales, pero hay unos más adecuados que otros. Uno muy utilizado, sobre todo en aplicaciones estadísticas, es el de Cholesky que describimos a continuación.

Teorema 2 Sea M ∈ ℝ^{n×n} simétrica y definida positiva. Entonces existe una única matriz triangular superior R con elementos positivos en la diagonal tal que M = RR^{T} .

Cuando se ha computado R, factor de Cholesky de A^{T} A, la solución de (1.1) se obtiene resolviendo dos sistemas triangulares:

R^{T} z = d Rx = z

La descomposición de Cholesky puede utilizarse también de la siguiente manera para resolver el problema que nos ocupa:

Teorema 3 : Sea el sistema aumentado

A= (A, b) ∈ ℝ^{m×(n+1)}

y el correspondiente producto

M=A^{T} A= [ \begin{matrix}M & d \\ d^{T} & b^{T} b \end{matrix} ] .

Sea R el factor de Cholesky de M,

R= [ \begin{matrix}R & z \\ 0 & ρ \end{matrix} ] .

Entonces la solución de (1.1) y la norma del residuo se obtienen de

Rx = z ∥b − Ax∥_{2} = ρ.

Cabe señalar que trabajar con el sistema aumentado en muchas ocasiones simplifica el trabajo y es más ventajoso que otros métodos.

El estudio del error al aplicar el método de Cholesky fue hecho con detalle por Wilkinson pero cabe señalar que el método puede fallar, pudiendo aparecer raíces cuadradas de números negativos cuando el número de condición de la matriz A, κ(A) es del orden de ε^{− \frac{1}{2}} , si ε denota el machine epsilon de la computadora.

Aunque las ecuaciones normales definen biunívocamente la solución cuando A es de rango completo, por lo que hemos mencionado, no son siempre recomendables para hallarla explícitamente y desde el punto de vista numérico hay otros métodos más satisfactorios.

Por ejemplo, si rank(A) = n, descomposición QR conduce a una ortonormalización numéricamente estable de las columnas de A.

El siguiente teorema caracteriza la solución de (1.1) cuando se aplica este método:

Teorema 4 Sea m ≥ n > 0, b ∈ ℝ^{m}, A ∈ ℝ^{m×n}de rango completo. Entonces existe una descomposición A = QR con Q ∈ ℝ^{m×m} matriz ortogonal, R ∈ ℝ^{m×n}, triangular superior y R_{u}, las n primeras filas de R, matriz no singular triangular superior. La única solución de (1.1) es

x_{∗} = R_{u}^{−1}(Q^{T} b)_{u},

donde (Q^{T} b)_{u} = ((Q^{T} b)_{1}, ..., (Q^{T} b)_{n}); y el mínimo de ∥Ax − b∥_{2}^{2} es

∑_{i=n+1}^{m}[(Q^{T} b)_{i}]^{2}.

La prueba de este teorema dada en [2] muestra como computar Q y R utilizando transformaciones ortogonales. El uso sistemático de transformaciones ortogonales para reducir matrices a formas más simples fue iniciada por Givens y Householder. Si bien Householder trabajó con cuadrados mínimos, fue Golub [20] quien aplicó las transformaciones al problema que nos ocupa.

Si A no es de rango completo, como ya mencionamos, la solución no es única y en este caso nos restringimos hallar el de mínima norma, y puede plantearse así:

Dada A ∈ ℝ^{m×n}, b ∈ ℝ^{m},resolver:

min {∥x^{∗}∥_{2} : ∥Ax^{∗} − b∥_{2} ≤ ∥Ax − b∥ , ∀x ∈ ℝ^{n}} .  (1.6)

Este problema es de interés para m mayor, igual o menor que n, y en todos los casos puede resolverse utilizando la descomposición en valores singulares (SVD) de la matriz A. Este método se define a continuación, los factores están caracterizados en el teorema que muestra cómo utilizar SVD para resolver el problema (1.6).

Definición 1 : Sean A ∈ ℝ^{m×n}, k = min{m, n}. La descomposición en valores singulares de A es A = UDV^{T} , donde U ∈ ℝ^{m×m} y V ∈ ℝ^{n×n} son matrices ortogonales, D ∈ ℝ^{m×n} es una matriz con d_{ii} = σ_{i} ≥ 0, i = 1, ..., k, d_{ij} = 0 si i ≠ j. Las cantidades σ_{1}, ..., σ_{k}se llaman valores singulares de A .

Teorema 5 Sea A ∈ ℝ^{m×n}. Entonces la descomposición SVD de A existe, los valores singulares son las raíces cuadradas no negativas de los autovalores de AA^{T} si m ≥ n , y las columnas de U y V son los autovectores de AA^{T} y A^{T} A respectivamente. El número de valores singulares no nulos es igual a rank(A).

A continuación veamos cómo caracterizar la solución de (1.6) utilizando descomposición en valores singulares .

Teorema 6 Sea A ∈ ℝ^{m×n} con una descomposición en valores singulares como definimos más arriba. Sea A^{+}la matriz pseudoinversa de A, definida por:

A^{+} = V D^{+}U^{T}  D^{+} = \{ \begin{matrix}d_{ii}^{+} = \{ \begin{matrix}\frac{1}{σ_{ii}} σ_{i} > 0 \\ 0 σ_{i} = 0 \end{matrix} \\ d_{ij}^{+} = 0  i ≠ j. \end{matrix}

Entonces la única solución de (1.6) es x_{∗} = A^{+}b.

Esta solución se caracteriza geométricamente por verificar las siguientes condiciones :

x_{∗} ⊥ N (A), Ax_{∗} = P_{R(A)}b

donde P_{R(A)} es la proyección ortogonal sobre R(A).

La descomposición en valores singulares se encuentra en un proceso iterativo similar al de hallar los autovalores y autovectores de una matriz simétrica, y es más costoso de obtener que cualquier otra descomposición. Este método es recomendable para resolver (1.6) cuando se desconoce el rango de A . Si m < n y A es de rango completo, puede probarse que es más eficiente la descomposición LQ. Esto es A = LQ , L ∈ ℝ^{n×n} triangular inferior, Q ∈ ℝ^{n×n}ortogonal.

La descomposición SVD es la técnica más recomendable para computar el rango de una matriz y puede también utilizarse para determinar la sensibilidad de x = A^{−1}b para una matriz cuadrada A con respecto a una variación en los datos △b. Sin embargo, la ventaja de QR sobre SV D reside en que esa descomposición puede actualizarse en forma más eficiente cuando A está sujeta a modificaciones de bajo rango.

1.3 Métodos directos para problemas sparse

A continuación haremos una descripción de métodos que se utilizan para resolver el problema (1.1) en el caso en que A es una matriz sparse.

Recordemos la definición de J.H. Wilkinson de matriz sparse: ‘una matriz con suficientes elementos nulos como para intentar tomar ventaja de ello’. A menudo la ganancia en operaciones y almacenamiento es tan importante que de otra manera sería imposible resolver el problema. Notar que si los problemas son muy grandes, es conveniente que sean sparse o tengan alguna estructura especial para que puedan ser resueltos por un método que explote eficientemente estas características. Estos casos suelen presentarse, por ejemplo, en problemas que estudian la estructura molecular, el campo gravitatorio de la tierra, estudio de fuerzas en análisis estructural, reconstrucción de imágenes, etc.

Un problema de cuadrados mínimos muy grande fue descripto por Kolata en 1978, y consiste en el ajuste de coordenadas de las estaciones geodésicas comprendidas en el North American Datum. Eran aproximadamente 6.5 millones de ecuaciones y 540.000 incógnitas. Las ecuaciones eran medianamente no lineales y por lo tanto fue suficiente resolver dos o tres sistemas lineales de este tamaño. En el futuro está previsto el ajuste de las medidas geodésicas de toda la tierra, lo que muestra la importancia de los problemas de gran tamaño.

Supongamos que trabajamos con una matriz A de rango completo. Los problemas de cuadrados mínimos lineales pueden resolverse por métodos directos o iterativos. Resolver un problema de cuadrados mínimos sparse usando el método de ecuaciones normales o descomposición QR es similar a resolver un sistema sparse definido positivo utilizando factorización de Cholesky. Una forma distinta de resolverlos esta basada en el sistema aumentado y utiliza métodos que computan descomposiciones sparse de matrices sparse, simétricas e indefinidas.

La elección del método para la resolución de un problema sparse depende, en parte, de las disponibilidades computacionales con que se cuente y la eficiencia de un algoritmo depende considerablemente de ellas.

1.3.1 Problemas de cuadrados mínimos de matrices banda.

La clase más simple de matrices sparse es la de tipo ‘banda’. Una ‘matriz banda’ A ∈ ℝ^{m×n} tiene la propiedad que en toda fila los elementos no nulos yacen en una banda fija. Definimos

f_{i}(A) = {j/a_{ij} ≠ 0} l_{i}(A) = {j/a_{ij} ≠ 0}

los subíndices de las primera y última columna no nula de la i-ésima fila.

Definición 2 Dada la matriz A ∈ ℝ^{m×n} se dice que A tiene filas de ancho de banda ω si

ω(A) = max_{1≤i≤m} (l_{i}(A) − f_{i}(A) + 1).

Para que la estructura tenga sentido práctico debe ser ω ≪ n. Notar que si bien el ancho de banda no depende del ordenamiento de las filas, es dependiente del ordenamiento de las columnas. Para matrices simétricas, el ancho de banda se define como el ancho de banda de su parte estrictamente triangular superior (o inferior).

Definición 3 El ancho de banda ( o media banda) de una matriz simétrica C ∈ ℝ^{n×n} se define como

p = max_{1≤i≤n} (l_{i}(C) − i).

Ahora veremos una relación entre el ancho de fila ω de la matriz A y el ancho de banda p de la correspondiente matriz de ecuaciones normales C = A^{T} A .

Teorema 7 Sea A ∈ ℝ^{m×n} de filas con ancho de banda ω . Entonces el ancho de banda de C = A^{T} A es a lo sumo ω − 1.

A menudo es ventajoso reordenar las filas de A de modo que los índices f_{i}(A) i = 1, ..., n formen una sucesión no decreciente, es decir,

i ≤ k ⟹ f_{i}(A) ≤ f_{k}(A).

En este caso, diremos que la matriz A tiene forma standard, y entonces puede escribirse

A = ( \begin{matrix}A_{1} \\ A_{2} \\ . \\ . \\ A_{q} \end{matrix} )  q ≤ n,  (1.7)

donde el bloque A_{k} está formado por todas las filas de A cuyo primer elemento no nulo está en la k-ésima columna. Puede probarse por inducción el siguiente teorema:

Teorema 8 Si ω_{k} es el ancho de banda del bloque A_{k}, entonces la estructura de banda de C = A^{T} A verifica:

l_{1}(C) = ω_{1} l_{k}(C) = max {ω_{k} + k − 1, l_{k−1}(C)} k = 2, ..., q.

A continuación veamos la estructura resultante del factor de Cholesky R de A^{T} A . Para ello es necesario una definición previa:

Definición 4 La cápsula de una matriz simétrica C es el conjunto de índices:

Env(C) = {(i, j)/i ≤ j ≤ l_{i}(C)} .

Teorema 9 Si R es el factor de Cholesky de la matriz simétrica y definida positiva C = A^{T} A, entonces

Env(R + R^{T}) = Env(C).

1.3.2 Almacenamiento de matrices banda.

La estructura de una matriz banda debe explotarse adecuadamente cuando se la almacena, con el objeto de evitar ocupar inútilmente espacios de memoria. Uno de los métodos para hacerlo, se conoce como ‘almacenamiento empaquetado por filas’, en el cual los elementos son almacenados en un vector en forma contigua fila por fila, comenzando con el primer elemento no nulo de cada fila y terminando con el último elemento no nulo de la misma.

De acuerdo al Teorema 7 la matriz de ecuaciones normales de una matriz banda, tiene estructura simétrica de banda, y en este caso es aconsejable almacenar de cada fila los elementos que van desde la diagonal hasta el último no nulo.

1.3.3 Ecuaciones normales para matrices banda.

Cuando se resuelve el problema de cuadrados mínimos utilizando ecuaciones normales, se comienza calculando

C = A^{T} A = ∑_{i=1}^{m} a_{i}a_{i}^{T} d = A^{T} b = ∑_{i=1}^{m} b_{i}a_{i}

donde a_{i}^{T} , i = 1, ..., m son las filas de A. Luego se calcula la factorización de Cholesky C = R^{T} R y finalmente se resuelven dos sistemas triangulares.

Si asumimos que un flop representa una operación de multiplicación y suma [2], puede probarse que si A tiene ancho de fila constante ω, la construcción de C y d requiere aproximadamente \frac{1}{2}mω(ω+s) flops, suponiendo que b no tiene elementos nulos.

De los Teoremas 7 y 9 se deduce que C y R + R^{T} tienen ancho de banda p = ω − 1 y cualquier algoritmo que calcule el factor C debe explotar este hecho.

Puede comprobarse que si p ≪ n la factorización de C requiere aproximadamente \frac{1}{2}n(p+1)(p+2) flops y n raíces cuadradas. En las sustituciones hacia atrás y adelante que se utilizan para resolver los sistemas triangulares resultantes se efectúan (2n − p − 1)(p + 2) multiplicaciones. Si se almacenan adecuadamente las matrices involucradas, la solución del problema de cuadrados mínimos donde A es una matriz banda de ancho ω, requiere \frac{1}{2}(mω(ω + 3) + n(ω − 1)(ω + 2)) + n(2ω − 1) + O(m + n) flops.

1.3.4 Descomposición QR para problemas con matrices banda.

En primer lugar estudiemos el caso en que la descomposición QR se obtiene utilizando transformaciones de Givens para resolver el problema de cuadrados mínimos min_{x} ∥Ax − b∥_{2} . De los Teoremas 7 y 9 resulta que R y

la porción triangular superior de A^{T} A son matrices banda con elementos no nulos sólo en las primeras p = ω −1 superdiagonales. Por lo tanto, R tiene ω elementos no nulos en cada fila. Explotando estos detalles, Cox [6] propuso un algoritmo basado en rotaciones de Givens que se aplica a matrices banda en forma standard. La ortonormalización requiere en este caso alrededor de 2nω^{2}flops y puede efectuarse utilizando \frac{1}{2}ω(ω + 3) lugares de memoria. Cabe destacar que si la matriz A no está ordenada en forma standard, la cantidad de flops necesarios se incrementa en un factor \frac{n}{ω} [2].

Este algoritmo puede modificarse adecuadamente para aplicarlo a matrices de ancho de banda variable, y en este caso conviene tener en cuenta la cápsula de A^{T} A para construir R.

En el caso en que m ≫ n, es más conveniente utilizar transformaciones de Householder. Fue Reid [39] quien las aplicó por primera vez a matrices banda. Si A está ordenada en la forma (1.7), el algoritmo procede en q etapas. Después de k − 1 pasos, los primeros k − 1 bloques han sido reducidos por transformaciones de Householder a una matriz banda trapezoidal superior R_{k−1}. En la etapa k, se trabaja con el bloque A_{k} y se computa

Q_{k}^{T} \binom{R_{k−1}}{A_{k}} = R_{k}

donde Q_{k} es producto de transformaciones de Householder y R_{k} es trapezoidal superior. Notar que como el bloque A_{k} tiene su primer elemento no nulo en la k-ésima columna, ésta y otras etapas posteriores no considerarán las primeras (k −1) filas y columnas de R_{k−1}, razón por la cual, las primeras filas (k − 1) de R_{k−1} son filas de la matriz final, R.

Utilizando este método son necesarios ω(ω + 1)(m + \frac{3n}{2}) flops, lo que representa aproximadamente la mitad de los necesarios si se aplicara rotaciones de Givens. Ambos métodos pueden aplicarse a varios términos independientes b para producir c = Q^{T} b y resolver en tal caso, otros tantos problemas de cuadrados mínimos por sustitución hacia atrás del sistema triangular Rx = c. Cabe destacar que es esencial subdividir las transformaciones de Householder tal como se indicó para evitar el crecimiento inmediato de la cantidad de operaciones.

Hay otros tipos de estructuras que aparecen frecuentemente en muchos problemas de cuadrados mínimos de gran tamaño, por ejemplo, la forma de ‘bloque angular’

A = ( \begin{matrix}A_{1} & & & & & B_{1} \\ & A_{2} & & & & B_{2} \\ & & . & & & . \\ & & & . & & . \\ & & & & . & . \\ & & & & A_{M} & B_{M} \end{matrix} ) , x = ( \begin{matrix}x_{1} \\ x_{2} \\ . \\ . \\ . \\ x_{M+1} \end{matrix} ) , b = ( \begin{matrix}b_{1} \\ b_{2} \\ . \\ . \\ . \\ b_{M} \end{matrix} ) ,

A_{i} ∈ ℝ^{m_{i}×n_{i}}  B_{i} ∈ ℝ^{m_{i}×n_{M+1}}  i = 1, 2, ..., M

m = m_{1} + m_{2} + ... + m_{M} n = n_{1} + n_{2} + ... + n_{M+1} aparece en aplicaciones de fotogrametría [21], posición de radar de Doppler

[30] y problemas en geodesia [22].

En estos casos la matriz C = A^{T} A de ecuaciones normales y el factor de Cholesky [26] correspondiente tienen una estructura previsible y el método que se utilice para resolver el problema de cuadrados mínimos debe aplicarse explotando este hecho.

1.4 Cuadrados mínimos lineales con restricciones

En esta sección trataremos el problema de cuadrados mínimos lineales con restricciones. Tales problemas se plantean en una amplia variedad de aplicaciones; que contienen restricciones de igualdad aparecen frecuentemente con la necesidad de interpolar datos o ajustar curvas que se ‘peguen’ con continuidad ellas y sus derivadas hasta cierto orden; las relaciones de desigualdad provienen de requerimientos tales como la positividad, monotonicidad y convexidad.

A los problemas de cuadrados mínimos lineales con restricciones de igualdad los denotamos CMI y a los que involucran desigualdades CMD.

A continuación describiremos brevemente estos problemas y los métodos más usuales que los resuelve.

1.4.1 Problemas de cuadrados mínimos con restricciones de igualdad (CMI).

Dadas las matrices A ∈ ℝ^{m×n}, B ∈ ℝ^{p×n} y los vectores b ∈ ℝ^{m}, d ∈ ℝ^{p}, el problema consiste en hallar x ∈ ℝ^{n} tal que

min_{x} ∥Ax − b∥  s.a. Bx = d.  (1.8)

Obviamente este problema tiene solución sí y sólo sí el sistema lineal Bx = d es consistente, y una condición para que esto se cumpla es que rank(B) = p.

Además, la solución de (1.8) es única sí y sólo sí

rank(B) = p y rank\binom{A}{B} = n.  (1.9)

Puede probarse que la segunda condición de (1.9), es equivalente a que N (A) ∩ N (B) = {0}.

Un algoritmo robusto para el problema CMI necesariamente debe chequear la posible inconsistencia de las restricciones.

Si las restricciones son consistentes, el problema (1.8) pueden reformularse como sugirieron Leringe y Wedin [27], como un problema de cuadrados mínimos secuencial:

min_{x∈S} ∥Ax − b∥_{2}  S = {x/ ∥Bx − d∥_{2} = min}  (1.10)

que tiene solución única de norma mínima. La mayoría de los algoritmos descriptos a continución pueden adaptarse para resolver (1.10) con una ligera modificación.

La forma más natural de resolver un CMI es transformarlo en un problema equivalente de CM sin restricciones de menor dimensión. Básicamente existen dos formas de llevar a cabo esta reducción: por eliminación directa y por el método del espacio nulo.

1.4.2 Método de eliminación directa.

Se basa en la descomposición QR de la matriz B. Existe una matriz de permutación P_{B} tal que si r = rank(B) ≤ p, se puede escribir:

Q_{B}^{T}BP_{B} = ( \begin{matrix}R_{11} & R_{12} \\ 0 & 0 \end{matrix} ) \begin{matrix}\}p \\ \}p − r \end{matrix}

aplicando Q_{B}^{T} al vector d resulta

(R_{11}, R_{12}) x=d_{1}, d= Q_{B}^{T}d = (\begin{matrix}d & 1 \\ d & 2 \end{matrix}),

con x= P_{B}^{T} x y d_{2}= 0 sí y sólo sí las restricciones son consistentes.

Aplicando la permutación P a la matriz A y utilizando las expresiones de arriba, puede probarse que

Ax − b =Ax −b = (A_{1}, A_{2}) (\begin{matrix}x_{1} \\ x_{2} \end{matrix}) − b.

Si definimos

A^{∼}_{2} = A_{2} − A_{1} R_{11}^{−1}R_{12}
b^{∼} = b− A_{1} R_{11}^{−1}d_{1}^{−1}

se obtiene el siguiente problema de cuadrados mínimos sin restricciones:

min_{x} ∥A^{∼}_{2}x_{2} − b^{∼}∥_{2} A^{∼}_{2}∈ ℝ^{m×(n−r)}

equivalente al problema CMI, y puede resolverse utilizando la descomposición QR de A^{∼}_{2}. Más detalles pueden consultarse en [2].

1.4.3 Método del espacio nulo.

Se aplica cuando la matriz B es de rango completo. Si QR es su descomposición y se particiona Q de la forma Q = (Q_{1}, Q_{2}), Q_{1} ∈ ℝ^{n×p}, Q_{2} ∈ ℝ^{n×(n−p)}, puede probarse que las columnas de Q_{2} constituyen una base ortonormal del subespacio nulo de B.

Usando este hecho se verifica que existe x_{1} ∈ ℝ^{n} tal que la solución de norma mínima del problema

min_{y} ∥(AQ_{2})y − (b − Ax_{1})∥_{2}

es tal que x = x_{1} +Q_{2}y es la solución de norma mínima del problema CMI.

Si bien los dos métodos descriptos anteriormente son los más naturales para resolver el problema, existen otros no menos exitosos como los que proponen Van Loan [43] y Anda y Park [1].

El caso de restricciones cuadráticas de igualdad

min_{x} ∥Ax − b∥_{2} s.a. ∥Cx − d∥ = α

fue tratado por Gander [15] y su propuesta se basa en la resolución de las ecuaciones normales que resultan al plantear la condición suficiente para existencia de mínimo para la función Lagrangiano:

L(x, λ) = ∥Ax − b∥^{2} + λ (∥Cx − d∥^{2} − α^{2}) .

1.4.4 El problema de cuadrados mínimos con restricciones de desigualdad.

A continuación haremos una breve reseña de los distintos tipos de restricciones de desigualdad (CMD) que pueden presentarse y mencionaremos algunos métodos para resolverlos

Problemas con restricciones lineales de desigualdad.

min_{x} ∥Ax − b∥_{2} s.a. l ≤ Cx ≤ u A ∈ ℝ^{m×n}, C ∈ ℝ^{p×n}. (1.11)

Observar que si además hay restricciones de igualdad, éstas pueden eliminarse utilizando alguno de los métodos propuestos para resolver problemas CMI.

Un caso importante es aquél en que las restricciones son simplemente cotas en las variables:

CMB  min_{x} ∥Ax − b∥_{2} s.a. l ≤ x ≤ u.  (1.12)

Esta problema aparece en muchas aplicaciones como problemas de reconstrucción en geodésica y tomografía, problemas de contacto para sistemas mecánicos y modelado de circulación oceánica.

Si A es una matriz de rango completo, (1.12) es un problema de optimización estrictamente convexo y por lo tanto, tiene solución única para cualquier vector b.

Otra posibilidad que aparece frecuentemente en la práctica son las restricciones de no negatividad,

CMNN  min_{x} ∥Ax − b∥_{2} s.a. x ≥ 0.  (1.13)

Otro caso particular de CMD es el problema de distancia mínima,

CMD  min_{x} ∥x∥_{2} s.a. g ≤ Gx ≤ h  (1.14)

o más generalmente,

min_{x} ∥x∥_{2} s.a. g ≤ G_{1}x_{1} + G_{2}x_{2} ≤ h  (1.15)
con x = \binom{x_{1}}{x_{2}} x_{1} ∈ ℝ^{k}, x_{2} ∈ ℝ^{n−k}.

Cabe señalar que el problema (1.11) es equivalente al problema de programación cuadrática

min_{x} (x^{T} Bx + c^{T} x) s.a. l ≤ Cx ≤ u con B = A^{T} A c = −2A^{T} b,

(1.16)

que aparece como subproblema en muchos algoritmos de programación no lineal. Por esta razón ha sido estudiado ampliamente y hay muchos algoritmos para resolverlo. En la aplicación al problema CMD, la matriz B resulta semidefinida positiva y por lo tanto (1.16) es un problema convexo.

Desde el punto de vista numérico no es aconsejable aplicar algoritmos de programación cuadrática a (1.11) ya que conviene evitar el cálculo explícito de la matriz B y por ello hay muchos algoritmos para problemas de programación cuadrática que se han adaptado a trabajar directamente con la matriz A. En algunas aplicaciones importantes el problema CMD puede definirse con una matriz triangular A, y en este caso se transforma en un problema de distancia mínima [40].

Lawson y Hanson [26], Haskell y Hanson [24], siguen la misma idea y usan aproximación dual para resolverlo. Estos últimos y Cline [5], describieron cómo modificar el problema CMD sólo con cotas superiores en un problema resoluble en dos etapas, en cada una de las cuales se resuelve un problema del tipo CMNN.

Mención aparte merecen los algoritmos que utilizan la estrategia de conjunto de restricciones activas [19],[26]. Es un método iterativo que genera una sucesión de puntos en una cara del politopo factible hasta alcanzar un punto que minimiza la función objetivo en la cara o un punto perteneciente a la frontera de la misma. En el primer caso se pasa a trabajar en una cara de dimensión mayor, y en el otro caso, a una cara de menor dimensión.

Estos algoritmos usualmente restringen los cambios en le conjunto de trabajo a una sóla dimensión, es decir, en cada iteración se incorpora o descarta una sola restricción, lo que implica que en los problemas de gran escala se requieren, en general, un gran número de iteraciones. Otra desventaja que cabe señalar es que antes de desechar restricciones, se debe hallar el mínimo exacto en cada cara, lo cual también puede demandar muchas iteraciones.

Para resolver el problema CMD Stoer [42] propone un algoritmo basado en este método para el caso de matrices de rango completo. Schittkowski y Stoer [41] lo implementaron utilizando descomposición de Gram Schmidt y la ventaja de esta propuesta reside en que la variación de los datos es tomada en cuenta con bastante facilidad.

Un conjunto muy completo de subrutinas en lenguaje Fortran para este problema fue desarrollado por Gill y otros [18].

Para el problema particular CMB, se utiliza también esta estrategia y hay exitosas implementaciones en Fortran como la de Lawson y Hanson [26] y una adaptación del método propuesto por Stoer para el caso general hecha por Zimmermann [46].

El problema de cuadrados mínimos no lineales

Sea R : ℝ^{n} → ℝ^{m}  R(x) = (r_{1}(x), r_{2}(x), ...., r_{m}(x))  una función no

lineal en x.

Nos proponemos resolver

(CMNL)  min_{x∈ℝ^{n}} f(x) = min_{x∈ℝ^{n}} \frac{1}{2}R^{T} R = min_{x∈ℝ^{n}} ∑_{i=1}^{m} r_{i}^{2}(x).  (2.1)

El problema surge generalmente cuando se intenta ajustar m datos (t_{i}, y_{i}) i = 1...m, por un modelo m(x, t) no lineal en x. En este caso r_{i}(x) = m(x, t_{i}) − y_{i} son las componentes de la función residual o simplemente residuo R(x).

Se elige la norma euclídea para minimizar el residuo porque de esta forma el problema está bien planteado matemáticamente y por su utilidad en aplicaciones a campo de la Estadística. Generalmente m >> n. En el caso particular en que m = n, este problema puede pensarse como un sistema no lineal y para todo m es un problema de minimización sin restricciones. Sin embargo, no es aconsejable tratarlo con técnicas tan generales que no tienen en cuenta la estructura particular del problema y esto justifica el desarrollo de técnicas específicas.

Los métodos básicos que resuelven el problema CMNL requieren información sobre las derivadas de las funciones r_{i}(x). Supondremos a partir de ahora que r_{i} ∈ C^{2}(ℝ).

Comencemos escribiendo las derivadas de R(x) y de f(x) = R^{T} (x)R(x). Denotemos con J(x) a la matriz Jacobiana de R(x), es decir,

J(x)_{i,j} = \frac{∂r_{i}}{∂x_{j}} , i, j = 1, ..., n.

(∇f)_{i} = \frac{∂}{∂x_{i}} f(x) = \frac{∂}{∂x_{i}} ∑_{j=1}^{m} \frac{1}{2}r_{j}^{2}(x) = ∑_{j=1}^{m}(\frac{∂}{∂x_{i}} r_{j}(x))r_{j}(x) =

∑_{j=1}^{m} J(x)_{i,j}r_{j}(x) = (J^{T} R)_{i}(x), i = 1, ..., n.

De manera análoga puede comprobarse que

∇^{2}f = ∑_{i=1}^{m}(∇r_{i}(x)∇r_{i}(x)^{T} + r_{i}(x)∇^{2}r_{i}(x)) = J(x)^{T} J (x) + S(x) (2.2)

donde

S(x) = ∑_{i=1}^{m} r_{i}(x)∇^{2}r_{i}(x).

Las formas especiales de ∇f(x) y de ∇^{2}f(x) pueden explotarse en los métodos que resuelven el problema que nos ocupa.

Se sabe que una condición necesaria para que x_{∗} sea un mínimo local de f(x) es que sea un punto crítico de f(x), es decir, ∇f(x_{∗}) = J(x_{∗})^{T} R(x_{∗}) = 0. Estableceremos ahora una condición suficiente para que un punto crítico x_{∗} sea un mínimo local de f(x), siguiendo la idea de aproximación geométrica de Wedin [44].

Si asumimos que J(x_{∗}) es de rango completo, entonces, J(x_{∗})^{+} J(x_{∗}) = I_{n} donde J(x_{∗})^{+} es la matriz pseudoinversa de J(x_{∗})y podemos reescribir (2.2) como

∇^{2}f = J^{T} J − γG_{ω} = J^{T} (I − γ(J^{+})^{T} G_{ω}J^{+})J  (2.3)
con G_{ω} = ∑_{i=1}^{m} ω_{i}∇^{2}r_{i},  ω = −\frac{R}{∥R∥_{2}} ,  γ = ∥R∥_{2} ,

siendo evaluadas las cantidades en x_{∗}.

La matriz simétrica K = (J^{+})^{T} G_{ω}J^{+} se llama matriz de curvatura de la superficie n-dimensional R(x) ∈ ℝ^{m} con respecto al vector normal ω.

Sean κ_{1} ≥ κ_{2≥}... ≥ κ_{n} los autovalores de K. Las cantidades ρ_{i} = \frac{1}{κ_{i}} , κ_{i} ≠ 0 son los radios principales de curvatura de la superficie con respecto al vector normal ω.

Notemos de (2.3) que ∇^{2}f(x_{∗}) es definida positiva y x_{∗} es un mínimo local sí y sólo sí 1 − γκ_{1} > 0.

Si 1 − γκ_{1} ≤ 0 , el problema de cuadrados mínimos tiene un punto de ensilladura en x_{∗} , y si 1 − γκ_{n} < 0, un máximo local en x_{∗}.

En el caso de ajuste de datos, el problema que nos ocupa se interpreta geométricamente, como el problema de hallar el punto de la superficie más cercano al vector de observación y ∈ ℝ^{m} [10].

Existen básicamente dos puntos de vista distintos para resolver el problema. Se puede pensar que proviene de un sistema de ecuaciones no lineales sobredeterminado

r_{i}(x) = 0 i = 1, 2, ...m

y en este caso es natural aproximar R(x) por un modelo afín M(x) en un entorno de x. El sistema

M(x) = 0

es sobredeterminado y puede resolverse en el sentido de los problemas de cuadrados mínimos no lineales. Esta idea es la base de los métodos de tipo Gauss-Newton que se caracterizan porque sólo usan información sobre las derivadas primeras.

El otro punto de vista consiste en pensar el problema como un caso especial de optimización en ℝ^{n} y resolverlo minimizando una aproximación cuadrática de f(x). De aquí provienen los métodos de tipo Newton.

Daremos a continuación una descripción de ambos .

2.1 Métodos de tipo Gauss-Newton

Como ya mencionamos, estos métodos iterativos se basan en una aproximación afín de R(x) alrededor del iterado actual x_{c}

M_{c}(x) = R(x_{c}) + J(x_{c})(x − x_{c})  (2.4)
M_{c} : ℝ^{n} → ℝ^{m} m > n.

Es claro que M(x) = 0 es un sistema sobredeterminado y por lo tanto, en general no es posible encontrar una solución no trivial, es por ello que al nuevo iterado x_{+}, se lo define como la solución del problema CML :

min_{x∈ℝ^{n}} \frac{1}{2} ∥M_{c}(x)∥_{2}^{2} .  (2.5)

Como mencionamos en el capítulo anterior, si J_{c}(x) tiene rango completo, la solución es

x_{+} = x_{c} − [J(x_{c})^{T} J(x_{c})]^{−1}J(x_{c})^{T} R(x_{c}).

Por cuestiones numéricas en general se evita el cálculo explícito de matrices inversas, razón por la cual para el cómputo de x_{c} no se utiliza esta fómula sino la descomposición QR de J(x_{c}).

La convergencia del método iterativo definido por este esquema está garantizada por el siguiente teorema [9] :

Teorema 10 Sea R : ℝ^{n} → ℝ^{m}, y f(x)=\frac{1}{2}R(x)^{T} R(x) ∈ C^{2}(D), con D ⊂ ℝ^{n}, abierto convexo. Si J(x) ∈ Lip_{γ}(D), con ∥J(x)∥_{2} ≤ α para todo x ∈ D, y existen x_{∗} ∈ D y λ, σ ≥ 0, tal que J(x_{∗})^{T} R (x_{∗}) = 0, λ es el menor autovalor de J(x_{∗})^{T} J(x_{∗}),

∥(J(x) − J(x_{∗}))^{T} R(x_{∗})∥_{2} ≤ σ ∥x − x_{∗}∥_{2}  (2.6)

para todo x ∈ D. Si σ < λ, entonces para todo c ∈ (1, \frac{λ}{σ}) existe ϵ > 0 tal que para todo x_{0} ∈ N(x_{∗}, ϵ), la sucesión generada por el método de Gauss Newton

x_{k+1} − x_{k} − (J(x_{k})^{T} j(x_{k}))^{−1}J(x_{k})^{T} R(x_{k})

está bien definida, converge a x_{∗}, y verifica:

∥x_{k+1} − x_{∗}∥_{2} ≤ \frac{cσ}{λ} ∥x_{k} − x_{∗}∥_{2} + \frac{cσγ}{2λ} ∥x_{k} − x_{∗}∥_{2}^{2}

y

∥x_{k+1} − x_{∗}∥_{2} ≤ \frac{cσ + λ}{2λ} ∥x_{k} − x_{∗}∥_{2} < ∥x_{k} − x_{∗}∥_{2} .

Corolario 11 Bajo las hipótesis el teorema anterior, si R(x_{∗}) = 0, entonces existe ε > 0 tal que para todo x_{0} ∈ N(x_{∗}, ε), la sucesión {x_{k}} generada por el método de Gauss Newton está bien definida y converge q-cuadráticamente a x_{∗}.

Notar que de acuerdo al teorema el método de Gauss Newton es menos satisfactorio que muchos métodos locales conocidos para resolver problemas de minimización sin restricciones, y esto se debe a que en muchos casos aún siendo J(x_{∗})^{T} J(x_{∗}) no singular, la convergencia es lenta o no hay convergencia y depende fuertemente de la medida de S(x_{∗}).

La constante σ juega un papel fundamental. Puede probarse que para x suficientemente próximo a x_{∗} resulta,

J(x_{∗})^{T} J(x_{∗})^{T} R(x_{∗}) ≈ S(x_{∗})(x − x_{∗})

por lo que σ puede interpretarse como ∥S(x_{∗})∥ . De esta interpretación y de (2.6) podemos decir que que σ es una medida absoluta de la no linealidad del problema y del tamaño del residuo, pues si R(x) es lineal o R(x_{∗}) = 0 de (2.6) se concluye que σ = 0. El hecho que debe ser \frac{σ}{λ} < 1 para garantizar convergencia, permite interpretar este cociente como una medida relativa de la no linealidad y del tamaño del residuo. Por la tanto, el teorema asegura que la rapidez del método de Gauss-Newton decrece con la relativa no linealidad o con el crecimiento relativo del residuo; si alguno de estos valores es muy grande el método puede no converger.

Aún cuando este tipo de métodos tiene varios inconvenientes, es la base de muchos algoritmos que resuelven el problema de cuadrados mínimos.

Más detalles sobre el teorema anterior se pueden encontrar en [9] y [2]. VENTAJAS Y DESVENTAJAS Entre las ventajas caben señalar:

• Convergencia local q-cuadrática en problemas de residuo nulo.

• Rápida convergencia local, q-lineal en problemas que no son ‘demasiado’ no lineales y con residuos ‘razonables’.

• Resuelven problemas de cuadrados mínimos lineales en una iteración.

Entre las desventajas de los métodos de tipo Gauss Newton, figuran:

• Convergen localmente en forma lenta y q-lineal en problemas ‘suficientemente’ no lineales y con residuos ‘razonablemente’ grandes.

• No convergen localmente en problemas ‘muy no lineales’ ni en aquellos que tienen residuos muy grandes.

• No están bien definidos si J(x_{c}) no es de rango completo.

• No son necesariamente globalmente convergentes.

Este método ha sido adaptado convenientemente para el caso en que la matriz jacobiana es sparse [31]. Las características principales de este algoritmo, son las siguientes: en cada iteración el problema

J(x_{k})z_{k} = −R(x_{k})

es resuelto parcialmente usando el método de gradientes conjugados, precondicionando el sistema y conservando el mismo precondicionador durante varias iteraciones. El nuevo iterado se calcula en el subespacio generado por z_{k} y ∇f(x_{k}) . Se consigue así un método globalmente convergente bajo ciertas hipótesis de compacidad.

Hasta ahora hemos trabajado con la hipótesis J(x_{c}) de rango completo, en el caso que no lo fuera, la actualización del iterado se define como sigue:

x_{+} = x_{c} − J(x_{c})^{+}R(x_{c}).

Es claro que para el caso en que la matriz Jacobiana es de rango completo, el método de Gauss Newton define una dirección de descenso, y esto sugiere básicamente dos modificaciones posibles para adaptarlo al caso general obteniéndose las siguientes estrategias:

1. El método de Gauss Newton amortiguado.

2. El método de Levenberg-Marquardt.

El primero de ellos combina el método de Gauss Newton con una estrategia de búsqueda lineal, definiendo

x_{k+1} = x_{k} + α_{k}p_{k}

donde p_{k} es el paso de Gauss-Newton y α_{k}, longitud del paso, es un valor a determinar.

Cuando J_{k} es de rango deficiente, p_{k} puede elegirse como la solución de norma mínima del problema (CML)

min_{p} ∥f(x_{k}) + J(x_{k})p∥_{2} , p ∈ ℝ^{n}

es decir,

p_{k} = −J^{+}(x_{k})f(x_{k}).

Para que el algoritmo resulte viable, la longitud del paso debe elegirse cuidadosamente, siendo los criterios más usuales los que sugieren examinar que el par (α_{k}, p_{k}) verifique una condición de tipo Armijo para asegurar suficiente reducción :

∥R(x_{k})∥_{2}^{2} − ∥R(x_{k} + α_{k}p_{k})∥_{2}^{2} ≥ \frac{1}{2}α_{k} ∥J(x_{k})p_{k}∥_{2}^{2} .

Otra opción es elegir α_{k} como una solución aproximada del problema unidimensional:

min_{α} ∥R(x_{k} + α_{k}p_{k})∥_{2}^{2} .

Lindström and Wedin [28] propusieron un algoritmo que define una aproximación p(α) a la curva R(x_{k} + αp_{k)} en ℝ^{m} y se minimiza ∥p(α)∥_{2} .

La dificultad de los métodos de tipo Gauss Newton amortiguado reside en el hecho que en muchas ocasiones, cuando la matriz Jacobiana no es de rango completo, en algunos pasos intermedios se presentan dificultades numéricas, en tales caso es aconsejable intentar resolver el problema con otra estrategia de globalización: la estrategia de región de confianza. Esta es la propuesta del algoritmo de Levenberg-Marquardt, en el cual la dirección p_{k} se computa como la solución de

min_{∥p_{k}∥≤δ_{k}} ∥R(x_{k}) + J(x_{k})p∥_{2}^{2} .

Puede probarse que la solución es

p_{k} = −(J_{k}^{T} J_{k} + μ_{k}I)^{−1}J_{k}^{T} R(x_{k})

con μ_{k} = 0 si δ_{k} ≥ ∥(J_{k}^{T} J_{k})^{−1}J_{k}^{T} R_{k}∥ y μ_{k} > 0 en otro caso.

Este método ha sido codificado utilizando diversas estrategias para la elección de μ_{k} y δ_{k}, siendo Moré [33] quien dió una implementación como un algoritmo de región de confianza.

El siguiente teorema garantiza la convergencia local del método de Levenberg-Marquardt [9] :

Teorema 12 Bajo las hipótesis del teorema anterior, si la sucesión de numeros {μ_{k}} reales no negativos, está acotada por b>0, y σ < λ, entonces, para cualquier c ∈ (1, (λ + b)/(σ + b)), existe ε > 0, tal que para todo x_{0} ∈ N(x_{∗}, ε), la sucesión generada por el método de Levenberg-Marquardt

x_{k+1} = x_{k} − (J(x_{k})^{T} J(x_{k}) + μ_{k}I)^{−1}J(x_{k})^{T} R(x_{k})

está bien definida y verifica:

∥x_{k+1} − x_{∗}∥_{2} ≤ \frac{c(σ + b)}{(λ + b)} ∥x_{k} − x_{∗}∥_{2} + \frac{cαγ}{2(λ + b)} ∥x_{k} − x_{∗}∥_{2}^{2}

y

∥x_{k+1} − x_{∗}∥_{2} ≤ \frac{c(σ + b) + (λ + b)}{2(λ + b)} ∥x_{k} − x_{∗}∥_{2} < ∥x_{k} − x_{∗}∥_{2}^{2} .

Si R(x_{∗}) = 0 y μ_{k} = O(∥J(x_{k})^{T} R(x_{k})∥_{2}), entonces {x_{k}} converge qcuadráticamente a x_{∗}.

Si bien el teorema sugiere que el algoritmo de Levenberg-Marquardt puede converger lentamente cuando el residuo es grande o en problemas ‘muy no lineales’, en la práctica la implementación de Moré es satisfactoria y una de las razones es que el algoritmo está bien definido aún en los casos de Jacobianos de rango deficiente. En los casos en que Gauss Newton provee un paso demasiado largo, este algoritmo encuentra una dirección próxima a la de máximo descenso.

Por otra parte Powell (1975) [37], Osborne (1976) y Moré [33] (1977) probaron que el método es globalmente convergente.

2.2 Métodos de tipo Newton

El análisis de la sección previa mostró que para problemas con residuos ‘grandes’ o ‘fuertemente no lineales’ la convergencia de los métodos de tipo Gauss-Newton puede ser muy lenta, además pueden tener problemas en puntos donde el Jacobiano tiene rango deficiente. Es por eso que cuando se dispone de información de las derivadas segundas de f(x), puede resultar conveniente la aplicación de Métodos del tipo Newton. Estos proponen minimizar el modelo cuadrático sugerido por la serie de Taylor de f(x) alrededor del iterado actual x_{k}

m_{k}(x) = \frac{1}{2}R(x_{k})^{T} R(x_{k}) + R(x_{k})^{T} J(x_{k})(x − x_{k})

+\frac{1}{2}(x − x_{k})^{T} (J(x_{k})^{T} J(x_{k}) + S(x_{k}))(x − x_{k}).

En este caso el punto crítico del modelo es un mínimo bajo las siguientes condiciones:

• ∇^{2}f(x) es Lipschitz continua en un entorno de x_{k}.

• ∇^{2}f(x_{∗}) definida positiva.

Las propiedades de convergencia local del método de Newton aplicado al problema CMNL son superiores a las de otros métodos, sin embargo el inconveniente de este método radica en el hecho que algunas veces no se dispone del término S(x_{k}) ya que raramente puede calcularse a un costo razonable.

Si no se dispone de la expresión analítica de J_{k}, puede aproximarse por diferencias finitas para obtener ∇f_{k} = J_{k}^{T} R_{k} . Esto requiere n evaluaciones del residuo por iteración, pero el costo de aproximar ∇^{2}f_{k} o S_{k}, es del orden de n^{2}, lo cual lo torna inaceptable.

Una aplicación directa de un método secante no es deseable, porque éstas aproximan ∇^{2}f_{k}, y en los problemas de cuadrados mínimos el término J_{k}^{T} J_{k} está siempre disponible pues se lo calcula para obtener ∇f_{k}. Por esta razón sólo es necesario aproximar S(x_{c}). Una alternativa es hallar una aproximación secante para cada uno de sus sumandos, ∇^{2}r_{i}(x_{k}), pero esto no es adecuado por el requerimiento de memoria que significa almacenar m matrices de orden n.

Hay propuestas que parcialmente tienen en cuenta este término de segundo orden. Por ejemplo, Gil y Murray sugieren aproximar la matriz Hessiana por J(x_{k})^{T} J(x_{k}) en el subespacio invariante correspondiente a los valores singulares más grandes de J(x_{k}) y en su complemento consideran el término S(x_{k}), como se indica a continuación.

Sea la descomposición en valores singulares de J(x_{k}) :

J(x_{k}) = U\binom{∑}{0}V^{T} U ∈ ℝ^{m×m} V ∈ ℝ^{m×n}

donde U y V son matrices ortogonales y ∑ = diag(σ_{1}, σ_{2}, ..., σ_{n}) es la matriz de los valores singulares ordenados en forma decreciente.

Entonces la ecuación que define la dirección de Newton puede escribirse de la forma:

(Σ^{2} + V^{T} S_{k}V)q = −Σ R  (2.7)

donde q = V^{T} (x_{+} − x_{c}) y R denota las n componentes del vector U^{T} R(x_{k}).

Si se dividen los valores singulares en dos grupos, se puede escribir ∑ = diag(∑_{1}, ∑_{2}), donde ∑_{1} = diag(σ_{1}, σ_{2}, ..., σ_{r}) son los valores singulares más grandes. Particionando V, q y R adecuadamente (2.7) puede escribirse como

(Σ_{1}^{2} + V_{1}^{T} S_{c}V_{1})q_{1} + V_{1}^{T} S_{c}V_{2}q_{2} = −Σ_{1} R_{1} .

Si los términos de S_{k} son despreciables frente a los de Σ_{1}^{2}q_{1}, se obtiene

q_{1} = −Σ_{1}^{−1} R_{1} .

Si se sustituye en las (n − r) ecuaciones restantes se puede resolver para q_{2}:

(Σ_{2}^{2} + V_{2}^{T} S_{c}V_{2})q_{2} = −Σ_{2} R_{2} −V_{2}^{T} S_{c}V_{1}q_{1}

y el paso de Newton viene dado por

p_{k} = V_{q} = V_{1}q_{1} + V_{2}q_{2}.

La descomposición en valores singulares se actualiza en cada iteración, y la idea es mantener el valor de r próximo a n.

Otra alternativa posible es reemplazar cada matriz ∇^{2}r_{i}(x) por una aproximación secante. Sin embargo, esto no es recomendable por el requerimiento de memoria que ello implica. La idea más acertada es aproximar ∇^{2}f(x_{k}) por J_{k}^{T} J_{k} + A_{k}, donde A_{k} es una aproximación secante de S(x_{k}), por ejemplo las actualizaciones BFGS o DFP [7].

Cuadrados mínimos no lineales con restricciones

En este capítulo nos dedicaremos a describir el problema general de cuadrados mínimos no lineales con restricciones y presentaremos un nuevo algoritmo para el caso en que las restricciones sean cotas en las variables.

La formulación general del problema es la siguiente:

Sean R : ℝ^{n} → ℝ^{m}, ḣ : ℝ^{n} → ℝ^{p}, p < n, se desea resolver:

min_{x∈ℝ^{n}} f(x) = min_{x∈ℝ^{n}} \frac{1}{2}R^{T} R = min_{x∈ℝ^{n}} ∑_{i=1}^{m} r_{i}^{2}(x) s.a. h(x) = 0.  (3.1)

Este problema puede resolverse linealizándolo y generalizar el método de Gauss Newton al caso con restricciones, es decir, en cada iteración se resuelve

min_{p} ∥R(x_{k}) + J(x_{k})p∥_{2}  s.a. h(x_{k}) + C(x_{k})p = 0  (3.2)

donde J y C son los jacobianos de R y h respectivamente, utilizando alguno de los métodos descriptos en el capítulo anterior para el caso lineal.

Puede probarse que la dirección p_{k} obtenida en (3.2) es una dirección de descenso de la función de mérito

ψ(x_{μ}) = ∥R(x)∥_{2}^{2} + μ ∥h(x)∥_{2}^{2} ,

en el punto x_{k} para μ suficientemente grande. Esto posibilita estabilizar el método de Gauss Newton con alguna estrategia de globalización.

Este problema ha sido tratado por diversos autores para distintos tipos de restricciones. Entre ellos cabe destacar el trabajo de Holt y Fletcher [25] en el que estudian el caso de variables acotadas y monótonas, es decir:

min_{x∈ℝ^{n}} f(x) = min_{x∈ℝ^{n}} \frac{1}{2}R^{T} (x)R(x)
s.a. x_{i} ≥ l_{i}  i ∈ I_{1}
x_{i} ≤ u_{i}  i ∈ I_{2}
x_{i} ≤ x_{i+1} i = 1, ..., k − 1

r_{i}(x) ∈ C^{2}, I_{1}, I_{2} ∈ {1, 2, ..., n} k ≤ n.

Para este caso, la propuesta es linealizar el residuo y en cada iteración identificar las restricciones activas en el iterado actual para resolver un problema de tipo CML sin restricciones en las variables libres. Es en este subespacio donde en general se busca la dirección de Levenberg-Marquardt, excepto cuando la predicción del modelo aconseje redefinir el conjunto activo. En este trabajo se adoptó una nueva elección del parámetro λ, que combina básicamente las propuestas de Marquardt [32] y Fletcher [11] para adecuarlas al caso con restricciones. Un problema clave en este contexto resulta de la posible variación de la dimensión del conjunto activo en cada iteración. Esto se manifiesta en el cambio sustancial del Jacobiano y en consecuencia, de la descomposición en valores singulares. Así, un adecuado valor del parámetro en una iteración puede no serlo en la siguiente. En este trabajo, se propone mantener en lo posible la medida de la región en la cual se busca el mínimo, independientemente de la variación del conjunto activo, y en función de este valor se estima λ.

Este algoritmo fue generalizado por Wright y Holt [45] para el caso de restricciones lineales de desigualdad y para este caso también propusieron una solución basada en el método de Bunch y Nielsen. La diferencia sustancial entre ambos algoritmos radica en que en el primero se encuentra una solución aproximada del subproblema mientras que en el segundo se halla la solución exacta. Para problemas de dimensiones pequeñas (n ≤ 10, m ≤ 40) el primero resulta ser significativamente más rápido. En problemas de grandes dimensiones, como la cantidad de operaciones en común que tienen ambos algoritmos es mucho más grande que las particulares de cada uno, se recomienda el segundo, ya que usualmente se encuentra el óptimo en menos iteraciones. Ambos tiene incorporada una técnica de ‘anti-zigzag’ para evitar la oscilaciones de los iterados en distintas variedades con el consiguiente consumo inútil de iteraciones.

Finalmente cabe señalar el trabajo de Lindstrom y Wedin [29] en el que tratan el problema

min_{x∈ℝ^{n}} f(x) = min_{x∈ℝ^{n}} \frac{1}{2}R^{T} (x)R(x)
s.a. C_{j}(x) = 0 i = 1, ..., q
C_{j}(x) ≥ 0 i = q + 1, ..., l.

La idea básica en esta propuesta es adecuar el método de Gauss Newton para el caso en que la matriz Jacobiana sea de rango deficiente o la aproximación lineal no sea adecuada, lo que puede ocurrir lejos de la solución o cuando el radio de curvatura de la superficie es muy pequeño. Sabemos que la dirección de Gauss Newton se computa mediante una descomposicón QR, en el primero de los casos mencionados, se aplica la técnica conocida como ‘Gauss Newton Truncado’ que consiste en resolver el problema eliminando las n − s últimas columnas de la matriz R en las cuales los elementos r_{ss} de la diagonal (s > 1) son muy pequeños en relación a r_{11}. En el segundo caso, se utiliza un método de tipo Newton.

El problema de cuadrados mínimos no lineales con restricciones de caja

4.1 Descripción del problema

Sean R(x) = (r_{1}(x), ...r_{m}(x)), l_{i} ≤ u_{i} (i = 1, .., n), m > n, el problema de cuadrados mínimos no lineales con restricciones de caja es enunciado de la siguiente manera:

min_{l≤x≤u} f(x) = min_{l≤x≤u} \frac{1}{2}R^{T} (x)R(x) = min_{l≤x≤u} ∑_{i=1}^{m} r_{i}^{2}(x)

donde r_{i}(x) ∈ C^{1}(ℝ^{n}).

Nuestro propósito es describir un algoritmo que aplica programación cuadrática sucesiva y en muchos aspectos generaliza la propuesta que para el caso sin restricciones desarrollaron Dennis Gay y Welsh [8], ya que el modelo cuadrático puede escogerse entre uno de tipo Gauss Newton y uno Cuasi Newton, según el decrecimento que proporcionen en la función objetivo, cuando el punto de prueba no es suficientemente satisfactorio.

Se utiliza región de confianza como estrategia de globalización. A pesar que en la definición se usa una norma arbitraria, en la implementación utilizamos la norma l_{∞}, por el tipo de restricciones que aparecen en el problema.

4.1.1 Definiciones y notación.

Por comodidad se resumen las principales definiciones necesarias para entender la descripción del algoritmo y la notación que utilizaremos.

Región de factibilidad del problema, el politopo Ω,

Ω = {x ∈ ℝ^{n}/ l ≤ x ≤ u}.

Si I_{l} e I_{u} son los conjuntos disjuntos de índices que corresponde a las restricciones activas en las cotas inferiores y superiores respectivamente, definimos I = I_{l} ∪ I_{u} y F_{I} la cara abierta de Ω asociada al conjunto I :

F_{I} = \{ \begin{matrix}x ∈ Ω/x_{i} = l_{i} si i ∈ I_{l}, x_{i} = u_{i} si i ∈ I_{u}, \\ l_{i} < x_{i} < u_{i} en caso contrario \end{matrix} \} .

Si q(x) es una función definida en Ω, para cada punto ω ∈ Ω se define el gradiente proyectado de ω al vector

∇_{p}q(ω)_{i} = \{ \begin{matrix}0 & si & ω_{i} = l_{i} y ∇q(ω)_{i} > 0 \\ 0 & si & ω_{i} = u_{i} y ∇q(ω)_{i} < 0 \\ ∇f(ω)_{i} & & en caso contrario. \end{matrix}

Para cada punto x ∈ F_{I} , se divide al gradiente proyectado en dos componentes, el gradiente interno, ∇_{i}q(ω) paralelo a la cara y el gradiente cortado, ∇_{c}q(ω) perpendicular a la misma, de modo que ∇_{p}q(ω) = ∇_{i}q(ω) + ∇_{c}q(ω).

Las condiciones necesarias y suficientes para que un punto x_{∗} sea un mínimo global de una función cuadrática q(x) en Ω, pueden sintetizarse en la condición citada por Friedlander y Martínez [12],

∇_{p}q(x_{∗}) = 0.

Con δ denotamos el radio de la región de confianza.

4.2 Algoritmo básico

El método que proponemos puede sintetizarse básicamente en el siguiente algoritmo.

Algoritmo Dados x_{k}, δ_{k} > 0, α > 1, β < 1, repetir hasta convergencia: Paso 1 Construir un modelo cuadrático q_{k}(x) de f(x) alrededor de x_{k}

q_{k}(x) = f(x_{k}) + (x − x_{k})^{T} g_{k} + \frac{1}{2}(x − x_{k})^{T} B_{k}(x − x_{k}).

Paso 2 Resolver aproximadamente

s_{k} = arg min_{l−x_{k}≤s≤u−x_{k}_{∥x∥_{∞}≤δ_{k}.}} q_{k}(s)

Paso 3 Evaluar el paso de prueba s_{k}.

Si el paso es adecuado definir: x_{k+1} = x_{k} + s_{k}

δ_{k+1} ∈ [δ_{k}, αδ_{k}]
k = k + 1

Ir a 1.

En caso contrario, δ_{k} = βδ_{k} Ir a 3.

Estos tres pasos sintetizan los puntos cruciales de nuestra propuesta, que pasamos a detallar a continuación.

4.2.1 Construcción del modelo cuadrático.

Se inicia el cálculo con un punto arbitrario, que en caso de no ser factible es proyectado dentro del politopo. En la primera iteración se construye un modelo afín de la función residual R(x) (método de Gauss Newton). En iteraciones posteriores, siguiendo las ideas de Dennis, Gay y Welsh, este modelo puede cambiarse por un modelo cuadrático de la función objetivo f(x) (método cuasi Newton), y en este caso, el término del matriz Hessiana que contiene información de segundo orden, al que se denomina generalmente término ‘no computable ’, se aproxima con una actualización de tipo BFGS. Una vez escogido un modelo lo llamamos preferencial y al otro, alternativo.

Puede utilizarse gradiente exacto o aproximado.

4.2.2 Resolución del subproblema.

En el Paso 3, se resuelve el subproblema cuadrático en un conjunto compacto, el cual resulta de la intersección del conjunto factible con la región de confianza.

Como en todos lo métodos modernos que aplican región de confianza, no es necesario resolverlo exactamente para garantizar la convergencia del algoritmo principal.

Como se verá más adelante es suficiente para conseguir un algoritmo bien definido, pedir que la solución aproximada satisfaga la condición

q_{k}(s_{k}) ≤ β_{1} min {q_{k}(s) : s = γg_{k}, ∥s∥ ≤ δ_{k}}  (4.1)

donde g_{k} es el gradiente de la función objetivo, o una aproximación del mismo, calculado en el punto x_{k}.

En el caso en que l = −∞, u = +∞, la norma con la que se trabaje en la definición de la región de confianza sea la norma l_{2}, y se tome como dirección de descenso la dirección opuesta del gradiente, esta condición requerida es la conocida como Fracción de Decrecimento en la dirección de Cauchy.

4.2.3 Aceptación del paso.

Una vez obtenido un punto de prueba es fundamental la decisión de aceptarlo o no. Para ello verificamos si es un punto de Khun Tucker de nuestro problema, si no lo es, se computa otro punto y se elige aquél que proporciona al menos, un decrecimiento igual al que minimiza el modelo preferencial a lo largo de cualquier dirección de descenso dentro de la región de factibilidad de esa etapa.

Estos pasos pueden resumirse en el algoritmo.que se describe a continuación.

4.2.4 El algoritmo.

Sea s la solución del subproblema cuadrático, para decidir si se acepta como nuevo iterado, es necesario estimar de alguna manera, cuán bueno es el modelo cuadrático propuesto en la región en la cual se trabaja. Para ello se comparan dos cantidades: la reducción actual en la función objetivo,

ared_{k} = f(x_{k}) − f(x_{k} + s)  (4.2)

y la reducción predicha por el modelo

pred_{k} = f(x_{k}) − q_{k}(s)  (4.3)

y básicamente se procede como se detalla a continuación:

Dados x_{k}, s, ared_{k}, pred_{k}, δ_{k}, η_{1} > 0, η_{2} > 0, α > 1, β < 1

Si ared ≥ η_{1} ∗ pred

( La iteración es satisfactoria y posiblemente se incremente el radio) Se define δ′ = αδ_{k}. Se calcula un nuevo paso de prueba s′. Si f(x_{k} + s′) < f(x_{k} + s)

x_{k+1} = x_{k} + s′
δ_{k+1} = δ′

sino

x_{k+1} = x_{k} + s
δ_{k+1} = δ_{k}
k = k + 1.

Si ared < η_{1} ∗ pred

(El paso no es suficientemente bueno.)

Si es la primera vez que se testea el paso en esta

iteración :

Se calcula un punto de prueba x′′

con el modelo alternativo

Se comparan las predicciones de ambos modelos.

Se escoge el modelo que mejor predice y

se acepta el paso proporcionado por dicho modelo. k = k + 1

δ_{k} = δ_{k}

Si no es la primera vez que se testea el paso

Si ared > η_{2} ∗ pred

(Se acepta el punto se reduce

el radio de la región de confianza).

x_{k+1} = x_{k} + s

δ_{k+1} = β ∗ δ_{k}

k = k + 1

Si no

(Se rechaza el punto se reduce el

radio de la región de confianza).

δ_{k} = β ∗ δ_{k}.

4.3 Convergencia del algoritmo

En esta sección se probará que bajo ciertas hipótesis, el algoritmo está bien definido y genera una sucesión de iterados que converge a un punto estacionario del problema.

Ambas afirmaciones pueden deducirse de la teoría general de algoritmos de minimización sin restricciones que utiliza región de confianza como estrategia de globalización. En efecto, por el tipo especial de restricciones del problema que nos ocupa, la solución aproximada del subproblema cuadrático se busca en la intersección del politopo y la región de confianza, por lo tanto, es una solución que verifica (4.1) en el conjunto {s/ ∥s∥ ≤ δ_{k}}, razón por la cual puede aplicarse la teoría antes mencionada.

A continuación se detallan resultados que garantizan que el algoritmo está bien definido, es decir, si x_{k} no es un punto estacionario, bajo ciertas hipótesis puede asegurarse que después de un número finito de iteraciones puede hallarse un paso satisfactorio que permite definir el iterado x_{k+1}.

Veamos alguno resultados preliminares [37], [34], [3], cuyas demostraciones repetimos con el fin de dar una presentación completa del tema.

El siguiente es un resultado clásico de la teoría de región de confianza, mencionado por Moré en [34].

Lema 13 Si s_{k} verifica 4.1 y si ∥.∥ = ∥∥_{2} ,entonces

pred_{k}(s_{k}) = −q_{k}(s_{k}) ≥ \frac{1}{2} ∥g_{k}∥ min \{δ_{k}, \frac{∥g_{k}∥}{∥B_{k}∥}\} .

Demostración: Sean φ : ℝ → ℝ, definida por

φ(t) = q_{k} (t\frac{g_{k}}{∥g_{k}∥}) = −t ∥g_{k}∥ + \frac{1}{2}t^{2}μ

con μ = \frac{g_{k}^{T} B_{k}g_{k}}{∥g_{k}∥} y t^{∗} el mínimo de φ en [0, δ_{k}].

Observar que si x_{k} no es punto estacionario del problema, no puede ser t^{∗} = 0.

Si t^{∗} ∈ (0, δ_{k}) , es t^{∗} = \frac{∥g_{k}∥}{μ}, y en este caso,

φ(t^{∗}) = −\frac{1}{2} \frac{∥g_{k}∥}{μ} ≤ −\frac{1}{2} \frac{∥g_{k}∥}{∥B_{k}∥}.  (4.4)

Si t^{∗} = δ_{k}, la función es decreciente en [0, δ_{k}] , por lo tanto,

φ^{,}(t) = − ∥g_{k}∥ + tμ ≤ 0,

de donde resulta δ_{k}μ ≤ ∥g_{k}∥ y

φ(t^{∗}) ≤ −δ_{k} ∥g_{k}∥ + \frac{1}{2}δ_{k}^{2}μ ≤ −δ_{k} ∥g_{k}∥ .  (4.5)

Como por hipótesis q_{k}(s_{k}) ≤ β_{1}φ(t^{∗}), de (4.4) y (4.5) resulta la tesis. Carter [3], estableció condiciones suficientes sobre la aproximación del

gradiente que permiten asegurar la convergencia del algoritmo.

Veamos un lema previo:

Lema 14 Sea f(x) : ℝ^{n} → ℝ una función continuamente diferenciable, con gradiente Lipschitziano con constante L. Sean ared_{k} y pred_{k} las funciones definidas en (4.2) y (4.3) respectivamente, λ_{k}^{min} y λ_{k}^{max} el menor y mayor autovalor respectivamente de B_{k}. Si se define el error

e_{k} = g_{k} − ∇f_{k}

entonces, para todo s ∈ ℝ^{n} tal que x_{k} + s es factible, se verifica

−\frac{1}{2} ∥s∥^{2} (L + λ_{k}^{max})−e_{k}^{T} s ≤ pred_{k}(s)−ared_{k}(s) ≤ −\frac{1}{2} ∥s∥^{2} (L + λ_{k}^{min})−e_{k}^{T} s.

Demostración: Utilizando la representación integral de ared_{k}(s), es:

pred_{k}(s) − ared_{k}(s) = −g_{k}^{T} s − \frac{1}{2}s^{T} B_{k}s + ∫_{0}^{1} ∇f(x_{k} + λs)^{T} sdλ =

−e_{k}^{T} s − \frac{1}{2}s^{T} B_{k}s + ∫_{0}^{1} [∇f(x_{k} + λs) − ∇f(x_{k})]^{T} sdλ.  (4.6)

Como λ_{k}^{min} ∥s∥^{2} ≤ s^{T} B_{k}s ≤ λ_{k}^{max} ∥s∥^{2} y

∫_{0}^{1} [∇f(x_{k} + λs) − ∇f(x_{k})]^{T} sdλ ≤ ∫_{0}^{1} L ∥λs∥ ∥s∥ dλ = \frac{1}{2}L ∥s∥^{2} ,

sustituyendo en (4.6), resulta la tesis.

A continuación, utilizando estos resultados probaremos que bajo ciertas hipótesis siempre puede encontrarse un paso satisfactorio.

Denotaremos con θ^{i} al ángulo entre s^{i} y −g_{k}, por lo tanto,

cos θ^{i} = \frac{−g_{k}^{T} s^{i}}{∥g_{k}∥ ∥s^{i}∥}.

Teorema 15 Bajo las hipótesis del lema anterior, si en la iteración k − ésima se genera una una sucesión {s^{i}}de pasos de prueba insatisfactorios asociados a radios δ^{i}, y lim_{δ^{i}→0} cos θ^{i} = 1, g_{k} ≠ 0 y \frac{∥e_{k}∥}{∥g_{k}∥} ≤ ζ,con ζ ∈ [0, 1 − η_{1}] ,

entonces para δ^{i} suficientemente pequeño se cumple que ρ^{i} = \frac{ared_{k}(s^{i})}{pred_{k}(s^{i})} > η_{1}.

Demostración: Es fácil comprobar que

1 − ρ^{i} = \frac{pred_{k}(s^{i}) − ared_{k}(s^{i})}{pred_{k}(s^{i})} ≤ \frac{\frac{1}{2}(L − λ_{k}^{min}) ∥s^{j}∥^{2} − e_{k}^{T} s^{i}}{−g_{k}^{T} s^{i′} − \frac{1}{2}s_{i}^{T} Bs^{i}} .  (4.7)

Utilizando la desigualdad de Schwartz y la definición de θ^{i}, (4.7) puede escribirse como:

1 − ρ^{i} ≤ \frac{∥e_{k}∥ ∥s^{i}∥ + \frac{1}{2} (L − λ_{k}^{min}) ∥s^{i}∥^{2}}{−g_{k}^{T} s^{i′} − \frac{1}{2}s_{i}^{T} Bs^{i}} = \frac{∥e_{k}∥ ∥s^{i}∥ + \frac{1}{2} (L − λ_{k}^{min}) ∥s^{i}∥^{2}}{− ∥g_{k}∥ ∥s^{i}∥ cos θ^{i} − \frac{1}{2}s_{i}^{T} Bs^{i}}

(4.8)

=\frac{∥e_{k}∥ + \frac{1}{2} (L − λ_{k}^{min}) ∥s^{i}∥}{− ∥g_{k}∥ cos θ^{i} − (\frac{1}{2}s_{i}^{T} Bs^{i}) / ∥s^{i}∥} ≤ \frac{1}{∥g_{k}∥} \frac{∥e_{k}∥ + \frac{1}{2} (L − λ_{k}^{min}) ∥s^{i}∥}{cos θ^{i} − \frac{\frac{1}{2} s_{i}^{T} Bs^{i}}{∥g_{k}∥∥s^{i}∥}}

Por lo tanto,

lim_{δ^{i}→0} 1 − ρ^{i} ≤ \frac{∥e_{k}∥}{∥g_{k}∥} ≤ ζ ≤ 1 − η_{1},

lo cual garantiza que para δ^{i} suficientemente pequeño, es ρ^{i} ≥ η_{1}.

Como consecuencia inmediata de este teorema, se puede asegurar que bajo las hipótesis mencionadas, el algoritmo genera una sucesión infinita de iterados {x_{k}} , o encuentra un punto con g_{k} = 0. Más aún, Carter [3] probó que si el subproblema es resuelto con las técnicas tradicionales de región de confianza [16], la condición de mantener el error relativo de las aproximaciones del gradiente uniformemente acotado, asegura que el algoritmo genera una sucesión infinita de puntos que satisfacen f(x_{k+1}) < f(x_{k}), o termina en un punto que verifica ∥∇f_{k}∥ = 0.

Si el algoritmo genera una sucesión infinita de iterados es deseable que ésta converja a un punto estacionario. Este hecho ha sido estudiado por diversos autores, y son destacables los trabajos de [37], [34] y [3].

A continuación se mencionan los resultados obtenidos por estos autores, en trabajos que se suceden cronológicamente y permiten comprobar cómo se ha logrado debilitar las hipótesis.

Powell probó los siguientes teoremas que se caracterizan fundamentalmente por utilizar gradiente exacto en el modelo cuadrático:

Teorema 16 Si f : ℝ^{n} → ℝ es una función acotada inferiormente, diferenciable y ∇f(x) uniformemente continuo, entonces lim inf_{k→∞} ∇f(x_{k}) = 0.

Teorema 17 Si x_{k} → x_{∗}, f(x) es de clase C^{2} en un entorno de x_{∗}, y B(x_{∗}) es definido positivo, entonces x_{k} → x_{∗} q-superlinealmente.

Por lo tanto, bajo las hipótesis de los dos últimos teoremas, se puede garantizar que si la sucesión converge, lo hace q-superlinealmente a un mínimo local.

Moré [34] generalizó estos resultados para el caso en que se trabaje con un problema escalado, es decir si el subproblema cuadrático es

min {φ_{k}(s); ∥Ds∥ ≤ δ_{k}}

con D matriz no singular, y φ_{k}(s) = g^{T} s + \frac{1}{2}s^{T} Bs modelo cuadrático de la función objetivo en un entorno del iterado x_{k}.

La diferencia sustancial con el trabajo de Powell radica en que muchos de los resultados son válidos aún cuando g es una aproximación del gradiente.

Si se supone que el paso s_{k} verifica la condición

φ_{k}(s_{k}) ≤ β_{1} min {φ_{k}(s) : D_{k}^{T} D_{k}s = νg_{k}, ∥Ds∥ ≤ δ_{k}} , ∥Ds_{k}∥ ≤ β_{2}δ_{k}

con β_{1} < 1, β_{2} ∈ [0, 2] , y las sucesiones ∥D_{k}∥ y ∥B_{k}∥ están uniformemente acotadas, el teorema 16 se generaliza afirmando, que lim inf_{k→∞} ∥g(x_{k})∥_{C_{k}^{−1}} = 0,

donde∥.∥_{C_{k}^{−1}} denota la norma inducida por la matriz C_{k}^{−1} = (D_{k}^{T} D_{k})^{−1} . En el caso que existan constantes α_{1}, α_{2}, tales que

∥x_{k} − x_{m}∥ ≤ α_{1}, k = m, ...l ⇒ ∥D_{l} − D_{m}∥ ≤ α_{2}

g_{k} = ∇f_{k}, f esté acotada inferiormente y ∇f sea uniformente continuo, puede afirmarse que

lim_{k→∞} ∥∇f(x_{k})∥_{C_{k}^{−1}}= 0.

Observar que en el caso particular que la sucesión {x_{k}} converja, puede asegurarse que lo hace a un punto estacionario.

La condición g_{k} = ∇f_{k}, es obviamente muy restrictiva, ya que en la práctica muchas veces se desconoce la expresión analítica del gradiente. Este hecho fue resuelto por Carter [3], quien reemplaza esta condición por una hipótesis que se refiere a acotar uniformemente el error relativo en la aproximación del gradiente, tal como lo establece el siguiente resultado:

Teorema 18 Sea f : ℝ^{n} → ℝ continuamente diferenciable en un abierto convexo D que contiene al punto x_{k}, sean {B_{k}} , {D_{k}^{T} D_{k}} y {(D_{k}^{T} D_{k})^{−1}} uniformemente acotadas, {x_{k}} la sucesión de iterados definida por el algoritmo. Si \frac{∥e_{k}∥}{∥g_{k}∥} ≤ ζ < 1 − η_{2}, entonces

lim inf_{k→∞} ∥g(x_{k})∥_{C_{k}^{−1}} =lim inf_{k→∞} ∥∇f(x_{k})∥_{C_{k}^{−1}} = 0

o

g(x_{k}) = ∇f(x_{k}) = 0

para algún iterado x_{k}.

4.4 Criterios de convergencia

Los criterios adoptados para detener el algoritmo son los sugeridos por Dennis, Gay y Welsh [8]:

1) Criterio de convergencia absoluta:

∥f(x_{k})∥ ≤ ε_{1}.

2) Criterio de gradiente proyectado:

∥∇_{p}f(x_{k})∥ ≤ ε_{2}.

Los siguientes criterios tienen en cuenta la confiablidad del modelo, pues sólo se aplican si se verifica la siguiente relación:

ared ≤ 2 ∗ pred

y ellos son:

3) (X-convergencia)

Reld(x_{k}, x_{k+1}) = \frac{∥x_{k} − x_{k+1}∥_{∞}}{∥x_{k}∥_{∞} + ∥x_{k+1}∥_{∞}} ≤ ε_{X}.

Observar que esto ocurre cuando x_{k} ≈ x^{∗} y éste es un mínimo fuerte, es decir, ∇f(x^{∗}) = 0 y ∇^{2}f(x^{∗}) es definida positiva.

4) (Convergencia relativa)

f(x_{k}) − q_{k}(x_{k+1}) ≤ ε_{r}f(x_{k}).

Sucede cuando f(x_{k}) ≈ f(x^{∗}) y x^{∗} es mínimo fuerte.

5) Si no se verifican ninguno de los criterios arriba mencionados pero,

Reld(x_{k}, x_{k+1}) ≤ ε_{F}

declaramos ‘falsa convergencia’. Esto sucede cuando las iteraciones convergen a un punto que no es crítico.

6) Finalmente cabe señalar que detenemos el algoritmo cuando el radio de la región de confianza es demasiado pequeño y en este caso no se declara convergencia.

4.5 Detalles de implementación

Implementamos este algoritmo en lenguaje FORTRAN con doble precisión, usando el compilador Microsoft FORTRAN Power Station (versión 1.0) en una computadora personal con procesador Pentium II.

El subproblema cuadrático fue resuelto utilizando la subrutina QUA-CON que implementa los resultados teóricos obtenidos por Friedlander y Martinez [12] y nuevas propuestas de minimización dentro de un cara basadas en alternativas diferentes descriptas en[13] [14].

Los algoritmos tradicionales para resolver este tipo de problemas se basan en el métodos de restricciones activas y de gradiente proyectado.

El primero se caracteriza por que en la k-ésima iteración se define el conjunto I de restricciones activas en x_{k} y se busca una dirección de descenso d_{k} en la cara definida por I y se ignoran las demás restricciones. Se determina λ tal que x^{−}= x_{k} + λd_{k} minimiza q(x) en la cara. Si una restricción pasó a ser activa en x^{−}, es decir se llegó al borde de la cara, se actualiza I y se pasa a la iteración siguiente.

Por otro lado, en caso en que d_{k} = 0, se ha encontrado el mínimo en la cara asociada a I , se calculan los multiplicadores de Lagrange correspondientes a las restricciones activas, se verifican ciertas condiciones de optimalidad y entonces, o se detiene el algoritmo o se remueve una de las restricciones de I y se pasa a otra cara.

El método del gradiente proyectado intenta corregir la lentitud con que el método anterior agrega o remueve restricciones, sobre todo en los problemas grandes, y la otra desventaja que es que en cada cara hay que resolver un problema de minimización.

Los métodos más modernos emplean gradiente conjugado para hallar direcciones de descenso dentro de una cara o gradiente proyectado para alterar varias componentes de I antes de encontrar el punto de mínimo en la cara.

El método propuesto en [12] combina las estrategias de restricciones activas con gradiente proyectado y se caracteriza básicamente por lo siguiente:

• Utiliza el gradiente cortado para cambiar de cara.

• Abandona una cara de acuerdo a una relación entre la norma de la componente del gradiente perpendicular a la cara y la componente tangencial a la misma.

• Tiene la posibilidad de utilizar distintos tipos de métodos de minimización de acuerdo a la dimensión de la cara en la cual se trabaja, incluyendo descomposición de Cholesky (recomendada para problemas de dimensiones pequeñas ), gradientes conjugados y gradientes con retardo [13].

Los valores escogidos para los parámetros en los criterios de convergencia son: ε_{r} = 10^{−10}, ε_{F} = 10^{−14}, ε_{x} = 10^{−8}.

Se inicia cada problema con δ_{0} =min_{i} {u_{i} − x_{0i}, x_{0i} − l_{i}}, y el valor mínimo admisible para el radio de la región de confianza se ha elegido igual a 10^{−12}. El valor máximo obviamente depende del conjunto factible. Los valores de α y β utilizados para redefinir el radio de la región de confianza, se han elegido iguales a 2 y 0.5 respectivamente.

Por el tipo de restricciones del problema, se ha escogido la norma l_{∞} para definir el radio de la región de confianza.

Debido a la dimensión de los problemas testeados se ha trabajado con descomposición de Cholesky como método para la resolución del subproblema cuadrático.

Cuando el modelo preferencial es del tipo Gauss Newton, se utiliza en su construcción la matriz Hessiana exacta; cuando se trata de un modelo Quasi Newton, el término de la Hessiana que involucra información de segundo orden se actualiza con una aproximación de tipo BFG, y se utilizó la subrutina CSECUNF del paquete UCMIN. Se han usado los escalamientos sugeridos en [38], pero no se ha notado una variación en la perfomance de los problemas que justifique su uso.

Se hicieron dos versiones de este algoritmo, en la primera de ellas el modelo cuadrático utiliza el gradiente analítico de la función objetivo y se comprobó la eficiencia del mismo resolviendo los problemas propuestos por Moré, Garbow y Hillstrom [35], con las cotas sugeridas por Gay [17]. La razón de esta elección se fundamenta en el deseo de comparar la perfomance de este algoritmo con el propuesto por Lindström y Wedin [29].

Los resultados están contenidos en la tabla que siguiente. En dicha tabla, la primera línea de un problema dado corresponde a los resultados obtenidos tomando como punto inicial aquél propuesto por [35], la segunda corresponde a 10x_{0}, y así sucesivamente.

La columna encabezada por Nuevo indica el número de iteraciones requeridas por el algoritmo aquí descripto para llegar a la solución y la encabezada por L&W, las que fueron necesarias para el algoritmo de Lindström y Wedin [29].

La variable C indica el criterio de convergencia verificado por el problema en consideración, con los números de los criterios arriba mencionados.

4.5.1 Tabla 1.

PROBLEMA  F  ∥Grad.P roy.∥  Iteraciones C
Dimensiones  Nuevo L&W

Linear-full rank (10,5)

2.624999999999999 {8.845400755139059E- \atop 016} 1 2 2

Linear-full rank (50,5)

22.625000000000010 {4.884731606239931E- \atop 015} 1 2 2

Linear rank 1 (10,5)

1.07142857142857 {5.253957563165094E- \atop 012} 1 1 2

Linear rank 1 (50,5)

6.064356435643565 {1.684928366288680E- \atop 010} 1 1 2

Linear-rank1,0’s (10,5)

1.823529411764706 {1.341942022882196E- \atop 011} 1 3 2

Linear-rank1,0’s (50,5)

97.520000000000680 5252.800000000019000 2 3

4

Rosenbrock (2,2)

{1.250000000000000E- \atop 001}

0.000000000000E+000 9 9

2

1.d1* x0  {1.250000000000000E- \atop 001} 0.000000000000E+000 12  6  2
1.d2* x0  {1.250000000000000E- \atop 001} 0.000000000000E+000 14  64  2

Helycal Vally (3,3)

{4.952118423816692E- \atop 001} {4.707941763538051E- \atop 007} 11 9 3

1.d1* x0  {4.952118423816692E- \atop 001} {4.707941763538051E- \atop 007} 11  11  3
1.d2*x0  {4.952118423816692E- \atop 001} {4.707941763538051E- \atop 007} 11  10  3

Powell Singular (4,4)

{9.390981620569023E- \atop 005} {7.464962023325440E- \atop 007} 18 8 3

1.d1*x0  {9.390981579966693E- \atop 005} {6.041596474724106E- \atop 007} 22  10  3
1.d2* x0  {9.390981594654762E- \atop 005} {6.591819268326061E- \atop 007} 21  12  3

Freudenstein, Roth (2,2)

24.492601000000000 {8.881784197001252E- \atop 016} 3 4 2

1.d1* x0  24.492601000000000 {8.881784197001252E- \atop 016} 8  8  2
1.d2* x0  24.492601000000000 {8.881784197001252E- \atop 016} 9  8  2
Bard (15,3)  {4.791142360576711E- \atop 003} {9.271987015073297E- \atop 010} 5  5  2
1.d1* x0  {4.791142382080187E- \atop 003} {5.147528986063535E- \atop 007} 9  11  3
1.d2* x0  {4.791142360576659E- \atop 003} {6.170008435797764E- \atop 010} 21  12  2

Kowalick Osborne (11,4)

{1.539153937497921E- \atop 004} {9.76514742369141E- \atop 007} 8 5 3

1.d1* x_{0}  {1.539153762361323E- \atop 004} {6.504988587035316E- \atop 007} 11  8  3
1.d2* x_{0}  {1.539153617333943E- \atop 004} {6.025898432654351E- \atop 009} 31  46  3
Meyer (16,3)  no converge  6
1.d1*x_{0}  63.725885605523770 {3.487120442641597E- \atop 007} 74  782  3
Watson (31,6)  {1.144803219125418E- \atop 003} {1.411872071196223E- \atop 007} 6  7  3
1.d1* x_{0}  {1.144803219127584E- \atop 003} {2.744291338568853E- \atop 007} 11  13  4
1.d2* x_{0}  {1.144803219137375E- \atop 003} {5.137874876166839E- \atop 007} 14  17  3
Watson (31,9)  {1.869696212899957E- \atop 002} {5.458265413318252E- \atop 010} 7  12  2
x(i)=1  {1.869696212899972E- \atop 002} {8.422089974079289E- \atop 008} 11  16  3
Watson (31,12)  {3.578573893095571E- \atop 002} {7.979597457005185E- \atop 007} 23  17  2
1.d1* x_{0}  {3.578573839133194E- \atop 002} {4.800309845796224E- \atop 007} 21  31  3
Box 3D (10,3)  {5.717409013214852E- \atop 005} {4.887055067965334E- \atop 009} 6  5  2

Jenrich, Sampson (10,2)

62.191556458188930 {1.174030576683860E- \atop 005} 11 10 4

Brown, Dennis (20,4)

44430.239883750250000 {1.303245380888788E- \atop 011} 11 11 2

1.d1* x_{0}  44430.239883750250000 {3.409095394208285E- \atop 007} 13  12  3
1.d2* x0  44430.239883750250000 {1.178458166897514E- \atop 011} 14  9  2

Brown alm. lin. (10,10)

{3.046589290326493E- \atop 003} {5.645928009777020E- \atop 007} 7 10 3

1.d1* x_{0}  {3.046589291786746E- \atop 003} {7.407235204921822E- \atop 007} 14  7  3
1.d2* x_{0}  {3.046589288737870E- \atop 003} {2.623760975161611E- \atop 007} 10  31  3

Brown alm. lin. (30,30)

{1.441793436461608E- \atop 002} {5.470749858922064E- \atop 007} 16 6 3

Brown alm. lin. (40,40)

{2.031892506793320E- \atop 002} {8.767545683932677E- \atop 007} 26 13 3

Osborne 1 (33,5) {2.735141671686873E- \atop 005} {7.630200683605875E- \atop 009} 6  8  2
Osborne 2 (66,11) {2.551152610242830E- \atop 001} {5.397308904986227E- \atop 007} 23  24  3

La segunda versión del algoritmo se implementó para el caso en que se desconoce la expresión analítica de la función objetivo, y por lo tanto se utiliza gradiente aproximado en la definición del modelo cuadrático. Se aplicó en un problema de identificación de parámetros en la solución de un sistema de ecuaciones algebraico-diferenciales, que se plantea al modelar procesos bioquímicos.

En el ejemplo considerado el modelo a ajustar es la solución de un sistema de ecuaciones diferenciales ordinarias

\{ \begin{matrix}\frac{dy}{dt} = f(t, y_{1}, y_{2}, ..., y_{p,}x_{1}, ..., x_{n}) \\ y(t_{0}, x_{1}, x_{2}, ..., x_{n}) = y_{0} \end{matrix}

donde y : ℝ → ℝ^{p} es una función no lineal en los parámetros x_{1}, x_{2}, ..., x_{n}, algunos de los cuales debido a la naturaleza del problema, satisfacen restricciones de canalización:

l_{i} ≤ x_{i} ≤ u_{i}.

Asumiendo que se han realizado observaciones ŷ(t_{i}, x), esto es, se conocen valores

ŷ(t_{i}, x) = ( \begin{matrix}ŷ_{1}(t_{i}, x_{1}, ..., x_{n}) \\ ŷ_{2}(t_{i}, x_{1}, ..., x_{n}) \\ ........................ \\ ŷ_{p}(t_{i}, x_{1}, ..., x_{n}) \end{matrix} )

en cada instante t_{i}, i = 1, ..., m, que dependen de los parámetros x_{i}, i = 1, ..., n; si y_{∗}(t_{i}, x_{1}, ..., x_{n}) denota la solución del sistema de ecuaciones calculada en t_{i} (i = 1, ..., m) el objetivo es determinar los valores de los parámetros (x_{1}, ...x_{n}) que resuelven el problema

min_{x} f(x) = ∑_{i=1}^{m}(ŷ(t_{i}, x) − y_{∗}(t_{i}, x))^{2}.

El proceso bioquímico aquí considerado describe el tratamiento de fermentación de células en un sustrato (levadura, penicilina, etc.) y en este caso el sistema de ecuaciones diferenciales ordinarias está formado por las ecuaciones que describen:

• la cinética microbiana:

\frac{dX}{dt} = mX − kdX

• el consumo de sustrato:

\frac{dS}{dt} = \frac{mX}{Y xs} + ms.X + \frac{qp.X}{Y ps}

• y la formación de producto:

\frac{dP}{dt} = qp.X − kdp.X

donde las incógnitas y parámetros representan: X= concentración de células en la fermentación. (g/L), m= velocidad específica de crecimento (1/h), kd= tasa específica de pérdida de viabilidad celular (1/h), S= concentración de sustrato de fermentación (g/L), Yxs= cociente de rendimiento de biomasa con respecto al sustrato (g(X)/g(S)), ms= cociente de mantenimiento celular con respecto al sustrato (g(X)/g(S).h), qp= velocidad específica de formación del producto (g(P)/g(X).h), Yps= cociente de rendimiento de producto con respecto al sustrato (g(P)/g(S)), kdp= constante de degradación del producto (1/h), P= concentración de producto en la fermentación (g/L).

Este problema fue proporcionado por la Dra. Longina Castellanos Noda [4], con datos aportados a partir de experiencias en la Universidad Autónoma de Méjico. Hasta la actualidad había sido resuelto con la rutina LBFGSB, que implementa un método de memoria limitada para problemas de gran escala, con los criterios de terminación propuestos en [36].

En la tabla siguiente se puede comprobar la eficiencia del algoritmo que proponemos frente al utilizado hasta ahora, comparando el número de iteraciones necesarias para lograr convergencia con ambos métodos.

En la primera columna se indican las dimensiones del problema, en la segunda el número de iteraciones requeridas por nuestro algoritmo y en la tercera las necesarias por el algoritmo de memoria limitada. La última columna indica el criterio de convergencia satisfecho.

4.5.2 Tabla 2.

Problema.  Iteraciones  C
Dimensiones Alg. nuevo Alg. LBFGSB.
1. (49,2)  7  24  4
2. (42,2)  7  21  4
3. (43,2)  4  15  4
4. (25,3)  15  92  4
5. (25,3)  17  87  1
6. (25,3)  160  81  4
7. (25,4)  21  46  4
8. (21,4)  85  42  4
9. (25,3)  15  35  4
10. (26,3)  9  310  4
11. (25,3)  3  83  4
12. (25,3)  6  74  4

4.6 Conclusiones

Se ha obtenido un nuevo algoritmo que utiliza técnicas desarrolladas por reconocidos matemáticos que permite resolver problemas de cuadrados mínimos no lineales con restricciones en las variables. Ha sido probado con problemas clásicos de la especialidad y en otros proporcionados por quienes trabajan en esos temas específicos. La versatilidad del mismo permite su aplicación en muchos problemas reales en los cuales se desconoce la expresión analítica de la función a minimizar.

Actualmente se está trabajando en su aplicación a un problema de estimación de la velocidad del agua en un medio poroso (acuíferos), que tienen la característica de ser constante por zonas y que lleva a plantear un problema de cuadrados mínimos con residuos grandes. La valores de la función objetivo se obtienen a partir de una ecuación diferencial de tipo parabólica, la que se resolvió con el método de diferencias finitas y el mismo se empleó para aproximar el gradiente.

Nos proponemos generalizarlo para el caso en que el problema además de requerir cotas en las varibles, presente restricciones de igualdad. En este caso se aplicará la teoría para el caso general de minimización con restricciones de igualdad desarrollada por Gomez, Maciel y Martinez [23].

Bibliografía

[1] A.A. Anda and H. Park. Self scaling fast rotations for stiff least squares a problems. Linear Algebra Appl., 234 (1996),pp.137-162.

[2] Å. Björck. Numerical Methods for Least Squares Problems. SIAM. 1996.

[3] R.G. Carter. On the Global Convergence of Trust Region Algorithms Using Inexact Gradient Information. SIAM J. Numer. Anal., Vol. 28, (1991) pp. 251-265.

[4] L. Castellanos Noda. Instituto de Cibernética, Física y Matemática de Cuba. Comunicación privada, 1998.

[5] A.K. Cline. The transformation of a quadratic programming problem into solvable form. Tech. Rep. ICASE 75-14, NASA, Langley Research Center, Hampton,VA., 1975.

[6] M.G. Cox. The least square solution of overdeterminated linear equations having band or augmented band structure. IMA J. Nume. Anal., 1 (1981) pp. 3-22.

[7] J.E. Dennis, D.M. Gay and R.E.Welsch. An adaptive nonlinear leastsquares algorithm. ACM Trans. Math Software, 7(1981), pp. 348-368.

[8] J.E. Dennis, D.M. Gay and R.E.Welsch. Algorithm 573 NL2SOL:An adaptive nonlinear least-squares algorithm. ACM Trans. Math Software, 7(1981), pp. 369-383.

[9] J.E. Dennis and R. Schnabel. Numerical Methods for Unconstrained optimization and Nonlinear Equations. Prentice Hall,Inc., 1983.

[10] N.R. Draper and H. Smith. Applied Regression Analysis, 2d.ed. John Wiley, New York, 1981.

[11] R. Fletcher. A modified Marquardt subroutine for non-linear least squares. Harwell Report AERE-R.6799. (1972).

[12] A. Friedlander and J.M. Martínez. On the maximization of a concave quadratic function with box constraints. SIAM Journal on Optimization, (1994) 4(1):177-92.

[13] A. Friedlander, J.M. Martínez, B. Molina and M. Raydan. Gradient method with retards and generalizations. SIAM J. Num. Anal., 36 (1999), pp. 275-289.

[14] A. Friedlander, J.M. Martínez and M. Raydan. A new method for largescale box constrained convex quadratic minimization problems. Optimization Methods and Software, (1995), Vol 5, pp. 57-74.

[15] W. Gander. Least squares with a quadratic constraint. Num. Math., 36 (1981), pp. 291-307.

[16] D. Gay. Computing optimal locally constrained steps. SIAM J. Sci. Statist. Comput., 2, (1981), pp.186-197.

[17] D. Gay. A trust region approach to linearly constrained optimization. Lecture Notes in Mathematics 1066.Springer Verlag, Berlin, Heidelberg, New York, Tokio, 1984.

[18] P.E. Gill, S.J. Hammarling, W. Murray, M.A. Saunders, and M.H. Wright. User ’s guide for LSSOL (versión 1.0): A Fortran Package for constrained linear least squares and convex quadratic programming. Report SOl, Department of operations Research, Stanford university, CA, 1986.

[19] P.E. Gill, W. Murray, and M.H. Wright, Practical Optimization, Academic Press, London and New York, 1981.

[20] G.H. Golub. Numerical Methods for solving least squares problems. Numer,Math.,7(1965), pp. 206-216.

[21] G.H. Golub, F.T. Luk and M. Pagano. A sparse least squares problem in photogrammetry. Proceedings of the Computer Science and Statistics 12th Annual Simposium on the Interface, J.F.Gentleman, ed., University of Waterloo, Canada, 1979, pp.26-30.

[22] G.H. Golub and R.J. Plemmons. Large-scale geodetics least-squares adjustment by dissection and ortogonal descomposition. Lin. Alg. its Appl., 34 (1980) pp. 3-28.

[23] F. Gomez, M.C. Maciel and J.M. Martínez. Non linear programming algorithms using trust regions and augmented Lagrangians with nonmonotone penalty parameters. To appear in Matehmatical Programming, 1998.

[24] K.H. Haskell and R.J. Hanson. An algoritm for linear least squares problems with equality and nonnegativity constraints. Math. Programming, 21 (1981), pp. 98-118.

[25] J.N. Holt and R. Fletcher. An algorithm for constrained non-linear least squares. J.Inst.Maths. Applic., 23 (1979), pp.449-463.

[26] C. Lawson and R.Hanson. Solving Least Squares Problems. Prentice Hall, Englewood Cliffs, NJ,1974.

[27] O. Leringe and P. Wedin. A comparison between different methods to compute a vector x wich minimize ∥Ax − b∥_{2} when Gx = h. Tech. Report, Department of Computer Science, Lund University, Lund, Sweden, 1970.

[28] P. Lindström and P. Wedin. A new linesearch algorithm for unconstrained nonlinear least squares problems. Math. Programming, 29 (1984), pp 268-296.

[29] P. Lindström and P. Wedin. Gauss-Newton based algorithms for constrained nonlinear least squares problems. Tech. Report, Department of Computer Science, University of Umeå, Sweden, 1987.

[30] P. Manneback, C. Murigande and P.L.Toint. A modification of an algorithm by Golub and Plemmons for large linear least squares in the context of Doppler positioning. IMA J.Numer. Anal., 5 (1985), pp. 221-234.

[31] J.M. Martinez. An Algorithm for Solving sparse Nonlinear Least Squares Problems. Computing 39, (1987 ) pp. 307-325.

[32] D.W. Marquardt. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Indust. Appl. Math., 11 (1963), pp. 431-441.

[33] J.J. Moré. The Levenberg-Marquandt algorithm: Implementation and theory. Numerical Analysis, G.A. Watson, ed., Lecture Notes in Math. 630, Springer Verlag, Berlín, (1978), pp.105-116.

[34] J.J. Moré. Recent Developments in Algorithms and Software for Trust Region Methods. Mathematical Programming. The State of the Art, Proceedings Bonn 1982, A. Bachem, M Grötchel, and B. Korte, eds., Springr-Verlag, Berlin, 1983, pp. 258-287.

[35] J.J.Moré, B. Garbow and K. Hillstrom. Testing unsconstrained optimization software. Technical Memorandum. Aragonne National Laboratorie. 1978.

[36] D.C. Liu and J. Nocedal. On the limited memory BFGS method for large scale optimization methods. Mathematical Programming 45 (1989), pp 503-528.

[37] M.J.D. Powell. Convergence properties of a class of minimization algorithms. Nonlinear Programming 2, O.L. Mangasarian, R.R. Meyer, and S.M. Robinson,eds., Ac. Press, 1975.

[38] S. Oren and D. Luenberger. Self-scaling variable metric (SSVM) algorithms. Management Science. Vol 20. N^{0}5, (1974), pp. 845-862.

[39] J.K. Reid A note on the least squares solutions of a band system of linear equations by Householder reductions. Comp.J, 10 (1967),pp. 188-189.

[40] K. Schittkowski. The numerical solution of constrained linear least squares problems. IMA Journal of Numerical Analysis,3 (1983), pp. 11-36

[41] K. Schittkowski and J. Stoer. A factorization method for constrained linear least squares problems allowing subsequent data changes. Numer. Math., 31 (1979), pp. 431-463.

[42] J. Stoer. On the numerical solution of constrained least squares problems.SIAM J. Num. Anal., 8 (1971), pp. 382-411.

[43] C.F. Van Loan. On the method of weighting for equality constrained least squares. SIAM J. Numer. Anal., 22 (1985), pp. 851-864.

[44] P.A. Wedin. On the Gauss Newton method for the nonlinear least squares problems, Working Paper 24, Intitute for Applied Mathematics, Stockholm, Sweden, 1974.

[45] S.J. Wright and J.N. Holt. Algorithms for nonlinear least squares with linear inequality constraints. SIAM J. Sci. Stat. Comput. Vol 6 N^{0}4, (1985), pp. 1033-1048.

[46] P. Zimmermann. Ein algorithmus zur Lösung lenearer Squares Probleme mit unteren und oberen Scharken als Nebebbedingungen. Diplomarbeit, Universität Würzburg, Germany, 1977. Introduction