Precondicionamiento de grandes sistemas lineales acoplados para descomposición de dominios

Luis Manuel Hernández Ramos

Escuela de Computación, Universidad Central de Venezuela Caracas-Venezuela ^{1}

1Trabajo realizado en parte durante una pasantia en el Laboratorio de Analisis Numérico de la Universidad ”Pierre et Marie Curie”(Paris VI) Paris, Francia con ayuda del CDCH-UCV y del Programa ALFA de la C.E.E.

A la Memoria de Rafael Anthony Underwood 1

Agradecimientos

Debo agradecer a un gran número de personas e instituciones que hicieron posible este trabajo. En primer lugar, a la Prof. Hilda López quien se empeñó desde hace mucho tiempo en sus labores de Madre Académica e inclusive de Madre de Verdad y me guío siempre por el camino correcto. A ella debo el haber cumplido con este sueño de viajar a Paris, Francia, una de las ciudades mas bellas del mundo a trabajar en un laboratorio de primer orden como lo es el Laboratorio de Analisis Numerico de la Universidad ”Pierre et Marie Curie”. Tambien debo agradecer al Profesor Pascal Joly, quien fue mi tutor durante mi estadía en Francia, le agradezco sobretodo su amabilidad, cortesía y orientacion académica. A mi familia por haberme apoyado siempre. A la profesora Vivette Girault quien fue un angel salvador para mi en todo momento, le agradezco haber compartido este proyecto, haberse ocupado de mi en Francia, haber hecho todas esas diligencias, haberme prestado dinero, etc, etc, etc, y sobretodo su inmensa calidad humana. A la profesora Brigida Molina quien me ayudo de las más variadas formas: Desde invitarme a su casa a cenar y haber sentido la calidez hogareña que me hacia falta despues de tantos días solos y frios, pasando por haberme prestado la bufanda que utilicé durante todo el otoño, invierno y parte de la primavera, hasta la orientación técnica, academica y personal decisiva y oportuna para el éxito de este trabajo. Al profesor Marcos Raydan, por sus orientaciones y apoyo. A mis companeros en el laboratorio de Analisis Numerico en Paris VI. A Claudia León y a su hijo Diego, quienes me brindaron su camaradería y me permitieron tener un refugio Venezolano en Paris. A mis vecinos y compañeros de la MINA, copartícipes de una de las experiencias mas hermosas que he tenido, como lo es vivir en hermandad, en una comunidad de personas de todos los paises. A todos los receptores de ”Voyage”en Internet, quienes me dieron palabras de estimulo y apoyo, y me hicieron sentir siempre acompañado. No puedo dejar de agradecer a la Cité Universitaire de Paris y a su RestaU, sin el cual no hubiese podido sobrevivir. A Meli por haberme acompañado y apoyado todo este tiempo. Y a todas aquellas personas que en algún momento pusieron su granito de arena en pro de la realizacion de este sueño.

GRACIAS A TODOS!!!.

Índice General

1 Introducción 5

2 Métodos Iterativos 7

2.1 Métodos de Proyección . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Definiciones Básicas y Algoritmos . . . . . . . . . . . . . . 7 2.1.2 Procesos de Proyección Unidimensional . . . . . . . . . . 10

2.2 Métodos de Krylov . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.2 Subespacios de Krylov . . . . . . . . . . . . . . . . . . . . 12 2.2.3 Método de Arnoldi . . . . . . . . . . . . . . . . . . . . . . 13 2.2.4 Método de Arnoldi para Sistemas Lineales (FOM) …. 15 2.2.5 Método GMRES . . . . . . . . . . . . . . . . . . . . . . . 17

3 Precondicionamiento 20 3.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Operaciones de Precondicionamiento . . . . . . . . . . . . . . . . 21

3.2.1 Precondicionamiento por la Izquierda . . . . . . . . . . . 22 3.2.2 Precondicionamiento por la Derecha . . . . . . . . . . . . 23 3.2.3 Precondicionamiento SPLIT o dividido . . . . . . . . . . . 24 3.2.4 Variantes Flexibles . . . . . . . . . . . . . . . . . . . . . . 24

3.3 Algunas Técnicas de Precondicionamiento . . . . . . . . . . . . . 25 3.3.1 Jacobi, SOR, SSOR . . . . . . . . . . . . . . . . . . . . . 25 3.3.2 Factorización Incompleta (ILU) . . . . . . . . . . . . . . . 27 3.3.3 Inversas Aproximadas . . . . . . . . . . . . . . . . . . . . 29

4 Método de Descomposición de Dominios 31 4.1 Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Descomposición de un Problema de S.E.L. . . . . . . . . . . . . . 31 4.3 Precondicionamiento para Problemas de Descomposición de Do-

minios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.3.1 Precondicionadores en Bloques Diagonales . . . . . . . . . 34 4.3.2 Precondicionamiento por una aproximación de la factor-

izacion de Gauss . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.3 Precondicionamiento por Pseudoinversa . . . . . . . . . . 36 4.3.4 El Precondicionador de Pseudoinversa desde el punto de

vista de Factorizaciones Aproximadas . . . . . . . . . . . 38 4.3.5 Escogencias para C . . . . . . . . . . . . . . . . . . . . . . 38 4.3.6 Implementacion práctica del precondicionador por pseu-

doinversa . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5 Experimentos 40

6 Conclusiones 41

Capítulo 1

Introducción

La resolución numérica de sistemas de ecuaciones lineales dispersos sigue siendo hoy en dia un tema activo de investigación; debido, en parte, a la enorme influencia que éste tiene sobre otras áreas de desarrollo: principalmente en la resolucion numérica de ecuaciones diferenciales en derivadas parciales, provenientes de la simulación de problemas de la fisica y la ingeniería, aplicados a la industria y a la investigación misma. Los problemas a los cuales se enfrentan la industria y la investigación, requieren la simulación de situaciones cada vez más precisas, y estas simulaciones requieren a menudo, la resolución de sistemas de ecuaciones lineales cada vez más grandes, y que exigen el uso de computadores cada vez más potentes. El advenimiento de las máquinas paralelas ha permitido afrontar muchos de esos nuevos problemas, y ha abierto nuevos horizontes en la investigación y el desarrollo, pero igualmente, ha multiplicado los retos. La creación de programas paralelos no ha quedado solo en adaptar los algoritmos existentes a trabajar sobre máquinas paralelas, sino que ha impulsado la tarea de replantear los problemas bajo una óptica diferente, dando nacimiento a una nueva disciplina de investigación que requiere el trabajo conjunto de especialistas en diferentes áreas. Los esfuerzos en la búsqueda de algoritmos paralelos que trabajen de manera general sobre los diferentes problemas e inclusive sobre las diferentes arquitecturas se han visto complicados por una gran dependencia entre dichos componentes, y esta circunstancia ha impedido igualmente su aplicación a niveles menos especializados, ya que es menester desarrollar algoritmos específicos que aprovechen las particularidades de cada problema y cada arquitectura en detrimento grave de su generalidad. Los métodos de descomposición de dominios representan una alternativa para aprovechar el paralelismo inherente en cierto tipo de problemas. Dichos métodos se refieren a descomponer un problema global en varios subproblemas locales independientes, mediante la introducción de un parámetro de relajación conocido como multiplicador de Lagrange, el cual modifica los datos de los subproblemas locales, a manera de garantizar su consistencia con el problema global.

La resolución, tanto de los subproblemas locales como del multiplicador de Lagrange, pueden realizarse mediante métodos directos o iterativos, segun la conveniencia del caso. En los últimos años ha habido un auge por el desarrollo de nuevos métodos iterativos para la resolución numérica de sistemas de ecuaciones lineales, basados en proyecciones sobre espacios de Krylov. Estos métodos se han implementado con éxito sobre máquinas paralelas. Pero la debilidad de los métodos iterativos radica en su falta de robustez y predictabilidad. Es por ello que es necesario desarrollar precondicionadores que no solo aseguren buenas caracteristícas de convergencia para los métodos, sino que tambien se adapten bien a su implementación en paralelo. En este trabajo, se aprovechan las caracteristícas inherentes a las matrices de sistemas acoplados provenientes de los métodos de descomposición de dominios, para proponer precondicionadores que permitan su resolución eficiente por medio de técnicas iterativas. El documento se divide de la siguiente manera: En el capítulo 2 se hará un breve resumen de algunas de las técnicas iterativas más importantes actualmente en la resolución de grandes sistemas de ecuaciones lineales, como son los métodos de proyección sobre subespacios de Krylov, y se mostrarán algunos casos importantes tales como FOM y GMRES. En el capítulo 3 se hablará un poco del significado y la utilidad del precondicionamiento, se mostrarán las operaciones de precondicionamiento más comunes y se ejemplificarán en el caso del método GMRES. Posteriormente se describirán brevemente algunas de las técnicas más importantes para conseguir el precondicionador, en particular, se hablará de los precondicionadores Jacobi y SOR, de factorizaciones aproximadas LU y de inversas aproximadas. Luego, en el capítulo 4 se desarrollará una técnica de descomposición de dominios, para el caso de 2 subdominios; y se mostrarán técnicas de precondicionamiento específicas para el sistema resultante, tales como: los precondicionadores en bloques diagonales, las factorizaciones de Gauss aproximadas y la propuesta de este trabajo que es el precondicionamiento por pseudoinversa. En el capítulo 6 se mostrarán los resultados de algunos experimentos de la aplicación del precondicionamiento por pseudoinversa, comparándolos con el precondicionamiento en bloques diagonales y finalmente, en el capítulo 7 se darán las conclusiones.

Métodos Iterativos

En este capítulo se tratarán brevemente los métodos iterativos basados en proyecciones sobre los subespacios de Krylov. En la primera parte, se explicarán de manera general los métodos de proyección, se mostrarán algunos casos sencillos de proyecciones sobre subespacios unidimensionales; y en la segunda parte se abordarán los importantes métodos de proyección en el caso específico, donde se proyecta sobre subespacios de Krylov. Este capítulo es un resumen breve del libro de la persona quien propuso el método GMRES, Yousef Saad [1] sobre los métodos iterativos.

2.1 Métodos de Proyección

2.1.1 Definiciones Básicas y Algoritmos

Muchas de las técnicas prácticas para resolver Sistemas de Ecuaciones lineales involucran de una u otra forma procesos de proyección. Un proceso de proyección representa una manera de obtener una aproximacion para la solución de un sistema lineal desde un subespacio. En este capítulo se presentaran dichas tecnicas desde un punto de vista muy general al igual que algunas nociones teóricas interesantes. Consideremos el Sistema de Ecuaciones Lineales:

Ax = b,  (2.1)

donde A es una matriz real de n×n. La idea de técnicas de proyeccion consiste en extraer la solución aproximada del problema anterior desde un subespacio de R^{n}. Si K es este subespacio de candidatos aproximantes o subespacio de búsqueda y m es su dimensión, entonces, en general, se imponen m restricciones para extraer tal aproximación. Una manera típica de describir dichas restricciones es imponer m (independientes) condiciones de ortogonalidad. Especificamente, el vector residual b − Ax se restringe a ser ortogonal a m vectores pertenecientes a R^{n} linealmente independientes. Estos definen otro subespacio L de dimensión m el cual denominaremos subespacio de restricciones o subespacio izquierdo.

Existen dos grandes clases de métodos de proyección: ortogonal y oblícua. En una técnica de proyeccion ortogonal, el subespacio L es el mismo K. En un método de proyección oblícua, L es diferente de K y puede estar totalmente desligado a él.

