El método de gradiente conjugado

proyectado precondicionado (PPGC) en

optimización con restricciones
* ** ***
Leobardo Valera Miguel Argáez Brígida Molina
6 de diciembre de 2011

Resumen
En este trabajo se propone utilizar el algoritmo gradiente conju-
gado proyectado precondidionado (PPCG) para resolver un problema
de minimización cuadrática con restricciones de igualdad. Para acele-
rar la velocidad de convergencia del método de gradiente conjugado
proyectado (PCG), se plantea precondicionar el sistema KKT usando
la factorización incompleta de Cholesky del Hessiano de la función
objetivo. Finalmente se muestran algunos resultados numéricos que
se obtienen al resolver varios problemas de la librería CUTEr, usan-
do el algoritmo PCG y su versión precondicionada PPCG, los cuales
muestran la efectividad del algoritmo.

1. Introducción
Problemas de optimización lineal con restricciones de igualdad surgen
en muchas y variadas aplicaciones entre otras en economía, en la solución de
*
Escuela de Matemáticas, Universidad Metropolitana, [email protected]
**
Department of Mathematical Sciences University of Texas at El Paso,
[email protected]
***
Centro de Cálculo Cienttífico y Tecnológico (CCCT), Escuela de Computación, Facul-
tad de Ciencias, Universidad Central de Venezuela. [email protected] Facultad
de Ciencias, Universidad Central de Venezuela. [email protected]

1
ecuaciones diferenciales parciales para materiales incomprensibles, en la re-
cuperación de señales de entrada. Un problema de optimización cuadrática
con restricciones de igualdad tiene la forma:
1 T
minimizar x Gx + gT x (1)
2
sujeto a Ax = b

donde G ∈ Rn×n es una matriz simétrica, A ∈ Rm×n con m < n, es una

matriz de rango completo , g ∈ Rn , b ∈ Rm son vectores dados, x ∈ Rn es el
vector a determinar. Las condiciones necesarias y suficientes para resolver
(1) son A sea de rango completo y G sea definida positiva sobre el núcleo de
A. Resolver (1) es equivalente a resolver el siguiente sistema de ecuaciones
lineales, llamado las condiciones Karush-Kuhn-Tucker (KKT)
! ! !
G AT x −g
= , (2)
A 0 y b
donde y ∈ Rm es el multiplicador de Lagrange [12] .
Este tipo de sistemas aparecen en muchas investigaciones en ciencia e
ingeniería. Resolver este problema eficientemente es la línea de investiga-
ción actual [5, 6, 7, 13, 14, 15, 16]. Benzi et al. en [5] presenta un estudio
exhaustivo de este problema.
Los métodos para resolver (2) se pueden clasificar en dos grupos: El
primero usa un método iterativo sobre todo el sistema como por ejemplo:
Gradiente Conjugado (CG), Mínimo Residuo (MINRES), entre otros. El
segundo reduce el sistema de una manera equivalente a un sistema de menor
dimensión, entre ellos está el Shur Complement, y el Gradiente Conju-
gado Proyectado, PCG [12] entre otros.
Uno de los propositos de este artículo es divulgar los trabajos de inves-
tigación del Profesor Argaez en el área PCG, y presentar nuevas líneas de
investigación donde se resuelva el sistema eficientemente sin la necesidad
del almacenarlo [1, 2, 3, 4]. La idea fundamental del PCG es conseguir una
solución particular para las restricciones lineales y luego buscar la solución
del problema moviendose en una variedad afín determinada por la solución
particular y el subespacio generado por el núcleo de A.

2
 La solución particular es determinada por la solución de norma mínima,
