Reformulación del VPC (A, B) como un problema de complementariedad no lineal

En general, dada F: ℝ ⁿ⁺¹ - ℝ ⁿ⁺¹, F (y) = (F ₁ (y),...,F _n+1 (y)) ^T , continuamente diferenciable, el Problema de Complementariedad No Lineal consiste en encontrar un vector y e ℝ ⁿ+¹ que satisfaga las tres condiciones siguientes,

Veamos a continuación que, en efecto, el problema VPC2(A,B) puede reformularse como un problema de complementariedad no lineal.

A partir de (3), VPC2(A,B) puede expresarse de la siguiente forma (^{Júdice y cols., 2007}): hallar λ > 0 y x ∈ ℝ ⁿ tal que

Haciendo un cambio de variable, se definen F: ℝⁿ+¹ - ℝⁿ+¹ y e Rⁿ+¹ tales que,

con lo cual, el problema VPC2(A, B), es de la forma (4); es decir, es un problema de complementariedad lineal.

Una estrategia de solución para problemas de complementariedad no lineal, tal vez, la más popular y ampliamente usada, consiste en reformularlo como un sistema de ecuaciones no lineales mediante ciertas funciones especiales, ϕ: ℝ² → ℝ, llamadas funciones de complementariedad (^{Yong, 2010}) que satisfacen la equivalencia, ϕ(a,b) = 0 ⇐⇒ a ≥ 0, b ≥ 0, ab = 0.

En este artículo usaremos la familia de funciones de complementariedad introducida en (^{Kanzow y Kleinmichel, 1998}), y definida por,

la cual ha resultado competitiva en la reformulación de problemas de complementariedad (^{Arenas, Martínez, y Perez, 2014}) (^{Perez, Arenas, Martínez, y Arias, 2019}). Así, la reformulación del problema de valores propios complementarios puede, a su vez, reformularse como el siguiente sistema de ecuaciones no lineales no diferenciable

Este sistema se puede resolver mediante un método tipo Newton generalizado (^{Ferris y Pang, 1997}; ^{Facchinei y Pang, 2003}; ^{Kanzow y Kleinmichel, 1998}) o cuasi-Newton (^{Arenas y cols., 2014}; ^{Arias y cols., 2017}), entre otros (^{Arenas, Martínez, y Perez, 2020}; ^{Sanchez, Perez, y Martínez, 2021}). En el primer caso, la iteración básica, es de la forma

donde Hk ∈ ∂Φ _τ (y _k ), es decir B _k es una matriz en el jacobiano generalizado de Φ _τ en y _k (^{Clarke, 1975}). En el caso, de un método cuasi-Newton, la iteración básica, es de la forma

donde B _k ≈ Hk ∈ ∂Φτ (y _k ) .

Por otra parte, con el fin de globalizar este tipo de algoritmos que resuelven sistemas de ecuaciones no lineales como el planteado en (7), se introduce una función de merito y se resuelve un problema de minimization (^{Dennis y Schnabel, 1996}). La función de mérito natural y ampliamente usada es Ψ _τ definida por Ψ _τ (y) = 12 Φ _τ (y)T Φ _τ (y) , y el problema de minimización es el siguiente:

En (^{Arias y cols., 2017}), los autores proponen un algoritmo global cuasi-Newton, concretamente un algoritmo secante de cambio mínimo, para resolver problemas de complementariedad no lineal, como el problema (8), el cual usa búsqueda lineal como estrategia de globalization. Este algoritmo mostro resultados excelentes al ser probado con funciones de difícil comportamiento numérico. Por ello, consideramos que es la mejor opción para resolver el problema VPC (A, B) vía la reformulation como un problema de complementariedad no lineal. En la experimentation numérica, presentamos más detalles de este algoritmo; en lo que sigue, nos referiremos a el como Algoritmo 1.

Verificación de hipótesis

El Algoritmo 1 converge con tasa superlineal, y esta convergencia esta condicionada a que satisfaga las siguientes hipotesis:

H.1 Existe una solution y* del sistema Φ _T (y) = 0.

H.2 Las funciones F _i son continuamente diferenciables y su gradiente es una funcion semi-suave.