Métodos de Proyección General

Sea A una matriz real de n×n y sean K y L dos subespacios de R^{n} de dimensión m. Una técnica de proyección sobre un subespacio K y ortogonal a L es un proceso el cual encuentra una solución aproximada x̃ en K y con un nuevo vector residual ortogonal a L,

Encontrar x̃ ∈ K tal que b − Ax̃ ⊥ L.

Si se utiliza un valor aproximado inicial x_{0} para la solución, entonces la aproximación se hará sobre el espacio afin x_{0} + K en vez del espacio vectorial homogéneo K. Esto requiere una ligera modificación de la formulación anterior. El problema de aproximación debe ser redefinido como:

Encontrar x̃ ∈ x_{0} + K tal que b − Ax̃ ⊥ L.

Notese que x̃ es escrito de la forma x̃ = x_{0} +δ, y el vector residual r_{0} es definido como,

r_{0} = b − Ax_{0},  (2.2)

entonces la ecuación anterior se escribe como b − A(x_{0} + δ) ⊥ L o también,

r_{0} − Aδ ⊥ L.  (2.3)

En otras palabras, la solución aproximada puede ser definida como:

x̃ = x_{0} + δ, δ ∈ K,  (2.4)
(r_{0} − Aδ, w) = 0, ∀w ∈ L.  (2.5)

Este es el paso básico de proyección en su forma más general. Muchas técnicas estándar utilizan una sucesión de tales proyecciones. Típicamente, un nuevo paso de proyección utiliza un nuevo par de subespacios K y L y un aproximado inicial x_{0} igual a la más reciente aproximación obtenida del paso previo de proyección. Los métodos de proyección comprenden muchos de los métodos mejor conocidos en computación científica, de hecho, virtualmente todas las técnicas iterativas básicas pueden ser consideradas técnicas de proyección.

Representación matricial

Sea V = [v_{1}, …, v_{m}], una matriz n × m cuyas vectores columna forman una base de K y, similarmente, W = [w_{1}, …, w_{m}], una matriz n×m cuyas columnas forman una base de L. Si la solución aproximada es escrita como

x = x_{0} + V y,

entonces la condición de ortogonalidad determina el siguiente sistema de ecuaciones para el vector y:

W^{T} AV y = W^{T} r_{0}.

Si se supuso que la matriz W^{T} AV de m×m es no singular, la siguiente expresión para la solución aproximada x̃ resulta,

x̃ = x_{0} + V (W^{T} AV)^{−1}W^{T} r_{0}.  (2.6)

Un método de proyección prototipo es representado por el siguiente algoritmo: Algoritmo: Método de Proyección Prototipo 1.- Mientras no-convergencia hacer 2.- Seleccionar los subespacios K y L 3.- Elegir bases V = [v_{1}, ..., v_{m}] y W = [w_{1}, ..., w_{m}] para K y L respectivamente. 4.- r := b − Ax 5.- y := (W^{T} AV)^{−1}W^{T} r 6.- x := x + V y 7.- fmientras La solución aproximada es definida sólo cuando la matriz W^{T} AV es no singular, condición que no es necesariamente garantizada cuando A es no singular. Sin embargo, existen dos casos particulares importantes en donde la no singularidad de la matriz W^{T} AV es garantizada, estos serán discutidos en la siguiente proposición.

Proposición 2.1 Si A, L, K satisfacen las siguientes condiciones, i. A es positivo definida y L = K ii. A es no singular y L = AK Entonces la matriz B = W^{T} AV es no singular cualquiera sean las bases V y W de K y L, respectivamente.

Demostración Consideremos el caso (i). Sea V cualquier base de K y W cualquier base de L. De hecho, dado que L y K son los mismos, W siempre puede ser expresado como W = V G, donde G es alguna matriz no singular m × m. Entonces,

B = W^{T} AV = G^{T} V^{T} AV.

Dado que A es positivo definida, también lo es V^{T} AV y esto muestra que B es no singular. Consideremos ahora el caso (ii). Sea V cualquier base de K y W cualquier base de L. Dado que L = AK, W puede ser expresado en este caso como W = AV G en donde G es alguna matriz no singular m × m. Entonces

B = W^{T} AV = G^{T} (AV)^{T} AV.  (2.7)

Dado que A es una matriz no singular m × m, AV tiene rango completo y como resultado (AV)^{T} AV es no singular. Por lo tanto B es no singular.

Ahora consideremos el caso donde A es simétrica(real) y una técnica de proyección ortogonal es utilizada. En esta situación, la misma base puede ser usada para L y K, los cuales son subespacios identicos, y la matriz proyectada, la cual es B = V^{T} AV, es simétrica. En adición, si la matriz A es simétrica y positivo definida, entonces B también lo es.

2.1.2 Procesos de Proyección Unidimensional

En esta sección examinaremos algunos ejemplos dados por los procesos de proyección más simples en una dimensión. En lo que sigue, el vector r denota el vector residual b − Ax para la aproximación actual x. Los procesos de proyección en una dimensión son definidos cuando,

K = span{v} y L = span{w},

donde v y w son dos vectores. En este caso, la nueva aproximación toma la forma x ← x + αv y la condición es: r − Aδ ⊥ w, lo que implica,

α = \frac{(r,w)}{(Av,w)}.

A continuación, las tres elecciones más populares serán consideradas.

Algoritmo del Descenso más Rápido

El algoritmo del descenso más rápido se define para el caso donde la matriz A es Simétrica Positivo Definida. Este consiste en tomar a cada paso v = r y w = r. Esto da como resultado una iteración que es descrita por el siguiente algoritmo. Algoritmo: Descenso más Rápido 1.- Mientras no-convergencia hacer 2.- r ← b − Ax 3.- α ← (r, r)/(Ar, r) 4.- x ← x + αr 5.- Fin Cada paso de la iteración anterior minimiza

f(x) = ∥x − x_{∗}∥_{A}^{2} = (A(x − x_{∗}), (x − x_{∗})),

sobre todos los vectores de la forma x + αd, donde d es el negativo de la dirección del gradiente −∇f. El negativo de la dirección del gradiente es localmente la dirección de decrecimiento mayor, dando la velocidad de convergencia más grande para f.

Algoritmo de Mínimo Residual

Ahora asumiremos que A no es necesariamente simétrica, sino solamente positivo definida, es decir, su parte simétrica A + A^{T} es simétrico positivo definida. Tomando a cada paso v = r y w = Ar, el siguiente procedimiento iterativo resulta. Algoritmo: Minimo Residual 1.- Hasta convergencia, hacer 2.- r ← b − Ax 3.- α ← (Ar, r)/(Ar, Ar) 4.- x ← x + αr 5.- Fin Aqui, cada paso minimiza f(x) = ∥b − Ax∥_{2}^{2} en la dirección de r.

Algoritmo de Descenso más Rápido de la Norma del Residual

En este algoritmo, no es necesario suponer que la matriz A sea positivo definida, de hecho, el único requerimiento es que A sea una matriz cuadrada no singular. A cada paso el algoritmo utiliza v = A^{T} r y w = Av, escrito directamente, el algoritmo da lugar a una secuencia de operaciones que requiere tres productos matriz-vector, el cual es tres veces más costoso que los algoritmos anteriores. Sin embargo, el número de operaciones matriz-vector puede ser reducido a dos por iteración calculando el residual de manera diferente; tal como la variante que viene a continuación: Algoritmo: Descenso más Rápido de la norma del Residual 1.- Calcular r := b − Ax 2.- Mientrasno-convergencia hacer 3.- v := A^{T} r 4.- Calcular Av y α := ∥v∥_{2}^{2}/∥Av∥_{2}^{2} 5.- x := x + αv 6.- r := r − αAv 7.- Fin Aqui, cada paso minimiza f(x) = ∥b−Ax∥_{2}^{2} en la dirección −∇f. Este algoritmo es equivalente al algoritmo de Descenso del Gradiente aplicado al Sistema de Ecuaciones Normales A^{T} Ax = A^{T} b. Dado que si A es no singular, entonces A^{T} A es positivo definida, se tiene que el método converge siempre que A sea no singular.

2.2 Métodos de Krylov

2.2.1 Introducción

En este capítulo se examinan algunos métodos que son considerados actualmente como las técnicas iterativas más importantes en la resolución de sistemas de ecuaciones lineales. Dichas técnicas están basadas en procesos de proyección, tanto ortogonales como oblícuos, sobre subespacios de Krylov, los cuales son los subespacios generados por vectores de la forma p(A)v donde p es un polinomial. En resumen, esas técnicas aproximan A^{−1}b por p(A)b, donde p es un ”buen”polinomial. Recordemos que un método de proyección general para resolver el sistema lineal,

Ax = b,  (2.8)

es un método el cual busca una solución aproximada x_{m} sobre un espacio afin x_{0} + K_{m} de dimensión m mediante la imposición de la condición de Petrov-Galerkin,

b − Ax_{m} ⊥ L_{m},

donde L_{m} es otro subespacio de dimensión m. Aqui, x_{0} representa una solución inicial aproximada arbitraria. Un método de subespacio de Krylov es un método para el cual el subespacio K_{m} es el subespacio,

K_{m}(A, r_{0}) = span{r_{0}, Ar_{0}, A^{2}r_{0}, ..., A^{m−1}r_{0}},

donde r_{0} = b − Ax_{0}. Cuando no exista ambiguedad, K_{m}(A, r_{0}) sera denotado simplemente por K_{m}. Las diferentes versiones de los métodos basados en subespacios de Krylov vendran dados por las diferentes elecciones del subespacio L_{m} y por la manera en como el sistema es precondicionado. Visto desde el ángulo de la teoría de aproximación, es claro que las aproximaciones obtenidas por los métodos basados en subespacios de Krylov son de la forma

A^{−1}b = x_{m} = x_{0} + q_{m−1}(A)r_{0}

en el cual q_{m−1} es algún polinomio de grado m − 1. Aunque todas las técnicas proveen el mismo tipo de aproximación polinomial, la elección de L_{m}, es decir, de las restricciones utilizadas para construir tal aproximación, tienen un efecto importante en la técnica iterativa. Para las técnicas mejor conocidas, existen dos fuertes elecciones. La primera es simple L_{m} = K_{m}, y la otra conocida como variación de minimo residual es L_{m} = AK_{m}.

2.2.2 Subespacios de Krylov

En esta sección consideraremos los métodos de proyección sobre los subespacios de Krylov, es decir, sobre los subespacios de la forma

K_{m} ≡ span{r_{0}, Ar_{0}, A^{2}r_{0}, ..., A^{m−1}r_{0}},  (2.9)

el cual denotaremos simplemente por K_{m} cuando no exista ambiguedad. La dimensión del subespacio de aproximantes se incrementa en uno en cada paso del proceso de aproximación. Unas pocas propiedades elementales de los subespacios de Krylov pueden ser establecidas. La primera propiedad, es que K_{m} es el subespacio de todos los vectores en R^{n} los cuales pueden ser escritos como x = p(A)v, en donde p es un polinomio de grado ≤ m − 1.

2.2.3 Método de Arnoldi