la cual esta en el espacio fila de la matriz A. Luego el problema se reduce a
encontrar la solución de un sistema lineal positivo definida sobre el núcleo
de A, este sistema lineal es llamada ecuación proyectada. Para resolver la
ecuación proyectada aplicamos el algoritmo CG el cual nos garantiza la
solución en n − m iteraciones trabajando en aritmética infinita.
El algoritmo tiene la propiedad que cada iteración satisface las restric-
ciones de igualdad del problema (1), la sucesión de aproximaciones a la
solución, y la sucesión de residuos son estrictamente crecientes y estric-
tamente decrecientes en norma euclideana, respectivamente. Estas propie-
dades serían las propiedades que tendríamos en el caso que el sistema (2)
fuera difinido positivo, desafortunadamente el sistema no lo es, por lo tanto
el uso del CG no garantiza una solución del problema.
En [4] se demuestra que un buen precondicionador para la información
de segundo orden del problema produce un buen precondicionador para el
algoritmo PCG, esto es, basta con encontrar una buena aproximación de
G.
En este trabajo proponemos seleccionar, para resolver el problema (2), al
método PPCG usando como precondicionador a la factorización incompleta
de Cholesky del Hessiano de la función objetivo.
El artículo esta organizado de la siguiente manera: en la sección 2 se
hace una breve revisión de como encontrar la solución de un problema de
optimización cuadrática con restricciones lineales de forma directa; poste-
riormente, en esa misma sección, se plantea el mismo problema pero ahora
en función de un operador proyección sobre el espacio nulo de A para
finalizar la sección con la descripción del método gradiente conjugado pro-
yectado y se muestra el algoritmo (PCG) [4]. En la sección 3 se presenta el
algoritmo gradiente conjugado proyectado precondicionado. En la sección
4 se presentan los resultados numéricos al resolver problemas de la librería
CUTEr usando el algoritmo PCG es sus versiones sin precondicionar y pre-
condicionado. En la sección 5 se proponen algunas líneas de investigación
que se pueden estudiar. Las conclusiones se hacen en la sección 6.

3
2. Gradiente Conjugado Proyectado (PCG)
El objetivo de esta sección es presentar el método Gradiente Conju-
gado Proyectado para resolver problemas de minimización cuadrática con
restricciones de igualdad propuesto en [4].
Consideremos el problema (2). Si A es de rango completo, entonces las
columnas de AT son linealmente independientes y generan un subespacio
de dimensión m de Rn . Si consideramos ahora Z, una base del espacio nulo
de A, dicho subespacio tendrá dimensión n − m y [AT Z] forma una base de
Rn .
Supongamos ahora que x es la solución de (1), entonces,

x = AT x∗ + Zw. (3)

De la ecuación (3) podemos concluir que x es la suma de dos vectores;

además, multiplicando a la expresión anterior por A, se tiene que:

b = Ax = A(AT x∗ ) + A(Zw) = A(AT x∗ ) + (AZ)w = A(AT )x∗ (4)

es decir, xp = AT x∗ es una solución particular de Ax = b y xh = Zw ∈

N(A), donde N(A) es el espacio nulo de A. Así, la solución x puede expre-
sarse como:

x = xp + xh . (5)
Por otra parte, como A es de rango completo con m < n, se tiene que
AAT es simétrica definida positiva, por lo tanto, es invertible y así de la
expresión (4) se obtiene que:

x∗ = (AAT )−1 b
xp = AT (AAT )−1 b
xh = (I − AT (AAT )−1 A)w para w ∈ Rn .

Observe que P = I − AT (AAT )−1 A, es una proyección ortogonal sobre el

espacio nulo de A.

4
 Una vez que se halla xp , se puede resolver el sistema aumentado (2)
en términos de xh e y. En efecto, sustituyendo (5), junto con las expre-
siones obtenidas para xp y xh en (2) y realizando algunas manipulaciones
algebraícas, el sistema aumentado se reduce a la siguiente ecuación:

Gxh + AT y = −(g + Gxp ). (6)

Resolviendo para y, se obtiene que:

y = −(AAT )−1 Ag − (AAT )−1 AGxp − (AAT )−1 AGxh . (7)