H.3 La matriz jacobiana de F es Lipschitz continua en una vecindad de y*.

A continuation, analizaremos hasta que punto la función F satisface las hipótesis que requiere el Algoritmo 1 para alcanzar la convergencia. Para ello, usamos (5) y deducimos que la matriz Jacobiana de F, esta dada por

Estas expresiones nos facilitaran la verification de las hipótesis que garantizan la viabilidad del Algoritmo 1 para este problema particular.

H.1 En general, para determinar si un problema de complementariedad tiene solución, bastaría verificar que la función que define el problema es estrictamente monótona (^{Ferris y Pang, 1997}); sin embargo, para el problema (4) con F definida en (5), este camino conduce a imponer condiciones muy fuertes sobre A y B para garantizar la monotonía de F y con ello, la solución al problema.

Una alternativa para garantizar la existencia de al menos una solución del problema VPC (A, B) es imponer condiciones menos restrictivas sobre A y B. Esto fue hecho en (^{Queiroz, Judice, y Humes, 2004}), donde los autores demuestran que, si todas las componentes de A son positivas y B = I _n , el problema VPC (A, B) tiene solución. Adicionalmente, si A es simétrica la solución es única. Este resultado, nos permite verificar la hipótesis H1 y diseñar experimentos numéricos, como veremos más adelante.

Vale la pena mencionar que la cantidad maxima de soluciones que puede tener el problema VPC (A, B) es n2^n-1 (^{Seeger, 1999}).

H.2 Las funciones F _i son diferenciables. En efecto, cuando 1 ≤ i ≤ n, su respectivo gradiente está dado por

y, para n+1, el ∇Fi(y) es el vector de unos de n+1 componentes. Cada uno diferenciable, con además para 1 ≤ i ≤ n

y, para i = n +1, la ∇ ² F _i (y) es la matriz nula de orden n +1. Esto demuestra que cada gradiente ∇F _i es continuamente diferenciable; por lo tanto, es también una función semisuave .

H.3 Probemos que F' es Lipschitz continua en todo su dominio. Sean y,z G ℝⁿ⁺¹ de tal forma que

Considérese la diferencia

Caso 1. Si máx

Caso 2. Si máx

por lo tanto, F' es Lipschitz continua con constante ∥A∥₁.

Al satisfacer de manera razonable las hipótesis del algoritmo, se procede á mostrar los resultados numéricos.

Experimentación numérica

En esta sección, resolvemos numericamente la reformulation del problema VPC2(A, B) como un problema de complementariedad no lineal, el cual a su vez, lo hemos reformulado como el problema de minimization (8). Para ello usaremos el algoritmo propuesto en (^{Arias y cols., 2017}), el cual hemos denominado en este documento Algoritmo 1. Este es un algoritmo cuasi-Newton, exactamente secante de cambio mínimo (^{Martínez, 1992}), globalizado mediante una estrategia de busqueda lineal (^{Dennis y Schnabel, 1996}). Para mayor claridad en la lectura de este documento lo incluimos a continuation.

Algoritmo 1 Cuasi Newton 

El Algoritmo 1, hace un llamado al Algoritmo 2, que selecciona la actualización secante, que actualizara la matriz B _k y que será usada en la siguiente iteración. Las actualizaciones que consideraremos son las siguientes: buena de ^{Broyden (Broyden, 1965)}, mala de Broyden (^{Brown y Saad, 1994}), ^{Schubert (Schubert, 1970)}, COLUM (^{Martínez, 2000}) e ICUM (^{Martínez, 1992}).

Algoritmo 2 Función Actualizar 

Se considero como algoritmo para hacer las comparaciones respectivas el método de Newton generalizado, el cual llamaremos Algoritmo 3, y lo incluimos a continuación para claridad en la lectura de la sección. Por supuesto, aprovechamos la estructura del jacobiano de F indicada en (9) ya que, no es muy costosa su evaluación, ni su implementación.

Algoritmo 3 Newton generalizado 

Construcción de los problemas