El procedimiento de Arnoldi es un algoritmo para construir una base ortogonal de un subespacio de Krylov K_{m}. En aritmética exacta, una variante del algoritmo es la siguiente: Algoritmo: Arnoldi 1.- Elegir un vector v_{1} de norma igual a 1. 2.- Para j = 1, .., m hacer 3.- Calcular h_{ij} = (Av_{j}, v_{i}) para i = 1, 2, .., j 4.- Calcular w_{j} := Av_{j} − ∑_{i=1}^{j} h_{ij}v_{i} 5.- h_{j+1,j} = ∥w_{j}∥_{2} 6.- Si h_{j+1,j} = 0 entonces Stop 7.- v_{j+1} = w_{j}/h_{j+1,j} 8.- fpara A cada paso, el algoritmo multiplica el vector de Arnoldi previo v_{j} por A y entonces ortonormaliza el vector resultante w_{j} con respecto a todos los anteriores v_{i}’s por un procedimiento de Gram-Schmidt standard. Este procedimiento se detiene si el vector w_{j} calculado en la linea 4 se hace cero.

Proposición 2.2 Supongamos que el algoritmo anterior no se detiene hasta el paso m-ésimo. Entonces los vectores v_{1}, v_{2}, ..., v_{m} forman una base ortonormal del espacio de Krylov

K_{m} = span{v_{1}, Av_{1}, ..., A^{m−1}v_{1}}

.

Demostración Los vectores v_{j}, i = 1, 2, …, m son ortonormales por construcción. El hecho de que ellos generen K_{m} sigue del hecho de que cada vector v_{j} es de la forma q_{j−1}v_{1} donde q_{j−1} es un polinomio de grado j − 1. Esto puede ser probado por inducción sobre j como haremos. El resultado es obvio para j = 1, dado que v_{1} = q_{0}(A)v_{1} con q_{0}(t) ≡ 1. Supongamos que el resultado es cierto para todos los enteros ≤ j y consideremos v_{j+1}. Tenemos,

h_{j+1}v_{j+1} = Av_{j} − ∑_{i=1}^{j} h_{ij}v_{i} = Aq_{j−1}(A)v_{1} − ∑_{i=1}^{j} h_{ij}q_{i−1}(A)v_{1}  (2.10)

Lo cual muestra que v_{j+1} puede ser expresado como q_{j}(A)v_{1} donde q_{j} es de grado j y esto completa la demostración.

Proposición 2.3 Denotemos por V_{m}, la matriz n × m cuyos vectores columna son los vectores v_{1}, …, v_{m}, por H̄_{m}, la matriz (m+1)×m de Hessemberg cuyos elementos no cero h_{ij} son definidos por el algoritmo de Arnoldi, y por H_{m} la matriz obtenida desde H̄_{m} eliminandole su última fila. Entonces el siguiente resultado es obtenido:

AV_{m} = V_{m}H_{m} + w_{m}e_{m}^{T}  (2.11)
= V_{m+1}H̄_{m},  (2.12)
V_{m}^{T} AV_{m} = H_{m}.  (2.13)

Demostración La relación 2.12 viene directamente de la siguiente igualdad generada por las lineas 4, 5 y 7 del Algoritmo de Arnoldi,

Av_{j} = ∑_{i=1}^{j+1} h_{ij}v_{i}, j = 1, 2, …, m  (2.14)

La relación 2.11 es la reformulación matricial de 2.14. La relación 2.13 proviene de multiplicar ambos lados de 2.11 por V_{m}^{T} y haciendo uso de la ortonormalidad de {v_{1}, …, v_{m}}. Observemos que el algoritmo se detiene en el caso en que la norma de w_{j} se hace cero a un cierto paso j. En este caso, el vector v_{j+1} no puede ser calculado y el algoritmo se detiene. Cuando esto sucede quiere decir que el subespacio K_{j} es invariante bajo A y por lo tanto no tiene sentido seguir iterando ya que estariamos proyectando siempre sobre el mismo subespacio K_{j}.

Implementaciones Prácticas

En la descripción previa del proceso de Arnoldi, se supuso por simplicidad una aritmética exacta. En la práctica, se puede ganar mucho utilizando los Algoritmos de Gram-Schmitd Modificado o Householder, en vez del algoritmo stándard de Gram-Schmidt. Con la alternativa de Gram-Schmitd modificado, el proceso de Arnoldi toma la siguiente forma: Algoritmo: Arnoldi Gram-Schmidt Modificado

1. Elegir un vector v_{1} de norma= 1

2. Para j = 1 Hasta j = m Hacer

(a) Calcular w_{j} := Av_{j}

(b) Para i = 1 Hasta i = j Hacer

i. h_{ij} = (w_{j}, v_{i}) ii. w_{j} := w_{j} − h_{ij}v_{i}

(c) FPara

(d) h_{j+1,j} = ∥w_{j}∥_{2}. Si h_{j+1,j} = 0 Stop

(e) v_{j+1} = w_{j}/h_{j+1,j}

3. FPara

En aritmética exacta, este algoritmo y el algoritmo Gram-Schmidt standard son matemáticamente equivalentes. En la presencia de errores de redondeo, la formulación anterior es mucho más estable. Sin embargo, existen casos en donde la cancelación puede ser severa en los pasos de ortogonalización, de tal manera que este proceso de Gram-Schmidt sea inadecuado.

2.2.4 Método de Arnoldi para Sistemas Lineales (FOM)

Dado un aproximado inicial x_{0} del sistema lineal original Ax = b, ahora consideraremos un método de proyección ortogonal, el cual toma L = K = K_{m}(A, r_{0}), con:

K_{m}(A, r_{0}) = span{r_{0}, Ar_{0}, A^{2}r_{0}, …, A^{m−1}r_{0}},  (2.15)

en el cual r_{0} = b−Ax_{0}. Este método busca una solución aproximada x_{m} desde el subespacio afín x_{0} + K_{m} de dimensión m imponiendo la condición de Galerkin,

b − Ax_{m} ⊥ K_{m}.  (2.16)

Si v_{1} = r_{0}/∥r_{0}∥_{2} en el método de Arnoldi, y pongamos β = ∥r_{0}∥_{2}, entonces

V_{m}^{T} AV_{m} = H_{m},

por 2.13 y

V_{m}^{T} r_{0} = V_{m}^{T} (βv_{1}) = βe_{1}.

Como resultado, la solución aproximada utilizando el subespacio m-dimensional anterior es dada por:

x_{m} = x_{0} + V_{m}y_{m},  (2.17)
y_{m} = H_{m}^{−1}(βe_{1}).  (2.18)

El método basado en esta aproximación llamado Método de Ortogonalización Completa (Full Orthogonalization Method - FOM) es descrito a continuación utilizando Gram-Schmidt Modificado como paso de Arnoldi. Algoritmo: Full Orthogonalization Method (FOM)

1. Calcular r_{0} = b − Ax_{0}, β := ∥r_{0}∥_{2} y v_{1} := r_{0}/β

2. Definir la matriz m × m; H_{m} = {h_{ij}}_{i,j=1,...,m} ; Hacer H_{m} = 0

3. Para j = 1, 2, …, m Hacer

(a) Calcular w_{j} := Av_{j}

(b) Para i=1,…,j Hacer

i. h_{ij} = (w_{j}, v_{i}) ii. w_{j} := w_{j} − h_{ij}v_{i}

(c) FPara

(d) Calcular h_{j+1,j} = ∥w_{j}∥_{2}. Si h_{j+1,j} = 0 Hacer m := j y Ir a paso 4

(e) v_{j+1} = w_{j}/h_{j+1,j}

4. Fpara

5. Calcular y_{m} = H_{m}^{−1}(βe_{1}) y x_{m} = x_{0} + V_{m}y_{m} El algoritmo anterior depende del parametro m el cual es la dimensión del subespacio de Krylov. En la práctica es deseable elegir m de una manera dinámica. Esto sería posible si se dispone de la norma del residual de la solución x_{m} sin haber calculado aun x_{m}. Entonces, el algoritmo pudiera ser detenido al paso apropiado utilizando esta información. La siguiente proposición da un resultado en esta dirección.

Proposición 2.4 El vector residual de la solución aproximada x_{m} calculados por el algoritmo FOM es tal que

b − Ax_{m} = −h_{m+1,m}e_{m}^{T}y_{m}v_{m+1}

y por lo tanto,

∥b − Ax_{m}∥_{2} = h_{m+1,m}|e_{m}^{T}y_{m}|.  (2.19)

Demostración Tenemos las relaciones

b − Ax_{m} = b − A(x_{0} + V_{m}y_{m})
= r_{0} − AV_{m}y_{m}
= βv_{1} − V_{m}H_{m}y_{m} − h_{m+1,m}e_{m}^{T}y_{m}v_{m+1}

Por la definición de y_{m}, H_{m}y_{m} = βe_{1} y por tanto βv_{1} − V_{m}H_{m}y_{m} = 0 de lo cual el resultado sigue inmediatamente. Un grueso estimado del costo de cada paso del algoritmo es determinado como sigue. Si N_{Z}(A) es el número de elementos no cero de A, entonces m pasos del procedimiento de Arnoldi requieren m productos matriz-vector a un costo de 2m×N_{z}(A) operaciones. Cada paso de Gram-Schmidt cuesta aproximadamente 4 × j × n operaciones, j = 1 … m, el cual da para m pasos un aproximado de 2m^{2}n. En total, en promedio, un paso de FOM cuesta aproximadamente,

2N_{z}(A) + 2mn.

En almacenamiento, m vectores de tamaño n son requeridos para guardar la base V_{m}. Adicionalmente hay que almacenar los vectores para el lado derecho y para la solución actual y un vector temporal para el producto matriz-vector. Además, la matriz de Hessemberg H_{m} debe ser almacenada. El total es a grosso modo,

(m + 3)n + \frac{m^{2}}{2} ,

operaciones. En muchas situaciones m es pequeño relativamente a n, por lo tanto el costo es dominado por el primer término.

Variante: FOM con Reinicio

Consideremos ahora el algoritmo desde un punto de vista práctico. Cuando m crece, el costo computacional se incrementa en O(m^{2})n a causa de la ortogonalización de Gram-Schmidt. El costo de memoria se incrementa en O(mn). Para valores grandes de n se límita el valor más grande de m que puede ser utilizado. Existen dos remedios. El primero es reinicializar el algoritmo periódicamente y el segundo es ”truncar”la ortogonalización en el algoritmo de Arnoldi. Veamos la primera opción. Algoritmo: FOM con Reinicio (FOM(m))

1. Calcular r_{0} = b − Ax_{0}, β = ∥r_{0}∥_{2} y v_{1} = r_{0}/β

2. Generar la base de Arnoldi y la matriz H_{m} utilizando el algoritmo de Arnoldi arrancando con v_{1}.

3. Calcular y_{m} = H_{m}^{−1}βe_{1} y x_{m} = x_{0} + V_{m}y_{m}. Si satisfecho Stop

4. Hacer x_{0} := x_{m} y goto paso 1

El Algoritmo de Gradiente Conjugado

El bien conocido algoritmo de Gradiente Conjugado, el cual es una de las mejores técnicas conocidas para resolver Sistemas Lineales Dispersos Positivo Definidos (Simétricos) es matemáticamente equivalente al Método FOM. Sin embargo, dado que A es simétrica, algunas simplificaciones se suceden provenientes de una recurrencia de Lanczos de tres términos, que da lugar a algoritmos más elegantes y eficientes.

2.2.5 Método GMRES

