Implementaciones secuenciales y paralelas del algoritmo Hessenberg-Schur

Resumen

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.

Introducción

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)

Bibliotecas de Software

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

Conclusiones

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.

Bibliografía

[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