Un algoritmo cuasi Newton global para problemas

de complementariedad no lineal

Carlos Andrés Arias Torres

Universidad del Cauca

Facultad de Ciencias Naturales, Exactas y de la Educación
Departamento de Matemáticas
Maestrı́a en Ciencias Matemáticas
Popayán
2014
 Un algoritmo cuasi Newton global para problemas
de complementariedad no lineal

Carlos Andrés Arias Torres

Trabajo de investigación presentado como requisito parcial para optar al

tı́tulo de Magı́ster en Ciencias Matemáticas

Directora
Dra. Rosana Pérez Mera
Profesora de la Universidad del Cauca

Codirector
Dr. Héctor Jairo Martı́nez Romero
Profesor de la Universidad del Valle

Universidad del Cauca

Facultad de Ciencias Naturales, Exactas y de la Educación
Departamento de Matemáticas
Maestrı́a en Ciencias Matemáticas
Popayán
2014
 Nota de aceptación

Director: Dra. Rosana Pérez Mera

Jurado: Dr. Jhon Jairo Duque Robles

Jurado: Mg. Favián Enrique Arenas Aparicio

Popayán, 24 de febrero de 2014

 AGRADECIMIENTOS

A mi familia por motivar cada uno de mis esfuerzos. A mi padre y a mi madre porque
han hecho posible cada uno de mis sueños. A Isabel por brindarme su amor y su apoyo
incondicional.

A mi directora de tesis, Dra. Rosana Pérez Mera, por el conocimiento que me compartió y
la motivación que me brindó en este arduo camino. Porque con su visión y su experiencia
me guió hacia este logro.

A mi codirector Dr. Héctor Jairo Martinez, por sus aportes e invaluables consejos que
consolidaron esta investigación.

Al Dr. Jhon Jairo Duque Robles y al Mg. Favián Enrique Arenas Aparicio por sus
acertadas sugerencias, fruto de sus conocimientos y su experiencia investigativa.

A mis profesores por enseñarme a ser un ser humano ı́ntegro.

Carlos Andrés Arias Torres

I
 TABLA DE CONTENIDOS

Índice de tablas IV

1. Introducción 1

2. Preliminares 6

3. Algoritmo y Teorı́a de Convergencia 12

3.1. Hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.2. Teorı́a de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4. Pruebas Numéricas 23

4.1. Pruebas con los Algoritmos 3.1 y 4.1 . . . . . . . . . . . . . . . . . . . . 27

4.2. Búsqueda lineal no monótona . . . . . . . . . . . . . . . . . . . . . . . . 30

4.3. Variación del parámetro λ . . . . . . . . . . . . . . . . . . . . . . . . . . 36

II
TABLA DE CONTENIDOS III

5. Comentarios finales 39

Bibliografı́a 41
 ÍNDICE DE TABLAS

4.1. Resultados numéricos para los Algoritmos 3.1 y 4.1. . . . . . . . . . . 29

4.2. Resultados del Algoritmo 4.2 con búsqueda no monótona. . . . . . . . 33

4.3. Resultados numéricos al variar el parámetro M . . . . . . . . . . . . . . . 34

4.4. Desempeño del Algoritmo 4.2 con M = 10. . . . . . . . . . . . . . . . 35

4.5. Desempeño del Algoritmo 3.1 al variar λ. . . . . . . . . . . . . . . . . 38

IV
 CAPÍTULO 1

INTRODUCCIÓN

Sea F : Rn −→ Rn , F (x) = (F1 (x), F2 (x), ..., Fn (x)) , una función no lineal y continua-
mente diferenciable. El Problema de Complementariedad No Lineal (PCNL), consiste en
hallar un vector x ∈ Rn tal que

x ≥ 0 ; F (x) ≥ 0 ; xT F (x) = 0, (1.1)

donde la expresión y ≥ 0 para un vector en Rn significa que todas sus componentes

son no negativas.

La tercera condición de (1.1) implica que xi Fi (x) = 0 para i = 1, 2, ..., n y por lo tanto,
xi = 0 o Fi (x) = 0 , siendo esta la razón por la cual el problema recibe el calificativo de
complementariedad. Esta condición trae consigo implı́citamente la búsqueda de un equi-
librio entre la variable del problema y el valor que toma la función que define el problema
en dicha variable. Ası́, el concepto de complementariedad es sinónimo del concepto de
sistema en equilibrio [17].

Por lo mencionado anteriormente, problemas de complementariedad no lineal surgen de

forma natural en campos como la Ingenierı́a y la Fı́sica, por ejemplo, cuando se presentan
situaciones de contacto mecánico [17], equilibrio del tráfico [35] y problemas de lubricación
elasto-hidrodinámicos [23], ası́ como en el campo de la Economı́a, particularmente cuando

1
1. Introducción 2

se busca un equilibrio en los sistemas económicos [17].

Es por ello que ha sido de gran interés para la comunidad matemática el estudio de
diversos métodos que permitan solucionar el PCNL, entre ellos, métodos de homotopı́a,
que son derivados de métodos de punto fijo [37] [14]; métodos basados en redes neuro-
nales [39][24][25] y métodos de ecuaciones no diferenciables [12][11]. En particular, estos
últimos, consisten en reformular el PCNL como un sistema de ecuaciones no lineales, no
diferenciable, usando para ello una función ϕ : R2 −→ R tal que