Luego, sustituyendo (7) en (6) se deduce la siguiente ecuación proyectada:

PGxh = −P(Gxp + g)

donde P = I − AT (AAT )−1 A. Finalmente, usando P en la expresión anterior,

ésta puede ser reescrita como:

PGPw = −P(Gxp + g) (8)

donde Pw = xh para algun w ∈ Rn .

En [4] se demuestra que la ecuación proyectada (8) tiene una única
solución en el espacio nulo de A y ésta es la solución de norma mínima;
además, como G es simétrica definida positiva sobre el espacio nulo de A, se
propone obtener dicha solución a usando el método de gradiente conjugado
[11] para el sistema (8) con iterado inicial cualquier vector en el espacio
nulo de A. Una de las ventajas de este procedimiento es que las soluciones
aproximadas, las direcciones y los vectores residuales permanecen en el
espacio nulo de A. Más aún, en aritmética exacta, el algoritmo encuentra
la solución en a lo sumo n − m iteraciones puesto que PGP tiene máximo
n − m autovalores diferentes de cero.
El algoritmo Gradiente Conjugado Proyectado (PCG) se presenta
en Algoritmo 1.

3. Precondicionamiento
Un procedimiento para acelerar la velocidad de convergencia del méto-
do iterativo PCG es mejorar la distribución de los autovalores de la matriz

5
Algoritmo 1 Algoritmo Gradiente Conjugado Proyectado. PCG
Entrada: G, A, b, g, ε, maxiter
Salida: x, y
1: x = AT (AAT )−1 b
2: f = Gx + g
3: y = −(AAT )−1 Af
4: r = −(f + AT y)
5: si ||r|| < ε entonces
6: devolver x, y
7: alto, terminar programa
8: fin si
9: d = r
10: βn = rT r
11: para i = 1 hasta maxiter hacer
12: αd = dT Gd
13: x = x + βαnd d
14: f = Gx + g
15: y = −(AAT )−1 Af
16: r = −(f + AT y)
17: si ||r|| < ε entonces
18: devolver x, y
19: alto, terminar programa
20: fin si
21: βd = βn
22: β n = rT r
23: d = r + ββnd d
24: fin para
25: devolver x, y

PGP o reducir su número de condición [8, 9]. Además, los precondiciona-

dores deben ser seleccionados de tal manera que se explote la estructura y
la esparcidad de las matrices involucradas en el problema.
Argaéz [4] demostró que PMP es un buen precondicionador para el siste-
ma (8) si M es una buena aproximación para la matriz G. La demostración

6
se basó en aproximar a la inversa Moore-Penrose de la matriz PGP del sis-
tema (8) por la inversa de Moore-Penrose para la matriz PMP. Definiendo
N = M − G, se demuestra que mientras mayor reducción se logre en kNk
se mejora la calidad del precondicionador PMP para el sistema (8) y esto
también implica una mejora de calidad de M como un preconditioner para
G; por lo tanto, las propiedades y estructura de G son importantes para la
elección de precondicionadores adecuados para el algoritmo de PCG.
El algoritmo Gradiente Conjugado Proyectado Precondicionado (PPCG)
se presenta en Algoritmo 2.
En [4], se proponen como aproximaciones para G las siguientes:
1. M = D + τ, donde D es la matriz diagonal de G y τ es un escalar
dado.
1
2. M = ω(2−ω) (D + ωL)D−1 (D + ωLT ), donde D y L son la matriz dia-
gonal y la parte triangular estricta de G, respectivamente. El escalar
ω es el paraámetro de sobrerelajación del método de SOR y debe ser
seleccionado de tal manera que 0 < ω < 2.
Aún cuando los resultados obtenidos con estos dos precondicionadores
en [4] fueron bastante satisfactorios, ellos son bastante básicos. En este
trabajo proponemos utilizar como precondicionador para G a la matriz
M = L̃L̃T , donde L̃ es el factor de Cholesky incompleto para la matriz G
con el mismo patrón de esparcidad que la parte inferior de G.