El Método de Minimo Residual Generalizado (GMRES, Generalized Minimum Residual) es un método de proyección basado en tomar K = K_{m} y L = AK_{m}, en el cual K_{m} es el m-ésimo subespacio de Krylov, y tomando v_{1} = r_{0}/∥r_{0}∥_{2}. Como vimos anteriormente, esta técnica minimiza la norma del residual sobre todos los vectores en x_{0} + K_{m}. La implementación de este algoritmo, basada en esta aproximación es similar a la del algoritmo FOM. A continuación describiremos la idea básica y posteriormente se discutirán algunas variaciones prácticas.

El Algoritmo Básico GMRES

Existen dos maneras de derivar el algoritmo. La primera explota la propiedad de optimalidad y la relación 2.12. Cualquier vector x en x_{0} + K_{m} puede ser escrito como,

x = x_{0} + V_{m}y,  (2.20)

donde y es un vector de tamaño m. Definiendo

J(y) = ∥b − Ax∥_{2} = ∥b − A(x_{0} + V_{m}y)∥_{2},  (2.21)

la relación 2.12 resulta en

b − Ax = b − A(x_{0} + V_{m}y),  (2.22)
= r_{0} − AV_{m}y,  (2.23)
= βv_{1} − V_{m+1}H̄_{m}y,  (2.24)
= V_{m+1}(βe_{1} − H̄_{m}y).  (2.25)

dado que los vectores columna de V_{m+1} conforman una base ortonormal, entonces

J(y) ≡ ∥b − A(x_{0} + V_{m}y)∥_{2} = ∥βe_{1} − H̄_{m}y∥_{2}.  (2.26)

La aproximación GMRES es el único vector de x_{0} + K_{m} que minimiza 2.21. Por 2.20 y 2.26, esta aproximación puede ser obtenida de manera simple por x_{m} = x_{0} + V_{m}y_{m} donde y_{m} minimiza la función J(y) = ∥βe_{1} − H̄_{m}y∥_{2}, es decir

x_{m} = x_{0} + V_{m}y_{m}, donde,  (2.27)
y_{m} = argmin_{y}∥βe_{1} − H̄_{m}y∥_{2}.  (2.28)

El minimizador y_{m} no es costoso de calcular dado que requiere la solución de un problema de minimos cuadrados (m + 1) × m donde m es típicamente pequeño. Esto resulta en el siguiente algoritmo. Algoritmo: GMRES

1. Calcular r_{0} = b − Ax_{0}, β = ∥r_{0}∥_{2}, y v_{1} = r_{0}/β

2. Definir la matriz (m+1)×m, H̄_{m} = {h_{ij}}_{i≤i≤m+1,1≤j≤m}. Hacer H̄_{m} = 0

3. Para j = 1, 2, …, m Hacer

(a) Calcular w_{j} := Av_{j}

(b) Para i = 1, …, j Hacer

i. h_{ij} := (w_{j}, v_{i}) ii. w_{j} := w_{j} − h_{ij}v_{i}

(c) FPara

(d) h_{j+1,j} = ∥w_{j}∥_{2}. Si h_{j+1,j} = 0 Hacer m := j y Goto 5

(e) v_{j+1} = w_{j}/h_{j+1,j}

4. FPara

5. Calcular y_{m} el minimizador de ∥βe_{1} − H̄_{m}y∥_{2} y x_{m} = x_{0} + V_{m}y_{m}.

Una segunda manera es derivar el Algoritmo GMRES utilizando la ecuacion 2.6 con W_{m} = AV_{m}.

Variante del Método GMRES: Reinicio

Similarmente al algoritmo FOM presentado anteriormente, el algoritmo GM-RES se hace impráctico cuando m es grande a causa del crecimiento de los requerimientos de memoria, entre otros. Existen algunas variantes del GMRES para afrontar dicho problema. Una de las opciones de la de establecer una versión del algoritmo GMRES(m) con reinicio, en la que se predetermina un valor m máximo para el número de iteraciones y se reinicia el algoritmo GMRES tomando como solución x_{0} inicial igual al valor x_{m} dado por el ciclo interno iterativo anterior. A pesar de que este algoritmo ha sido exitoso en muchos casos, es bien conocido que la versión de GMRES(m) no garantiza la convergencia para m < n.

Precondicionamiento

3.1 Introducción

La falta de robustez es una debilidad ampliamente reconocida de los resolvedores iterativos, en relación con los resolvedores directos. Este problema limita la aceptación de los métodos iterativos en aplicaciones industriales, a pesar de su conveniencia para sistemas lineales muy grandes. Tanto la eficiencia como la robustez de las técnicas iterativas pueden ser mejoradas por el uso de precondicionadores. El precondicionamiento es simplemente una manera de transformar el problema lineal original en uno que tenga la misma solución, pero que sea mejor de resolver iterativamente. Encontrar un buen precondicionador para resolver un problema lineal disperso dado es una combinación de arte y ciencia. Los resultados teóricos son raros y algunos métodos trabajan sorpresivamente bien y otros no cubren las expectativas. Un precondicionador puede ser definido como un resolvedor aproximado auxiliar, el cual es combinado con otra técnica de iteración externa, típicamente una de las iteraciones de subespacios de Krylov. Los precondicionadores son frecuentemente construidos de la matriz de coeficientes original del problema, pero también frecuentemente son derivados del conocimiento del problema físico original. Hablando a grosso modo, un precondicionador es cualquier forma de modificación, implicita o explicita, al sistema lineal original que haga más ”fácil”su resolución. Por ejemplo, el escalamiento de todas las filas del sistema lineal para hacer que los elementos de la diagonal sean iguales a uno es una manera explícita de precondicionamiento. El sistema resultante puede ser resuelto por un método basado en subespacios de Krylov y debería requerir menos pasos para converger que el sistema original (aunque esto puede no estar garantizado). Como otro ejemplo, resolver el sistema lineal,

M^{−1}Ax = M^{−1}b,

donde M^{−1} es alguna función complicada que envuelva el cálculo de transformadas FFT, cálculos integrales y solución de sistemas lineales auxiliares podría ser otra forma de precondicionamiento. Aqui, no es deseable que la matriz M y M^{−1}A sean calculadas explícitamente. Afortunadamente, los procesos iterativos casi nunca requieren operar con esos valores directamente. En la práctica, la operación de precondicionamiento M^{−1} se aplica a un vector arbirtrario y dicha operación debe procurarse que no sea muy costosa. Una de las maneras más simples de definir un precondicionador es la de desarrollar una factorización incompleta de la matriz original A. Esto implica una descomposición de la forma A = LU − R donde L y U tienen la misma estructura de no ceros que las triangulares inferior y superior de A respectivamente, y R es el residual o error de la factorización. Esta factorización incompleta, conocida como ILU(0) es relativamente fácil y barata de calcular. Por otra parte, frecuentemente este precondicionamiento no es suficiente y esto redunda en que el método iterativo emplee muchas iteraciones para converger. Para remediar esto, varias alternativas de factorización incompleta han sido desarrolladas que permitan más elementos de ”relleno”en L y U. En general, la factorización más exacta del tipo ILU requiere pocas iteraciones para converger, pero el costo de preprocesamiento para calcular los factores es muy alto. Entonces, existe un compromiso entre costo de precondicionamiento y el costo de las iteraciones de los resolvedores iterativos, con una inclinación hacia la robustez de los métodos, especialmente cuando se trata de resolver varios sistemas con la misma matriz, en que el precondicionamiento se usa en el mismo sentido de las factorizaciones LU. Existen otra serie de consideraciones a ser tomadas a la hora de precondicionar un sistema, y es que la elección del precondicionador depende del método iterativo que va a ser utilizado y muchas veces de las características de la matriz. Por ejemplo, si el método de gradiente conjugado va a ser utilizado para resolver un problema con una matriz simétrica, el precondicionamiento debe preservar la simetría de la matriz. También la elección del precondicionador depende del esquema de almacenamiento utilizado para la matriz, sobretodo cuando se trata de resolución en ambientes distribuidos, y de la distribución de la data entre los diferentes procesadores. Toda esta serie de consideraciones dificultan la elección del esquema de precondicionamiento en paquetes de uso general. La convergencia de los métodos iterativos depende generalmente de las propiedades espectrales de la matriz, es por eso que los diferentes esquemas de precondicionamiento buscan mejorar dichas características.

3.2 Operaciones de Precondicionamiento

Las tres operaciones de precondicionamiento más comunes en los métodos iterativos son las siguientes:

3.2.1 Precondicionamiento por la Izquierda

El precondicionamiento por la izquierda consiste en sustituir la resolución del sistema lineal Ax = b por la resolución de un nuevo sistema precondicionado:

M^{−1}Ax = M^{−1}b,  (3.1)

que consiste en multiplicar por la izquierda a ambos lados de la ecuación lineal por la matriz M^{−1}. Veamos como ejemplo en que consiste la aplicación de la operación de precondicionamiento por la izquierda al método GMRES, dando origen a una versión precondicionada de dicho método. Algoritmo: GMRES con Precondicionamiento por la Izquierda.

1. Calcular r_{0} = M^{−1}(b − Ax_{0}), β = ∥r_{0}∥_{2}, y v_{1} = r_{0}/β

2. Definir la matriz (m+1)×m H̄_{m} = {h_{ij}}_{i≤i≤m+1,1≤j≤m}. Hacer H̄_{m} = 0

3. Para j = 1, 2, …, m Hacer

(a) Calcular w_{j} := M^{−1}Av_{j}

(b) Para i = 1, …, j Hacer

i. h_{ij} := (w_{j}, v_{i}) ii. w_{j} := w_{j} − h_{ij}v_{i}

(c) FPara

(d) h_{j+1,j} = ∥w_{j}∥_{2}. Si h_{j+1,j} = 0 Hacer m := j y Goto 5

(e) v_{j+1} = w_{j}/h_{j+1,j}

4. FPara

5. Calcular y_{m} el minimizador de ∥βe_{1} − H̄_{m}y∥_{2} y x_{m} = x_{0} + V_{m}y_{m}.

En este caso el ciclo de Arnoldi construye una base ortogonal del subespacio de Krylov precondicionado por la izquierda

Span{r_{0}, M^{−1}Ar_{0}, …, (M^{−1}A)^{m−1}r_{0}}.

El ciclo de Arnoldi utiliza un proceso de Gram-Schmidt modificado, en el cual el nuevo vector a ser ortogonalizado es obtenido del vector previo en el proceso. Todos los vectores residuales y sus normas que son calculados por el algoritmo corresponden a los residuales precondicionados, llamados z_{m} = M^{−1}(b − Ax_{m}), en vez del residual original (sin precondicionar) b − Ax_{m}. por lo tanto, el acceso a los residuales sin precondicionar no es fácil, dado que hay que calcularlos explicitamente, por ejemplo, multiplicando los residuales precondicionados por M. Esto puede causar algunas dificultades si el criterio de parada se basa en los residuales sin precondicionar.

3.2.2 Precondicionamiento por la Derecha

El precondicionamiento por la derecha se basa en resolver:

AM^{−1}u = b, u = Mx.  (3.2)

En el GMRES y en los otros métodos de proyección, por ejemplo, es posible realizar este proceso sin necesidad de invocar explícitamente la variable u. Una vez que el residual inicial b−Ax_{0} = b−AM^{−1}u_{0} es calculado, todos los vectores subsecuentes del espacio de Krylov pueden ser obtenidos sin hacer referencia a las variables u. Notese que u_{0} no es necesario calcularlo. El residual inicial para el sistema precondicionado puede ser calculado desde r_{0} = b − Ax_{0}, el cual es el mismo que b − AM^{−1}u_{0}. En la práctica, es usual disponer de x_{0} y no de u_{0}. Al final, la solución aproximada en la variable u al sistema 3.2 es dada por,