ϕ(a, b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0 , (1.2)

conocida como función de complementariedad [21][38].

Geométricamente, a partir de (1.2), se infiere que la traza de la función ϕ obtenida por

la intersección con el plano xy, es la curva formada por los semiejes positivos x y y ,
la cual no es diferenciable en (0, 0) , razón por la que ϕ no es diferenciable [1].

Dos de las funciones de complementariedad más utilizadas en la solución de PCNL son

la función mı́nimo[34][29] y la función de Fischer [19][18], definidas respectivamente por:
√
ϕ(a, b) = mı́n{a, b} ϕ(a, b) = a2 + b2 − a − b .

En 1998, Kanzow y Kleinmichel [21] presentaron una familia de funciones de comple-

mentariedad ϕλ dependientes de un parámetro λ ∈ (0, 4), definida por
p
ϕλ (a, b) = (a − b)2 + λab − a − b . (1.3)

Observemos que esta familia de funciones incluye la función de Fischer como un caso
particular, cuando λ = 2, y converge a un múltiplo de la función mı́nimo cuando λ
tiende a cero.

Ası́, si definimos la función Φλ : Rn −→ Rn por

 
ϕλ (x1 , F1 (x))
Φλ (x ) = 
 .. 
(1.4)
. 
ϕλ (xn , Fn (x))

entonces un vector x∗ es solución del PCNL, si y solo si, x∗ es solución del sistema de
ecuaciones
Φλ (x) = 0. (1.5)
Como consecuencia de la no diferenciabilidad de ϕλ , el sistema de ecuaciones no lineales
1. Introducción 3

(1.5) no es diferenciable. De esta forma, el PCNL es reformulado como un sistema de

ecuaciones no lineales, no diferenciable.

Entre los métodos más populares para resolver sistemas de ecuaciones no lineales dife-
renciables de la forma G(x) = 0, se encuentran el método de Newton [27] y métodos
Cuasi Newton [13]. La principal diferencia entre ellos radica en que en el primero es
necesario calcular, en cada iteración, la matriz jacobiana de G, lo que computacional-
mente resulta costoso, mientras que en el segundo tipo de métodos basta con tomar una
aproximación a dicha matriz, pagando a cambio un cierto precio que se ve reflejado en
una disminución de la velocidad de convergencia del algoritmo respectivo [13].

Sin embargo, cuando el sistema de ecuaciones G(x) = 0 no es diferenciable, como es

el caso del sistema de ecuaciones (1.5), no tiene sentido hablar de la “matriz jacobiana
de G”. En 1973, Frank H. Clarke [7][8][9] presentó el concepto de jacobiano generalizado,
definido para G : Rn −→ Rn Lipschitz continua 1 , como el conjunto de matrices
n o
0 n×n
∂G(x) = conv lı́m G (xk ) ∈ R : lı́m xk = x, xk ∈ DG , (1.6)
k→∞ k→∞

donde DG es el conjunto de todos los vectores de Rn en los que G es diferenciable y

conv {A} representa la envolvente convexa del conjunto A . En particular, el conjunto
n o
∂B G(x) = lı́m G0 (xk ) ∈ Rn×n : lı́m xk = x, xk ∈ DG (1.7)
k→∞ k→∞

se conoce como B-jacobiano Generalizado de G en x [34]. Es decir, de (1.6) y (1.7), el

Jacobiano generalizado de G en x puede verse como el conjunto

∂G(x) = conv {∂B G(x)}. (1.8)

Dado que cada función de complementariedad ϕλ es Lipschitz continua, entonces la

función Φλ también lo es y por tanto, su jacobiano generalizado existe.

Usando matrices del conjunto (1.8), Qi [33] propone el llamado método de Newton ge-
neralizado para resolver sistemas de ecuaciones no lineales no diferenciables, el cual, en
cada iteración, resuelve un sistema lineal cuya matriz de coeficientes es una de las matri-
ces del conjunto ∂B G(x). Por otra parte, para problemas de complementariedad no lineal
también se han propuesto métodos cuasi Newton [31][1][6], los cuales han mostrado un
1
Una función G : Rn −→ Rn es Lipschitz continua, si existe una constante positiva K tal que

kG(x) − G(y)k2 ≤ K kx − yk2 ,

para todo x y y en Rn
1. Introducción 4

buen desempeño numérico.

En general, el método de Newton para resolver sistemas de ecuaciones no lineales, bajo

ciertas hipótesis, presenta una tasa de convergencia q-cuadrática [13] siempre y cuando
el punto inicial del algoritmo esté suficientemente cerca de la solución del sistema (es
decir, es un método local ). Preguntas que surgen al implementar algoritmos basados en
el método de Newton y que pueden llevar a problemas en su convergencia son ¿qué tan
cerca debe estar el punto inicial de la solución del sistema? ¿cómo hacer una buena
elección del punto inicial para implementar el algoritmo?

Con el fin de evitar inconvenientes en la elección del punto inicial, los algoritmos locales
son globalizados, es decir, son dotados de un criterio que permite ejecutarlos cuando su
punto inicial está lejos de la solución del problema.

Algunos de esos criterios se basan en ciertas funciones conocidas como funciones de méri-
to [27]. La función de mérito más utilizada en el caso de sistemas de ecuaciones no lineales
de la forma G(x) = 0 , es la función f : Rn −→ R definida por f (x) = kG(x)k22 . Te-
niendo como base esta función, una nueva iteración en un algoritmo global será aceptada
siempre y cuando el valor de la función f decrezca con respecto a la iteración actual del
mismo; este requerimiento conduce a resolver el problema de minimización
1
M inimizar kG(x)k22 (1.9)
x∈R n 2

1
donde el factor 2
es adicionado por conveniencia algebraica.

Cabe destacar que toda solución del sistema de ecuaciones G(x) = 0 es un minimizador
de (1.9), sin embargo, pueden existir soluciones locales de (1.9) que no son soluciones del
sistema de ecuaciones no lineales asociado.

Teniendo en cuenta lo anterior y con el fin de globalizar algún método tipo Newton que
busque hallar una solución del sistema de ecuaciones no lineales (1.5), recurriendo al
jacobiano generalizado de Φλ , se puede hacer uso de una función de mérito asociada a
este sistema. En efecto, para tal fin basta considerar Ψλ : Rn −→ R , definida por
1 1
Ψλ (x) = Φλ (x)T Φλ (x) = kΦλ (x)k22 . (1.10)
2 2
Observemos que Ψλ es continuamente diferenciable (Teorema 3.1 en [21] y Proposición
3.4 en [15]). Con esta función de mérito, podemos reformular nuevamente el sistema de
ecuaciones no lineales, no diferenciable (1.5) como el siguiente problema de minimización
1. Introducción 5

diferenciable

M inimizar Ψλ (x) . (1.11)

x ∈ Rn

En [21], los autores proponen un algoritmo global tipo Newton generalizado para resolver
el problema de minimización (1.11). Recientemente, en [1], el autor propone un algoritmo
cuasi Newton local que permite resolver el PCNL recurriendo al sistema de ecuaciones
(1.5) y deja abierta la posibilidad de introducir estrategias de globalización para mejorar
el desempeño de su algoritmo.

Teniendo en cuenta lo anterior y pensando en las limitaciones que puede tener un al-
goritmo local, en este trabajo de investigación, proponemos un algoritmo cuasi Newton
global para resolver el PCNL a través de la minimización de la función de mérito (1.10),
es decir, proponemos una globalización del algoritmo presentado en [1], usando ideas
del algoritmo expuesto en [21]. Desarrollamos su teorı́a de convergencia global y rea-
lizamos pruebas numéricas preliminares que muestran un buen desempeño del mismo.
Adicionalmente, analizamos numéricamente el comportamiento del algoritmo propuesto
cuando cambiamos la búsqueda lineal estándar, por una búsqueda lineal no monótona
[20]. También incluimos una exploración preliminar de la variación del parámetro λ y
su impacto en la convergencia global de nuestro algoritmo

El documento está organizado de la siguiente forma: En el Capı́tulo 2, presentamos

algunos resultados teóricos que serán útiles en el desarrollo de la teorı́a de convergencia
del algoritmo que propondremos. En el Capı́tulo 3, proponemos un Algoritmo cuasi
Newton global para resolver el PCNL mediante la minimización de la función de mérito
definida en (1.10) y bajo ciertas hipótesis, desarrollamos su teorı́a de convergencia global.
En la primera parte del Capı́tulo 4, presentamos los resultados numéricos obtenidos
al resolver con el algoritmo cuasi Newton global propuesto, algunos problemas de com-
plementariedad no lineal y comparamos estos resultados con los obtenidos al resolver
los mismos problemas con el algoritmo tipo Newton generalizado presentado en [21]. En
la segunda parte, analizamos numéricamente el desempeño de nuestro algoritmo cuan-
do cambiamos la búsqueda lineal monótona por una búsqueda lineal no monótona [20].
En la última parte de este capı́tulo, presentamos los resultados obtenidos al explorar
el comportamiento del algoritmo cuando variamos el parámetro λ . Finalmente, en el
Capı́tulo 5, hacemos algunos comentarios finales y una propuesta para futuros trabajos
en el tema.
 CAPÍTULO 2

PRELIMINARES

En este capı́tulo, presentamos algunas definiciones y resultados teóricos que han sido
demostrados previamente por algunos autores y que serán útiles en el desarrollo de
nuestra investigación.

De (1.1), una solución x∗ del PCNL se caracteriza porque

x∗ ≥ 0 F (x∗ ) ≥ 0 F (x∗ )T x∗ = 0.

Asociados a esta solución, tenemos los siguientes conjuntos de ı́ndices:

α ={i ∈ I : x∗i > 0, Fi (x∗ ) = 0}

β ={i ∈ I : x∗i = 0 = Fi (x∗ )}
γ ={i ∈ I : x∗i = 0, Fi (x∗ ) > 0}.

Cabe mencionar que si β 6= ∅, la solución x∗ se denomina solución degenerada.

Dado que la función Φλ dada en (1.4) es no diferenciable y Lipschitz continua, su jaco-

6
2. Preliminares 7

biano generalizado existe y está definido por

n o
∂Φλ (x) = conv lı́m Φ0λ (xk ) ∈ Rn×n : lı́m xk = x, xk ∈ DΦλ = conv {∂B Φλ (x)}
k→∞ k→∞

donde DΦλ es el conjunto de todos los vectores de Rn en los que Φλ es diferenciable.

En [1], el autor construye matrices H ∈ ∂Φλ (x) de la forma

 
[H]1
H =  ...  (2.1)
 
[H]n
con

T T
(χ(xi , Fi (x)) − 1)ei + (ψ(xi , Fi (x)) − 1)∇Fi (x) ,
 i 6∈ β
[H]i =

(χ(zi , ∇Fi (x)T z) − 1)eTi + (ψ(zi , ∇Fi (x)T z) − 1)∇Fi (x)T , i∈β


donde z es un vector tal que zi 6= 0 cuando i ∈ β, {e1 , . . . , en } es la base canónica de

Rn y las funciones χ y ψ son las derivadas parciales de la función

p
Gλ (a, b) = (a − b)2 + λab , (2.2)

es decir,

2(a − b) + λb −2(a − b) + λa
χ(a, b) = y ψ(a, b) = · (2.3)
2Gλ (a, b) 2Gλ (a, b)

Además, teniendo en cuenta la forma de las matrices H ∈ ∂Φλ (x) dada en (2.1), el
autor de [1] propone aproximaciones cuasi Newton de dichas matrices de la forma
 
[B]1
B =  ...  (2.4)
 
[B]n
2. Preliminares 8

con 
T
(χ(xi , Fi (x)) − 1)ei + (ψ(xi , Fi (x)) − 1)[A]i ,
 i 6∈ β
[B]i =

(χ(zi , [A]i z) − 1)eTi + (ψ(zi , [A]i z) − 1)[A]i , i∈β


donde la matriz  
[A]1
A =  ...  (2.5)
 
[A]n

es una aproximación a F 0 (x), la matriz jacobiana de F en x.

Con respecto a las derivadas parciales dadas en (2.3), cabe destacar que, para cualquier
vector no nulo (a, b ), estas derivadas están acotadas de la siguiente forma [1][21]:

√
|χ(a, b)| ≤ 1 y |ψ(a, b)| ≤ 2. (2.6)

Con base en estas consideraciones, en [1], el autor propone un algoritmo cuasi Newton
local para resolver el PCNL a través su reformulación como el sistema de ecuaciones
no lineales, no diferenciable (1.5), cuya forma genérica presentamos a continuación ya
que nuestro interés es introducir en dicho algoritmo una estrategia de globalización que
permita obtener un nuevo algoritmo cuasi Newton global para resolver el PCNL.

Algoritmo 2.1. Tipo cuasi Newton local para Φλ (x) = 0.

Dados x0 y λ ∈ (0, 4), para k = 1, 2, . . . , se genera la iteración siguiente por

xk+1 = xk − Bk−1 Φλ (xk );

donde Bk es de la forma (2.4).

Por otra parte, en [1], se demuestra que la distancia entre una matriz H y su respectiva
aproximación B está acotada por un valor constante y, que bajo ciertas condiciones, la
matriz B de la forma (2.4) es no singular. Formalmente, los resultados son los siguientes

Teorema 2.1. [1] Sean x∗ una solución del PCNL; F : Rn −→ Rn , una función de
clase C 1 cuya matriz jacobiana es Lipschitz continua,  y δ constantes positivas dadas;
2. Preliminares 9

H y B definidas por (2.1) y (2.4), respectivamente. Entonces, para cada x ∈ B(x∗ , )

y para cada A ∈ B(F 0 (x∗ ), δ), existe una constante θ > 0 tal que

kH − Bk∞ ≤ θ.

Teorema 2.2. [1] Sean x∗ una solución del PCNL y B la matriz definida por (2.4).
Existen constantes positivas 0 y δ0 tales que si

kx − x∗ k∞ ≤ 0 y kA − F 0 (x∗ )k∞ ≤ δ0 ,

la función Q definida por

Q(x, A) = x − B −1 Φλ (x)
está bien definida.

En particular, este teorema permite garantizar que, si una sucesión {x k } converge a x ∗

entonces existe un entero positivo k tal que, para todo k > k, la matriz Bk−1 existe.

En el desarrollo de la teorı́a de convergencia para algoritmos que resuelvan el PCNL

es importante establecer condiciones suficientes para la no singularidad de las matrices
en el jacobiano generalizado en una solución del problema. Por ello, introducimos a
continuación el concepto de regularidad que nos ayuda para tal fin.

Definición 2.1. Sea x∗ una solución del PCNL.

1. Si todas las matrices H ∈ ∂B Φλ (x∗ ) son no singulares, x∗ es llamada una solución

BD-regular, .

2. Si la submatriz F 0 (x∗ )αα 1 es no singular y el complemento de Schur

F 0 (x∗ )ββ − F 0 (x∗ )βα F 0 (x∗ )−1 0 ∗

αα F (x )αβ

es una P-matriz 2 , x∗ es llamada una solución R-regular.

1
Dada una matriz A = (aij ) de tamaño m × n y conjuntos de ı́ndices η y τ , la matriz Aητ es
aquella con componentes aij tal que i ∈ η y j ∈ τ .
2
Una matriz M ∈ Rn×n es una P-matriz, si para cualquier vector no nulo y ∈ Rn , existe un ı́ndice
i0 = i0 (y) ∈ {1, 2, . . . , n} , tal que yi0 [M y]i0 > 0.
2. Preliminares 10

El siguiente teorema da una condición suficiente para garantizar la no singularidad de

los elementos del jacobiano generalizado de Φλ en una solución del PCNL [21].

Teorema 2.3. [21] Si x∗ ∈ Rn es una solución R-regular del PCNL, entonces las
matrices en el jacobiano generalizado ∂Φλ (x∗ ) son no singulares.

Una consecuencia inmediata de este resultado es dada en el siguiente corolario.

Corolario 2.1. Si x∗ es una solución R-regular del PCNL, entonces x∗ es una solución
BD-regular.

Demostración. Supongamos que x∗ es una solución R-regular del PCNL. Por hipótesis
y el Teorema 2.3, todas la matrices de ∂Φλ (x∗ ) son no singulares. Por otra parte, de
(1.8) tenemos que
∂Φλ (x∗ ) = conv {∂B Φλ (x∗ )},
con lo cual
∂B Φλ (x∗ ) ⊆ ∂Φλ (x∗ ).
Ası́, si H ∈ ∂B Φλ (x∗ ), entonces H ∈ ∂Φλ (x∗ ) y por lo tanto, es no singular. En
consecuencia, x∗ es una solución BD-regular del PCNL. 

Las dos definiciones siguientes son útiles para la presentación del Teorema 2.4, cuya
demostración se puede consultar en [16].

Definición 2.2. [21] Una función G : Rn −→ Rm localmente Lipschitz continua y di-

reccionalmente diferenciable es llamada semisuave en un punto x, si para cada d que
converge al vector cero y para cada matriz H ∈ ∂G(x + d) se tiene que

Hd − G0 (x ; d) = o(kdk).

Definición 2.3. [12] Una función f : Rn −→ R es de clase SC 1 , si es continuamente

diferenciable y su gradiente es una función semisuave.

Para el enunciado del siguiente teorema, supongamos que x∗ es un punto de acumulación

de la sucesión {xk }.
2. Preliminares 11

Teorema 2.4. [16] Sean f : Rn −→ R una función de clase SC 1 y δ > 0 . Si Df (x ; d)

es una aproximación de la derivada direccional de f en x en la dirección d, tal que
para todo xk ∈ B(x∗ ; δ)

∇f (xk )T dk = Df (xk ; dk ) + o(kdk k2 ),

entonces para cualquier σ ∈ 0, 21 , existe un k tal que para todo k > k

f (xk + d) ≤ f (xk ) + σDf (xk ; dk ).

El siguiente resultado presenta una caracterı́stica particular de la función de mérito que

queremos minimizar y cuya demostración se puede consultar en [12].

Teorema 2.5. [12] Sean F : Rn −→ Rn , F (x) = (F1 (x), F2 (x), . . . , Fn (x)) una función
no lineal y continuamente diferenciable, Φλ la función dada en (1.4) asociada a F y
Ψλ la función de mérito (1.10) asociada a Φλ . Si para todo i = 1, 2, . . . , n, la función
Fi es de clase SC 1 , entonces Ψλ es también de clase SC 1 .

Para finalizar, recordemos que la función de mérito dada en (1.10) es diferenciable. En

[21] prueban que el gradiente de esta función se puede calcular usando las matrices
H ∈ ∂Φλ (x) de la siguiente manera

Teorema 2.6. [21] Para cualquier matriz H ∈ ∂Φλ (x), la función de mérito Ψλ dada
en (1.10) es continuamente diferenciable y

∇Ψλ (x) = H T Φλ (x). (2.7)

 CAPÍTULO 3

ALGORITMO Y TEORÍA DE
CONVERGENCIA

En este capı́tulo, proponemos un nuevo algoritmo cuasi Newton global para resolver
el PCNL indirectamente, a través de la solución del problema de minimización (1.11).
Además, para este algoritmo, desarrollamos su teorı́a de convergencia global.

El algoritmo que proponemos es básicamente un algoritmo de búsqueda direccional, cuya

idea geométrica es tomar pasos que conduzcan a un descenso de la función Ψλ . Para
cumplir este objetivo, el paso o dirección utilizada en el algoritmo, siempre que sea
posible, será un paso cuasi Newton como el sugerido en [1] y donde, para seleccionar
el tamaño del paso, usamos una estrategia de búsqueda lineal (monótona) [13] como la
usada en el algoritmo tipo Newton generalizado propuesto en [21].

De acuerdo con lo mencionado en el párrafo anterior, el algoritmo que proponemos a

continuación es una globalización del algoritmo cuasi Newton local propuesto en [1].

12
3. Algoritmo y Teorı́a de Convergencia 13

Algoritmo 3.1. Tipo cuasi Newton global para Φλ (x) = 0

P.0 Inicialización

Escoja λ ∈ (0, 4), x0 ∈ Rn , ρ > 0, µ ∈ (0, 1), σ ∈ (0, 21 ), p > 2, ε ≥ 0, N > 0 y haga
k := 0 y A0 = F 0 (x0 ).

P.1 Criterio de parada

Si kBkT Φλ (xk )k ≤ ε o k > N pare.

P.2 Búsqueda direccional

Calcule Bk como en (2.4) y halle dk ∈ Rn tal que

Bk dk = −Φλ (xk ). (3.1)

Si el sistema no tiene solución o si Φλ (xk )T Bk dk > −ρkdk kp haga dk = −BkT Φλ (xk ).

P.3 Búsqueda lineal

Haga tk := máx{µl | l = 0, 1, 2, ...} tal que

Ψλ (xk + tk dk ) ≤ Ψλ (xk ) + σtk BkT Φλ (xk )dk . (3.2)

P.4 Actualización

Haga xk+1 = xk + tk dk , k = k + 1 , actualize Ak como en (3.4), (3.5) o (3.6) y vuelva

a P.1.

Observemos que la matriz aproximación Bk dada por (2.4) depende a su vez de la

aproximación Ak de la matriz jacobiana de F en xk . Ası́, una pregunta natural que
surge es ¿cómo actualizar la matriz Ak ? Es decir, ¿quién es Ak+1 ?

La forma como se actualice Ak+1 da lugar a un variado abanico de métodos cuasi

Newton, entre los que se encuentran los llamados métodos secantes de cambio mı́nimo,
en los que la actualización de Ak , la matriz Ak+1 , debe satisfacer la llamada ecuación
secante [27]
Ak+1 (xk+1 − xk ) = F (xk+1 ) − F (xk ) (3.3)
3. Algoritmo y Teorı́a de Convergencia 14

y su cambio, medido en alguna norma, con respecto a la matriz Ak , debe ser mı́nimo.
Ejemplos de este tipo de actualizaciones son las siguientes:

1. Actualización “buena” de Broyden [5]

(yk − Bk sk )sTk
Ak+1 = Ak + · (3.4)
sTk sk

2. Actualización “mala” de Broyden [4]

(yk − Bk sk )eTjk Bk
Ak+1 = Ak + (3.5)
eTjk Bk sK

donde jk es definido por kyk k∞ = |eTjk yk | .

3. Actualización de Schubert [10]

F (xk )sTk
Ak+1 = Ak + (3.6)
sTk sk

donde yk = F (xk+1 ) − F (xk ) y sk = xk+1 − xk .

Las tres actualizaciones anteriores las usaremos en el Capı́tulo 4, cuando analicemos el

desempeño numérico de nuestro algoritmo.

Teniendo en cuenta que nuestro deseo es no calcular las matrices Hk en el jacobiano

generalizado de la función Φλ , proponemos algunas variaciones en las condiciones clásicas
de búsqueda direccional y en la búsqueda lineal (condición de Armijo [27]).

Por el Teorema 2.6, para calcular el gradiente de la función Ψλ en cada iteración nece-
sitarı́amos las matrices Hk en el jacobiano generalizado de Φλ en xk , como alternativa
a dicho cálculo proponemos realizar una aproximación a dicho gradiente, esta aproxima-
ción, que denotaremos ∇Ψλ (xk ), requiere de las matrices Bk definidas en (2.4) y de la
evaluación de la función Φλ en xk , explı́citamente, proponemos la siguiente aproxima-
ción a ∇Ψλ (xk ),
∇Ψλ (xk ) = BkT Φλ (xk ). (3.7)
A continuación presentamos las hipótesis bajo las cuales desarrollamos la teorı́a de con-
vergencia global del Algoritmo 3.1.
3.1. Hipótesis 15

3.1. Hipótesis

H1. Existe x∗ ∈ Rn tal que Φλ (x∗ ) = 0.

H2. Las funciones Fi son de clase SC 1 , es decir, son continuamente diferenciables y su

gradiente es una función semisuave.

H3. La matriz jacobiana de F es Lipschitz continua en una vecindad de x∗ .

3.2. Teorı́a de convergencia

A continuación, presentamos dos lemas que serán útiles en la demostración de los teo-
remas de convergencia del Algoritmo 3.1. En el primero de ellos, probamos que bajo
ciertas condiciones, la norma de las matrices que aproximan a las respectivas matrices
del jacobiano generalizado de Φλ , es acotada. Cabe mencionar que, en nuestra demostra-
ción, hacemos uso de la norma infinito para vectores y su respectiva norma inducida para
matrices [32]; sin embargo, teniendo en cuenta que las normas en Rn son equivalentes,
el resultado probado es independiente de la norma usada. En el segundo lema, probamos
una equivalencia entre el gradiente de la función de mérito Ψλ y la aproximación a este
que proponemos en (3.7).

Lema 3.1. Sean B como en (2.4) y δ > 0. Para toda matriz A ∈ B(F 0 (x∗ ); δ), existe
una constante η tal que
kBk ≤ η. (3.8)

Demostración. Consideremos la matriz B definida en (2.4). De la definición de norma

matricial infinito [32],
kBk∞ = máx {k[B]i k1 }. (3.9)
1≤i≤n

Supongamos que el máximo en (3.9) se alcanza en la fila j. Si j 6∈ β entonces xj 6= 0 o

Fj (x) 6= 0, luego por (2.6)
√
|χ(xj , Fj (x))| ≤ 1 y |ψ(xj , Fj (x))| ≤ 2,

además,
3.2 Resultados de Convergencia 16

kBk∞ = k(χ(xj , Fj (x)) − 1)eTj + (ψ(xj , Fj (x)) − 1)[A]j k∞

≤ |χ(xj , Fj (x)) − 1| + |ψ(xj , Fj (x)) − 1|k[A]j k∞

 √ 
≤ 2 + 1 + 2 k[A]j k∞
 √ 
≤ 2 + 1 + 2 (kF 0 (x∗ )k∞ + δ).

Si j ∈ β, siguiendo un procedimiento idéntico, obtenemos la misma cota anterior. Por

lo tanto, existe  √ 
η = 2 + 1 + 2 (kF 0 (x∗ )k∞ + δ)

tal que kBk∞ ≤ η. 

Lema 3.2. Si ∇Ψλ (xk ) es la aproximación a ∇Ψλ (xk ) dada en (3.7), entonces

∇Ψλ (xk )T dk = ∇Ψλ (xk )T dk + o(kdk k2∞ ). (3.10)

Demostración. De (2.7) y (3.7),

|∇Ψλ (xk )T dk − ∇Ψλ (xk )T dk | = |Φλ (xk )T Hk dk − Φλ (xk )T Bk dk |

(3.11)
= |Φλ (xk )T (Hk − Bk )dk |.

Usando la desigualdad de Cauchy-Schwarz [36] y la relación entre las normas vectoriales

dos e infinito en (3.11), obtenemos la siguiente desigualdad

|∇Ψλ (xk )T dk − ∇Ψλ (xk )T dk | ≤ n kΦλ (x k )k∞ kHk − Bk k∞ kdk k∞ .

Por el Teorema 2.1, para un k suficientemente grande, existe una constante θ tal que
kHk − Bk k∞ ≤ θ , en consecuencia

|∇Ψλ (xk )T dk − ∇Ψλ (xk )T dk | ≤ n θ kΦλ (x k )k∞ kdk k∞

de donde,
|∇Ψλ (xk )T dk − ∇Ψλ (xk )T dk | n θ kΦλ (x k )k∞
0≤ 2
≤ ,
kdk k∞ kdk k∞
3.2 Resultados de Convergencia 17

por lo tanto, teniendo en cuenta que con el Algoritmo 3.1 se busca minimizar la
función de mérito Ψλ y que dk es una dirección de descenso, podemos inferir que
aunque la sucesión {kdk k} converja a cero, la sucesión {kΦλ (xk )k} lo hará más rápido,
en consecuencia
|∇Ψλ (xk )T dk − ∇Ψλ (xk )T dk |
lı́m =0
k→∞ kdk k2∞
lo que demuestra (3.10). 

Los siguientes teoremas de convergencia son análogos a los presentados para un método
tipo Newton (no suave) que busque resolver PCNL.

Teorema 3.1. Todo punto de acumulación de la sucesión {xk } generada por el Algo-
ritmo 3.1 es un punto estacionario de Ψλ .

Demostración. Sea x∗ un punto de acumulación de la sucesión {xk } generada por el

Algoritmo 3.1.

Si para un conjunto infinito de ı́ndices J, la dirección dada en el paso 2 del Algoritmo

3.1 es siempre dk = −∇Ψλ (xk ), entonces cualquier punto lı́mite de la subsucesión
{xk }J es un punto estacionario de Ψλ [2]. Luego, sin pérdida de generalidad, podemos
suponer que la dirección siempre está dada por

Bk dk = −Φλ (xk ). (3.12)

Supongamos que x∗ no es un punto estacionario de Ψλ , es decir, ∇Ψλ (x∗ ) no es el

vector cero. Demostraremos que existen constantes positivas m y M tales que

m ≤ kdk k ≤ M. (3.13)

En efecto, de (3.12)
kΦλ (xk )k ≤ kBk kkdk k
de donde
kΦλ (xk )k
≤ kdk k
kBk k
por tanto, si la sucesión {dk } converge al vector cero, entonces la sucesión {kΦλ (xk )k}
converge a cero ya que, por el Lema 3.1, kBk k es acotada para todo k , en consecuencia,
la sucesión {Φλ (xk )} converge al vector cero, luego ∇Ψλ (xk ) también converge al vector
3.2 Resultados de Convergencia 18

cero y dado que la función Φλ es continua, entonces

∇Ψλ (x∗ ) = 0,

lo que contradice nuestro supuesto. Es decir, existe una constante positiva m tal que
m ≤ kdk k.

Por otra parte, si kdk k no es acotada superiormente, teniendo en cuenta que ∇Ψλ (xk )
es acotado y que p ≥ 2, entonces

k∇Ψλ (xk )k cos θ

lı́m =0
k→∞ kdk kp−1

donde θ es el ángulo formado por los vectores ∇Ψλ (xk ) y dk . Equivalentemente,

∇Ψλ (xk )T dk
lı́m =0
k→∞ kdk kp

lo que contradice el hecho que

∇Ψλ (xk )T dk ≤ −ρkdk kp . (3.14)

Esto significa que existe una constante positiva M tal que kdk k ≤ M .

Ahora, dado que la sucesión {xk } converge a x∗ y que la iteración siguiente en el

Algoritmo 3.1 es generada por xk+1 = xk + tk dk entonces

tk dk → 0. (3.15)

Por (3.13), podemos asumir que si k → ∞, dk → d 6= 0, luego tk = µlk → 0. Además,

de (3.2)
Ψλ (xk + µlk −1 dk ) − Ψλ (xk )
l −1
> σ∇Ψλ (xk )T dk (3.16)
µ k

Ası́, de (3.16), si k → ∞

∇Ψλ (x∗ )T d ≥ σ ∇Ψλ (x∗ )T d

luego, por el Lema (3.2) podemos inferir que

σ ≥ 1,

lo que contradice el hecho que σ ∈ 0, 12 . Ası́, de (3.15) tenemos que, necesariamente

3.2 Resultados de Convergencia 19

la sucesión {dk } converge al vector cero y por lo tanto, la sucesión {kdk k} converge a
cero, lo cual contradice el hecho que ∇Ψλ (x∗ ) 6= 0. En conclusión, ∇Ψλ (x∗ ) = 0. 

En el Teorema 3.2, hacemos referencia a un término aparentemente contradictorio, por

ello y para evitar ambigüedades realizamos la siguiente definición.

Definición 3.1. Sea Ω el conjunto de puntos de acumulación de la sucesión {xk }.

Decimos que x∗ ∈ Ω es un punto de acumulación aislado, si x∗ es un punto aislado en
Ω, es decir, si existe δ > 0 tal que

B(x∗ ; δ) ∩ Ω = {x∗ }.

Teorema 3.2. Si x∗ es un punto de acumulación aislado de la sucesión generada por

el Algoritmo 3.1, entonces la sucesión converge a x∗ .

Demostración. Sean {xk } la sucesión generada por el Algoritmo 3.1 y x∗ un punto

de acumulación aislado de esta sucesión. Como Ψ(x∗ ) = 0 y Ψ es una función convexa,
entonces x∗ es un minimizador global de Ψ.

Sea Ω el conjunto de puntos de acumulación de {xk }, entonces Ω 6= ∅, pues x∗ ∈ Ω .

Definamos
(
dist(x∗ ; Ω\{x∗ }) si Ω\{x∗ } =6 ∅
δ=
1 si Ω = {x∗ }
donde
dist(a, A) = ı́nf ka − xk
x∈A

indica la distancia del punto a al conjunto A .

Teniendo en cuenta que x∗ es único, localmente, entonces δ > 0. Si

 
n δ
Ω1 = y ∈ R : dist(y; Ω) ≤ ,
4

entonces existe k tal que xk ∈ Ω1 para todo k ≥ k.

Ahora, si  
∗ δ
Λ= k ∈ N : kxk − x k ≤ ,
4
3.2 Resultados de Convergencia 20

entonces {xk }Λ ⊂ B x∗ ; 4δ . Como x∗ es el único punto de acumulación de {xk } en

B x∗ ; 4δ , podemos concluir que {xk }Λ converge a x∗ y como Ψ(x∗ ) = 0, entonces

{kΨ(xk )k}Λ converge a cero, lo que a su vez, por (3.14) y (3.10) implica que la sucesión
{dk } converge al vector cero. Ası́, existe e
k tal que si k ∈ Λ y k ≥ e
k ≥ k entonces

δ
kdk k ≤ ·
4

k ∈ Λ tal que b
Sea b k≥e
k , observemos que

xbk+1 = xbk + tbk dbk (3.17)

con tbk ∈ (0, 1]. Luego

δ
kxbk+1 − xbk k ≤ kdbk k ≤ ,
4
por lo tanto,

dist(x∗ ; Ω\{x∗ }) ≤ dist(xbk+1 ; Ω\{x∗ }) + kx∗ − xbk+1 k

≤ dist(xbk+1 ; Ω\{x∗ }) + kx∗ − xbk k + kxbk − xbk+1 k
δ δ
≤ dist(xbk+1 ; Ω\{x∗ }) + +
4 4
es decir,
δ δ δ
dist(xbk+1 ; Ω\{x∗ }) ≥ dist(x∗ ; Ω\{x∗ }) − ≥ δ − =
2 2 2
   
∗ δ ∗ δ
en consecuencia, xbk+1 6∈ Ω1 \B x ; y como xbk+1 ∈ Ω1 entonces xbk+1 ∈ B x ; ,
4 4
k + 1 ∈ Λ y dado queb
es decir, b k + 1>e k, entonces por inducción, bk ∈ Λ para todo
δ
k>e
b k, en consecuencia, xk ∈ B x∗ ; k, ası́, {xk } converge a x∗ . 
para todo k ≥ e
4

El siguiente resultado muestra que existe un entero positivo k tal que para todo k > k, el
Algoritmo 3.1 tiene un comportamiento local idéntico al del Algoritmo 4.1 propuesto
en [1] lo que nos permitirá garantizar una tasa de convergencia similar a la obtenida
para dicho algoritmo. Cabe mencionar que si x∗ es R-regular, la hipótesis H3 en [1] es
inmediatamente satisfecha (teorema 2.3).
3.2 Resultados de Convergencia 21

Teorema 3.3. Suponga que x∗ es un punto de acumulación de la sucesión {xk } gene-

rada por el Algoritmo 3.1 tal que x∗ es una solución R-regular del PCNL, entonces
la sucesión {xk } converge a x∗ , la dirección de descenso es eventualmente dada por la
solución del sistema de ecuaciones Bk dk = −Φλ (xk ) y el tamaño del paso tk = 1 es
aceptado para todo k suficientemente grande.

Demostración. Si x∗ es una solución R-regular del PCNL, x∗ es también una solución

BD-regular (Corolario 2.1), luego por la Proposición 3 en [28], x∗ es un punto de acu-
mulación aislado de la sucesión {xk } generada por el Algoritmo 3.1 y en consecuencia,
del Teorema 3.2, se sigue que la sucesión {xk } converge a x∗ .

Por otra parte, por el Teorema 2.2, para algún k suficientemente grande, la matriz
Bk−1 existe y por lo tanto, el sistema de ecuaciones

Bk dk = −Φλ (xk )

tiene solución. Ası́, si dk es una solución de dicho sistema,

dk = −Bk−1 Φλ (xk ), (3.18)

luego, para cualquier norma vectorial y su respectiva norma matricial inducida

kdk k ≤ kΦλ (xk )kkBk−1 k. (3.19)

De (3.7), (3.18) y (3.19)

∇Ψλ (xk )T dk = Φλ (xk )T Bk dk

= −Φλ (xk )T Bk Bk−1 Φλ (xk )
= −kΦλ (xk )k2
kdk k2
≤ − −1 2
kBk k

luego, podemos tomar la constante ρ = η 2 , donde η es la cota superior de kBk k dada

en el Teorema 2.2. Es decir, teniendo en cuenta que la sucesión {kdk k} converge a
cero entonces
∇Ψλ (xk )dk ≤ −ρkdk kp
para todo p ≥ 2. Esto prueba que para un k suficientemente grande, la dirección de
descenso es siempre dada por (3.13).

Ahora, por el Lema 3.2, el Teorema 2.4 y teniendo en cuenta que Ψλ es una función
3.2 Resultados de Convergencia 22

1


de clase SC 1 , se tiene que para cualquier σ ∈ 0, 2
existe k tal que para todo k > k

Ψλ (xk + dk ) ≤ Ψλ (xk ) + σBkT Φλ (xk )dk

es decir, el tamaño del paso tk = 1 es aceptado para k suficientemente grande. Esto

completa la demostración. 
 CAPÍTULO 4

PRUEBAS NUMÉRICAS

En la primera parte de este capı́tulo, analizamos numéricamente el algoritmo propuesto

en el capı́tulo anterior (Algoritmo 3.1). Para ello, comparamos su desempeño con el
algoritmo tipo Newton generalizado presentado en [21] que llamaremos Algoritmo 4.1
el cual, para mayor claridad en la lectura de este documento, incluimos a continuación.

Algoritmo 4.1. Newton generalizado para Φλ (x) = 0.

P.0 Inicialización:

Escoja λ ∈ (0, 4), x0 ∈ Rn , ρ > 0, µ ∈ (0, 1), σ ∈ (0, 21 ), p > 2, ε ≥ 0, N > 0 y haga
k := 0.

P.1 Criterio de parada:

Si k∇Ψλ (xk )k ≤ ε o k > N pare.

23
4. Pruebas numéricas 24

P.2 Búsqueda direccional:

Calcule Hk ∈ ∂Φλ (xk ) y halle dk ∈ Rn tal que

Hk dk = −Φλ (xk ).

Si este sistema no tiene solución o si ∇Ψλ (xk )T dk > −ρkdk kp haga dk := −∇Ψλ (xk ).

P.3 Búsqueda lineal:

Haga tk := máx{µl | l = 0, 1, 2, ...} tal que

Ψλ (xk + tk dk ) ≤ Ψλ (xk ) + σtk ∇Ψλ (xk )T dk

P.4 Actualización:

Haga xk+1 = xk + tk dk , k = k + 1 y vuelva a P.1.

Como se puede ver, este algoritmo usa en cada iteración las matrices Hk ∈ ∂Φλ (x) . Es
básicamente en esta etapa donde radica su diferencia con el Algoritmo 3.1 que propo-
nemos en el capı́tulo anterior, ya que en lugar de calcular dichas matrices, calculamos
aproximaciones a ellas como la dada en (2.4). De (2.7) y por lo mencionado anteriormen-
te, es imposible calcular el gradiente de la función de mérito Ψλ en cada iteración con
el Algoritmo 3.1, motivo por el cual precisamente, proponemos la aproximación (3.7)
a dicho gradiente para ser usada en nuestro algoritmo.

Motivados por el buen desempeño numérico de la búsqueda lineal no monótona [20] en

algoritmos de complementariedad no lineal y problemas relacionados [21][30], en la se-
gunda parte de este capı́tulo, presentamos el Algoritmo 4.2 en el cual, modificamos el
paso P.3 del Algoritmo 3.1 e introducimos una búsqueda lineal no monótona. Anali-
zamos numéricamente el Algoritmo 4.2 comparando su desempeño con el respectivo
Algoritmo 4.1 tras hacerle la misma modificación y además realizamos pruebas numéri-
cas con el Algoritmo 4.2 al variar los parámetros propios del tipo de búsqueda utilizada
en el mismo.

Finalmente, teniendo en cuenta el buen desempeño local que han mostrado los algoritmos
de complementariedad que usan la función mı́nimo y el buen desempeño global que han
mostrado los algoritmos de complementariedad que usan la función de Fischer [31] [21],
en la tercera parte de este capı́tulo realizamos una exploración numérica preliminar del
comportamiento del Algoritmo 3.1 al variar el parámetro inicial λ que define la función
4. Pruebas numéricas 25

de complementariedad dada en (1.3).

Para las pruebas numéricas, consideramos cuatro problemas de complementariedad no

lineal asociados respectivamente a las siguientes funciones:

1. Kojima Shindo (Koj-Shi): F : R4 → R4 definido por

 2 
3x1 + 2x1 x2 + 2x22 + x3 + 3x4 − 6
 2x21 + x22 + x1 + 10x3 + 2x4 − 2 
F (x) =  ·
 3x21 + x1 x2 + 2x22 + 2x3 + 9x4 − 9 
x21 + 3x22 + 2x3 + 3x4 − 3

2. Kojima Josephy (Koj-Jo): F : R4 → R4 definido por

 2 
3x1 + 2x1 x2 + 2x22 + x3 + 3x4 − 6
 2x21 + x22 + x1 + 3x3 + 2x4 − 2 
F (x) =  ·
 3x21 + x1 x2 + 2x22 + 2x3 + 3x4 − 9 
x21 + 3x22 + 2x3 + 3x4 − 3

3. Mathiesen Modificado (Mathiesen): F : R4 → R4 definido por

 
−x2 + x3 + x4
 
4.5x3 + 2.7x4 
 
 x1 −


 x 2 + 1 
F (x) =  ·
 

 5−x − 0.5x 3 + 0.3x 
4 
1
x3 + 1
 
 
 
3 − x1

Cabe mencionar que este problema tiene infinitas soluciones de la forma

x∗ = (a 0 0 0 ),

con a ∈ [0, 3].

4. Billups (Billups): F : R → R definido por

F (x) = (x − 1)2 − 1.1.

4. Pruebas numéricas 26

Estos problemas, tomados de [28] y [3], son considerados problemas difı́ciles en cuanto a
su convergencia y por ello, de uso obligado para probar algoritmos de complementariedad
no lineal. Cabe destacar que los problemas de Kojima Shindo (aplicación en Economı́a al
problema de equilibrio económico [22]), Kojima Josephy y Mathiesen Modificado (apli-
cación al problema de equilibrio económico walrasiano [26]) son analizados localmente
en [1]. Además, es de subrayar que el problema de Billups fue construido por Billups [3]
con el objetivo de hacer que los métodos usados para resolverlo fallen [21], es decir, es
un problema de alta exigencia en cuanto a su convergencia.

Para escribir los códigos de los algoritmos y de las funciones de prueba, utilizamos el soft-
ware Matlabr y realizamos las pruebas numéricas en un computador con procesador
Intel(R) Atom(TM) de 1.67 GHz.

En la implementación de los Algoritmos 3.1 y 4.1, para cada problema, usamos como
puntos iniciales los utilizados en [21] para sus pruebas. Además, usamos los siguientes
valores para los parámetros requeridos

ρ = 10−8 , σ = 10−4 , p = 2.1, ε = 10−12 , N = 200, µ = 0.5.

De igual forma, teniendo en cuenta los resultados presentados en [21], el buen desempeño
que ha mostrado la función de Fischer cuando se ejecutan algoritmos iterativos con
puntos iniciales alejados de la solución y la eficiencia local de la función mı́nimo para
este tipo de algoritmos, implementamos una escogencia dinámica para el parámetro λ,
el cual se actualiza en cada iteración de la siguiente manera:

1. Escoja λ = 2

2. Si Ψλ (xk ) ≤ γ1
λ := Ψλ (xk ),
si no,
λ := mı́n{c1 Ψ(xk ), λ}

3. Si Ψ(xk ) ≤ γ2 ,
λ := mı́n{c2 , λ}

donde
γ1 = 10−2 , γ2 = 10−4 , c1 = 10, c2 = 10−8 .
Esta elección del parámetro λ permite que los algoritmos realicen sus primeras iteracio-
nes usando la función de complementariedad de Fischer [19][18], la cual ha mostrado un
4.1. Pruebas con los Algoritmos 3.1 y 4.1 27

mejor comportamiento cuando el punto inicial en los algoritmos está alejado de la solu-
ción y culmine su ciclo con funciones de complementariedad que tienden a ser múltiplos
de la función mı́nimo [34][29], la cual ha mostrado un buen desempeño local [34].

Presentamos los resultados de las pruebas numéricas en tablas cuyas columnas contienen
la siguiente información:

Problema: Nombre del problema.

t: Tiempo (en segundos) de ejecución de los Algoritmos.
x0 : Punto inicial.
k: Número de iteraciones.
x∗ : Solución.

De igual manera, incluimos una columna en la que hacemos alusión al algoritmo y la

actualización usada:

NG: Algoritmo 4.1

BB: Algoritmo 3.1 y actualización (3.4) de Broyden “bueno”.
BM: Algoritmo 3.1 y actualización (3.5) de Broyden “malo”.
SN: Algoritmo 3.1 y actualización (3.6) de Shubert.
-: indica que hubo divergencia.

4.1. Pruebas con los Algoritmos 3.1 y 4.1

En la siguiente tabla, encontramos los resultados obtenidos al resolver los problemas de

complementariedad anteriormente expuestos con los Algoritmos 4.1 y 3.1, este último
con las actualizaciones dadas en (3.4), (3.5) y (3.6), respectivamente.

Problema Algoritmo x0 t (seg) k x∗

Koj-Jo NG (0 0 0 0)T 0.152 21 (1 0 3 0)T
Koj-Jo BB (0 0 0 0)T 0.141 21 (1 0 3 0)T
Koj-Jo BM (0 0 0 0)T - - -
Koj-Jo SN (0 0 0 0)T - - -
Koj-Jo NG (1 1 1 1)T 0.060 8 (1 0 3 0)T
Koj-Jo BB (1 1 1 1)T 0.086 13 (1 0 3 0)T
Koj-Jo BM (1 1 1 1)T 0.074 12 (1 0 3 0)T
Koj-Jo SN (1 1 1 1)T - - -
4.1 Algoritmos 3.1 y 4.1 28

Problema Algoritmo x0 t (seg) k x∗

Koj-Jo NG (100 100 100 100)T - - -
Koj-Jo BB (100 100 100 100)T - - -
Koj-Jo BM (100 100 100 100)T - - -
Koj-Jo SN (100 100 100 100)T - - -
Koj-Jo NG (1 0 1 0)T 0.038 6 (1 0 3 0)T
Koj-Jo BB (1 0 1 0)T 0.051 8 (1 0 3 0)T
Koj-Jo BM (1 0 1 0)T 0.056 8 (1 0 3 0)T
Koj-Jo SN (1 0 1 0)T - - -
Koj-Jo NG (1 0 0 0)T 0.082 8 (1 0 3 0)T
Koj-Jo BB (1 0 0 0)T 0.067 10 (1 0 3 0)T
Koj-Jo BM (1 0 0 0)T 0.087 12 (1 0 3 0)T
Koj-Jo SN (1 0 0 0)T - - -
Koj-Jo NG (0 1 1 0)T 0.066 9 (1 0 3 0)T
Koj-Jo BB (0 1 1 0)T 0.089 14 (1 0 3 0)T
Koj-Jo BM (0 1 1 0)T 0.229 17 (1 0 3 0)T
Koj-Jo SN (0 1 1 0)T - - -
Koj-Shi NG (0 0 0 0)T 0.212 15 (1 0 3 0)T
Koj-Shi BB (0 0 0 0)T 0.131 21 (1 0 3 0)T
Koj-Shi BM (0 0 0 0)T - - -
Koj-Shi SN (0 0 0 0)T - - -
Koj-Shi GN (1 1 1 1)T 0.067 10 (1 0 3 0)T
Koj-Shi BB (1 1 1 1)T 0.141 16 (1 0 3 0)T
Koj-Shi BM (1 1 1 1)T 0.102 16 (1 0 3 0)T
Koj-Shi SN (1 1 1 1)T - - -
√ T
6 1
Koj-Shi GN (100 100 100 100)T 0.076 12 2
0 0 2
Koj-Shi BB (100 100 100 100)T - - -
Koj-Shi BM (100 100 100 100)T - - -
Koj-Shi SN (100 100 100 100)T - - -
Koj-Shi GN (1 0 1 0)T 0.039 6 (1 0 3 0)T
Koj-Shi BB (1 0 1 0)T 0.046 7 (1 0 3 0)T
Koj-Shi BM (1 0 1 0)T 0.044 7 (1 0 3 0)T
√ T
6 1
Koj-Shi GN (1 0 0 0)T 0.051 8 2
00 2
√ T
6 1
Koj-Shi BB (1 0 0 0)T 0.060 9 2
00 2
√ T
6 1
Koj-Shi BM (1 0 0 0)T 0.057 9 2
00 2
T
Koj-Shi SN (1 0 0 0) - - -
4.1 Algoritmos 3.1 y 4.1 29

Problema Algoritmo x0 t (seg) k x∗

Koj-Shi GN (0 1 1 0)T 0.070 10 (1 0 3 0)T
Koj-Shi BB (0 1 1 0)T 0.093 15 (1 0 3 0)T
Koj-Shi BM (0 1 1 0)T 0.125 15 (1 0 3 0)T
Koj-Shi SN (0 1 1 0)T - - -
Mathiesen GN (1 1 1 1)T 0.186 11 (0.0302 0 0 0)T
Mathiesen BB (1 1 1 1)T 0.178 11 (0.0302 0 0 0)T
Mathiesen BM (1 1 1 1)T 0.091 11 (0.0302 0 0 0)T
Mathiesen SN (1 1 1 1)T 0.070 5 (1.497 0 0 0)T
Mathiesen GN (100 100 100 100)T 0.102 13 (2.999 0 0 0)T
Mathiesen BB (100 100 100 100)T 0.132 7 (3 0 0 0)T
Mathiesen BM (100 100 100 100)T 0.054 7 (3 0 0 0)T
Mathiesen SN (100 100 100 100)T 0.094 9 (2.320 0 0 0)T
Mathiesen GN (1 0 1 0)T 0.091 4 (0 0 0 0)T
Mathiesen BB (1 0 1 0)T 0.046 6 (0 0 0 0)T
Mathiesen BM (1 0 1 0)T 0.043 6 (0 0 0 0)T
Mathiesen SN (1 0 1 0)T 0.071 6 (0.018 0 0 0)T
Mathiesen GN (0 1 1 0)T 0.037 5 (0.752 0 0 0)T
Mathiesen BB (0 1 1 0)T 0.036 5 (0.608 0 0 0)T
Mathiesen BM (0 1 1 0)T 0.036 5 (0.647 0 0 0)T
Mathiesen SN (0 1 1 0)T 0.037 5 (0.778 0 0 0)T
Billups NG 0 - - -
Billups BB 0 2.843 16 2.0488
Billups BM 0 2.827 16 2.0488
Billups SN 0 - - -

Tabla 4.1: Resultados numéricos para los Algoritmos 3.1 y 4.1.

Como se puede ver, estos resultados muestran la importancia de tener en cuenta el

tiempo de ejecución de los algoritmos ya que, evaluar solo el número de iteraciones
realizadas por cada uno de ellos no es un buen comparativo pues, si bien en general el
número de iteraciones realizadas por el Algoritmo 3.1 es mayor que las realizadas por
el Algoritmo 4.1, el número de operaciones que realiza en cada iteración el Algoritmo
3.1 es en general menor que el número de operaciones que realiza en cada iteración el
Algoritmo 4.1.

Prueba de lo mencionado anteriormente son los resultados obtenidos para el problema de

Kojima-Shindo con punto inicial (0 0 0 0), ya que el Algoritmo 4.1 resuelve el problema
4.2. Búsqueda lineal no monótona 30

en 15 iteraciones, mientras que el Algoritmo 3.1 lo resuelve en 21 iteraciones, sin

embargo, si comparamos los tiempos de ejecución vemos que el algoritmo propuesto
en esta investigación fue un poco más rápido que el Algoritmo 4.1. Esto fue algo
que sucedió en varias ocasiones en las pruebas pero que sobresalió con el problema
de Mathiesen pues en todas las ocasiones el Algoritmo 3.1 fue más rápido que el
Algoritmo 4.1.

Por otra parte, los resultados muestran el buen desempeño que tiene el Algoritmo 3.1
cuando se utiliza la actualización de Broyden “buena” pues convergió para un 93 % de los
problemas y puntos iniciales, sin contar los resultados obtenidos al resolver el problema
de Billups, pues en este caso, el Algoritmo 4.1 divergió, entre tanto, el Algoritmo 3.1
convergió.

4.2. Búsqueda lineal no monótona

Como lo mencionamos en la primera parte de este capı́tulo, motivados en el buen compor-

tamiento de los algoritmos de complementariedad que usan búsqueda lineal no monótona,
en esta parte, analizamos el efecto, en cuanto a convergencia y velocidad de convergen-
cia, de introducir una búsqueda lineal no monótona en los Algoritmos 3.1 y 4.1. Para
mayor claridad en la lectura, mostramos a continuación la estructura del nuevo algoritmo.

Algoritmo 4.2. Tipo cuasi Newton globalizado con búsqueda no monótona

P.0 Inicialización

Escoja λ ∈ (0, 4), x0 ∈ Rn , ρ > 0, µ ∈ (0, 1), σ ∈ (0, 12 ), p > 2, ε ≥ 0, N > 0 s > 0,
M > 0 y haga k = 0 mk = 0 y A0 = F 0 (x0 )

P.1 Criterio de parada

Si kBkT Φλ (xk )k ≤ ε o k > N pare.

P.2 Búsqueda direccional

Calcule Bk como en (2.4) y halle dk ∈ Rn tal que

Bk dk = −Φλ (xk ).

Si el sistema no tiene solución o si Φλ (xk )T Bk dk > −ρkdk kp , haga dk = −BkT Φλ (xk ).

4.2 Búsqueda lineal no monótona 31

P.3 Búsqueda lineal no monótona

Si k ≤ s o dk = −BkT Φλ (xk ),

haga mk = 0

si no

haga mk = mı́n{mk−1 + 1, M }

Calcule
Rk = máx Ψλ (xj )
k−mk ≤j≤k

Haga tk = máx{µl | l = 0, 1, 2, ...} tal que

Ψλ (xk + tk dk ) ≤ Rk + σtk BkT Φλ (xk )T dk .

P.4 Actualización

Haga xk+1 = xk + tk dk , k = k + 1 , actualize Ak como en (3.4), (3.5) o (3.6) y vuelva

a P.1.

Cabe destacar que con la búsqueda lineal no monótona usada en el Algoritmo 4.2, es
posible comparar el valor que toma la función a minimizar en la iteración actual con el
máximo valor que toma la función en alguna iteraciones anteriores con el fin de seleccionar
el tamaño del paso, es decir, es posible que la función no decrezca entre dos iteraciones
consecutivas. El número de iteraciones anteriores con la cuales se puede comparar el valor
de la función en la iteración actual depende del parámetro M, por lo tanto, si M = 0
tendremos una búsqueda lineal monótona. De igual manera, el parámetro s indica el
número de iteraciones iniciales con las que el Algoritmo 4.2 realizará una búsqueda
lineal monótona.

Siguiendo la propuesta hecha en [21], decidimos comparar el valor que toma la función
de prueba en la iteración actual hasta con ocho de las iteraciones inmediatamente an-
teriores, es decir, tomamos M = 8, además, pedimos al Algoritmo 4.2 que realizara
una búsqueda lineal monótona sólo en la primera iteración, es decir, tomamos s = 1.
Los resultados obtenidos se muestran en la siguiente tabla.
4.2 Búsqueda lineal no monótona 32

Problema Algoritmo x0 t (seg) k x∗

Koj-Jo NG (0 0 0 0)T 1.57 134 (1 0 3 0)T
Koj-Jo BB (0 0 0 0)T 0.98 76 (1 0 3 0)T
Koj-Jo BM (0 0 0 0)T - - -
Koj-Jo SN (0 0 0 0)T - - -
Koj-Jo NG (1 1 1 1)T 0.096 10 (1 0 3 0)T
Koj-Jo BB (1 1 1 1)T 0.129 13 (1 0 3 0)T
Koj-Jo BM (1 1 1 1)T 0.122 12 (1 0 3 0)T
Koj-Jo SN (1 1 1 1)T - - -
Koj-Jo NG (100 100 100 100)T 1.53 131 (1 0 3 0)T
Koj-Jo BB (100 100 100 100)T - - -
Koj-Jo BM (100 100 100 100)T - - -
Koj-Jo SN (100 100 100 100)T - - -
Koj-Jo NG (1 0 1 0)T 0.052 6 (1 0 3 0)T
Koj-Jo BB (1 0 1 0)T 0.068 8 (1 0 3 0)T
Koj-Jo BM (1 0 1 0)T 0.086 8 (1 0 3 0)T
Koj-Jo SN (1 0 1 0)T - - -
Koj-Jo NG (1 0 0 0)T 0.082 8 (1 0 3 0)T
Koj-Jo BB (1 0 0 0)T 0.098 10 (1 0 3 0)T
Koj-Jo BM (1 0 0 0)T 0.110 11 (1 0 3 0)T
Koj-Jo SN (1 0 0 0)T - - -
Koj-Jo NG (0 1 1 0)T 0.072 8 (1 0 3 0)T
Koj-Jo BB (0 1 1 0)T 0.132 13 (1 0 3 0)T
Koj-Jo BM (0 1 1 0)T 0.154 15 (1 0 3 0)T
Koj-Jo SN (0 1 1 0)T - - -
√ T
6 1
Koj-Shi NG (0 0 0 0)T 0.161 15 2
0 0 2
√ T
6 1
Koj-Shi BB (0 0 0 0)T 0.494 44 2
0 0 2
T
Koj-Shi BM (0 0 0 0) - - -
Koj-Shi SN (0 0 0 0)T - - -
Koj-Shi GN (1 1 1 1)T 0.080 9 (1 0 3 0)T
Koj-Shi BB (1 1 1 1)T 0.149 14 (1 0 3 0)T
Koj-Shi BM (1 1 1 1)T - - -
Koj-Shi SN (1 1 1 1)T - - -
√ T
6 1
Koj-Shi GN (100 100 100 100)T 0.141 13 2
0 0 2
√ T
6 1
Koj-Shi BB (100 100 100 100)T 1.23 102 2
0 0 2
4.2 Búsqueda lineal no monótona 33

Problema Algoritmo x0 t (seg) k x∗

Koj-Shi BM (100 100 100 100)T - - -
Koj-Shi SN (100 100 100 100)T - - -
Koj-Shi GN (1 0 1 0)T 0.054 6 (1 0 3 0)T
Koj-Shi BB (1 0 1 0)T 0.061 7 (1 0 3 0)T
Koj-Shi BM (1 0 1 0)T 0.065 7 (1 0 3 0)T
√ T
6 1
Koj-Shi GN (1 0 0 0)T 0.054 6 2
00 2
√ T
6 1
Koj-Shi BB (1 0 0 0)T 0.080 9 2
00 2
√ T
6 1
Koj-Shi BM (1 0 0 0)T 0.071 8 2
0 0 2
T
Koj-Shi GN (0 1 1 0) 0.091 10 (1 0 3 0)T
Koj-Shi BB (0 1 1 0)T 0.155 15 (1 0 3 0)T
Koj-Shi BM (0 1 1 0)T 0.173 17 (1 0 3 0)T
Koj-Shi SN (0 1 1 0)T - - -
Mathiesen GN (1 1 1 1)T 0.094 11 (0.0302 0 0 0)T
Mathiesen BB (1 1 1 1)T 0.125 11 (0.0302 0 0 0)T
Mathiesen BM (1 1 1 1)T 0.105 11 (0.0302 0 0 0)T
Mathiesen SN (1 1 1 1)T 0.053 5 (1.497 0 0 0)T
Mathiesen GN (100 100 100 100)T 0.090 9 (3 0 0 0)T
Mathiesen BB (100 100 100 100)T 0.077 8 (2.999 0 0 0)T
Mathiesen BM (100 100 100 100)T 0.073 8 (3 0 0 0)T
Mathiesen SN (100 100 100 100)T 0.639 9 (2.320 0 0 0)T
Mathiesen GN (1 0 1 0)T 0.089 11 (0 0 0 0)T
Mathiesen BB (1 0 1 0)T 0.054 6 (0 0 0 0)T
Mathiesen BM (1 0 1 0)T 0.054 6 (0 0 0 0)T
Mathiesen SN (1 0 1 0)T - - -
Mathiesen GN (0 1 1 0)T 0.041 5 (0.752 0 0 0)T
Mathiesen BB (0 1 1 0)T 0.029 4 (2.087 0 0 0)T
Mathiesen BM (0 1 1 0)T 0.065 7 (0 0 0 0)T
Mathiesen SN (0 1 1 0)T 0.042 5 (0.778 0 0 0)T
Billups NG 0 - - -
Billups BB 0 - - -
Billups BM 0 - - -
Billups SN 0 - - -

Tabla 4.2: Resultados del Algoritmo 4.2 con búsqueda no monótona.

4.2 Búsqueda lineal no monótona 34

Observemos que el Algoritmo 4.2 en general tuvo un desempeño similar al Algoritmo

3.1 que usa búsqueda lineal monótona, sin embargo, hay que destacar que aunque la
diferencia en tiempo de ejecución entre un algoritmo y otro es relativamente pequeña,
esta siempre es a favor del Algoritmo 3.1, algo que no ocurre cuando se implementa el
Algoritmo 4.1 haciendo búsqueda lineal no monótona en el paso P.3, ya que en este
caso la diferencia en tiempo de ejecución es a favor del algoritmo con búsqueda lineal no
monótona.

Un caso que llama la atención es el problema de Kojima-Shindo con punto inicial

(100 100 100 100)T , pues el Algoritmo 3.1 diverge cuando se implementa la actua-
lización de Broyden “buena” dada en (3.4), sin embargo, el Algoritmo 4.2 converge
bajo estas mismas condiciones. Por el contrario, el problema de Billups no fue resuelto
por el Algoritmo 4.2 mientras que el Algoritmo 3.1 si lo hizo.

Estos comportamientos nos motivaron a explorar el desempeño del Algoritmo 4.2,

cuando se varı́a el valor del parámetro M.

A continuación, presentamos algunos resultados relevantes, para los cuales cabe mencio-
nar que tomamos s = 1.

Problema Algoritmo M x0 t (seg) k x∗

Mathiesen SN 0 (500 − 250 500 − 250)T - - -
Mathiesen SN 2 (500 − 250 500 − 250)T 0.075 8 (1.100 0 0 0)T
Mathiesen SN 5 (500 − 250 500 − 250)T 0.068 8 (1.100 0 0 0)T
Mathiesen SN 10 (500 − 250 500 − 250)T 0.070 8 (1.100 0 0 0)T
Koj-shi BB 0 (100 100 100 100)T - - -
Koj-shi BB 2 (100 100 100 100)T - - -
Koj-shi BB 5 (100 100 100 100)T - - -
√ T
6 1
Koj-shi BB 10 (100 100 100 100)T 1.63 121 2
00 2
T
Koj-jo BB 0 (0 0 0 ) 0.145 21 (1 0 3 0)T
Koj-jo BB 2 (0 0 0 )T - - -
Koj-jo BB 5 (0 0 0 )T 0.309 32 (1 0 3 0)T
Koj-jo BB 10 (0 0 0 )T 1.543 117 (1 0 3 0)T
Koj-jo BB 0 (0 1 1 0)T 0.093 15 (1 0 3 0)T
Koj-jo BB 2 (0 1 1 0)T 0.152 15 (1 0 3 0)T
Koj-jo BB 5 (0 1 1 0)T 0.138 15 (1 0 3 0)T
Koj-jo BB 10 (0 1 1 0)T 0.158 15 (1 0 3 0)T

Tabla 4.3: Resultados numéricos al variar el parámetro M .

4.2 Búsqueda lineal no monótona 35

Esto en primera instancia nos permite, hipotéticamente, deducir que usar la búsqueda
lineal no monótona con 5 ≤ M ≤ 10 es una buena opción cuando el Algoritmo 3.1 no
converja, sin embargo, en caso que este algoritmo converja, en general lo hará un poco
más rápido que si se usa la búsqueda lineal no monótona, ası́ el número de iteraciones
realizadas por los dos algoritmos sea el mismo.

Similarmente, estudiamos el comportamiento del Algoritmo 4.2 al variar el número

de iteraciones iniciales para las cuales realiza una búsqueda lineal monótona, es decir,
el desempeño del algoritmo al variar el parámetro s . En la Tabla 4.4, reportamos los
resultados más relevantes de nuestras pruebas.

Problema Algoritmo s x0 t (seg) k x∗

Koj-Jo BB 1 (0 0 0 0)T 1.525 117 (1 0 3 0)T
Koj-Jo BB 2 (0 0 0 0)T 0.713 57 (1 0 3 0)T
Koj-Jo BB 3 (0 0 0 0)T 0.248 23 (1 0 3 0)T
Koj-Jo BB 5 (0 0 0 0)T 0..264 25 (1 0 3 0)T
√ T
6 1
Koj-Shi BB 1 (0 0 0 0)T 0.531 44 2
0 0 2
√ T
6 1
Koj-Shi BB 2 (0 0 0 0)T 0.255 22 2
0 0 2
√ T
6 1
Koj-Shi BB 3 (0 0 0 0)T 0.198 19 2
0 0 2
√ T
6 1
Koj-Shi BB 5 (0 0 0 0)T 0.186 19 2
0 0 2
Mathiesen SN 1 (100 100 100 100)T - - -
Mathiesen SN 2 (100 100 100 100)T 0.087 9 (2.322 0 0 0)T
Mathiesen SN 3 (100 100 100 100)T 0.070 9 (2.322 0 0 0)T
Mathiesen SN 5 (100 100 100 100)T 0.065 9 (2.322 0 0 0)T

Tabla 4.4: Desempeño del Algoritmo 4.2 con M = 10.

Estos resultados nos permiten ratificar las conclusiones realizadas por Grippo, Lampa-
riello y Lucidi en [20], pues los mejores resultados se obtienen para s entre 3 y 5, es
decir, el Algoritmo 4.2 presenta un mejor desempeño si se realizan entre tres y cinco
iteraciones iniciales usando una búsqueda lineal monótona.
4.3. Variación del parámetro λ 36

4.3. Variación del parámetro λ

Teniendo en cuenta el buen desempeño local que presentan los algoritmos de complemen-
tariedad que hacen uso de la función mı́nimo y el buen desempeño global que presentan
los algoritmos de complementariedad que hacen uso de la función de Fischer, realiza-
mos una exploración numérica al variar el parámetro inicial λ en el Algoritmo 3.1.
Recordemos que la función de complementariedad dada en (1.3) es idéntica a la función
de Fischer si su parámetro λ es igual a dos y converge a un múltiplo de la función
mı́nimo si λ converge a cero, por lo tanto, hemos decidido tomar valores aleatorios para
el parámetro inicial λ iniciando en uno y haciendo incrementos de 0.1 unidades hasta
tres.

En la siguiente tabla mostramos los resultados obtenidos para algunos de los problemas
analizados anteriormente.

Problema Algoritmo λ x0 t (seg) k x∗

Koj-jo BB 1 (0 0 0 0)T 0.232 23 (1 0 3 0)T
Koj-jo BB 1.1 (0 0 0 0)T - - -
Koj-jo BB 1.2 (0 0 0 0)T 0.314 29 (1 0 3 0)T
Koj-jo BB 1.3 (0 0 0 0)T 0.192 26 (1 0 3 0)T
Koj-jo BB 1.4 (0 0 0 0)T - - -
Koj-jo BB 1.5 (0 0 0 0)T 0.265 32 (1 0 3 0)T
Koj-jo BB 1.6 (0 0 0 0)T 0.150 23 (1 0 3 0)T
Koj-jo BB 1.7 (0 0 0 0)T - - -
Koj-jo BB 1.8 (0 0 0 0)T - - -
Koj-jo BB 1.9 (0 0 0 0)T - - -
Koj-jo BB 2 (0 0 0 0)T 0.141 21 (1 0 3 0)T
Koj-jo BB 2.1 (0 0 0 0)T 0.140 21 (1 0 3 0)T
Koj-jo BB 2.2 (0 0 0 0)T 0.156 23 (1 0 3 0)T
Koj-jo BB 2.3 (0 0 0 0)T 0.185 22 (1 0 3 0)T
Koj-jo BB 2.4 (0 0 0 0)T 0.129 18 (1 0 3 0)T
Koj-jo BB 2.5 (0 0 0 0)T 0.139 21 (1 0 3 0)T
Koj-jo BB 2.6 (0 0 0 0)T 0.152 21 (1 0 3 0)T
Koj-jo BB 2.7 (0 0 0 0)T 0.137 19 (1 0 3 0)T
Koj-jo BB 2.8 (0 0 0 0)T 0.198 29 (1 0 3 0)T
Koj-jo BB 2.9 (0 0 0 0)T 0.162 24 (1 0 3 0)T
Koj-jo BB 3 (0 0 0 0)T - - -
Mathiesen SN 1 (550 − 250 500 − 250)T - - -
4.3 Variación del parámetro λ 37

Problema Algoritmo λ x0 t (seg) k x∗

Mathiesen SN 1.1 (550 − 250 500 − 250)T - - -
Mathiesen SN 1.2 (550 − 250 500 − 250)T 0.088 8 (2.6727 0 0 0)T
Mathiesen SN 1.3 (550 − 250 500 − 250)T 0.081 8 (2.2882 0 0 0)T
Mathiesen SN 1.4 (550 − 250 500 − 250)T 0.087 8 (1.9697 0 0 0)T
Mathiesen SN 1.5 (550 − 250 500 − 250)T - - -
Mathiesen SN 1.6 (550 − 250 500 − 250)T - - -
Mathiesen SN 1.7 (550 − 250 500 − 250)T - - -
Mathiesen SN 1.8 (550 − 250 500 − 250)T - - -
Mathiesen SN 1.9 (550 − 250 500 − 250)T - - -
Mathiesen SN 2 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.1 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.2 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.3 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.4 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.5 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.6 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.7 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.8 (550 − 250 500 − 250)T - - -
Mathiesen SN 2.9 (550 − 250 500 − 250)T - - -
Mathiesen SN 3 (550 − 250 500 − 250)T - - -
Mathiesen BB 1 (1 1 1 1)T 0.187 11 (0.0089 0 0 0)T
Mathiesen BB 1.1 (1 1 1 1)T 0.087 11 (0.0107 0 0 0)T
Mathiesen BB 1.2 (1 1 1 1)T 0.081 11 (0.0128 0 0 0)T
Mathiesen BB 1.3 (1 1 1 1)T 0.109 11 (0.0149 0 0 0)T
Mathiesen BB 1.4 (1 1 1 1)T 0.079 11 (0.0170 0 0 0)T
Mathiesen BB 1.5 (1 1 1 1)T 0.084 11 (0.0191 0 0 0)T
Mathiesen BB 1.6 (1 1 1 1)T 0.078 11 (0.0213 0 0 0)T
Mathiesen BB 1.7 (1 1 1 1)T 0.091 11 (0.0235 0 0 0)T
Mathiesen BB 1.8 (1 1 1 1)T 0.095 11 (0.0258 0 0 0)T
Mathiesen BB 1.9 (1 1 1 1)T 0.094 11 (0.0280 0 0 0)T
Mathiesen BB 2 (1 1 1 1)T 0.178 11 (0.0302 0 0 0)T
Mathiesen BB 2.1 (1 1 1 1)T 0.085 11 (0.0346 0 0 0)T
Mathiesen BB 2.2 (1 1 1 1)T 0.167 11 (0.0346 0 0 0)T
Mathiesen BB 2.3 (1 1 1 1)T 0.064 7 (0.0421 0 0 0)T
Mathiesen BB 2.4 (1 1 1 1)T 0.139 11 (0.0389 0 0 0)T
Mathiesen BB 2.5 (1 1 1 1)T 0.079 11 (0.0410 0 0 0)T
Mathiesen BB 2.6 (1 1 1 1)T 0.169 11 (0.0430 0 0 0)T
4.3 Variación del parámetro λ 38

Problema Algoritmo λ x0 t (seg) k x∗

Mathiesen BB 2.7 (1 1 1 1)T 0.086 11 (0.0449 0 0 0)T
Mathiesen BB 2.8 (1 1 1 1)T 0.099 11 (0.0468 0 0 0)T
Mathiesen BB 2.9 (1 1 1 1)T 0.144 11 (0.0486 0 0 0)T
Mathiesen BB 3 (1 1 1 1)T 0.089 11 (0.0502 0 0 0)T

Tabla 4.5: Desempeño del Algoritmo 3.1 al variar λ.

Estos resultados sugieren una elección del parámetro λ entre uno y dos, pues en este
rango convergió el algoritmo con mayor frecuencia, aunque en general no lo hizo más
rápido que para valores de λ mayores que dos. Un dato que ratifica nuestras observa-
ciones y que no está condensado en la Tabla 4.5 es el comportamiento del Algoritmo
3.1 al resolver problema de Kojima Shindo cuando se tomó como punto inicial el punto
(100 100 100 100), ya que bajo estas condiciones, convergió solo para λ = 1.1, para los
demás valores el algoritmo divergió.

Cabe destacar que la variación del parámetro λ es sólo una experiencia numérica y
por tanto, deja las puertas abiertas para futuras investigaciones en algoritmos que usen
funciones de complementariedad.
 CAPÍTULO 5

COMENTARIOS FINALES

Resolver el PCNL a través de la minimización de la función de mérito (1.10) es impor-

tante en el sentido que da robustés a los algoritmos cuyo objetivo es resolver el sistema
de ecuaciones no lineales (1.5), los cuales por si solos ya son bastante significativos en la
solución de problemas de este tipo.

En este trabajo, propusimos un algoritmo cuasi Newton global que permite resolver el
PCNL y que puede ser útil, en primer lugar, cuando la matriz jacobiana de la función que
define el problema es costosa de calcular y en segundo lugar, cuando el usuario no tiene
información de la posible solución del problema. Para dicho algoritmo desarrollamos la
teorı́a de convergencia global y demostramos que después de algunas iteraciones, este
tiene una tasa de convergencia idéntica a la del algoritmo local propuesto en [1].

En un algoritmo de minimización, el gradiente de la función a minimizar (la función de

mérito, en nuestro caso) desempeña un papel fundamental. En nuestro problema, dicho
gradiente depende de las matrices del jacobiano generalizado de Φλ , las cuales, siguiendo
la filosofı́a cuasi Newton, no estábamos interesados en calcular, por ello, propusimos una
aproximación al gradiente y analizamos, qué tan buena es esta aproximación.

Algunos experimentos numéricos preliminares realizados en esta investigación nos permi-

39
5. Comentarios Finales 40

tieron ver la efectividad y competitividad del algoritmo propuesto, cuando comparamos

su rendimiento con el de un algoritmo tipo Newton generalizado.

Además, las variaciones realizadas a nuestro algoritmo nos permitieron formular conclu-
siones a priori sobre el desempeño del mismo, por lo tanto, creemos conveniente realizar
más pruebas numéricas que consideren dichas variaciones. De igual forma, creemos ne-
cesario incorporar estrategias que permitan elegir el parámetro inicial λ con el objetivo
de mejorar la efectividad en cuanto a convergencia del algoritmo.
 BIBLIOGRAFÍA

[1] Favian E. Arenas, Métodos secantes de cambio mı́nimo para el problema de comple-
mentariedad no lineal, Tesis de Maestrı́a, Universidad del Cauca, 2013.

[2] Dimitri Bertzekas, Constrained optimization and lagrange multiplier methods, At-
hena Scientific, 1996.

[3] Stephen C. Billups, Algorithms for complementarity problems and Generalized Equa-
tions, Ph.D. thesis, University of Wisconsin, 1995.

[4] Peter N. Brown and Youcef Saad, Convergence theory of nonlinear newton-krylov
algorithms, SIAM Journal on Optimization 4 (1994), no. 2, 297–330.

[5] Charles G. Broyden, A Class of Methods for Solving Nonlinear Simultaneous Equa-
tions, Mathematics of Computation 19 (1965), no. 92, 577–593.

[6] Sandra Buhmiler and Nataža Krejić, A new smoothing quasi-newton method for
nonlinear complementarity problems, Journal of Computational and Applied Mat-
hematics 211 (2008), no. 2, 141 – 155.

[7] Frank H. Clarke, Necessary Conditions for Nonsmooth Problems in Optimal Control
and the Calculus of Variations, Ph.D. thesis, University of Washington, 1973.

[8] Frank H. Clarke, Generalized gradients and applications, Transactions of the Ame-
rican Society 205 (1975), 247–262.

[9] Frank H. Clarke, Optimization and nonsmooth analisis, SIAM, 1990.

41
BIBLIOGRAFÍA 42

[10] Thomas F. Coleman, Burton S. Garbow, and Jorge J. More, Software for estimating
sparse jacobian matrices, ACM Trans. Math. Softw. 10 (1984), no. 3, 329–345.

[11] Tecla De luca, Francisco Facchinei, and Kanzow Christian, A theorical and nu-
merical comparison of some semismooth algorithms for complementarity problems,
Mathematical Programming (1997), 1 – 34.

[12] Tecla De Luca, Francisco Facchinei, and Christian Kanzow, A semismooth equa-
tion approach to the solution of nonlinear complementarity problems, Mathematical
Programming 75 (1996), 407–439.

[13] John E. Dennis and Robert B. Schnabel, Numerical methods for unconstrained opti-
mization and nonlinear equations, Society for Industrial and Applied Mathematics,
1996.

[14] Jundi Ding and Hongyou Yin, A new homotopy method for solving non-linear com-
plementarity problems, Numerical Mathematics 16 (2007), 155–163.

[15] F. Facchinei and J. Soarez, A new merit function for nonlinear complementarity
problems and a related algorithm, SIAM 7 (1997), 225–247.

[16] Francisco Facchinei, Minimization of SC 1 functions and the maratos effect, Opera-
tions research letters 23 (1994), no. 17, 131–137.

[17] Michael C. Ferris and Jong-Shi Pang, Engineering and economic applications of
complementarity problems, SIAM Review 39 (1997), 669–713.

[18] Andreas Fischer, A special newton-type optimization method, Optimization (1992),

269 – 284.

[19] Andreas Fischer and Christian Kanzow, On finite termination of an iterative method
for linear complementarity problems, Math. Program. 74 (1996), no. 3, 279–292.

[20] L. Grippo, F. Lampariello, and S Lucidi, A nonmonotone line search technique for
newton’s method, SIAM 23.

[21] Christian Kanzow and Helmut Kleinmichel, A new class of semismooth newton-type
methods for nonlinear complementarity problems, Comput. Optim. Appl. 11 (1998),
no. 3, 227–251.

[22] Masakazu Kojima and Susumu Shindoh, Extensions of newton and quasi-newton
methods to systems of P C 1 equations, Research reports on information sciences,
Inst. of Technology, Department of Information Sciences, 1986.
BIBLIOGRAFÍA 43

[23] Michael Kostreva, Elasto-hydrodynamic lubrication: A non-linear complementarity

problem, International Journal for Numerical Methods in Fluids 4 (1984), no. 4,
377–397.

[24] Li-Zhi Liao and Houduo Qi, A neural network for the linear complementarity pro-
blem, Math. Comput. Model. 29 (1999), 9 – 18.

[25] Li-Zhi Liao, Houduo Qi, and Liqun Qi, Solving nonlinear complementarity problems
with neural networks, Comput. Appl. Math. 131 (2001), 343 – 359.

[26] Lars Mathiesen, An algorithm based on a sequence of linear complementarity pro-

blems applied to a walrasian equilibrium model: An example, Mathematical Pro-
gramming 37 (1987), 1–18.

[27] Jorge Nocedal and Wright Stephen J., Numerical optimization, Springer Series in
Operations Research, 2000.

[28] Jong-Shi Pang and Steven A. Gabriel, A robust algorithm for the nonlinear comple-
mentarity problem, Mathematical Programming 60 (1993), 295–337.

[29] Jong-Shi Pang and Liqun Qi, Nonsmooth equations: Motivation and algorithms,
SIAM Journal on Optimization 3 (1993), no. 3, 443–465.

[30] Rosana Pérez and Vera Lucı́a Rocha, Recent applications and numerical implementa-
tion of quasi-newton methods for solving nonlinear systems of equations, Numerical
Algorithms 35 (2004), 261–285.

[31] Rosana Pérez M., Algoritmos para complementariedade não lineal e problemas re-
lacionados, Ph.D. thesis, IMECC, UNICAMP, Campinas, Brazil, 1997.

[32] Rosana Pérez M. and Diaz V. Tomás H., Minimización sin restricciones, Universi-
dad del Cauca, 2010.

[33] Liqun Qi, Convergence analysis of some algorithms for solving nonsmooth equations,
Math. Oper. Res. 18 (1993), no. 1, 227–244.

[34] Vera Lucia Rocha, Jose M. Martı́nez, and Rosana Pérez, On the local convergence
of quasi-newton methods for nonlinear complementary problems, Applied Numerical
Mathematics 30 (1999), 3–22.

[35] Jhon Glen Wardrop, Some theoretical aspects of road traffic research, Proceeding of
the Institute of Civil Engineers (1952), 325 – 378.

[36] David Watkins, Fundamentals of matrix computations, John Wiley and Sons, 2002.
BIBLIOGRAFÍA 44

[37] Qing Xu and Chuangyin Dang, A new homotopy method for solving non-linear
complementarity problems, Optimization 57 (2008), no. 5, 681–689.

[38] Longquan Yong, Nonlinear complementarity problem and solution methods, Procee-
dings of the 2010 International Conference on Artificial Intelligence and Compu-
tational Intelligence: Part I, Springer-Verlag, 2010, pp. 461–469.

[39] Xia Youshen and Feng Gang, A neural network for solving nonlinear projection
equations, Neural Networks 5 (2007), 577 – 589.