4. Resultados Numericos
Los experimentos numéricos fueron implementados en una computadora
Intel(R) Core(TM) i3-2330M CPU @ 2.20GHz 2.20 GHz, con 2.70 GB de
memoria Ram. Los algoritmos fueron implementados en MatLab R2008b.
Se usó como criterio de parada |Gx + AT y + g| < 1e − 6 y el tiempo es
medido en segundos. Los problemas seleccionados para las pruebas fueron
obtenidos de la colección CUTEr [10].
En la tabla (1) se comparan los tiempos de convergencia y el número de
iteraciones para algunos problemas de la colección CUTEr. El precondicio-
nador utilizado para el algoritmo (2) fue M = L̃L̃T , donde L̃ es el factor de

7
Algoritmo 2 Algoritmo Gradiente Conjugado Proyectado Precondiciona-
do. PPCG
Entrada: G, A, b, g, ε, maxiter, M
Salida: x, y
1: x = AT (AAT )−1 b
2: f = Gx + g
3: y = −(AAT )−1 Af
4: r = −(f + AT y)
5: si ||r|| < ε entonces
6: devolver x, y
7: alto, terminar programa
8: fin si
9: Resolver para w, Mw = r
10: rprec = Pw
11: d = rprec
12: βn = rT rprec
13: para i = 1 hasta maxiter hacer
14: αd = dT Gd
15: x = x + βαnd d
16: f = Gx + g
17: y = −(AAT )−1 Af
18: r = −(f + AT y)
19: si ||r|| < ε entonces
20: devolver x, y
21: alto, terminar programa
22: fin si
23: Resolver para w, Mw = r
24: rprec = Pw
25: βd = βn
26: βn = rT rprec
27: d = r + ββnd d
28: fin para
29: devolver x, y

8
 PCG PPCG
Problema m n it t seg. it t seg. ||Gx + AT y + g||
DUAL1 1 85 141 0.044 2 0.004 8.24e-07
DUAL2 1 96 66 0.024 2 0.004 6.55e-07
DUAL3 1 111 56 0.022 2 0.005 7.06e-07
DUAL4 1 75 29 0.011 2 0.003 7.52e-07
GENHS28 4900 5000 100 0.434 3 0.025 2.97e-11
MOSARQP1 700 2500 12 0.013 5 0.010 4.50e-07
MOSARQP2 600 900 17 0.018 6 0.011 7.38e-07
STCQP2 2052 4097 111 0.306 44 0.169 9.39e-07

Cuadro 1: Tamaño de los bloques, número de iteraciones requeridas por

PCG y PPCG para satisfacer el criterio de parada, tiempo de CPU y norma
de los residuales ||Gx + AT y + g|| para los diferentes problemas de prueba

Cholesky incompleto para G con con el mismo patron de esparcidad que la

parte inferior de G. En los resultados mostrados en la tabla (1) se observa
que se reduce de CPU desde un 24 % para el problema MOSARQP1 hasta
un 91 % para el problema DUAL1, cuando se resuelve el sistema usando el
algoritmo PPCG en lugar de PCG.

5. Líneas de Investigación
A pesar de que el PCG tiene una buena rata de convergencia, n − m
iteraciones, tiene el inconveniente de almacenar la matriz asociada a las res-
tricciones de igualdad. En algunas aplicaciones el tamaño de la matriz A es
tan grande que supera la capacidad de almacenamiento del hardware o la
matriz A no se conoce explícitamente, como en el caso de procesamiento de
señales. Por lo tanto proponemos resolver los problemas de mínimos cua-
drados asociados al problema usando métodos iterativos, como por ejemplo,
CG. Como el almacenamiento de A ya no es un obstaculo, podemos usar
(PMP)† como precondicionador con una matriz M los más cercana a G.