Para construir los problemas VPC (A, B) usamos el lineamiento descrito en (^{Júdice y cols., 2007}); es decir, dos categorías de problemas: simétricos y no simétricos; en todos, la matriz B es la identidad de tamaño n x n. La matriz A se construye como se indica en la Tabla 1, donde Escala significa el valor más grande que puede tomar cada componente de A, randn(n) es una matriz aleatoria de tamaño n x n y p esta dado por (3), generando un total de ocho configuraciones de parámetros.

Para escribir los códigos de los algoritmos y de las funciones objetivo de cada problema usamos el software Matlab® . Realizamos las pruebas numéricas en un computador Intel (R) Core (TM) i7-9700F CPU de 3.00 GHz.

Experimento preliminar

Uno de los interrogantes planteados antes de comenzar la experimentation fue el siguiente, ¿cual es el valor de x adecuado para la función ϕ _τ definida en (6)? En (^{Kanzow y Klein-michel, 1998}) se presento un algoritmo dinámico para variar x con buenos resultados, sin embargo, en el caso particular del problema VPC (A, B), en pruebas preliminares no se obtuvo un decrecimiento eficaz, al compararlo con el valor de x = 2 estático, el cual es un valor muy utilizado en funciones de complementariedad, esto llevo a diseñar el siguiente experimento.

Encontrar en cada paso el valor de x en el cual se pueda obtener el mayor decrecimiento de la función objetivo Ψ _t . Para ello, se utilizaron diferentes valores de t desde 0.01 hasta 3.99 y se seleccionó para cada paso aquel valor de t en que se obtuvo el máximo descenso de Ψ _t utilizando el método de Newton generalizado. Luego se repitió esta experiencia 10000 veces, los promedios de los 9110 ejecuciones exitosas se muestran en la Figura 1.

Figura 1 Valores de t que producen mayor decrecimiento y su respectivo Ψ _t (x _k ) en las primeras 7 iteraciones. 

Podemos observar que el mayor decrecimiento se obtiene en las primeras 3 iteraciones, en los cuales λ toma valores muy pequeños, es por ello que hemos decidido realizar todos los experimentos con un valor de λ = 10^-2.

Experimento principal

Realizamos para cada valor de n, cada valor de Escala y cada valor de p, 1000 ejecuciones; en cada ejecución se generó aleatoriamente un problema. El punto inicial y ₀ también se genero aleatoriamente con el comando randn(n + 1,1). Cada uno de estos problemas, se resolvió con los algoritmos 1 y 2 respectivamente.

Los resultados más relevantes obtenidos en esta experimentación, se presentan en las Tablas 2 a 4, donde n, es el tamaño de las matrices A y B, p, es el valor de la constante dada en (3); escala, es el valor mas grande que puede tener cada componente de la matriz A; en cada una de las actualizaciones E, indica el porcentaje de éxitos con cada una de las actualizaciones; es decir, el número de veces que se obtuvo convergencia entre 1000 expresado como porcentaje; Iter, indica el promedio de iteraciones redondeado a un entero de las ejecuciones que tuvieron éxito; finalmente, tiempo indica el promedio de tiempo CPU redondeado a cuatro decimales en las ejecuciones que tuvieron éxito.

En los resultados presentados en la Tabla 2 se destacan las actualizaciones Broyden M.e ICUM por su rapidez de convergencia, sin embargo, tienen poco éxito, por lo cual afirmamos que no son muy efectivas al tomar valores iniciales aleatorios, si el problema tuviese una forma de iniciar el algoritmo en una vecindad pequeña de la solución esta falencia pudiera ser superada. Por otro lado, los métodos de Newton generalizado y Broyden B son los métodos más eficaces para este caso no simétrico con valores pequenos de las entradas de la matriz A.

Los resultados de la Tabla 3 muestran que el método más eficiente en el caso no simétrico con entradas menores o iguales que 50, sigue siendo el método de Newton generalizado, seguido por la actualización de Broyden B también la actualización ICUM es la que converge mas rapidamente, aunque tiene poco exito.

De la Tabla 4 deducimos que para el caso simétrico, el método de Newton generalizado es el más indicado, por efectividad y eficiencia, en contraste, el método de Shubert no parece ser recomendable para problemas de tamaño mayor a 5. También llama la atención que en muchos casos no se alcanza la convergencia en algunos métodos cuasi-Newton, salvo el método Broyden Bueno, que es muy confiable.