Implementaciones secuenciales y paralelas del algoritmo Hessenberg-Schur
- Resumen
- Descripción
teórica del algoritmo - Algoritmo Hessenberg-Schur para
la resolución de la ecuación de
Sylvester - Bibliotecas de
Software - Implementación de la
resolución de la ecuación de
Sylvester - Algoritmo secuencial
utilizando LAPACK - Algoritmo paralelo
utilizando ScaLAPACK - Validación de las
soluciones obtenidas - Bibliotecas de
Software - Análisis de los
resultados experimentales - Prestaciones del algoritmo
secuencial - Prestaciones del algoritmo
paralelo - Conclusiones
- Bibliografía
La ecuación de Silvestre 
tiene numerosas aplicaciones en teoría
de control,
procesamiento de señales, filtrado, reducción de
modelos,
restauración de imágenes,
etc. Para la resolución de dicha ecuación existen
múltiples algoritmos
tanto por métodos
directos como iterativos. En este trabajo se
aborda el algoritmo de
Hessenberg-Schur y su implementación en forma secuencial y
paralela, haciéndose un análisis comparativo de los resultados
obtenidos de la versión paralela y las ventajas de dicha
variante sobre la secuencial dado el elevado costo
computacional de dicho algoritmo.
El presente trabajo tiene como objetivo
presentar el algoritmo de Hessenberg-Schur para la
resolución de la ecuación de Silvestre la cual
tiene la siguiente expresión:
(1)
donde
,
y 
.
La ecuación (1) tiene solución
única si y solo si 
[LT85], donde 
(Z) denota el espectro de la matriz
Z
Esta ecuación tiene numerosas aplicaciones en
teoría de control, procesamiento de señales,
filtrado, reducción de modelos, restauración de
imágenes, implementación de métodos
numéricos implícitos para ecuaciones
diferenciales ordinarias y la diagonalización por
bloques de matrices.
Ejemplos de ello pueden encontrarse en: [AL98], [CR96],
[GVL96].
Para la resolución de la ecuación de
Silvestre existen múltiples algoritmos. En problemas de
pequeña dimensión los métodos más
empleados son los de Bartels-Stewart [BS72] y el de
Hessenberg-Schur [GNVL79]..
Para problemas de gran dimensión varios
algoritmos iterativos han sido propuestos para la solución
aproximada de la ecuación de Sylvester ver [Jbi98]. Estos
métodos están basados en el subespacio de Krylov,
así tendremos por ejemplos métodos que se basan en
el algoritmo de Arnoldi [Sad93] o en el de Lanzos [JR], o
utilizando el algoritmo GMRES.
También podemos mencionar como métodos
iterativos para la resolución de la ecuación de
Sylvester el de Newton, el de
Newton-Schulz y el de Halley´s [BQOQO04].
En este trabajo nos centraremos en el algoritmo de
Hessenberg-Schur el cual se describe en la siguiente
sección.
Descripción
teórica del algoritmo
El algoritmo de Hessenberg-Schur para la
resolución de la ecuación de Sylvester
AX + XB = C, constituye una mejora significativa en
relación al algoritmo de Bartels-Stewart. La matriz A es
reducida por transformaciones ortogonales a la forma Hessenberg
superior.
La matriz B es reducida a la forma real superior de
Schur
Obtenemos ahora la ecuación
transformada
Donde
Para su resolución, si asumimos que 
tenemos que:
(2)
Siendo:
De ese modo podemos calcular 
a partir de los previamente calculados
resolviendo un
sistema de
ecuaciones
lineales de nxn con la matriz de los coeficientes de
Hessenberg. Este sistema puede resolverse fácilmente
empleando la eliminación Gaussiana con pivoteo parcial en
n2 flops
Si por el contrario 
tenemos entonces
(3)
Esto es un sistema lineal de 2nx2n para 
. Ordenando las variables
convenientemente la matriz de coeficiente es obtenida como una
matriz triangular superior con dos subdiagonales distintas de
cero. Utilizando la eliminación Gaussiana con pivoteo
parcial el sistema puede ser resuelto en 6n2 flops. .
Para realizar los cálculos se necesita un espacio de
trabajo de 2n2.
El incremento en el espacio de almacenamiento
necesario para este algoritmo es en parte compensado por el hecho
de que la matriz ortogonal U puede almacenarse en forma
factorizada debajo de la diagonal de la matriz 
.Esto implica que ya no
será necesario un arreglo de n x n para la matriz U
como sucede en el algoritmo de Bartels-Stewart.
Algoritmo
Hessenberg-Schur para la resolución de la ecuación
de Sylvester.
Reducir la matriz A a la forma Hessenberg
superior
(4)
Reducir BT a la forma real superior de
Schur
(5)
Calcular 
(6)
Para k=m,m-1,….,1
Si .
Calcular k a partir de (2) (7)
Si .
Calcular k-1, k a partir de(3) (8)
Calcular 
(9)
En esta sección haremos una breve descripción de las principales bibliotecas de
software empleadas en la programación del algoritmo para la
resolución de la ecuación de Sylvester
MPI(Message Passing Interface) es una biblioteca
estándar de paso de mensajes que facilita el desarrollo
portable de aplicaciones paralelas.).
BLAS (Basic Linear Algebra Subprograms) es una
biblioteca compuesta por una serie de rutinas de alta calidad
orientadas a bloques para la realización de operaciones
básicas con matrices y vectores.
.
LAPACK (Linear Algebra PACKage) es una librería
de rutinas escritas en Fortran 77 que implementan un gran
número de algoritmos que resuelven problemas
estándar de álgebra
lineal. LAPACK incluye rutinas para resolver sistemas de
ecuaciones lineales, problemas de mínimos cuadrados y
problemas de valores
propios y singulares
ScaLAPACK (Scalable LAPACK) es una biblioteca que
contiene un subconjunto de rutinas del LAPACK rediseñadas
para ejecutarse eficientemente en una arquitectura de
memoria
distribuida en computadoras
paralelas MIMD.
Implementación de la resolución de la
ecuación de Sylvester
En esta sección se presentan un conjunto de
rutinas tanto secuenciales como paralelas
que resuelven la ecuación matricial lineal de
Sylvester mediante el método de
Hessenberg-Schur..
Las rutinas del algoritmo secuencial se desarrollaron en
C utilizando en lo posible llamadas a subrutinas de las
librerías LAPACK, BLAS y ScaLAPACK_ Esto proporciona
portabilidad y al mismo tiempo
eficiencia.
En cuanto a la implementación paralela se ha
tratado de transportar el código
secuencial escrito en LAPACK a código paralelo sobre
ScaLAPACK, siguiendo la filosofía de la mayor parte de las
rutinas de ScaLAPACK.
Algoritmo secuencial utilizando
LAPACK
El algoritmo secuencial fue escrito en el lenguaje C
(gcc-2.96), aprovechando las potencialidades de las bibliotecas
BLAS (v3.0) y LAPACK (v3.0), de donde se tomaron funciones tanto
para la transformación de la matriz A a Hessenberg
superior (4), como para la reducción a forma real de Schur
de la matriz B (5), o para realizar todas las operaciones
matriz-vector, incluyendo la resolución de los sistemas de
ecuaciones (6).
Como se explicó anteriormente el método de
Hessenberg-Schur permite reducir la complejidad espacial del
mismo almacenando la matriz ortogonal U en forma compacta en la
parte inferior de la descomposición Hessenber superior de
la matriz A, para lo cual se utilizó la función
dgehrd_ combinada con la función dormhr_ lo
cual permite operar sobre la matriz U sin la necesidad de
reservar memoria adicional para almacenarla.
En la resolución del sistema de ecuaciones (2) se
aprovecha el hecho de que la matriz de coeficientes se encuentra
en forma casi triangular superior y la subdiagonal inferior
distinta de cero se anula empleando rotaciones de Jacobi para lo
cual se reutilizó una función previamente
implementada.
Algoritmo paralelo utilizando
ScaLAPACK
La filosofía de diseño
seguida consiste en convertir las llamadas a BLAS en llamadas a
la librería paralela PBLAS. Esta librería contiene
versiones paralelas de los núcleos computacionales de
BLAS. La transformación a Hessenberg superior de la matriz
A (4), se realizó convenientemente utilizando la
función pdgehrd_, sin embargo el uso de la
función pdormhr_ impone algunas restricciones con
respecto a las matrices A y C (6).
Validación
de las soluciones
obtenidas
La validación tanto del algoritmo secuencial como
el paralelo se realizó sustituyendo los valores de
la matriz X obtenida como solución por los algoritmos
implementados en la parte izquierda de la ecuación de
Sylvester AX + XB = C , restándole luego la matriz C y
calculando la norma de Frobenius de este resultado. 
Análisis
de los resultados experimentales
El análisis experimental de los algoritmos se
realiza a nivel de prestaciones.
Las prestaciones de los algoritmos secuenciales se evalúan
mediante el tiempo de ejecución secuencial
(TS). Los algoritmos paralelos se evalúan
mediante el tiempo de ejecución paralelo (Tp) en p
procesadores,
la ganancia de velocidad, Sp=
TS /Tp y la eficiencia, Ep= TS /(Tp*
p).
La evaluación
de prestaciones realizada pretende determinar el comportamiento
en tiempo de ejecución y eficiencia de los algoritmos
implementados a medida que varía la dimensión del
problema o el número de procesadores en el algoritmo
paralelo. Además se estudió el comportamiento del
algoritmo frente a variaciones de determinados parámetros
como el tamaño de bloque o las dimensiones de la malla de
procesadores en el algoritmo paralelo.
La plataforma empleada para la ejecución de los
experimentos
consistió en un cluster de 6 PC´s Pentium III a
833Mhz y 256MBytes de memoria interconectadas por un switch
FastEthernet a 100Mbs. Para esta red se estimaron los valores
que caracterizan al tiempo de comunicación
τ(latencia) y
β(tiempo de arranque). Para el cálculo de
los mismos se enviaron paquetes desde n=1 hasta n=1000, y luego
por el método de ajuste por mínimos cuadrados se
determinaron los valores de estos parámetros para el
modelo lineal
.
τ ≈ 1.80246 *
10-7s
β ≈ 1.4 *
10-4s
siendo el tiempo de comunicaciones. Para el valor de un
Flop utilizaremos la aproximación de 1.6935E-08 s. El
sistema operativo
del cluster es Linux Red Hat
7.3.
Para la corrida de los experimentos configuramos varias
topologías teniendo como máximo 6
procesadores por lo que obtuvimos las siguientes
configuraciones:
• Prueba secuencial
• 2 procesadores: 2×1, 1×2
• 3 procesadores: 1×3, 3×1
• 4 procesadores(doble de 2): 1×4, 4×1 y
2×2
• 6 procesadores (doble de 3): 1×6, 6×1, 2×3 y
3×2
En el algoritmo paralelo en la distribución de los datos se
consideraron bloques cuadrados de tamaño 16, 32 y
64.
Para los datos se tomaron matrices cuadradas de
magnitudes 128, 256, 512, 1024 y 2048. Las mismas fueron
generadas con la subrutina DLATM5 la cual genera las
matrices de prueba involucradas en la ecuación
generalizada de Sylvester. La misma puede obtenerse del
repositorio de matrices de prueba del LAPACK, ver
http://www.netlib.org/lapack/testing/matgen/
Prestaciones del
algoritmo secuencial
Como se explicó con anterioridad el algoritmo
secuencial fue escrito en lenguaje C y
haciendo uso de las potencialidades de las bibliotecas BLAS y
LAPACK. En la gráfica de la figura 1 se muestran los
tiempos de ejecución del algoritmo. Donde se puede
apreciar el alto costo computacional del mismo.
Figura 1Tiempo de ejecución del algoritmo
secuencial Hessenberg-Schur para la resolución de la
ecuación de Sylvester
Prestaciones del algoritmo
paralelo
En el gráfico de la figura 2 se muestran los
tiempos de ejecución del algoritmo paralelo usando 2
procesadores. Se consideró el mejor tamaño de
bloque para cada una de las matrices utilizadas. Para dimensiones
pequeñas de las matrices la malla de procesadores que
mejores resultados alcanza es la de 1×2, sin embargo al aumentar
las dimensiones del problema se invierten las prestaciones,
resultando superior la topología 2×1.
Figura 2: Tiempo de ejecución del algoritmo
paralelo
Hessenberg-Schur usando 2 procesadores
La Eficiencia alcanzada por el algoritmo paralelo para 4
procesadores se detalla en la figura 3. Se observa como para
valores pequeños de las matrices sin importar el tipo de
topología usada los valores de eficiencia son muy bajos
debido a que el parámetro β
de la red de interconexión
(β ≈ 1.4 * 10-4s)
es muy superior al tiempo de un flop (F≈1.69E-08 s) y por
tanto el alto costo de las comunicaciones degrada su rendimiento.
Se aprecia como lógica
consecuencia que con el aumento de la dimensión de las
matrices aumenta la eficiencia alcanzada.
Figura 3Eficiencia del algoritmo paralelo con 4
procesadores Hessenberg-Schur
En el gráfico de la figura 4 se muestra la
ganancia de velocidad obtenida por el algoritmo Hessenber-Schur
paralelo ejecutándose sobre 6 procesadores y distintos
valores de n. En el mismo podemos apreciar como a partir de 512
las distribuciones bidimensionales logran un mayor desempeño.
Figura 4 Ganancia de velocidad del algoritmo
Hessenberg-Schur paralelo para 6 procesadores
El objetivo de este trabajo consistía en la
implementación tanto secuencial como paralela del
algoritmo de Hessenberg-Schur para la resolución de la
ecuación de Sylvester. Este es uno de los métodos
directos más utilizados para la resolución de dicha
ecuación.
La implementación debería ser eficiente y
portable. Lo cual fue logrado mediante la utilización de
librerías estándar de algebra lineal
que incluyen un conjunto de núcleos computacionales que
son muy utilizados en la mayoría de los algoritmos de
algebra matricial.
Como se pudo comprobar debido al alto costo
computacional de este algoritmo es necesario abordar la
paralelización del mismo mediante algoritmos escalables,
portables y eficientes aspectos que se tratan de garantizar con
la utilización de las bibliotecas BLACS, PBLAS y
ScaLAPACK.
Para mantener una mayor compatibilidad con las
bibliotecas antes mencionadas se escogió para el
desarrollo del algoritmo un particionado cíclico por
bloques de los datos, y un mallado de procesadores de una y dos
dimensiones. A través del análisis de los
resultados experimentales se puede apreciar la influencia del
tiempo de comunicación sobre el rendimiento del algoritmo
paralelo, siendo más notable esto para valores
pequeños de n.
Se pudo verificar la factibilidad de
realizar una versión paralela del algoritmo de
Hessenberg-Scur para la resolución de la ecuación
de Sylvester logrando una eficiencia mayor del 70 % para
tamaños de matriz de 1024 y empleando hasta 6
procesadores.
También se pudo observar una tendencia a obtener
mejores resultados cuando la topología de procesadores es
bidimensional, lo cual posibilita una mejor interacción con otros algoritmos
construidos sobre ScaLAPACK.
[AL98] F.A. Aliev and V.B. Larin. "Optimization of
Linear Control Systems: Analytical Methods and Computational
Algorithms". Volume 8 of Stability and Control: Theory, Methods
and Applications. Gordon and Breach, 1998.
[BS72] R. Bartels and G. Stewart. "Solution of the
matrix
equation AX + XB = C: Algorithm 432". Comm. ACM,
15:820–826, 1972.
[BQOQO04]P. Benner, E. Quintana-Ortí, and G.
Quintana-Ortí. "Solving Stable Sylvester Equations via
Rational Iterative Schemes". Sonderforschungsbereich 393,
2004
[CR96]D. Calvetti and L. Reichel. "Application of ADI
iterative methods to the restoration of noisy images". SIAM J.
Matrix Anal. Appl., 17:165–186, 1996.
[Dat03]B. Datta. "Numerical Methods for Linear Control
Systems Design and Analysis". Elsevier Press, 2003.
[DOR88] L. Dieci, M.R. Osborne, and R.D. Russell. "A
Riccati transformation method for solving linear bvps. I:
Theoretical aspects." SIAM J. Numer. Anal.,
25(5):1055–1073,
1988.
[GGKK03] A. Grama, A. Gupta, G. Karypis, V. Kumar.
"Introduction to Parallel Computing".Second Edition. Addison
Wesley, ISBN : 0-201-64865-2, 2003.
[GNVL79] G. H. Golub, S. Nash, and C. F. Van Loan ."A
Hessenberg–Schur method for the problem AX + XB = C". IEEE
Trans. Automat. Control, AC-24:909–913, 1979.
[GVL96] G.H. Golub and C.F. Van Loan. "Matrix
Computations". Johns Hopkins University Press, Baltimore, third
edition, 1996.
[Jbi98] K. Jbilou "Krylov subspace methods for solving
the matrix Sylvester equation AX-XB=C "Technical report LMPA (71)
Universit e du Littoral Calais France
(1998).
[JR] K. Jbilou and J. Riquet "A nonsymmetric block
Lanczos algorithm for eigenvalue Problems".
[LT85] P. Lancaster and M. Tismenetsky. "The Theory of
Matrices". Academic Press, Orlando, 2nd edition, 1985.
[Sad93] M. Sadkane "Block Arnoldi and Davidson methods
for unsymmetric large eigenvalue problems Nume. Math. 64(1993) pp
687-706.
Autor:
Ernesto Díaz López
Jorge Luis Díaz Flores
Universidad central "Martha Abreu" de Las
Villas