u_{m} = u_{0} + ∑_{i=1}^{m} v_{i}η_{i}

con u_{0} = Mx_{0}. Multiplicando entonces por M^{−1} se obtiene la aproximación deseada en términos de la variable x.

x_{m} = x_{0} + M^{−1} [∑_{i=1}^{m} v_{i}η_{i}]

Por lo tanto, la operación de precondicionamiento es necesaria al final del ciclo exterior en vez de al principio como en el caso del precondicionamiento a la izquierda. Como ejemplo, a continuación daremos la versión del GMRES con precondicionamiento a la derecha: Algoritmo: GMRES con Precondicionamiento por la Derecha

1. Calcular r_{0} = b − Ax_{0}, β = ∥r_{0}∥_{2}, y v_{1} = r_{0}/β

2. Definir la matriz (m+1)×m H̄_{m} = {h_{ij}}_{i≤i≤m+1,1≤j≤m}. Hacer H̄_{m} = 0

3. Para j = 1, 2, …, m Hacer

(a) Calcular w_{j} := AM^{−1}v_{j}

(b) Para i = 1, …, j Hacer

i. h_{ij} := (w_{j}, v_{i}) ii. w_{j} := w_{j} − h_{ij}v_{i}

(c) FPara

(d) h_{j+1,j} = ∥w_{j}∥_{2}. Si h_{j+1,j} = 0 Hacer m := j y Goto 5

(e) v_{j+1} = w_{j}/h_{j+1,j}

4. FPara

5. Calcular y_{m} el minimizador de ∥βe_{1} − H̄_{m}y∥_{2} y x_{m} = x_{0} + M^{−1}V_{m}y_{m}.

Esta vez, el ciclo de Arnoldi construye una base ortogonal del subespacio de Krylov precondicionado por la derecha

Span{r_{0}, AM^{−1}r_{0}, …, (AM^{−1})^{m−1}r_{0}}

nótese que la norma del residual es ahora relativa al sistema inicial b−Ax_{m}, dado que el algoritmo obtiene el residual b−Ax_{m} = b−AM^{−1}u_{m} implícitamente. Esta es una diferencia esencial con respecto al precondicionamiento por la izquierda.

3.2.3 Precondicionamiento SPLIT o dividido

En muchos casos, M es el resultado de una factorización de la forma:

M = LU.

Entonces, una opción para el precondicionamiento es utilizar el sistema precondicionado dividido

L^{−1}AU^{−1}u = L^{−1}b, x = U^{−1}u.

En esta situación, es claro que es necesario, para los métodos como el GMRES, operar el residual inicial por L^{−1} al inicio del algoritmo y por U^{−1} en la combinación lineal V_{m}y_{m} para formar la solución aproximada. La norma del residual disponible es L^{−1}(b − Ax_{m}). Una pregunta a responderse es la diferencia entre las opciones de precondicionamiento a la derecha, a la izquierda y dividido. El hecho de disponer de diferentes versiones de los residuales en cada caso puede afectar los criterios de parada y puede causar que el algoritmo se detenga prematuramente o con retardo. Esto es particularmente dañino en el caso de que M sea mal condicionado. El grado de simetría, y por ende, el rendimiento, pueden ser afectados también por la manera de aplicar el precondicionador. Por ejemplo, un precondicionador dividido puede ser mucho mejor si A es cercana a ser simétrica.Fuera de esas dos situaciones, hay poca diferencia entre las tres opciones. Al comparar dichas opciones de precondicionamiento, una primera observación es que el espectro de los tres operadores asociados M^{−1}A, AM^{−1} y L^{−1}AU^{−1} son idénticos. Por lo tanto, en principio se espera que la convergencia sea similar, a pesar de que como es sabido, los autovalores no gobiernan siempre la convergencia. En la mayoría de las situaciones prácticas, la diferencia de convergencia entre las opciones de precondicionamiento no es significante, excepto cuando M está mal condicionada.

3.2.4 Variantes Flexibles

En la discusión sobre las técnicas de precondicionamiento asumió implícitamente que la matriz de precondicionamiento M era fija, es decir, no cambia entre las iteraciones. Sin embargo, en algunos casos, dicha matriz M no está disponible. Por ejemplo, la operación M^{−1}x es el resultado de algun cálculo no especificado, posiblemente un proceso iterativo. En tales casos, podría esperarse que M^{−1} no sea un operador constante. Los procesos iterativos precondicionados vistos podrían no converger si M no es constante. Existen algunas variantes de los procesos iterativos hallados en la literatura que tratan con precondicionadores variables, es decir, que el precondicionador puede variar paso a paso. Tales procedimientos iterativos son llamados iteraciones flexibles.

3.3 Algunas Técnicas de Precondicionamiento

En esta seccion resumiremos a manera de ejemplo algunas de las técnicas más importantes existentes para conseguir una matriz de precondicionamiento.

3.3.1 Jacobi, SOR, SSOR

Una iteración de punto fijo para resolver el sistema lineal:

Ax = b,

toma la forma general,

x_{k+1} = M^{−1}Nx_{k} + M^{−1}b,  (3.3)

donde M y N realizan una partición de A en,

A = M − N.  (3.4)

La iteración anterior es de la forma

x_{k+1} = Gx_{k} + f,  (3.5)

donde f = M^{−1}b y,

G = M^{−1}N = M^{−1}(M − A),  (3.6)
= I − M^{−1}A.  (3.7)

Para los métodos de Jacobi y Gauss-Seidel ha sido demostrado que:

G_{JA}(A) = I − D^{−1}A,  (3.8)
G_{GS}(A) = I − (D − E)^{−1}A.  (3.9)

donde A = D−E −F es una partición de A en donde D corresponde a la matriz diagonal de A, −E a la triangular inferior y −F a la triangular superior. La iteración 3.5 intenta resolver

(I − G)x = f,  (3.10)

lo cual, a causa de la expresión 3.7 para G, puede ser reescrita como,

M^{−1}Ax = M^{−1}b.  (3.11)

El sistema anterior es el sistema precondicionado asociado con la partición A = M − N y la iteración 3.5 no es más que una iteración de punto fijo sobre el sistema precondicionado. Similarmente, un método basado en subespacios de Krylov, por ejemplo, GMRES, puede ser utilizado para resolver 3.11, dando lugar a una versión precondicionada del método basado en el subespacio de Krylov, por ejemplo, GMRES precondicionado. En teoría, cualquier partición en la cual M sea no singular puede ser utilizada. Idealmente, M debería ser ”cercana^{\underline{a}} A en algún sentido. Notese que un sistema lineal con la matriz M debe ser resuelto en cada iteración del procedimiento iterativo. por lo tanto, un requerimiento práctico es que esa solución este disponible a bajo costo. El precondicionador SSOR está definido por:

M_{SSOR} = (D − ωE)D^{−1}(D − ωF).

Tomando ω = 1 obtenemos la iteración de Gauss-Seidel Simétrico (SGS),

M_{SGS} = (D − E)D^{−1}(D − F).  (3.12)

Una observación interesante es que D − E es la parte inferior de A, incluyendo la diagonal, y D − F es, similarmente, la parte superior de A. Por lo tanto

M_{SGS} = LU,

con

L ≡ (D − E)D^{−1} = I − ED^{−1}, U = D − F.

La matriz L es triangular superior unitaria y U es triangular superior. Con respecto a la implementación de la operación de precondicionamiento, para calcular w = M_{SGS}^{−1}x, se procede como sigue:

1. Resolver (I − ED^{−1})z = x,

2. Resolver (D − F)w = z.

Como vimos anteriormente, las matrices de precondicionamiento de SSOR o SGS son de la forma M = LU, donde L y U tienen el mismo patron que la parte-L y la parte-U de A, respectivamente. Aqui, la parte-L significa la parte triangular inferior y la parte-U es la parte triangular superior. Si se calcula la matriz error A − LU, entonces para SGS, por ejemplo, podemos encontrar,

A − LU = D − E − F − (I − ED^{−1})(D − F) = −ED^{−1}F.

Supongamos que restringimos a L a tener la misma estructura que la parte-L de A y a U a tener la misma estructura de la parte-U de A, la pregunta ahora es si es posible o no encontrar L y U tales que tengan un error más pequeño en algún sentido que el anterior. Podríamos, por ejemplo, intentar encontrar una factorización incompleta en la cual la matriz residual A − LU tenga elementos iguales a cero en las mismas posiciones en donde A los tiene. Esto es posible en general mediante el proceso de factorización incompleta ILU(0) que discutiremos a continuación.

3.3.2 Factorización Incompleta (ILU)

Consideremos una matriz dispersa general A cuyos elementos son a_{i,j}, i, j = 1, …, n. Una proceso general de Factorización LU Incompleta (ILU) calcula una matriz triangular inferior dispersa L y una matriz triangular superior dispersa U tales que la matriz residual R = LU − A satisfaga ciertas restricciones, tales como tener valores iguales a cero en ciertas posiciones. Un algoritmo general para construir una factorización LU incompleta puede ser derivado desarrollando la eliminación gaussiana y eliminando algunos elementos en posiciones predeterminadas fuera de la diagonal. Los elementos a ser eliminados podrían ser especificados a priori, mediante un patron determinado, asi como tambien pueden ser elegidos dinámicamente. Supongamos el primer caso en donde tenemos un patron de ceros dado por el conjunto P

P ⊂ {(i, j)/i ≠ j; i ≤ i, j ≤ n}.  (3.13)

Fijese que la única restricción al patrón de ceros es excluir los elementos de la diagonal. Dado este patrón, una factorización LU incompleta ILU_{P}, puede ser calculada como sigue: Algotitmo: ILU con un patron general estático

1. Para k = 1 … n − 1 Hacer

(a) Para i = k + 1, n y si (i, k) ∉ P Hacer

i. a_{ik} := a_{ik}/a_{kk}

ii. Para j = k + 1, …, n y (i, j) ∉ P Hacer

iii. a_{ij} = a_{ij} − a_{ik} ∗ a_{kj} iv. FPara

(b) FPara

2. FPara

ILU sin relleno o ILU(0)

La técnica de factorización incompleta sin rellenos, denotada por ILU(0), consiste en tomar el patron de ceros P precisamente como el patron de ceros de A He aqui el algoritmo para realizar dicha factorización incompleta sin rellenos. Algotitmo: ILU(0)

1. Para i = 2 … n Hacer

(a) Para k = 1, …, i − 1 y para (i, k)NZ(A) Hacer

i. Calcule a_{ik} := a_{ik}/a_{kk} ii. Para j = k + 1, …, n y para (i, j) ∉ NZ(A) Hacer iii. a_{ij} = a_{ij} − a_{ik} ∗ a_{kj} iv. FPara

(b) FPara

2. FPara

Aqui NZ(A) representa el patron de elementos no cero de la matriz A.

ILU(p)

La precisión de la factorización incompleta ILU(0) puede ser insuficiente para obtener una rata de convergencia adecuada. Es posible obtener factorizaciones incompletas ILU más precisas y frecuentemente más eficientes. Dichas factorizaciones más precisas difieren de ILU(0) por permitir algunos rellenos. Una manera de elegir cuales rellenos van a ser permitidos, es mediante el cálculo del Nivel de Relleno de un elemento a_{ij} de la matriz A y no permitir rellenos con un nivel de relleno superior a un valor p determinado. La determinación del nivel de relleno de un elemento a_{ij} se realiza según la siguiente definición:

Definición 3.1 El nivel inicial de relleno de un elemento a_{ij} de una matriz dispersa A se define por

lev_{ij} = \{ \begin{matrix}0 & si a_{ij} ≠ 0 o i = j \\ ∞ & en caso contrario \end{matrix}

Cada vez que hay una actualización del elemento a_{ij}, su nivel debe ser actualizado por:

lev_{ij} = min{lev_{ij}, lev_{ik} + lev_{kj} + 1}.  (3.14)

Obsérvese que el nivel de relleno de un elemento nunca se incrementa durante la eliminación. Si a_{ij} ≠ 0 en la matriz original A, entonces este tendrá nivel de relleno igual a cero durante todo el proceso de eliminación. Se provee entonces una manera natural de descartar elementos. En ILU(p), todos los elementos de relleno cuyo nivel de relleno no excede p son guardados y por lo tanto el patrón de ceros de ILU(p) viene dado por:

P_{p} = {(i, j)|lev_{ij} > p},

donde lev_{ij} es el nivel de relleno despues de todas las actualizaciones. El caso de p = 0 coincide por supuesto con la factorización ILU(0) y es por tanto consistente con la definición anterior. En las implementaciones prácticas de la factorización ILU(p) es común separar la fase simbólica (donde la estructura de los factores L y U son determinados) de la factorización numérica, cuando los valores numéricos son calculados. la variante que presentaremos aqui no separa esas dos fases. Algotitmo: ILU(p)

1. Para a_{ij} ≠ 0 definir lev_{ij} = 0 FPara

2. Para i = 2 … n Hacer

(a) Para k = 1, …, i − 1 y para lev(a_{ik}) ≤ p Hacer

i. Calcule a_{ik} := a_{ik}/a_{kk} ii. Calcule a_{i∗} = a_{i∗} − a_{ik}a_{k∗} iii. Actualice el nivel de relleno segun la ecuación 3.14

(b) FPara

(c) Reemplace cualquier elemento de la fila i con lev_{ij} > 0 por cero

3. FPara

3.3.3 Inversas Aproximadas

Una manera de precondicionamiento de un sistema lineal, para acelerar la convergencia de un método iterativo consiste en intentar, al menor costo, calcular una inversa aproximada M de A. Dicha inversa aproximada puede ser una inversa aproximada a la izquierda, es decir MA ≈ I o de una inversa aproximada a la derecha, es decir AM ≈ I dentro de un sentido a definir. Dentro de la óptica de un precondicionador paralelo, es interesante obtener un precondicionador en donde el producto matriz-vector asociado My sea fácil a paralelizar y poco costoso en términos de número de operaciones (+ y *). Es también interesante elegir una matriz M dispersa, a pesar de que una inversa dispersa de A no existe siempre se encuentra generalmente que bastantes coeficientes de A^{−1} son despreciables. Por ejemplo, si el problema proviene de una discretización de un problema en derivadas parciales, es interesante buscar una inversa aproximada dispersa. Los precondicionadores polinomiales (M es de la forma M = ρ(A)) son implícitamente paralelos, pero no son más eficientes para la convergencia de un método iterativo que las factorizaciones incompletas LU. A continuación vamos a interesarnos en una manera diferente de hallar inversas aproximadas: se busca la minimización de ∥AM − I∥ para un grafo fijo de M y una norma que se deberá definir.

Presentación del Método SPAI

La elección de la norma para la minimización de ∥AM − I∥ se realizará sobre la norma de Frobenius, la cual permitirá dar algunos grados de libertad dentro del cálculo de la matriz M. Además, dicha elección permitirá obtener un paralelismo implícito debido a que las columnas m_{k} de M pueden ser calculadas de manera independientes una de las otras. En efecto, dado

∥AM − I∥_{F}^{2} = ∑_{k=1}^{n} ∥(AM − I)e_{k}∥_{2}^{2}  (3.15)

la solución de la minimización de la norma de Frobenius se transforma en n problemas de mínimos cuadrados.

min_{m_{k}} ∥Am_{k} − e_{k}∥_{2}, k = 1, …, n.  (3.16)

donde e_{k} = (0, …, 0, 1, 0, …, 0)^{T}. De hecho, los n problemas de mínimos cuadrados se resuelven rápidamente explotando el hecho de que M es dispersa quedando n problemas pequeños de mínimos cuadrados. Sin embargo, fijando a priori el grafo de M, la dificultad viene de elegir dicho grafo para obtener una inversa aproximada dispersa eficaz, dado que con un grafo elegido arbitrariamente, la solución obtenida no resulta siempre en un precondicionador eficaz. Existen diversas heurísticas en la bibliografía que permiten capturar los elementos dominantes de A^{−1} automáticamente. En dichas heurísticas, el algoritmo comienza con un grafo fijo para M, por ejemplo, una estructura diagonal. Se rellena progresivamente la matriz M hasta que un cierto criterio de parada sea satisfecho.

Método de Descomposición de Dominios

4.1 Introducción

Los métodos de descomposición de dominios se refieren a una colección de técnicas las cuales utilizan el principio de ”divide y vencerás”. Tales métodos fueron desarrollados inicialmente para resolver ecuaciones diferenciales parciales sobre regiones en dos o tres dimensiones. Sin embargo, principios similares han sido explotados en otros contextos de la ciencia y la ingeniería. En este capitulo se verá una de tales técnicas desde el punto de vista puramente de álgebra lineal. Los métodos de descomposición de dominios son probablemente uno de los más promisiorios en cierto tipo de problemas, sobre todo dado el progreso de la tecnología multiprocesamiento, la cual es grandemente aprovechada.

4.2 Descomposición de un Problema de S.E.L.

Supongamos que tenemos el sistema de ecuaciones lineales Ax = f con una matriz A que se reordenó de la siguiente forma:

A = ( \begin{matrix}A_{11} & 0 & A_{13} \\ 0 & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{matrix} ),  (4.1)

donde las dimensiones de las submatrices A_{ij} son n_{i} × n_{j}, i, j = 1, 2, 3. Suponemos que en este reordenamiento el conjunto de índices de A tiene el siguiente orden Υ = (Ω_{1}, Ω_{2}, Γ_{3}), con dim(Ω_{1}) = n_{1}, dim(Ω_{2}) = n_{2} y dim(Γ_{3}) = n_{3}. Llamaremos a Ω_{i} el subdominio i : i = 1, 2 y a Γ_{3} la interfaz entre los subdominios 1, 2. Supongamos también que los vectores del lado derecho y el vector solucion han sido reordenados de la misma manera y divididos como:

f = ( \begin{matrix}f_{1} \\ f_{2} \\ f_{3} \end{matrix} ),  (4.2)

y

x = ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ).  (4.3)

Definimos

A_{i} = ( \begin{matrix}A_{ii} & A_{i3} \\ A_{3i} & A_{33}^{(i)} \end{matrix} ) i = 1, 2,  (4.4)

de tal manera que A_{33}^{(1)} + A_{33}^{(2)} = A_{33} y llamaremos A_{i} la matriz local sobre Ω_{i}. Igualmente definimos los vectores del lado derecho locales como:

F_{i} = ( \begin{matrix}f_{i} \\ f_{3}^{(i)} \end{matrix} ) i = 1, 2,  (4.5)

de manera que f_{3} = f_{3}^{(1)} + f_{3}^{(2)} y los vectores solución locales como:

X_{i} = ( \begin{matrix}x_{i} \\ x_{3}^{(i)} \end{matrix} ) i = 1, 2,  (4.6)

de manera que x_{3} = x_{3}^{(1)} + x_{3}^{(2)}. La resolución del sistema lineal Ax = f es entonces equivalente a la resolución del sistema de ecuaciones lineales siguiente:

\{ \begin{matrix}( \begin{matrix}A_{11} & A_{13} \\ A_{31} & A_{33}^{(1)} \end{matrix} )( \begin{matrix}x_{1} \\ x_{3} \end{matrix} ) = ( \begin{matrix}f_{1} \\ f_{3} − A_{32}x_{2} − A_{33}^{(2)}x_{3} \end{matrix} ) \\ ( \begin{matrix}A_{22} & A_{23} \\ A_{32} & A_{33}^{(2)} \end{matrix} )( \begin{matrix}x_{2} \\ x_{3} \end{matrix} ) = ( \begin{matrix}f_{2} \\ f_{3} − A_{31}x_{1} − A_{33}^{(1)}x_{3} \end{matrix} ) \end{matrix}  (4.7)

De la relación A_{31}x_{1} + A_{33}^{(1)}x_{3} + A_{32}x_{2} + A_{33}^{(2)}x_{3} = f_{3} y de f_{3}^{(1)} + f_{3}^{(2)} = f_{3}, resulta la relación:

f_{3}^{(2)} − A_{31}x_{1} − A_{33}^{(1)}x_{3} = −(f_{3}^{(1)} − A_{32}x_{2} − A_{33}^{(2)}x_{3}) ≡^{def} λ

. Hemos definido entonces el término λ llamado multiplicador de Lagrange. Es posible probar entonces que para toda descomposición de f_{3} de la forma f_{3} = f_{3}^{(1)} + f_{3}^{(2)} existe un solo λ tal que:

Ax = f ⟺ \{ \begin{matrix}( \begin{matrix}A_{11} & A_{13} \\ A_{31} & A_{33}^{(1)} \end{matrix} )( \begin{matrix}x_{1} \\ x_{3}^{(1)} \end{matrix} ) = ( \begin{matrix}f_{1} \\ f_{3}^{(1)} + λ \end{matrix} ) \\ ( \begin{matrix}A_{22} & A_{23} \\ A_{32} & A_{33}^{(2)} \end{matrix} )( \begin{matrix}x_{2} \\ x_{3}^{(2)} \end{matrix} ) = ( \begin{matrix}f_{2} \\ f_{3}^{(2)} − λ \end{matrix} ) \\ x_{3}^{(1)} − x_{3}^{(2)} = 0 \end{matrix}  (4.8)

y el sistema lineal 4.7 es por lo tanto equivalente a 4.8. Podemos reescribir 4.8 de una manera más compacta como:

\{ \begin{matrix}A_{1}X_{1} = F_{1} + ( \begin{matrix}0 \\ λ \end{matrix} ) \\ A_{2}X_{2} = F_{2} + ( \begin{matrix}0 \\ −λ \end{matrix} ) \\ (0, I_{Γ_{3}})X_{1} − (0, I_{Γ_{3}})X_{2} = 0 \end{matrix}  (4.9)

con la notación que habíamos definido anteriormente,

A_{i} = ( \begin{matrix}A_{ii} & A_{i3} \\ A_{3i} & A_{33}^{(i)} \end{matrix} ) X_{i} = ( \begin{matrix}x_{i} \\ x_{3}^{(i)} \end{matrix} ) F_{i} = ( \begin{matrix}f_{i} \\ f_{3}^{(i)} \end{matrix} ),

y donde X_{1}, X_{2}, λ son las incógnitas a buscar. Hemos descompuesto el problema original Ax = f en la resolución de dos problemas ”locales”modificados en el vector del lado derecho y de un valor λ denominado multiplicador de Lagrange, además de una condición que deben cumplir estas soluciones locales que se denomina condición de salto nulo en la interfaz. Podemos escribir el problema 4.2 de una manera global como:

( \begin{matrix}A_{1} & 0 & B_{1}^{T} \\ 0 & A_{2} & B_{2}^{T} \\ B_{1} & B_{2} & 0 \end{matrix} )( \begin{matrix}X_{1} \\ X_{2} \\ λ \end{matrix} ) = ( \begin{matrix}F_{1} \\ F_{2} \\ 0 \end{matrix} )  (4.10)

donde B_{1} = −B_{2} = (0, I_{Γ_{3}}) y de manera más compacta este sistema lineal puede escribirse bajo la forma:

( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} )( \begin{matrix}X \\ λ \end{matrix} ) = ( \begin{matrix}F \\ 0 \end{matrix} )  (4.11)

En 4.11,

A = ( \begin{matrix}A_{1} & 0 \\ 0 & A_{2} \end{matrix} )  (4.12)

es una matriz diagonal en bloques y la matriz

B = (B_{1}, B_{2})  (4.13)

permite imponer un salto nulo en la interfaz BX = 0.

4.3 Precondicionamiento para Problemas de Descomposición de Dominios

Se quiere precondicionar el sistema:

A = ( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} )( \begin{matrix}X \\ λ \end{matrix} ) = ( \begin{matrix}F \\ 0 \end{matrix} )  (4.14)

proveniente del problema de descomposición de dominios. Para ello existen diferentes alternativas, tales como: la familia de precondicionadores diagonales en Bloques M_{α} propuestos por Kuznetsov [18] [19] y los precondicionadores por aproximación de la factorización de Gauss.

4.3.1 Precondicionadores en Bloques Diagonales

Estos métodos de precondicionamiento están inspirados en precondicionadores utilizados para resolver problemas de Navier-Stokes que presentan una estructura parecida a los problemas de Descomposición de Dominios. Los diversos precondicionadores presentados a continuación son motivados por el precondicionador presentado por Kunetsov en [18] [19]. El siguiente teorema presenta y justifica la elección del precondicionador propuesto por Kunetsov.

Teorema 4.1 Sea

A = ( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} ),