9
6. Conclusiones
Basados en los resultados numéricos presentados en la tabla (1) ob-
servamos que precondicionando con una proyección ortogonal usando una
factorización incompleta de Choleski de G, del algoritmo (PPCG) mejora
de manera significativa los resultados obtenidos del algoritmo (PCG).

Referencias
[1] M. Argaez, H. Klie, M. Rame, and R. A. Tapia, An interior point
Krylov-orthogonal projection method for nonlinear program-
ming. Talk presented at the Fifth SIAM Conference on Optimi-
zation, Victoria, Canada, 1996.

[2] M. Argaez, H. Klie, M. Rame, and R. A. Tapia, An interior point

Krylov-orthogonal projection method for nonlinear program-
ming. TR97-16, Computational and Applied Mathematics, Rice
University, Houston Texas, 1997.

[3] M. Argaez, Exact and Inexact Newton Linesearch Interior-

Point Algorithms for Nonlinear Programming Problems.
PHD Thesis, Computational and Applied Mathematics, Rice
University, Houston Texas, 1997.

[4] M. Argáez, A projected conjugate gradient algorithm for KKT

systems, Tech. Rep.TR01-2010. Mathematical Sciences Dept.
University of Texas at El Paso.

[5] M. Benzi, G. H. Golub and J. Liesen, Numerical solution of

saddle point problems, Acta Numerica, Vol. 14, pp. 1-137, 2005.

[6] D. P. Bertsekas, Nonlinear Programming, Athena Scientific,

2th, 1999.

[7] Z. Cao, Fast Uzawa algorithm for generalized saddle point

problems, Appl. Numer. Math., Vol. 46, No. 2, pp. 157-171, 2003.

10
 [8] D. De Cecchis, H. López and B. Molina, Técnicas de Precondi-
cionamiento de los Métodos de Krylov para la Discretización
del Problema Stokes Generalizado, Proceedings VIII Interna-
tional Congress on Numerical Methods in Engineering and Ap-
plied Sciences, pp. TM1-TM7, 2006.

[9] D. De Cecchis, H. López and B. Molina, FGMRES Preconditio-

ning by Symmetric/skew-symmetric Decomposition of Gene-
ralized Stokes Problems, Mathematics and Computers in Simu-
lation, Vol. 79, No. 6, pp. 1862-1877, 2009.

[10] N. I. M. Gould, D. Orban and Ph. L. Toint, CUTEr (and Sif-

Dec), a Constrained and Unconstrained Testing Environ-
ment, Revisited, Tech. Report TR/PA/01/04, CERFACS, Tou-
louse, France, 2001.

[11] M.R. Hestenes and E.L. Stiefel, Methods of conjugate gradients

for solving linear systems, J. Res. Nat. Bur. Stand., Vol. 49, pp.
409-436, 1952.

[12] N.I.M. Gould, M.E. Hribar, and Nocedal. On the solution of

equality constrained quadratic programing problems arising
in optimization. SIAM J. Sci Computing, Vol 23, 4, pp. 1375-
1394, 2001.

[13] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM,

2nd ed., Philadelphia, PA, 2003.

[14] H. Uzawa, Iterative methods for concave programming, Stu-

dies in Linear and Nonlinear Programming, K. J. Arrow and
L. Hurwicz and H. Uzawa (Eds), Stanford University Press, pp.
154-165, Stanford, CA, 1958.

[15] M. H. Wright, Interior point methods for constrained optimi-

zation, Acta Numerica, Vol. 1, Cambridge University Press, pp.
341-407, 1992.

11
[16] S. J. Wright, Primal Dual Interior Point Methods, SIAM, Phi-
ladelphia, PA, 1997. N. I. M. GOULD, D. ORBAN, AND PH. L.
TOINT, CUTEr (and SifDec), a Constrained and Unconstrained
Testing Environment, Revisited, Tech. Report TR/PA/01/04,
CERFACS, Toulouse, France, 2001.

12