una matriz invertible, se supone que A es una matriz invertible. Entonces la matriz

MA = ( \begin{matrix}A & 0 \\ 0 & BA^{−1}B^{T} \end{matrix} )^{−1}( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} ),  (4.15)

admite a lo sumo tres valores propios distintos:

Λ_{1} = 1; Λ_{2} = \frac{1 + \sqrt{5}}{2}; Λ_{3} = \frac{1 − \sqrt{5}}{2};

Demostración Sea Λ un valor propio de MA. Sea

U = ( \begin{matrix}X \\ Y \end{matrix} )  (4.16)

un vector propio asociado al valor propio Λ. De la definición de vector propio tenemos que:

MAU = ΛU,

por 4.15 tenemos

( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} )( \begin{matrix}X \\ Y \end{matrix} ) = Λ( \begin{matrix}A & 0 \\ 0 & BA^{−1}B^{T} \end{matrix} )( \begin{matrix}X \\ Y \end{matrix} )  (4.17)

Supongamos Y = 0, entonces de 4.17 viene que AX = ΛAX. Dado que A es invertible y X ≠ 0 se tiene que si Y = 0 entonces Λ = 1 Supongamos ahora que Y ≠ 0. Entonces se tiene que:

\frac{1}{(Λ − 1)}BA^{−1}B^{T} Y = Λ(BA^{−1}B^{T})Y.  (4.18)

Como A es invertible, entonces el complemento de Schur BA^{−1}B^{T} es igualmente invertible. Además Y ≠ 0 por hipótesis. Por lo tanto Λ verifica la ecuación siguiente:

Λ(Λ − 1) = 1,  (4.19)

ecuación que admite dos soluciones reales:

Λ_{2} = \frac{1 + \sqrt{5}}{2}; Λ_{3} = \frac{1 − \sqrt{5}}{2}.

De esa manera, si U es un vector propio de M, el valor propio asociado a U solo puede ser uno de los tres valores propios encontrados en el teorema:

Λ_{1} = 1; Λ_{2} = \frac{1 + \sqrt{5}}{2}; Λ_{3} = \frac{1 − \sqrt{5}}{2};

El cálculo exacto de este precondicionador aplicado a una matriz cualquiera garantiza entonces que algunos métodos iterativos convergen a lo sumo en tres iteraciones [14]. Es posible mejorar este precondicionador de manera significativa sin que sea mucho más costoso. Para ello, se introduce un parámetro α ≠ 0 en el precondicionador M de la siguiente manera:

M_{α} = ( \begin{matrix}A & 0 \\ 0 & αBA^{−1}B^{T} \end{matrix} )^{−1}

Con una demostración similar a la anterior se obtienen los tres valores propios siguientes (no necesariamente reales):

Λ_{1} = 1; Λ_{2} = \frac{1}{2} + \frac{1}{2}\sqrt{1 + \frac{4}{α}}; Λ_{3} = \frac{1}{2} − \frac{1}{2}\sqrt{1 + \frac{4}{α}}.

Eligiendo convenientemente el parámetro α se obtienen diversos precondicionadores. Juvigny [14] ha seleccionado algunos casos interesantes tomando en cuenta los siguientes criterios:

1. Cantidad de autovalores diferentes en el sistema precondicionado.

2. Relación entre el autovalor más grande y el más pequeño en el sistema precondicionado

3. Características de Simetría y Positividad del Sistema Precondicionado.

Eligiendo:

• M_{−4}: Interesante en el sentido que el cálculo exacto del precondicionador permite obtener un espectro de M_{−4}A que posee a lo sumo dos valores propios reales positivos:

Λ_{1} = 1; Λ_{2} = \frac{1}{2};

• M_{\frac{1}{2}} que permite mejorar el número de condición respecto a M_{1} en el caso donde A es simétrica y BA^{−1}B^{T} positivo definida.

• M_{−1} quien es un precondicionador interesante en el caso general ya que puede ser vista como una aproximación de la factorización de Crout por Bloques.

La utilización de estos precondicionadores calculados exactamente requiere el cálculo de la inversa exacta de A y del complemento de Schur S = BA^{−1}B^{T} cuestión que puede resultar tan costosa como la resolución directa del sistema lineal. Sin embargo, una aproximación suficientemente buena de este precondicionador puede garantizar una buena convergencia de los métodos iterativos. Es suficiente para ello intentar aproximar la inversa de A y del complemento de Schur. Para aproximar A^{−1} se puede utilizar el método SPAI. Aún resta aproximar la inversa del complemento de Schur para lo cual hay diferentes técnicas explicadas en [14] y en [8].

4.3.2 Precondicionamiento por una aproximación de la factorizacion de Gauss

Este precondicionador utiliza una aproximación de la factorización de Gauss por bloques de la matriz del problema:

A = ( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} ) = LU = ( \begin{matrix}A & 0 \\ B & −S \end{matrix} )( \begin{matrix}I & A^{−1}B^{T} \\ 0 & I \end{matrix} )  (4.20)

en donde S = BA^{−1}B^{T} es el complemento de Schur. Este precondicionador realiza una aproximación a Ã^{−1} = Ũ^{−1}L̃^{−1} de la siguiente manera:

M_{G} = ( \begin{matrix}I & −Ã^{−1}B^{T} \\ 0 & I \end{matrix} )( \begin{matrix}Ã^{−1} & 0 \\ ˜ (BA^{−1}B^{T})^{−1}BÃ^{−1} & −  ˜ (BA^{−1}B^{T})^{−1} \end{matrix} ) (4.21)

El precondicionador de Gauss es más natural, desde el punto de vista algebraico, y más consistente que los precondicionadores anteriores. Sin embargo hay una gran diferencia en el costo computacional que radica sobretodo en el cálculo del valor ˜ (BA^{−1}B^{T})^{−1}BÃ^{−1}.

4.3.3 Precondicionamiento por Pseudoinversa

El precondicionamiento por pseudoinversa es la proposición principal de este trabajo y constituye una solución intermedia entre los precondicionadores propuestos por Kunetsov y los de aproximación de la factorización de Gauss en bloques. Supongamos el sistema resultante de descomposición de dominios:

A = ( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} )( \begin{matrix}X \\ λ \end{matrix} ) = ( \begin{matrix}F \\ 0 \end{matrix} ),  (4.22)

este sistema es equivalente al siguiente sistema,

( \begin{matrix}A & B^{T} \\ 0 & S \end{matrix} )( \begin{matrix}X \\ λ \end{matrix} ) = ( \begin{matrix}F \\ BA^{−1}F \end{matrix} )  (4.23)

en donde S = BA^{−1}B^{T} es el complemento de Schur asociado. Sin embargo, el sistema 4.23 presenta la ventaja de ser triangular superior en bloques, y por lo tanto es más sencillo de precondicionar. Estudiaremos la familia de precondicionadores dada por:

M_{C} = ( \begin{matrix}A^{−1} & C \\ 0 & S^{−1} \end{matrix} ).  (4.24)

Supongamos que tenemos los valores exactos de A^{−1} y para S^{−1}. Multiplicando a la derecha a ambos lados del sistema 4.23 por M_{C} obtendremos el sistema precondicionado A^{}˜ x = f̃ donde:

A^{}˜ = ( \begin{matrix}A^{−1} & C \\ 0 & S^{−1} \end{matrix} )( \begin{matrix}A & B^{T} \\ 0 & S \end{matrix} ) = ( \begin{matrix}I & (I + CB)A^{−1}B^{T} \\ 0 & I \end{matrix} ),(4.25)
f̃ = ( \begin{matrix}A^{−1} & C \\ 0 & S^{−1} \end{matrix} )( \begin{matrix}F \\ BA^{−1}F \end{matrix} ) = ( \begin{matrix}(I + CB)A^{−1}F \\ S^{−1}(BA^{−1}F) \end{matrix} ).  (4.26)

La matriz A^{}˜ tiene todos los autovalores asociados iguales a 1, lo cual es ideal para su resolución por métodos iterativos. Supongamos I +CB ≠ 0, está suposición es válida en el caso de descomposición de dominios puesto que la matriz B es invertible por la izquierda solo en el caso en que los subdominios no tienen punto de interfaz común. Bajo está suposición y por sustitución regresiva en bloques tenemos el sistema final precondicionado a resolver dado por:

λ = S^{−1}BA^{−1}F,  (4.27)
(I + CB)^{−1}x = A^{−1}(F − B^{T} λ).  (4.28)

Este sistema precondicionado es equivalente en la teoría a resolver:

Sλ = BA^{−1}F,  (4.29)
A(I + CB)^{−1}x = (F − B^{T} λ).  (4.30)

En la práctica al construir el precondicionador, se utilizan aproximaciones de las inversas de A^{−1} y de S^{−1}. Dichas inversas aproximadas pueden utilizarse como precondicionadores en el problema anterior de la manera siguiente: para la resolución del multiplicador de Lagrange precondicionamos por la izquierda con S̃^{−1} ≈ S^{−1}; y para la resolución de x precondicionamos por la izquierda con Ã^{−1} ≈ A^{−1} y tomamos la matriz (I + CB)^{−1} como un precondicionador a la derecha, haciendo:

y = (I + CB)^{−1}x,  (4.31)
x = (I + CB)y.  (4.32)

4.3.4 El Precondicionador de Pseudoinversa desde el punto de vista de Factorizaciones Aproximadas

La matriz original:

( \begin{matrix}A & B^{T} \\ B & 0 \end{matrix} )  (4.33)

admite la siguiente factorización:

A = LU = ( \begin{matrix}I & 0 \\ BA^{−1} & I \end{matrix} )( \begin{matrix}A & B^{T} \\ 0 & −S \end{matrix} ).  (4.34)

El cálculo de la inversa de A se podria realizar mediante el cálculo de la inversa de los factores L y U. El factor L^{−1} es claramente calculable teniendose:

L^{−1} = ( \begin{matrix}I & 0 \\ −BA^{−1} & I \end{matrix} ),  (4.35)

no siendo lo mismo el valor U^{−1}. Dado que U es una matriz triangular en bloques su inversa tiene la forma

U^{−1} = ( \begin{matrix}A^{−1} & X \\ 0 & −S^{−1} \end{matrix} ),  (4.36)

de manera que el valor U^{−1}U = I. Tenemos que

U^{−1}U = ( \begin{matrix}I & (I + XB)A^{−1}B^{T} \\ 0 & I \end{matrix} ),  (4.37)

por lo que el valor X = A^{−1}B^{T} (BA^{−1}B^{T})^{−1} es el valor exacto. Sin embargo calcular este valor exacto es una opción muy costosa por lo que intentaremos luego dar aproximaciones C ≈ X. Utilizando dicha aproximación tendriamos una factorización aproximada de la inversa de A como:

A^{−1} ≈ Ũ^{−1}L^{−1} = ( \begin{matrix}A^{−1} & C \\ 0 & −S^{−1} \end{matrix} )( \begin{matrix}I & 0 \\ −BA^{−1} & I \end{matrix} ).  (4.38)

4.3.5 Escogencias para C

Queremos encontrar una matriz C a bajo costo de tal manera que (I+CB)A^{−1}B^{T} ≈ 0. En este trabajo se propone la alternativa de elegir C = −pseudoinv(B). Como en el caso de los problemas de descomposición de dominios la matriz B tiene una estructura muy particular, se conoce que es posible determinar a muy bajo costo −pseudoinv(B). Además el valor x = (I +CB)y con C = −pseudoinv(B) es un vector que se mantiene exactamente igual a y en las componentes asociadas a nodos que no pertenecen a la interfaz, mientras que promedia los valores obtenidos en las componentes asociadas a nodos de interfaces comunes, lo que quiere decir que es muy sencillo calcular su aplicación en el precondicionador. Analizando esa escogencia, la matriz (I + CB)^{−1} como precondicionador a la derecha tiene el efecto de garantizar de cierta manera la consistencia de las resoluciones independientes hechas en cada subdominio, al forzar, promediando, un salto nulo sobre la interfaz Γ_{3}. En el caso de resolución sobre máquinas paralelas, la aplicación de este precondicionamiento a la derecha no requiere mayores cálculos, sino poca comunicación entre los nodos de la interfaz que se correspondan.

4.3.6 Implementacion práctica del precondicionador por pseudoinversa

Existen muchas maneras de implementar el precondicionamiento por pseudoinversa sobre máquinas paralelas. Una alternativa es la siguiente:

1. Cálculo del multiplicador de Lagrange λ, se requiere:

(a) Una inversa aproximada del Complemento de Schur S̃^{−1} = ˜ (BA^{−1}B^{T})^{−1} como precondicionador

(b) El cálculo de BA^{−1}F resolviendo el sistema Aw = F y luego haciendo Bw

(c) Resolución del Sistema Precondicionado relativo a λ: Sλ = BA^{−1}F.

2. Para el calculo de x, dado λ se requiere:

(a) Cálculo de una inversa aproximada de Ã^{−1} de A como precondicionador a la izquierda.

(b) Resolución del Sistema Ay = F − B^{T} λ con el precondicionador a la izquierda.

(c) Aplicación del Precondicionador a la Derecha x = (I + CB)y promediando los elementos de y pertenecientes al mismo nodo de interfaz.

Dado que la matriz A es una matriz diagonal en bloques, tanto la resolución de los sistemas en donde está involucrada dicha matriz, como el cálculo de la inversa aproximada pueden hacerse paralelamente con cada bloque de manera independiente. Para el cálculo de la inversa aproximada, es conveniente hacerlo de tal manera que se preserven las caracteristicas inherentes a está, como por ejemplo simetría, de tal manera que el sistema precondicionado guarde dichas características y puedan ser aprovechadas. Tanto el cálculo del complemento de Schur, como el de su inversa aproximada y la resolución del sistema precondicionado destinado a calcular λ constituye la parte más delicada del proceso, dado que es alli donde intervienen la mayor parte de las comunicaciones. El tamaño de la interfáz es vital en este punto, puesto que la magnitud de estos sistemas y aproximaciones está en función directa de dicha interfáz. La viabilidad de la aplicación de estas técnicas está directamente relacionada con la manera en que la partición del problema interviene en el tamaño de la interfáz y en el balance de la carga entre los procesadores.

Experimentos

Se llevaron a cabo experimentos con matrices no simétricas de la colleción Harwell Boeing. Las etapas realizadas para cada matriz consistieron en:

1. Hacer la división del conjunto de indices en Ω_{1}, Ω_{2} y Γ_{3}.

2. Crear las matrices A_{ii}, i = 1, 2, 3, tales que A_{12} = 0 y A_{21} = 0.

3. Ensamblar los sistemas de ecuaciones concernientes al problema de descomposición de dominios, calculando las matrices A_{1}, A_{2}, B y el complemento de Shur.

4. Calcular aproximaciones a la inversa del complemento de Shur y a la submatriz diagonal en bloques A.

5. Calcular los diferentes precondicionadores de Kunetsov y M_{C}.

6. Resolver los sistemas lineales precondicionados utilizando el acelerador GMRES.

Los cálculos fueron realizados en MATLAB 5.2, utilizando las rutinas MATLAB del GMRES y de aproximación de inversas. El propósito del experimento fue: comparar la velocidad de convergencia utilizando GMRES, entre los sistemas precondicionados por Kunetsov y los sistemas precondicionados por M_{C} A continuación se anexa una descripción de las matrices utilizadas y los resultados obtenidos.

Conclusiones

El precondicionador por pseudoinversa constituye una solución intermedia viable entre los precondicionadores de Kunetsov y la factorización de Gauss aproximada. Su implementación es mucho menos costosa que las aproximaciones de la factorizacion de Gauss y difiere muy poco del costo de los precondicionadores de Kunetsov, en cambio, presenta propiedades de convergencia mucho más satisfactorias El precondicionador por pseudoinversa tiene características teóricas que lo avalan: tanto por ser una aproximación a bajo costo de una factorización de la inversa de A, como por sus características espectrales. Además, desde el punto de vista intuitivo, es muy natural y consistente. El estudio del precondicionador por pseudoinversa abre nuevas perspectivas en la busqueda de precondicionadores para el problema de descomposición de dominios y para problemas con una estructura similar. Se recomienda continuar con las siguientes lineas:

• Consolidar el precondicionador haciendo experimentos en: Divisiones en más de dos subdominios y con interfaces variadas. Además de su aplicación en problemas no académicos sobre máquinas paralelas y con aproximaciones menos precisas de las inversas de A, del complemento de Schur y de la inversa del Complemento de Schur.

• Investigar la gama de precondicionadores que existen tomando otras aproximaciones similares a la Factorización de la inversa de A.

• Investigar la aplicación de este tipo de ideas a otro tipo de problemas con estructuras similares, por ejemplo, los provenientes de formulaciones de Navier-Stokes.

Bibliografía

[1] Y.Saad. Iterative Methods for Sparse Linear Systems. PSW Publishing Company, 1996.

[2] M. Benzi, M.Tuma. A Sparse Approximate Inverse Preconditioner for Nonsymetric Linear Systems. SIAM J. Sci Compt, Vol. 19, N. 3, pp. 968-994, 1998.

[3] M. Benzi, C. Meyer, M. Tuma. A Sparse Approximate Inverse Preconditioner for the Conjugate Gradient Method. SIAM J. Sci Compt, Vol. 17, N. 5, pp. 1135-1149, 1996.

[4] M. Benzi, W. Joubert, G. Mateescu. Numerical Experiments with Parallel ordering for ILU Preconditioners. Electronic Transactions on Numerical Analysis, Vol. 8, pp. 88-114, 1999

[5] M. Casarin, O. Widlund. A Hierarchical Preconditioner for the Mortar Finite Element Method Electronic Transactions on Numerical Analysis, Vol. 4, pp. 75-88, 1996.

[6] E. Correa, A. Egea, C. León, J. Lista. Una Introducción a MPI: Message Passing Interface, 1994.

[7] E. Chow, Y.Saad. Approximate Inverse Preconditioners via Sparse-Sparse Iterations. SIAM J. Sci Compt, Vol. 19, N. 3, pp. 995-1023, 1998.

[8] E. Chow, Y. Saad. Approximate Inverse Techniques for Block-Partitioned matrices. SIAM J. Sci Compt, Vol. 18, N. 6, pp. 1657-1675, 1997.

[9] E. Chow, Y. Saad. Approximate Inverse Preconditioners for General Sparse Matrices. 1994.

[10] J. Erhel. A Parallel GMRES version For General Sparse Matrices. Electronic Transactions on Numerical Analysis, Vol. 3, pp. 160-176, 1995

[11] J. Hontand. Simulation d’Ecoulements a Grands Nombres de Reynolds par la Methode des Elements avec Joints sur un Reseau de Stations, These de Doctorat de L’Universite ”Pierre et Marie Curie”, Paris, 1997 [12] T. Huckle, M. Grote. A New Approach to parallel Preconditioning with

Sparse Approximate Inverses. 1994

[13] P. Joly. Compléments d’Analyse Numerique Matricielle. Publications du Laboratoire d’Analyse Numerique Paris VI, Nro. A97003, Paris, 1997

[14] X.Juvigny. Resolution de Grands Systemes Lineares sur Machines Massivement Paralleles. These de Doctorat de L’Universite ”Pierre et Marie Curie”, Paris, 1997

[15] N. Gould, J.Scott. Sparse Approximate-Inverse Preconditioners using Norm-Minimization techniques. SIAM J. Sci Compt, Vol. 19, N. 2, pp. 605-625, 1998.

[16] W. Gropp, E. Lusk. Users Guide for mpich, a Portable Implementation of MPI. ANL/MCS-TM-000 Argonne National Laboratory

[17] C. Le Calvez. Acceleration de Methodes de Krylov pour la Resolution de Systemes Lineaires Creux sur Machines Paralleles. These Doctoral de L’Universite de Lille, 1998

[18] Y.Kunetsov. Efficient iterative solvers for elliptic finite element problems on nonmatching grids.Russian Journal of Numerical Analyis Mathematic Modelling, vol. 10. 1995. pp. 187-211

[19] G. Abdoulaev, Y. Kunetsov, O. Pironneau. The numerical implementation of the domain descomposition method with mortar finite elements for a 3d problem, por aparecer.

[20] Y. Saad, A. Malevsky. P-SPARSLIB: A Portable Library of Distributed Memory Sparse Iterative Solvers, 1995.

[21] Y. Saad. Krylov Subspace Methods in Distributed Computing Enviroments

[22] V. Sarin, A. Sameh. An Efficient Iterative Method for the Generalizated Stokes problem SIAM J. Sci Compt, Vol. 19, N. 1, pp. 206-226, 1998.