Construcción de inversas aproximadas tipo sparse basadas en la proyección ortogonal de Frobenius para el precondicionamiento de sistemas de ecuaciones no simétricas

'<U6* BIBLIOTECA UNIVERSITARIA LAS PALMAS DE G. CAIVARIA N." Doeumento JLSJrñS^ N°Copia '^^\2 ^1 Universidad de Las Palmas de Gran Canaria Departamento de Matemáticas Tesis Doctoral Construcción de inversas aproximadas tipo sparse basada en la proyección ortogonal de Frobenius para el precondicionamiento de sistemas de ecuaciones no simétricos Autor: Elizabeth Flórez Vázquez Directores: Gustavo Montero García Luis González Sánchez La Doctorando El Director El Director Fdo.: Ilizabeth Flórez Vázquez o.: Gustavo Montero García Fdo.: Luis González Sánchez Las Palmas de Gran Canaria, Marzo de 2003 A mi hijo Agradecimientos Quisiera agradecer primeramente a Gustavo Montero y Luis González sin cuya ayuda hubiera sido imposible realizar de este trabajo. Además a Eduardo Rodríguez Barrera que pacientemente me ha ayudado a resolver los imprevistos informáticos y a resolver todas las dudas en el tema. A todos los que de una forma u otra han contribuido a la terminación de esta tesis. A todos, muchas gracias. Esta tesis ha sido desarrollada en el marco del proyecto subvencionado por el Ministerio de Ciencia y Tecnología, REN2001-0925-C03-02/CLI, titulado Modelización numérica de transporte de contaminantes en la atmósfera. índice general Introducción 1 1.1. Métodos basados en la miriimización de la norma de Frobenius . 3 1.2. Inversas aproximadas sparse factorizadas 6 1.2.1. Método FSAI 7 1.2.2. Método AINV 8 1.2.3. Método de orlado 8 1.3. Técnicas de ILU inversa 10 1.3.1. Técnicas de factorización ILU inversa basadas en la expansión de Neumann truncada 11 Conceptos básicos 13 2.1. Normas matriciales 13 2.1.1. Notación 13 2.1.2. El Lema de Banach y las inversas aproximadas 15 2.1.3. El radio espectral 17 2.2. Métodos basados en subespacios de Krylov 18 2.2.1. Subespacios de Krylov 18 2.2.2. Métodos basados en los subespacios de Krylov 19 2.2.2.1. Método CGS (Conjugated Gradient Squared) . . 20 2.2.2.2. Método Bi-CGSTAB (Biconjugated Gradient Stabilized) 22 2.2.2.3. Método de QMRCGSTAB 26 2.2.3. Método GMBES (Generalized Minimal Residual) 28 2.3. Precondicionamiento 29 2.3.1. Precondicionador diagonal 30 2.3.2. Precondicionador SSOR 31 2.3.3. Precondicionador ILU(O) 32 2.3.4. Precondicionador Diagonal Óptimo 32 2.3.5. Algvmos métodos de Krylov precondicionados 33 2.3.5.1. Algoritmo Bi-CGSTAB 33 2.3.5.2. Algoritmo QMRCGSTAB 34 2.3.5.3. Algoritmo VGMRES 35 2.4. Esquemas de almacenamiento 36 índice general 2.4.1. Almacenamiento de la matriz del sistema 36 2.5. Reordenación 37 2.6. Precondicionadores explícitos e implícitos 39 2.6.1. Dos m.étodos para el cálculo de aproximadas inversas de matrices en banda por bloques 40 2.6.1.1. Un método explícito de aproximación 40 2.6.1.2. Método implícito de Aproximación 42 2.6.1.3. Matrices definidas positivas 44 2.6.2. Una Clase de Métodos para calcular las inversas aproximadas de matrices 45 2.6.3. Comparación de resultados 52 Inversa aproximada 59 3.1. Resultados teóricos 59 3.1.1. Mejor aproximación y producto escalar de Frobenius ... 61 3.2. <S-Inversa generalizada y complemento ortogonal de Frobenius . 63 3.2.1. 5-Inversa generalizada 63 3.2.2. Complemento ortogonal de Frobenius 64 3.2.2.1. Subespacio de las matrices cuadradas 64 3.2.2.2. Subespacio de las matrices sparse 65 3.2.2.3. - Subespacio de las matrices simétricas 65 3.2.2.4. Subespacio de las matrices hemisimétricas .... 65 3.2.2.5. Subespacio de las matrices que simetrizan aA.. 65 3.2.2.6. Subespacio de las matrices que antisimetrizan a A 66 3.2.2.7. Subespacio de las matrices simétricas que simetrizan a A 66 3.3. Expresiones explícitas para los mejores precondicionadores .... 66 3.4. Aplicación a algunos precondicionadores usuales 69 3.4.1. Precondicionador con patrón de sparsidad dado 69 3.4.2. Precondicionador diagonal 71 3.4.3. Precondicionador simétrico y hemisitnétrico 72 3.4.4. Precondicionador M tal que AM sea simétrica . 74 Valores propios y valores singulares 75 4.1. Menor valor propio y menor valor singular 75 4.2. Análisis de convergencia 78 Inversa aproximada sparse 81 5.1. Cálculo de inversas aproximadas spflrse 81 5.2. Método de la inversa aproximada mejorada 83 5.2.1. Efectividad teórica de la inversa aproximada mejorada . . 85 5.3. Experimentos numéricos 86 índice general xi 6. Efecto de la reordenación 93 6.1. Algoritmo de la inversa aproximada sparse 94 6.2. Algimos comentarios sobre reordenación 94 6.3. Experimentos numéricos 96 7. Conclusiones y líneas futuras 109 Bibliografía 111 índice de figuras 5.1. EstructaTa sparse de la matríz convdifhor 88 5.2. Estructura sparse de la inversa aproximada de la matriz convdifhor con Ck — 0,5 y máx nz(mofc) =50 88 5.3. Estructura sparse de la inversa aproximada de la matriz convdifhor con Sk = 0,05 y máx nz(mok) = 50 89 5.4. Estructura sparse de la inversa aproximada de la matriz convdifhor con Sk = 0,05 y máx nz{mok) = 200 89 5.5. Comportamiento de los prencondicionadores con BiCGSTAB para convdifhor 91 5.6. Comportamiento de los prencondicionadores con BiCGSTAB para isla 91 6.1. Comparación del comportamiento de BiCGSTAB-ILU(O) y BiCGSTAB-SPAI con reordenación para orsreg^l 98 6.2. Comparación del comportamiento de BiCGSTAB-SPAI con reordenación para convdifhor 100 6.3. Patrón de sparsidad de la matriz SPAI(0.3) original para convdifhor. 101 6.4. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con grado mínim.o para convdifhor 101 6.5. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con Reverse Cuthill-Mckee para coni;díf72or 102 6.6. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con mínimo vecino para convdifhor 102 6.7. Comparación del comportamiento de BiCGSTAB-SPAI con reordenación para cnaref. 104 6.8. Patrón de sparszdflíí de la matriz SPAI(0.3) original para CMcre/. . . 105 6.9. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con Grado M'inimo para CMflre/. 105 6.10. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con Reverse Cuthill-Mckee para cuaref. 106 6.11. Patrón de sparsidad de la matriz SPAI(0.3) reordenada con Mínimo Vecino para cuaref. 106 índice de tablas 5.1. Resultados de convergencia para convdifhor con BiCGSTAB precondicionado por la izquierda 87 5.2. Resultados de convergencia para oilgen con GMRES precondicionado por la izquierda 90 5.3. Resultados de convergencia para sherman con QMRCGSTAB precondicionado por la izquierda 90 5.4. Resultados de convergencia para isla con BiCGSTAB precondicionado por la izquierda 90 5.5. Comparación de los resultados de convergencia para pores con métodos de Bóylov e lAI 90 6.1. Resultados de convergencia para orsregl con el orden original y BiCGSTAB precondicionado por la izquierda 96 6.2. Resultados de convergencia para orsregl con Grado Mínimo y BiCGSTAB precondicionado por la izquierda 97 6.3. Resultados de convergencia para orsregl con Reverse CutHill McKee y BiCGSTAB precondicionado por la izquierda 97 6.4. Resultados de convergencia para convdifhor con el orden originial y BiCGSTAB precondicionado por la izquierda 99 6.5. Resultados de convergencia convdifhor con Grado Mínimo y BiCGSTAB precondicionado por la izquierda 99 6.6. Resultados de convergencia para convdifhor con Reverse CutHül McKee y BiCGSTAB precondicionado por la izquierda 99 6.7. Resultados de convergencia para convdifhor con Mínimo vecino y BiCGSTAB precondicionado por la izquierda 100 6.8. Resultados de convergencia para cuaref con ordenación original y BiCGSTAB precondicionado por la izquierda 103 6.9. Resultados de convergencia para cuaref con Grado Mínimo y BiCGSTAB precondicionado por la izquierda 103 6.10. Resultados de convergencia para cuaref con Reverse CutHill McKee y BiCGSTAB precondicionado por la izquierda 103 6.11. Resultados de convergencia para cuaref con Mínimo vecino y BiCGSTAB precondicionado por la izquierda 104 Capítulo 1 Introducción Uno de los problemas más importantes en la ciencia de la computación es el desarrollo de Métodos iterativos paralelizables y eficientes para resolver sistemas de ecuaciones lineales Ax = b con matriz de coeficientes A de orden elevado y sparse. Los Métodos basados en los subespacios de Krylov efectúan productos matriz-vector en cada iteración y pocas operaciones con vectores (producto escalar y actualizaciones de vectores). Estos métodos pueden ser implementados eficientemente en ordenadores de gran capacidad pero necesitan del precondicionamiento para ser efectivos. La mayoría de los precondicionadores de propósitos generales tales como los basados en la factorización incompleta de A, son suficientemente robustos y permiten obtener una buena velocidad de convergencia, pero son altamente secuenciales lo que dificulta su implementación eficiente en ordenadores en paralelo, especialmente para problemas no estructurados. Así, el precondicionamiento es normalmente la mayor dificultad en la resolución de grandes sistem.as sparse. Se ha realizado una gran cantidad de trabajos desde los iiúcios del los procesadores vectoriales y paralelos, enfocados a paralelizar todo lo posible los mejores precondicionadores tales como SSOR y los métodos de factorización incompleta. Como resultado de esto, es posible lograr tina buena ejecución para algimos tipos de problemas con matrices altamente estructuradas como los que se obtienen de la discretización de ecuaciones en derivadas parciales en nnallas regtdares. Por otra parte, todavía es muy difícil resolver eficientemente sistemas Lineales generales con un patrón de sparsidad irregular en ordenadores vectoriales. Otra línea de investigación consiste en el desarrollo alternativo de métodos de precondicionamiento que son paralelizables de forma natural. Dentro de las primeras técnicas de este tipo podemos mencionar los precondicionadores polinomiales, que se basan en aproximar la inversa de la matriz de coeficientes A con un poHnomio de grado pequeño en la matriz. Estos métodos tienen una larga historia, [27,30] pero se poptilarizaron sólo después que los primeros procesadores vectoriales aparecieran en el mercado, [42, 68]. Los precondicionadores polinomiales únicamente requieren productos matriz-vector con A y. Introducción 1.2.2. Método AINV Otro método de calcular una inversa aproximada factorizada es el basado en la biconjugación incompleta, propuesto primeramente en [14]. Este camino denominado método AINV se describe detalladamente en [20] y [22]. El método AI>JV no necesita el conocimiento previo del patrón de sparsidad y puede aplicarse a matrices con patrones de sparsidad generales. La construcción del precondicionador AINV se basa en un algoritmo que calcula dos conjuntos de vectores {-zj^^i y {^Wi}"=i que son ^—biconjugados tal que wfAzj = O si y solo si Í 7^ jf. Dada una matriz no singular A e M", existe una relación estrecha entre el problema de invertir A y él problema de calcular los dos conjuntos de vectores yl—biconjugados {-z¿}f=i y {•«¿'¿HLiSi ^ — {^ii ^2, •••, Zn} es la matriz cuya i-ésima columna es Zj y W = {wi,W2,...,Wn} es la matriz cuya ¿-ésima columna es Wi, entonces pi O --f O W^AZ = V O • • • o pn J donde pi = wfAzi 7^^ 0. Se sigue que W y Z son necesariamente no singulares y ^ ZiWf A-' = zD-'w^ = y ^-^ Pi De aquí que la inversa de A se conoce si los dos conjuntos de vectores A—biconjugados son conocidos. Existen infinitos conjuntos de estos vectores. Las matrices W y Z cuyas columnas son A—biconjugadas pueden calcularse explícitamente por medio de un proceso de biconjugación aplicado a las columnas de cualquier par de matrices no singulares Vl^^°\ Z^°^ € R"^". Una selección conveniente desde el punto de vista computacional es hacer W^^^ = Z^°^ — I, aplicándose el proceso de biconjugación a tina base canónica de vectores. 1.2.3. ]VIétodo de orlado. Un tercer camino que se utiHza para calcular el precondicionador del tipo inversa aproximada f actorizado directamente de la matriz de entrada A se basa en el orlado. En la actuaHdad son posibles muchos esquemas. El método que se describe aquí es una modificación propuesta por Saad en [94], (pag.308). Sea Ak Inversas aproximadas sparse factorizadas la submatriz principal de A de orden k x k . Considérese el siguente esquema de orlado: WT 0\ f Ak Vk \ f Zk Zk\ f Dk O donde Zjt = -ZkDl^Wlvk,Wk = -WkD^'^ZlykySk+i = ak+i+wjAkZk+ykZk + w^VkAquí Wk y Zk son matrices trangulares superiores de orden A;. Por otra parte, Wk,Vk,yky Zk son vectores de orden k, mientras que 6k+i y ctk+i son escalares, siendo a^+i = 0^+1,^+1. Comeiizando en fc = 1, este esquema sugiere un algoritmo obvio para el cálculo de los factores inversos de A (asumiendo que A admite una factorización LU). Cuando este esquema se ejecuta de form.a incompleta, se obtiene una factorización aproximada de A~^. Puede preservarse la sparsidad eliminando elementos de los vectores Wk y Zk una vez calculados, por ejemplo, usando cierta tolerancia al igual que en el proceso AINV. El precondicionador resultante de la factorización de la inversa aproximada sparse se denominará AIB (Approximate Inverse via Bordering), es decir, inversa aproximada sparse mediante orlado. Aparte de xxn producto matriz-vector con Ak, la construcción del precondicionador AIB requiere cuatro productos sparse matriz-vector que involucran a W^, Zfc y sus traspuestas en cada paso k, lo cual constituye la mayor parte del trabajo a realizar. Es importante que estas operaciones se efectúen en modo sparse-sparse. Nótese que los cálculos para los factores inversos Z yW están fuertemente acoplados, en contraste con el algoritmo de biconjugación. Como siempre, si A es simétrica, W = Z y el trabajo se redtice a la mitad. Más aún, si A es sim.étrica y definida positiva, tal como se muestra en [94], en aritmética exacta, 6k > O para todo k. Así, el precondicionador AIB está siempre bien definido en el caso simétrico y definido positivo. En el caso general, se requieren modificaciones en la diagonal para completar el proceso. Las inversas aproximadas sparse factorizadas no tienen otros problemas que limiten la efectividad como ocurre con otros métodos. Al contrario que los métodos explicados en las secciones previas, las inversas aproximadas sparse factorizadas pueden utiHzarse como precondicionadores para el método del gradiente conjugado para resolver problemas simétricos y definidos postitivos. Así, si A es simétrica y definida positiva, entonces Z = W y él precondicionador M = ZD'^Z'^ es simétrico y definido positivo según las entradas de la diagonal D sean todas positivas. Está claro que la no singularidad de M es triviaknente controlable cuando M se expresa de foma factorizada. Siguiendo [32] puede decirse que las formas factorizadas nos dan vina mejor aproximación de A~'^ para la misma capacidad de almacenamiento que las no factorizada porque pueden ser matrices más densas comparando con el número total de entradas no nulas de sus factores. Tal y como muestran los experimentos en [24], esta observación intuitiva casi se traduce en una mejora en la convergencia para el mismo número de entradas no nulas. Más aún, las formas factorizadas son más baratas de calcular y en la mayoría de los casos requieren definir menos para- W Introducción metros que las técnicas mencionadas en las secciones previas. Finalmente, las formas factorizadas son sensibles al reordenamiento de la matriz de coeficientes, propiedad que puede explotarse para reducir el efecto de llenado en los factores inversos y/o mejorar la velocidad de convergencia, ver [22], [14], [28]. Por otra parte, las formas factorizadas tiene un problema. Al ser métodos de factorización (inversa) incompleta pueden fallar debido a interrupciones en el proceso de factorización incompleta como el ILU. Mientras que las estrategias de desviación de la diagonal [80] o las de compensación de la diagonal pueden usarse como protección de los cálculos, no existe garantía de que el precondicionador resultante sea efectivo, especialmente cuando se necesitan grandes y/o numerosas desviaciones. El método FSAI requiere prescribir la estructura de los factores no nulos por adelantado, lo que hace difícil su uso en problemas con patrón de sparsidad generales. Otros métodos como AINV, ofrecen oportunidades limitadas de factorización en la fase de construcción del precondicionador. En su formulación corriente, AINV parece difícil de im^plementar en máquinas de memoria distribuida. Tam^bién el paralelismo en la aplicación de inversas aproximadas factorizadas es algo menos que en las formas no factorizadas, puesto que en las primeras es necesario ejecutar dos productos matrizvector con los factores secuencialmente. 1.3. Técnicas de ILU inversa Muchos autores han propuesto construir precondicionadores del tipo inversas aproxim.adas sparse factorizadas basados en el siguiente proceso en dos etapas: primero una factorización LU incompleta A ^ LÜ que se calcula usando técnicas clásicas y después una aproximaxión de la inversa de los factores L y Ü, ver [2, 40, 99,100, 44]. Existen varias posibilidades para calcular inversas aproximadas de L y t7, y cada tma nos permite obtener diferentes precondicionadores. Asumimos que los factores incom^pletos LyÜ existen. Entonces los factores inversos aproximados pueden calcularse resolviendo de forma inexacta los 2n sistemas lineales triangulares Lxi = e¿, Üyi = 6,: (1 < ¿ < n) Obsérvese que todos estos sistemas lineales pueden resolverse independientemente, por lo que podemos contar con realizar tma buena paralelización al menos en principio. Estos sistemas lineales se resuelven aproximadamente para reducir el tiempo de computación y porque debe conservarse la sparsidad en las columnas de los factores inversos aproximados. Una posibilidad es prefijar los patrones de sparsidad SL y Su para las inversas aproximadas de Z y 17. Las aproximaciones sparse pueden calcularse usando la norma Frobenius. El método adaptativo SPAI puede usarse como alternativa Técnicas de ILU inversa H para invertir aproximadamente L y ¿7, sin necesidad de prefijar los patrones de sparsidad. A partir de algimos experimentos realizados en [44] se concluyó que este camino no es recomendable. Un camino m.ejor consiste en resolver los 2n sistemas triangulares para las columnas de X~^y í/~^ por descenso y remonte respectivamente. Se conserva la sparsidad introduciendo las entradas en los vectores soluciones terüendo en cuenta tanto las posiciones (de forma más general, usando un esquema de nivel de llenado) como la tolerancia. Se han descrito con detalle muchos de estos esquemas [40,100,44]. Algunos autores han propuesto introducir en Xi = L'^EÍ y y- = Ü'^Ci, una vez que se calcula la solución exacta, pero tm esquema más práctico es introducir dtirante el proceso de sustitución mejor que después [44]. Los precondicionadores de esta clase comparten algunas ventajas de los métodos de inversas aproximadas factorizadas pero tienen ciertas desventajas que los precondicionadores descritos en las secciones previas no tienen. Estas desventajas radican en asumir que se ha calculado la factorización ILU. Esto implica que estos métodos no son aplicables si la factorización ILU no existe o si es inestable como es el caso, algimas veces, de problemas altamente no simétricos e indefinidos [35, 46]. Asumir esto también Umita la eficiencia de la paralelización de esta clase de métodos puesto que la fase de construcción del precondicionador no es enteramente paralelizable (calcular una factorización ILU es tm proceso altamente secuencial). Otra desventaja es que el cálculo del precondicionador involucra dos niveles de incompletitud, este es el caso de otros métodos de inversas aproximadas considerados anteriormente. Para algunos problemas, esto puede llevarnos a una degradación significativa en la calidad del precondicionador. Quizás , la presencia de dos niveles de incompletitud hace que estos métodos sean difíciles de aplicar en la práctica debido a la necesidad de seleccionar un gran número de parámetros definidos por el usuario. En general, los métodos de ILU inversa son mucho más difícües de usar que los precondicionadores descritos en las subsecciónes 1.1 y 1.2. 1.3.1. Técnicas de factorización ILU inversa basadas en la expansión de Neumann truncada Las técnicas de factorización ILU inversa basadas en la expansión de Neumann truncada no son métodos de inversa aproximada estrictamente hablando. Pueden considerarse un híbrido de las técnicas de precondicionamiento ILU y polinomial. Sin embargo son sitnüares a otros métodos ya descritos que se basan en una factorización ILU donde los factores incompletos se invierten de forma inexacta aplicando algún tipo de truncamiento. En particular, la sustitución por descenso y remonte se reemplazan por productos matriz-vector con matiices triangulares sparse. Esta idea proviene de van der Vorst [103] y se ha aplicado recientemente a los precondicionadores SSOR que pueden verse 12 Introducción como un tipo de factorización incompleta en Gustafsson y Lindskog [62]. El precondicionador SSOR de Neumann truncado para una matriz simétrica A se define como sigue. Sea A — E-¥D-\- E'^ la descomposición de A en sus partes estrictamente triangular inferior, diagonal y triangular superior. Considérese el precondicionador SSOR, M = [w-^D + E) [w-^D]'^ {w-^D + E^) donde O < w < 2 es un parámetro de relajación. Sea L = wED'^ y D = w~^D, entonces M = (7 + L) ¿> (/ + rf tal que M-^ = (/ + L)-^ D-^ (/ + L)-^ Note que A;=0 Para algunas matrices, por ejemplo, aquellas con diagonal dominante, \\L^\\ disminuye rápidamente cuando k aumenta; y la suma puede aproximarse con un polinomio de grado pequeño en L, esto es/^1 o 2. Por ejemplo, para primer orden, G={IL^) D-'^ (I-L)^ M-^ « A-^ puede considerarse como una inversa aproximada factorizada de A. Debido a que el precondicionador SSOR no requiere cálculos, excepto para la estimación de w, posiblemente, este es lin precondicionador libre desde el punto de vista virtual. También es muy fácil de implementar. Por otra parte, la efectividad de este método se restringe a problemas para los cuales el precondicionador SSOR estándar trabaja eficientemente, lo cual no ocurre en la generalidad de los problemas. Más aún, la expansión de Neumann truncada a un polinomio de grado pequeño puede provocar una degradación seria de la velocidad de convergencia, particularmente en problemas donde la diagonal no es dominante. La idea de la expansión de Neumann truncada puede también aplicarse a la factorización de Cholesky incompleta y, de forma más general, a los precondicionadores de factorización LU incompleta. Capítulo 2 Conceptos básicos 2.1. Normas matriciales 2.1.1. Notación Comenzaremos por revisar algunas ideas del álgebra lineal numérica. Una referencia excelente sobre las ideas básicas tanto del álgebra lineal numérica como de los métodos directos de resolución de sistemas de ecuaciones Lineales puede encontrarse en [71, 64,65]. Escribimos un sistema de ecuaciones Hneales como, Ax = h (2.1) donde A es vina matriz no singular de orden n x n y b G M" es conocido y definiremos la solución exacta de (2.1) como X* = A-^h e M" (2.2) En lo que sigue, denotaremos por x la solución aproximada de la ecuación (2.1) y {xk}k>o será la secuencia de iteraciones. Por otra parte, llamaremos {x)i a la ¿-ésima componente de x, y {xk)i a la i-ésima componente de Xk, aunque usualmente nos referiremos al vector y no a sus componentes. La norma en E" se simbolizará mediante ||• || al igual que la norma matricial inducida. Definición 1 Sea ||-|| una norma en R". La norma matricial inducida de una matriz A de orden nx nse define como, \\A\\ = máx \\Ax\\. Ikll=i Las normas inducidas cumplen ima importante propiedad, \\Ax\\ < \\A\\ \\x\\ (2.3) El número de condición de A con respecto a la norma ||- |i es, K{A) = \\A\\\\A-m (2.4) 14 Conceptos básicos donde K{A) es infinita si A es singular. Si || • || es la norma p N p ii^iip = !Ei('^)^n (2.5) escribiremos el número de condición como Kp. La mayoría de los métodos iterativos finalizan cuando el vector residuo, r = bAx (2.6) es lo suficientemente pequeño. Otro criterio de parada es, M < r {1.7) que puede relacionarse con el error mediante, e = x-x* (2.8) en términos del número de condición. Lema 2 Sean b, x, XQ E M", A una matriz no singular y x* = A~^b, entonces "^" <K{A)P\: (2.9) iieoll \\ro\ Demostración. Como, r = b — Ax = —Ae tenemos, ||e|| = \\A-'Ae\\ < \\A-^\\ \\Ae\\ = ¡¡A'^ \\r\\ llroll = Peoll < \\A\\ ||eo|| De aquí, llell ^ Ik'^lllIHI _ ^(A\M. INII - IIAirMlroll ~ '^^^^^ llroll como se afirmó. • El criterio de parada (2.7) depende de la iteración inicial y, cuando la iteración inicial es buena, puede resultar im trabajo innecesario, sin embargo si la iteración inicial está muy lejos de la solución inicial el resultado puede ser de calidad insuficiente. Por esta razón es preferible concltiir las iteraciones cuando "'''=" < T (2.10) \\b Las dos condiciones (2.7) y (2.10) son las mismas cuando partimos de XQ = O, lo cual se efectúa com.únmente, en especial cuando la iteración lineal se usa como parte de un método no lineal. Normas matriciales 25 2.1.2. El Lema de Banach y las inversas aproximadas. El camino más directo para la solucióii iterativa de un sistema lineal es reescribir la ecuación (2.1) como una iteración lineal de pimto fijo. Una forma es escribir Ax = b como x = {I-A)x + b (2.11) y definir la iteración de Richardson Xk+i = (/ - A)xk + b (2.12) Mostraremos métodos más generales donde {xk} está dado por Xk+i = Mxk + c (2.13) donde M es una matriz cuadrada de orden n denominada matriz de iteración. Los métodos iterativos de esta forma se denominan métodos iterativos estacionarios porque la forma de transición de Xk a Xk+i no depende del desarrollo de la iteración. Los métodos de Krylov, que se discutirán en la sección 2.2, son métodos iterativos no estacionarios. Todos estos restiltados están basados en el siguiente Lema, Lema 3 Si M es una matriz nx n con \\M\\ < 1 entonces I — M es no singular y Demostración. Probaremos que / — M es no singular y que (2.14) se cumple demostrando que oo £M' = (/ - M)-i Las sumas parciales Sk = Í:M' 1=0 forman una sucesión de Cauchy en ]R"^". Nótese que para todo m> k m \\Sk-Sm\\< E 11^11 l=k+l Ahora, ||M'|| < ||M||' porque ||-|| es ima norma matricial inducida por una norma vectorial. De aquí 115. - 5.11 < J^^ iiMii^ = mr^ (i^SF^) - o cuando m,k —> oo. Por tanto, la secuencia Sk converge a S. Como MSk +1 = Sk+i, tendremos MS + I = S j (I — M)S = I. Esto prueba que J - M es no singular y que S = {I — M)~^. Observando que, ||(/-M)-i||<EI|M||' = (l-||M||)-i 1=0 se demuestra (2.14) y completa la prueba • El siguiente corolario es ima consecuencia directa del Lema 3. 16 Conceptos básicos Corolario 4 Si \\M\\ < 1 entonces la iteración (2.13) converge a x = {I - M)~'^c para toda iteración inicial x^. Una consecuencia del corolario 4 es que la iteración de Richardson (2.12) convergerá si ||/ — ^|| < 1. Algunas veces es posible precondicionar el sistema lineal multiplicando ambos miembros de (2.1) por una matriz B BAx = Bh (2.15) de tal forma que se mejore la convergencia del método iterativo. En la sección (2.3) se profundiza en los métodos de precondicionamiento. En el contexto de la iteración de Richardson, las matrices B que permiten aplicar el lema de Banach y su corolario se denominan inversas aproximadas.. Definición S B es una inversa aproximada de Asi\\I — BA\\ < 1. El siguiente teorema se denomina comúnmente Lema de Banach. Teorema 6 Si Ay B son dos matrices de orden nxny B es una inversa aproximada de A. Entonces, Ay B son ambas no singulares y, II ^-111 ^ ll-^ll II p-iil ^ ll^ll n^c^ il^ ll-l-||J-5^||'ll^ ll-l-||J-5^||' ^^-^^^ y .-1 ^^^^\\B\\\\I-BA\\ 11, ^-r.^\\A\\\\I-BA ^11 ^ \-¡I-BA¡^ 11^ - ^11 ^ 1-117-^^11 (^-^^ Demostración. Sea M = I-BA. Por el Lema 3,1-M = I-{I-BA) = BA es no singular. Por tanto, ambas matrices Ay B son no singulares. Por la ecuación (2.14) ll^-.B-.|l = 11(7 - M)-i < ^^^^ = T^jr^m <'•''' como A"'^ = {I — M)~^ B, enconces la inecuación (2.18) implica la primera parte de la ecuación (2.16). La segianda parte se sigue de forma simüar de B~^ = A{I — M)~\ Para completar la prueba se debe tener en cuenta que A-^-B = {IBA) A-\AB-^ = B'^ (/ - BA) y (2.16). • La iteración de Richardson, precondicionada con inversa aproximada tiene la forma, Xfc+i = (J - BA) Xk + Bb (2.19) Si la norma de / — BA es pequeña, entonces, no solo convergerá la iteración rápidamente sino que, según indica el Lema 2, el criterio de parada basado en el residuo precondicionado Bb — BAx reflejará mejor el error real. Este método es una técnica muy efectiva para resolver ecuaciones diferenciales, ecuaciones integrales y problemas relacionados. También pueden interpretarse bajo esta luz los métodos multimalla. Normas matriciales 17^ 2.1.3. El radio espectral El análisis en la sección anterior relaciona la convergencia de la ecuación (2.13) con la norma de la matriz A. Sin embargo, la norma de M puede ser pequeña en algunos tipos de normas y muy grande en otros tipos. Aquí la iteración no está completamente descrita por ||M||. El concepto de radio espectral nos permite hacer una descripción completa. Denotando por a{A) el conjxmto de los autovalores de A. El radio espectral de vma matriz A de orden n x n es p(^)= máx|A|=lím||A'=||^ (2.20) Xe<7{A) fc—•oo El término de la derecha de la segunda igualdad en (2.20) es el límite que resulta al utilizar la caracterización por radicales de convergencia de la serie de potencias Yl ^• El radio espectral de M es independiente de cualquier norma matricial particular de M, de hecho p{A) < \\A\\ (2.21) para cualquier norma matricial inducida. La inecuación (2.21) concluye que el radio espectral es tma cota inferior del conjunto de los valores de las normas matriciales sobre A. Además, se cumple el siguiente teorema. Teorema 7 Sea A una matriz de orden n x n. Entonces para cualquier e > O existe una norma ||-|| en R" íaZ que p{A) > \\A\\ - e (2.22) En otras palabras, el radio espectral es el ínfimo del conjunto todas las normas matriciales. Corolario 8 El menor de los valores singulares de A nunca excede del menor de sus valores propios en módulo. Demostración. Denotemos los autovalores y valores singulares de A en orden decreciente, |Ai| > IA2I > ... > |An| o-i>a2> ... > (7„ como p{A) = mí\\A\\ 11-11 entonces P{A) < \\A\\ En particular, tomando como norma matricial la norma espectral o de Hübert, PÍA) < \\A\\, y |Ai(^)| < \Xi{AA^)\ es decir. 24 Conceptos básicos incorporando la ortogonalidad de (f)j{A)ro respecto de todos los vectores {Á^YrQ, con i < j y teniendo en cuenta que sólo los coeficientes principales de los polinomios son relevantes en el desarrollo del anterior producto, si j( es el principal coeficiente para el polinomio 4>j{A), entonces. Examinando las relaciones de recurrencia para 0j+i(^) y ipj+i{A), los coeficientes principales para estos polinomios fiaeron hallados para satisfacer las relaciones. Por tanto, Pj+i ^ Wj pj+i Pj Oij Pj que nos permite encontrar la siguiente relación para (3j, Pj ^j De forma similar, y por una simple fórmula de recurrencia podemos encontrar Qj como, _ {cPj{A)ro,MA'')ro) "^ {A7rjiA)ro,7rj{ATy,) De la misma manera que anteriormente, en los productos de polinomios, tanto en el numerador como en el denominador, sólo se consideran los términos correpondientes a sus respectivos coeficientes principales y como éstos son idénticos para (pj (A^)ro y TTJ (Á^)ro, podemos escribir. a, = _ (0,-(^)ro,0,(A^)rg) ' {A7rj{A)ro,MAT)r*o) _ {<Pj{A)ro,i;j{A^)r*o) {A7rj{A)ro,i^M''K) _ {i;j{A)^j{A)ro,r*) {A^-¡{A)7rj{A)ro,r*) Y como Pj = ipj{A)-Kj{A)ro, entonces, % = 7T^ (2-55) A partir de la ecuación (2.53), si hacem.os, Sj = rj — CíjApj (2.56) Métodos basados en subespacios de Krylov 25_ el valor óptimo para el parámetiro Wj que interviene en la constinicción del polinomio reductor ipj{A) y que figura en las relaciones de recurrencia del algoritmo lo obtendremos con la condición de minimizar la norma del vector residuo, Tjj^\ ^ Sj iJüj/iiSj \\rj+if = {sj - WjAsj,Sj - WjAsj) = {sj, Sj) - 2wj {sj, Asj) + w] {Asj, Asj) dwj con lo que. -2 {sj, Asj) + 2wj {Asj, Asj) = O w • = ^'^•^J' ^^^ (2.57) La ecuación (2.53) ahora resulta, Tj+i = Sj — WjAsj = Tj — ajApj — WjAsj (2.58) y la relación de recurrencia para el vector solución viene dada por, Xj+i = Xj + ajPj + WjSj (2.59) Una vez expresados todos los vectores en función del nuevo residuo y determinadas sus relaciones de recurrencia, así como los escalares que intervienen en las mismas, se puede escribir el algoritmo. ALGORITMO BI-CGSTAB Aproximación inicial XQ. ro = b — AXQ} TQ arbitrario. Fijar po = ro Desde j = 1,2,... hasta converger, hacer ' {Apj,r*) Sj = rj - ajApj UJ'i — y^/\Sj, /\.Sj j Xj+\ = Xj + ajPj + WjSj fj+i = Sj — WjAsj _ {rj+i,r*o) aj Pj — (fj^ro) w 3 Pj+i = ^j+i + Pj {Pj - WjApj) Fin En el algoritmo figuran dos productos matriz por vector y cuatro productos escalares, mientras que el CGS exige los mismos productos matriz por vector y sólo dos productos escalares. Sin embargo en la mayoría de los casos la convergencia del Bi-CGSTAB es más rápida y uniforme, necesitando menor carga computacional para alcanzar una determinada tolerancia, pues la reducción del número de iteraciones compensa su mayor coste. 26 Conceptos básicos 2.2.2.3. Método de QMRCGSTAB Desarrollado por Chan y otros [31] está basado en la aplicación del principio de minimización usado en el algoritmo Bi-CGSTAB al método QMR, de la misma forma que el TFQMR es derivado del CGS. Sea, ^fc = [yi,y2,---yk] , Wk+i = [wQ,Wi,...Wk] tal que. Y' y2i-i=Pi para Z = l,...,[(A; + l)/2] y2i = si para¿= l,...,[k/2] { W21-1 = Si para l = l,...,[{k + l)/2} W2i = ri para ¿ = 0,1,..., [k/2] donde [k/2] y [{k + l)/2] son la parte entera de k/2 y {k +1)/2 respectivamente. Definiendo [6i,52,... Sk] tal que, 021 = üJi para l = l,...,[ik + l)/2] S21-1 = oci para Z = 1,..., [(A; + l)/2] Entonces, para cada columna de H4+i y Yk, las expresiones (2.56) y (2.58) pueden escribirse como, Ay^ = {wm-i - Wm) d;^, m=l,...,k (2.60) o usando notación matricial, AYk = Wk+iEk+i donde Ek+i es una matriz bidiagonal (k + l) x k con elementos diagonales 6:;^^ y en la diagonal ii\ferior —S:^^. Esto puede ser fácilmente comprobado hasta que el grado de los polinomios correspondientes a los vectores r^, Sj y pj sean 2j, 2j — 1, y 2j — 2, respectivamente. Entonces, Yk y Wk generarán el mismo subespacio de Krylov generado por ro pero de grado fc — 1. La idea principal del QMRCGSTAB es encontrar ima aproximación a la solución del sistema (2.1) usando el subespacio de Krylov ICk-i en la forma, Xk^xo + YkQk con Qk e E" La expresión para el vector residuo queda, Tfc = ro - AYkQk = ro~ Wk+iEk+i9k Métodos basados en subespacios de Krylov 27 Teniendo en cuenta el hecho de que el primer vector de W^+i es justamente ro, entonces, Tfc = Wk+i (ei - Ek+igk) Como las columnas de Wk+\ no están normalizadas, se usa una matriz de escalado Sfe+i = diag (CTI, ..., Uk+i) con GJ = \\wj \\ para hacer unitarias las columnas de Wfc+i. Entonces, Tfc = Wk^illll^llk+i ifii - Ek+i9k) = Wk+iT.'l^li (cTiei - Hk+ign) con Hk+\ = '^k+iEk+iLa aproximación QMR consiste en la minimización de ||o"iei — ií^+i^fcH para algún g G E*^, donde este problema de mínimos cuadrados es resuelto usando la descomposición QR de la matriz Hk+i de forma incremental utilizando las rotaciones de Givens. Como -ff^+i es bidiagonal inferior, solamente es necesaria la rotación en los pasos previos. ALGORITMO QMRCGSTAB Aproxim.ación inicial XQ. ; ro = bAxo; Elegir ro tal que po = r^ro y^ O Po = vo = do = 0; Po = «o = ^0 = 1; T = ||ro||, 9o = 0,77o = 0; Desde A; = 1, 2,..., hacer: Pk = Tork-ü Í3k = (Pk^k-i) /Pk-i<^k-ú Pk = Tk-i + Pk {pk-i - uJkVk-i); Vk = Apk) «fe = Pk/r^Vk, Sk = rk-i - akVk; Primera cuasi-minimización ^k = \\sk\\ ¡r; c = ; T = T Qkc; rjk = c^ak) «fe = Pfc H "fc-i; ^~ Xk = Xk-i + rjkdk; HaUar: tk = Ask) {tk,tk) rk = Sk~ ujktk) Segunda cuasi-minimización dk = Ikfcll /r; c = • T = T dkc; Vl + ^fc 28 Concqytos básicos 'ü2~ _ dk = Sk -\ dk) Xk = Xk+ r¡kdk; Si Xk converge parar Fin 2.2.3. Método GMRES (Generalized Minimal Residual) El GMRES [92, 93, 94] es un método de proyección sobre un subespacio de Krylov Kk de dimensión k, basado en minimizar la norma del residuo. Así, el desarrollo del algoritmo consiste en encontrar un vector xáexQ + Kk tal que, x = xo + VkV Imponiendo la condición de mínimo para, J{y) = \\b-Ax\\ (2.61) como, b-Ax = b-A{xo + Vkv) = ro - AYuy y teniendo en cuenta que AVk = VkHk + Wkel = Vk+iHk (2.62) y que vi = ro/ || TQ ||, llamando /? =|| TQ ||, entonces, b-Ax = /3viVk+iHkV Pero, Vi = Vk+i&ir con ei € R''"^^ por tanto, b-Ax = Vk+i{l3ei-líky) (2.63) y la ecuación (2.61) quedará, J{y) = \\Vk+i{(3ei-líky) Como las columnas de la matriz 14+1 son ortonormales por construcción, podemos simplificar la expresión anterior, J{y) = \\{pei-Hky)\\ (2.64) El siguiente algoritmo del GMRES busca el único vector de xo + Kk que minimiza la ftmcional J{y). ALGORITMO GMRES Aproximación inicial XQ.rQ = b — AXQ; Precondicionamiento 29^ Definir la (A; + 1) x A; matriz Hk = {H}i<i<k+i,i<j<kPorier Hk = 0. Desde j = 1,..., k hacer Wj = Avj Desde i — I,..., j hacer {H}ij = {wj,Vi); Wj = Wj - {H}^. Vi) Fin {ií},+ij =li -^j II; Si {H}j+ij = o poner k = JY parar 1 Fin Hallar yk que niinintLza || (^ ei — iífcí/) ||; Determinar Xk — xo + VkVk siendo 14 = [vi,V2,...,Vk\; Calcular Vk — b — Axk El algoritmo GMRES, resulta impracticable cuando k es muy grande, ya que conlleva un elevado coste computacional y de almacenamiento. Por esta razón, se suele usar la técnica de restart y de truncamiento. ALGORITMO Restarted GMRES 1) Aproximación inicial XQ. TQ = 6 — Axo, (3 =|] ro ||, y wi = ro/P; 2) Generar la base de Amoldi y la matriz Hk usando el algoritmo de Amoldi; 3) Iniciar con vi; 4) Hallar yk que minimiza || {peí - Hky) \\yxk = xo + Vkyk, 5) Si se satisface entonces parar. En caso contrario poner XQ-.— Xkj volver al paso 1. 2.3. Precondicionamiento La convergencia de los métodos basados en los subespacios de Krylov mejora con el uso de las técnicas de precondicionamiento. Estas consisten generalmente en cambiar el sistema original (2.1) por otro de idéntica solución, de forma que el número de condicionamiento de la matriz del nuevo sistema sea menor que el de A, o bien que tenga una mejor distribución de autovalores. Para efectuar el precondicionamiento, se introduce una matriz M, llamada matriz de precondicionamiento. Para ello multiplicaremos ambos miembros del sistema (2.1) por la matriz M, MAx = Mb (2.65) tal que, K {MA) < K [A] 30 Conceptos básicos El menor valor de K{A) corresponde a M = A~^, de forma que K (AA'^) — 1, que es el caso ideal, para el cual el sistema convergería en tma sola iteración, pero el coste computacional del cálculo de A~^ equivaldría a resolver el sistema por un método directo. Por ello se sugiere que M sea una matriz lo más próxima a A~^ sin que su determinación suponga un coste elevado. Generalmente, se considera como matriz de precondicionamiento a M~^ y obtener M como una aproximación de A, esto es, M-^Ax = M'^h Por tanto, la matriz M debe ser fácilmente invertible para poder efectuar los productos M~^ por vector que aparecen en los algoritmos precondicionados sin excesivo coste adicional. Por ejemplo, en el caso en que M es una matriz diagonal o está factorizada adecuadamente, se pueden efectuar dichos productos mediante remonte sin necesidad de calcular M~^. Dependiendo de la forma de plantear el producto de M~^ por la matriz del sistema obtendremos distintas formas de precondicionamiento. Estas son, M~^Ax = M~^b (Precondicionamiento por la izquierda) AM~''-Mx = b (Precondicionamiento por la derecha) (2.66) M¡^^AM2^M2X = M~^b (Precondicionamiento por ambos lados) si M puede ser factorizada como M = M1M2. Los precondicionadores pues, han de cumplir dos requisitos fundamentales, fácil implementación, evitando un coste computacional excesivo del producto de M"^ por cualquier vector y mejorar la convergencia del método. Por tanto una matriz que sea ima aproximación más o menos cercana de A, obtenida con estos criterios puede dar lugar a un buen precondicionador. El campo de posibles precondicionadores es, así, muy ampHo. Algunos de los más usados son los bien conocidos Diagonal o de Jacobi, SSOR e ILU. Consideraremos estos, además de la Inversa Aproximada de estructura Diagonal que hemos denomimado Diagonal Óptimo, para comparar con los precondicionadores basados en la inversa aproximada de estructura sparse general propuestos en esta tesis. 2.3.1. Precondicionador diagonal Surge comparando la fórmula de recurrencia para la solución que resulta de aplicar el método de Richardson, cuya relación de recurrencia viene dada por Xi+i = Xi + a{b — Axi), con a > O, al sistema precondicionado con la fórmula correspondiente que se obtiene apHcando el método de Jacobi al sistema sin precondicionar. De la aplicación del método de Richardson al sistema precondicionado, M-^Ax = M-^b (2.67) Precondicionamiento 31^ se obtiene, para el cálculo de los sucesivos valores de la solución, Xi+i = Xi+a (M~^fe - M~'^Axi) Multiplicando por la matiiz de precondicionamiento M y haciendo a = 1, resulta Mxi+i = Mxi + a{bAxi) (2.68) Por otro lado, descomponiendo la matriz del sistema en A = D — E — F, (siendo D la matriz diagonal formada por los elementos de la diagonal de A y E y F matrices triangulares inferior y superior, respectivamente), y utilizando el método de Jacobi para la resolución del sistema, se obtiene, Dxi+i = Dxi + {bAxi) (2.69) Comparando las expresiones de recurrencia finales de ambos métodos, se observa que el método de Jacobi aplicado al sistema sin precondicionar, equivale al de Richardson, menos robusto y más simple, cuando este se aplica al sistema precondicionado con la matriz diagonal D. Resulta así vm precondicionador elem.ental, fácü de implementar y con matriz inversa que se determina con muy bajo coste computacional, ya que la matriz de precondicionamiento es diagonal y sus entradas son las de la diagonal de A. 2.3.2. Precondicionador SSOR Si aplicamos el método SSOR al sistema sin precondicionar, considerando la descomposición de la matriz AenA = D — E — F como en el caso anterior siendo u el parámetio de relajación, se obtiene para la solución, D \"Vl-í^„ „\ ÍD ^-^ u^-FY'-^(^-Ey\ operando para expresar esta relación de forma que se pueda comparar con la solución que resulta de aplicar el método de Richardson al sistenia precondicionado, -j^—. (D - uE) D-' {D - uF) Xi+^ UJ y¿ — üJ) = -Tx^^ T {D - LüE) D-^ {D - cüF) Xi + {bAxi) con lo que resulta como matiiz de precondicionamiento, M = —- ^ (D - üüE) D-^ {D - uF) (2.70) 32 Conceptos básicos que en el caso de sistemas simétricos, como se cumple que, {D ~ uF) = {D - üüEf podem.os expresarla como un producto de dos matrices triangulares traspuestas. M = jD - uE) D-^l^ y/uj (2 - Lü) ^üJ (2 - LO) TT (2.71) para el caso de sistemas no simétricos también lo podrem^os expresar como un producto de matrices triangulares, inferior y superior respectivamente. M= (J - uED-^) 2.3.3. Precondicionador ILU(O) D-ÜJF Uj{2-Lü) (2.72) Resulta de la aproximación de A por una factorización incompleta LU, guardando las misma entradas nulas en las matrices triangulares LyU [81], A^LU^ ILU{0) = M donde las entradas de M, rriij, son tales que. rriij = O if üij = O {A-LU}.j = 0 if aij^O (2.73) (2.74) (2.75) Es decir que lo elementos nulos de la matriz del sistema siguen siendo nulos en las posiciones respectivas de las matrices triangulares para no incrementar el coste computacional. 2.3.4. Precondicionador Diagonal Óptimo Resulta de resolver un problema de minimización [8], mín||MA-J| Mes" ' lA^^-JI donde S es el subespacio de las matrices diagonales de orden n. La solución a este problema, que puede verse en [56], es: N = diag Olí «22 |2' 2 J • • • > eÍAWtWe^Ar, ' \KA\\l n \\NA-I\\l = n-J2 \e7A\ñ (2.76) (2.77) Precondicionamiento 33^ 2.3.5. Algunos métodos de Krylov precondicionados 2.3.5.1. Algoritmo Bi-CGSTAB El método BICGSTAB introduce un nuevo parámetro cu en cada iteración para minimizar el residuo (ver [104]). Para cada forma de precondicionainiento, el valor del parámetro uj resulta: (Asfs {AsfAs (M-Ms)^(M-^5) {M-^Asf{M-^As) {L-'Asf{L-h) (L-i^s)^(L-Ms) LO — ——f— precondicionamiento por la derecha ) üj = ,V^^ , ^"T Vr T y^ precondicionamiento por la izquierda (2.78) u = — IjA — precondicionamiento por ambos lados De este modo, el cálculo de tó incluye un proceso de sustitución por iteración para el precondicionamiento por la izquierda y dos para el precondicionamiento por ambos lados. No obstante, si obtenemos u a partir de la minimización del residuo del sistema original sin precondicionar, todos los valores dados en la ecuación (2.78) coinciden con el del precondicionamiento por la derecha. En este caso se obtiene también un único algoritmo. ALGORITMO BICGSTAB PRECONDICIONADO Aproximación inicial XQ.rQ = h — AXQ; Elegir TQ arbitrario, po = TQ Desde j = 1,2,... hasta converger, hacer Resolver Mzj = yj = Apj Resolver Mvj = "j = Sj = Uj = tj = Uj - Xj+i ^j+l 01 = {zj^rD {vj,r^) Tj - ajVj : Zj - OijVj Auj \tj, Sj) V'j' ^j) = Xj + ajPj = Sj — U)jtj ^ó Vi + ^i'^j Oj "^3 Pj+i = Zj+i + f3j (pj - UjVj) Fin Cualquier otra elección del residuo inicial conducirá a una nueva forma de precondicionamiento (ver [98]). 40 Conceptos básicos Podemos también construir un método híbrido, por ejemplo, combinando los dos métodos en un camino con dos etapas donde se use un método de factorización incompleta (implícito) para la matriz particionada en bloques y, en la segunda etapa se use tm método directo para aproximar la inversa del bloque pivote que aparece en cada etapa de la factorización. En esta sección presentamos varios métodos para calcular aproximadas inversas de matrices. También discutimos cómo pueden aproximarse matrices que no se conocen explícitamente pero de las que se dispone de su producto matriz por vector Una estructura común de las matrices inversas aproximadas es la estructura en banda. Los rangos para las entradas de ima matriz dada disminuyen según la distancia a la diagonal aumenta, la exactitud de dichas aproximaciones es de primordial importancia y discutiremos este tópico al final de la sección. 2.6.1. Dos métodos para el cálculo de aproximadas inversas de matrices en banda por bloques Para ilustrar los métodos de aproximación explícita e implícita consideramos primeramente el cálculo de una aproximada inversa de ima matriz en banda. Considérese una matriz B = [6^], donde bij. = O si j < i — p o j > i + q. Por tanto, B es ima matriz en banda con un ancho de semibanda izquierdo p y un ancho de semibanda derecho q. De hecho, podemos considerar el caso más general donde las entradas de B son matrices en bloque. Podemos asumir que los bloques de la diagonal B^j son regulares y posteriomente hacer las suposiciones que se requieran. Se dice que B tiene un ancho de semibanda izquierdo del bloque igual apy vn ancho de semibanda derecho del bloque igual a. q si Bij =0,sij <i — po j > i + q. 2.6.1.1. Un método explícito de aproximación Consideraremos un método de cálculo de una inversa aproximada G de una matriz en banda B particionada por bloques. Sea G = [Gij] la partición correspondiente a la inversa aproximada de B a calcular, y asumimos que G tiene im ancho de semibanda del bloque pi y 9¿ respectivamente, donde Pi > p Y qi> q. Para calcular Gij i —pi < j <i + qi hacemos, {GB).^ = Y^ Gi,kBk,j = Aij, i-pi<j<i + qi (2.82) k=i—pi donde. ^..^ j ^i paraz = j '•^ ' O para el resto siendo /¿ la matriz identidad del mismo orden que Bij . La ecuación (2.82) da Pi + 9i + 1 ecuaciones con matrices en bloque que pueden usarse para calcular Precondicionadores explícitos e implícitos 42 el mismo número* de bloques de G^. Sin embargo, com.o Bkj = 0 si k<j — qo k > j +ptenemos, mín{i+gi ,j+p} Y2 Gi,kBk,j = \j, j = i-pi,i-pi + l,...,i + qi (2.83) iaáxk={i—pi,j—g} de forma similar, pero a menor escala las ecuaciones son válidas para i < 2pi +1 y i > m — qi, donde m es el orden de B tenemos. i—pi+p paraj = i-pi: j^ Gi,kBk,i-p^ = O k=i—pi i—pi+p+1 para j = i-pi + l: J2 Gi,kBk,i-p,+\ = O k=i—pi+l i—pi+p para j = i : j^ Gi,kBk,i = O k=i—p\ para j = i + qy: ^ Gi,kBk,i+q, = O k—i+qi Estas ecuaciones tienen solución si hacemos algunas suposiciones adicionales sobre B. Ahora mostramos el método para el caso especial, pero importante, donde pi = p = qi = q = 1, esto es, el caso donde B y G son bloques triangulares donde es suficiente asumir para poder obtener solución, que Bij y la matriz de Schur complementaria en (2.84) son no singulares. Entonces (2.83) muestra que Gt,i_iSi_i^j_i + GijBi^i-i = O Gi^i-\Bi-i^i -\- Gi^iBi^i + Gi^i+iBi+i^i = I, Gi,iBi^i+i + Gj^j+i^í+i^j+i = O así. Gi^i-i — —Gi^iBi^i_i[Bi-i^i-i) Gi,i+i = —Gi^iBi^i[Bi+i^i+i) Gi^i-i — [Bi^i — 5j,i_i(B¿_i^¿_i) Bi-i^i — Bi^i+i{Bi^.i^i+i) Bi+i^i\ (2.84) donde debemos asumir que Bu es no singular y que la inversa de la última matriz existe (complemento de Schur). Por ejemplo, si B es una íí-matriz en bloques, mostraremos que las condiciones de solubilidad se cumplen. Para i = 42 Conceptos básicos 1 y para i = n ninguno de los términos entre paréntesis en (2.84) está ausente. Se puede ver que GÍ,Í-I, GÍ,Í, GÍ^Í+I es la i-ésitna fila de la inversa exacta de /I B = o\ \o // Es interesante resaltar que la matriz entre paréntesis en (2.84) es la misma que la matriz obtenida (en la i-ésima fila) cuando hacemos un reordenamiento par-impar de la nnatriz tridiagonal por bloques dada respectivamente. Nótese también que G es no simétrica en general, porque las matrices Gj,¿_i y Gf_i.i no tienen que ser iguales necesariamente, aún cuando B sea simétrica. Las entradas de los bloques en cada fila de G pueden calcularse consecuentemente. Finalmente puede verse que GB contiene una diagonal principal unitaria y una sub y super diagonal a la distancia dos de la diagonal principal. El resto de las entradas son nulas. De aquí que sea sparse. Para una matriz tridiagonal tenemos gi,i = (bi biA-ik i—XA bi,i+ibi 1 + 1,5 bi-iA-i bi i+l,i+l -1 91,1-1 9i,i^i,i—'i. bi-i^i-i 9i,i+'i. 9i,ib 'i,i+l Ji+lA+l i = 1,2, ...,n 2.6.1.2. Método implícito de Aproximación Considérese otra vez una matriz con bloques en banda B. Existe xin algoritmo para calcular los bloques de la inversa exacta de B dentro de una banda consistente en pi >p bloques a la izquierda y qi> q bloques a la derecha de la diagonal, en cada fila del bloque. El algoritmo fue presentado en [9] y en [5] y se basa en tma idea presentada en [101] (ver también [47]). El algoritmo requiere la factorización previa de B. Por tanto, este método es prácticamente apHcable solo cuando los anchos de banda p,qde B sean pequeños y las entradas S^ sean escalares o, en caso de ser matrices en bloque, tener órdenes pequeños. Si B se factoriza en la forma B = {IL)D-\I - Ü) (2.85) Precondicionadores explícitos e implícitos 43 donde D es una matriz en bloque diagonal y L,U son bloques estrictamente triangulares inferiores y superiores respectivamente, entonces no es necesario calcular mas inversas de las matrices en bloques. Debe resaltarse que el cálculo de la matriz D en (2.85) requiere el cálculo de las inversas de las matrices en bloque pivotes. Una segunda ventaja importante del algoritmo es que no necesitamos calcular ninguna entrada de la inversa fuera de la parte de la banda requerida por B~^.El algoritmo se basa en las siguientes identidades: Lema 9 Sea B = LDU~^, donde L = I —L,U = I — ÜyLyÜ son estrictamente triangular inferior y superior respectivamente. Entonces; (a) B-^ = DL-^ + ÜB-^ (b) B-^ = U-^D + B-''L Demostración. Tenemos B'^ = U-'^DL así, (/ - Ü)B-^ = DL'^ y B-\I - L) = U~^D, los cual es (a) respectivamente (h) • Las relaciones (a) y (b) pueden usarse para el cálculo de los triángulos superiores e inferiores de B~^respectivamente. Como L"^ es triangular inferior con matriz en bloque diagonal igual a la identidad, L"^ no se considerará en el cálculo del triángulo superior de J5~^ de igual forma que Í7""^ no se considerará en el cálculo del triángulo inferior. El algoritmo para el cálculo de la banda de B~^ con bloques inferiores y superiores con ancho de semibanda pi y qi tienen la forma siguiente. Asumimos que B = [Bij] tiene un bloque de orden n, tm bloque con ancho de semibanda inferior y superior pyq respectivamente y ha sido factorizada en la forma (2.85) ALGORITMO BBI (INVERSA EN BLOQUES POR BANDA) Desde r = n, n — 1,..., 1, hacer míií(q,n—T) iB~ )r,r = Dr^r + /^ f^r,r+s(-B~ )r,r+s (2.86) Desde k = 1,2,..., gi, hacer niín(g,n—r+fc) {B~ )r-k,r = 2_^ Ur-k,r-k+s{B~ )r-k+s,r (2.87) Desde k = 1,2,..., gi, hacer mm(p,n—r+fc) -k+tLr-k+t,r-k (2.88) t=i Cabe destacar que sólo se realizan los productos matriz-matriz. Además si B es simétrica, solamente se necesita ima de las ecuaciones (2.87) o (2.88). 44 Conceptos básicos Comentario 10 La ecuación (2.86) muestra en particular que, donde Dn,n ^s el último bloque en D. Esto se convertirá en una propiedad muy útil para los métodos de estimación de valores singulares y números de condición de matrices precondicionadas. Comentario 11 Patrón de sparsidad de la envoltura: El algoritmo BBI puede extenderse al caso donde B tiene un patrón de sparsidad en la envoltura, esto es, cuando existen p, q tales que Bij = O para i — j > pty para j — i > QÍ. En particular, el método puede utilizarse para calcular una parte de la envoltura de la inversa de la matriz semibanda, para lo cual Bij = O, |¿ — j| > p, pero Bi^n ¥" OTales matrices aparecen, por ejemplo, para las ecuaciones en diferencias elípticas con condiciones de contorno periódicas. En este caso, el algoritmo toma la forma de la ecuación (2.86), pero debe añadirse el término Ür-k,n{B~^)n,r ci la ecuación (2.87) y el término (B~^)r,nLn,r-k a la ecuación (2.88) porque la última fila de L y la última columna de Ü están llenas generalmente. Más aún, para r = n debe calcularse {B'^)n-k,n y iB~^)n,n-k P^^^ todo k,k = 1,2,...,n— leu las ecuaciones (2.87) y (2.88) 2.6.1.3. Matrices definidas positivas El método anterior tiene una desventaja: aún cuando B sea definida positiva, la banda [5"^]^'^^ de B~^ puede no ser definida postiva como se muestra en el siguiente ejemplo. Ejemplo 12 Sea G = banda [G]^'^ de G, donde 1 -2 1 -2 5 -3 1 -3 4 Entonces G es definida positiva, pero la [or = 1 -2 O -2 O 5 -3 4 es indefinida. Más generalmente, esto puede observarse deG^^ = G — R. Aquí, G^^ es la banda simétrica, esto es, p = q, de G. R es indefinida puesto que tiene la diagonal nula. Por tanto, puede ocurrir que el menor autovalor de G^' sea negativo y entonces G^^'l será indefinida. Sin embargo, para una clase importante de matrices B, denominadas matrices monótonas, podemos modificar la matriz G^^ con una compensación diagonal para las entradas que se eliminan en R, tal que la matriz modificada se convierta en definida positiva. Definición 13 Una matriz real A se dice que es monótona si Ax>^ implica x>0. Precondicionadores explícitos e implícitos 45^ Lema 14 A es monótona si y sólo si A es no singular con A~^ > 0. Teorema 15 (Compensación diagonal). Sea B monótona, simétrica y definida positiva. Sea G = B~^ y R = G — G^', donde G^^ es la banda simétrica de G con ancho de semibanda p. Entonces G = G^^ + D es definida positiva, donde D es una matriz diagonal tal que Du = Ru, para algún vector positivo u. Demostración. Como G = B~^ > O, tenemos R>0. Por tanto, £> > 0. Sea y = diag(ui,U2,-,'"n) Entonces la matriz {D — R)V tiene diagonal dominante, y, por el teorema de los círculos de Gershgorin, los autovalores son no negativos. Lo mismo ocurre para V~'^{D - R)V y también para D - R, porque D - Res equivalente a esta última matriz. Por tanto, D — Res semidefiíüda positiva y la relación G-G = G^^+D-G = D-R muestra que G — G es semidefínida positiva. Entonces, G es definida positiva debido a que G también lo es, al ser la inversa de una matriz definida positiva. • Dentro de todos los posibles vectores u, debe seleccionarse aquel que mejore el número de condición de G^^B lo más posible. Esto es similar al uso de un vector para la modificación en los métodos de factorización incompleta. Tanto los m.étodos explícitos como los implícitos tratados anteriormente fueron usados en Axelsson, Brinkkemper e Il'in [9] para cacular aproximaciones tridiagonales de las inversas de m.atrices tridiagonales que se obtienen durante la factorización en bloque de matrices en las ecuaciones en diferencias para ecuaciones diferenciales elípticas de segundo orden en dos dimensiones. Concus, Golub y Meurant [37] utilizaron también un método implícito para calcular la parte de la banda tridiagonal de la inversa de una matriz tridiagonal, basados en un método que usa dos vectores para generar la matriz, el cual fue sugerido primeramente por Asplund [4]. Sin embargo ocurre que este último método no puede extenderse de forma estable para calcular la banda con anchos de semibanda m.ayores o iguales que 2, por lo que en la práctica está limitado al caso p = q= 1.' 2.6.2. Una Clase de Métodos para calcular las inversas aproximadas de matrices Consideraremos ahora tina rama general de la clase de métodos para construir inversas aproximadas, basados en los trabajos de [75] y [73]. Veremos que los métodos explícitos e implícitos presentados anteriormente se relacionan y además, para matrices simétricas son equivalentes a las dos versiones de esta clase de métodos. 46 Conceptos básicos La idea básica es la siguiente: Dado im patrón de sparsidad, se calcula la inversa aproximada G con dicho patrón de sparsidad para una matriz no singular A dada de orden n, esta será la mejor aproximación en alguna norma. Las normas basadas en las trazas de la matriz (J — GA)W(I — GA)'^, esto es, para el cuadrado de la norma de Frobenius con pesos de la matriz de error (J — GA), proporcionan métodos prácticos y eficientes para ciertas selecciones de la matriz de pesos W. Considérese la funcional, Fw{G) = 11/ - GAWl, = tr{{IGA)W{I - GAf} donde W es tma matriz simétrica y definida positiva y se asume que G depende de algunos parámetros libres, ai,..., ctp. Estos parám.etros pueden ser las entradas de G en posiciones definidas por un conjunto de p pares de índices (i, j) e S, que es un subconjunto del conjunto total de pares, {{i,j);l<i<n; 1 < i < n}, siendo n el orden de A. Más aún, gij = O para todo par de índices fuera del patrón de sparsidad S. Entonces, S define el patrón de sparsidad de G = {QÍJ} y gij = O para todo {i, j) e S'^ que es el complementario de S definido por S^ = {(i,j) ^ S;l < i < n;l < j < n} Al igual que en el caso anterior, S contiene al menos los pares de índices {{i, i); 1 <i < n}. Así, el menor patrón de sparsidad de G es el de la matriz diagonal. Obsérvese que FwiG) > O y Fw{G) = O si G = A~'^. Es más, el lema 16, nos muestra que ||-||^^ es una norma tal que ||7 — G^Hi^ nos da una medida del error cuando G se aproxima aA~^. Queremos calcular los parámetios {ai} para minimizar Fw{G) = Fw{ai,...,ap) La solución a este problema debe satisfacer las relaciones estacionarias dFw doi = 0,i = l,...,p Sea \\B\\y^^ = {ír(5WS^)}^. Para demostrar que ||-||^;,^ es ima norma, definimos el siguiente lema del cual se también se derivan otras propiedades útiles. Lema 16 Sean A, B dos matrices cuadradas de orden n. Sea W una matriz definida Precondicionadores explícitos e implícitos 47 positiva. Entonces, (a) tr{A) = ír(A^), tr{A + B) = tr{A) + tr{B) (b) triA) = f:Xi(A) ¿=i (c) tr{AA'^) = ¿ üij y tr{AA^) > tr{A^) (d) \\AB\\^r < \\A\\j\\B\\^,, donde \\-\\j es la norma Frobenius, m\j <{ijaT,lj}'= {HAA^)}' (e) \\A\\,<\\A\\^y\\AB\\^<\\A\\^\\B\\^,siysólosiW-Ies definida positiva (fl W^Ww — {^'''i^^^^)}^ ss una norma aditiva Demostración. El apartado (a) se sigue de la definición del operador traza, n tr{A) = ^üij y (b) se muestra en [8]. La primera propiedad de (c) se obtiene por el cálculo directo. Para demostrar la desigualdad, nótese que, i j i j i j Para demostrar id) partimos de \\AB\\p, < \\A\\p \\B\\p (Apéndice A de [8]). De esta forma, \\AB\\^ = {tr(ABWB^Á^)y = {ír(^M^A^)}' AB <\\A\ B =Uh\\B\ w donde AB = BW^. Ahora (e) se sigue de (d), y \\Al = {tr{AA^)y^ = [£HAA^)}' < {Y:Xi{AWA^)}' = {triAWA'')}'' = \\A\\^ (2.89) donde la desigualdad se cumple si y sólo si 14^ — 7 es semidefinida positiva. Finalmente, \\A + B\\^ = ¡tr UAW^ + BW^){AW^ + BW^f] V = \\AW^ + BW^ < AW^ + I BW^ < \\A\\y^ + \\B\\w 45 Conceptos básicos Considérese ahora el problema: hallar G e S tal que ||/-G'>1||^< I-GA paratodoGe5 w donde I-GA es la menor para G = G dentro de todas las matrices con w ^ patrón de sparsidad S. Para simplificar la notación usamos la misma que para el conjunto de matrices que tienen este patrón de sparsidad como el correspondiente conjunto de índices. Nótese que, Fw{G) = tr {(/ - GA)W{I - GAf} = trW - tr{GAW) - tr{GAWf + tr(GAWÁ^G'^) = trW - Y,9ij [{AW)ji + {AW^)ji] + tr{GAWA'^G^) Por lo tanto, las relaciones estacionarias, muestran que, -iAW)ji - {AW^)ji + {AWA^G'')ji + iGAWA'^)ij = O o, como W = W^, iGAWA%j^iWÁ^)ij, {i,j)eS (2.90) La ecuación (2.90) define el conjunto de ecuaciones que deben satisfacer las entradas de G € 5. En dependencia de la selección de S y/o A, estas ecuaciones pueden o no tener solución. Por otro lado, la ecuación (2.90) es un caso particular de {GAV)ij^Vij, {i,j)eS (2.91) donde V = WA'^ en la ecuación (2.90). Puede verse que (2.91) tiene solución única si todos los menores de AV, restringidos al conjunto S, son no singulares. Entonces, el problema es encontrar las matrices apropiadas W (oV) y los conjuntos 5. Comentario 17 (Cálculos en paralelo) Si asumimos que V es conocida deforma explícita, las ecuaciones para el cálculo de las entradas de G pueden paralelizarse, porque las entradas en cualquier fila de G se calculan independientemente de las entradas de otras filas, es decir, pueden ser calculadas en paralelo entre las filas. Una matriz G definida de esta forma será en general no simétrica aún si A es simétrica. Establecemos los resultados anteriores en un teorema y consideramos algunas selecciones particulares e importantes en la práctica de la matriz de pesos W. Precondicionadores explícitos e implícitos 49^ Teorema 18 Sea A no singular yseaGeS tal que, \\I - GA\\^ < I - GA para cualquier G G S Entonces G debe ser solución del sistema lineal {GAWA'')ij = iWA^)ij, {i,j)eS Para la siguiente selección de W, este sistema y Fw{G) son como sigue, (a) Para W = A~'^, donde A es simétrica y definida positiva, entonces, Fw{G) = tr{A-^) - tr{G) y G satisface, {GA)ij = 5ij, {i,j)eS (b) Para W = J„, entonces, Fw{g) =ntr{GA) y G satisface, {GAA^)ij = {A'^)ij, ii,j)eS (c) Para W = A'', donde k es un entero positivo y A es simétrica y definida positiva, entonces, Fw{G) = tr{A'' - GA''^^) y G satisface, iGA'^%j = {A'^^, {i,j)eS (d) Para W — (A^Ay^ entonces, Fw{G) = tr [{A-^ - G)A-'^] y G satisface, Gij = {A-%j, ii,j)eS Demostración. La parte general establecida es (2.90) la cual se ha demostrado anteriormente. Para la matriz G óptima, el lema 16 (a) muestra que, FwiG) = trW - tr{2GAW - GAWÁ^G'^) = trW - 2tr{GAW) + tr{WA'^G'^) = trW - tr(GAW) = triW - GAW) = tr{{I - GA)W) Con lo cual se obtienen las diferentes selecciones de W. m Note que para WA~^ tenemos el método explícito (2.82) discutido anteriormente para A simétrica y para W(A'^A)~'^ usamos el método implícito. 56 Conceptos básicos Teorema 26 Sea A una H-matriz y G, definida por (GA)ij = Sij, {i,j) e S. Entonces G está únicamente determinada y p{I — GA) < 1 Demostración. De forma similar a la demostración de la parte (b) del teorema 23, notamos primeramente que la i-ésima fila de G es la inversa exacta de a^^) 0 A para un producto matricial de Hadamard. Claramente, a^*^ © yl es una ií-matriz. Por tanto,, su inversa existe. Como esto es cierto para cualquier i, G está determinada de forma única. Sea M (A) la matriz de comparación y sea G — [^¿j] la inversa aproximada de M{A) para el conjunto S, esto es, {GM{A)}i,^^5i,j, ii,j)eS Primero, el teorema de Ostrovsky muestra que, para inversas exactas Gi y Gi de cn^') Q Ay a^^^ 0M{A), respectivamente, tenemos |G,| < Gi. Esto es, |G| < G Demostraremos a continuación que \GA\ < GM{A) (2.93) Esta inecuación es trivial para {i,j) € S , por lo que necesitamos considerar solamente (i,j) ^ S. Para tales posiciones-tenemos \{GA)ij\ = k;k^j{i,k)GS < E \9i,k\\{M{A)},j\ k;k^j{i,k)eS < E g,,i,[-{M(A)K,] = -{GM(A)),j k;k^j(i,k)eS Por tanto, \I-GA\< IGM{A) y, usando el teorema de Perron-Frobenius, p{I - GA) < p{\I - GA\) <p{IGM{A) ) = p{I-GM{A))< 1 donde la últim.a desigualdad se sigue del corolario 24 • Comentario 27 (H-Matrices en bloques) Tomando normas en vez de valores absolutos encontramos que el teorema 26 también se cumple para H-matrices en bloques y que \\GiA<h3 y II(/-G^)MII< {I-GMM))Í (2.94) donde gij son las entradas de la matriz G que satisfacen {GMb{A)}ij = 5ij, {i, j) e S, y Mb{A) es la matriz de comparación en bloques Precondicionadores explícitos e implícitos 57 Corolario 28 Sea A una H-matriz y M{A) su matriz de comparación, y sea x > O tal que M{A)x > 0. Entonces M(GA)x > O, donde G se define como en el teroema 26, esto es, la dominancia diagonal generalizada se mantiene. GM{A) Como Demostración. Por la ecuación (2.93) tenemos \GA\ < M{A)x > O para algún x positivo y G > O, tenemos GM{A)x > 0. Las entradas de la diagonal de GA son unos y {GÁ)ij <0,ij^j (ver la demostración del teorema 23). Entonces, M{GA)x ^x-ilGA)x >x-\I-GA\x > X I-GM{A) X X I-GM(A) X = GMiA)x > O donde hemos usado la propiedad de que GM{A) es tina Z-matriz. • Los métodos explícitos de precondicionamento fueron propuestos primeramente por Benson [12], Frederickson [50], Benson y Fredericson [13] y Ong [89], siendo las mejores aproximaciones en el sentido de la norma de Frobenius. Fueron extendidos poteriormente con un análisis cuidadoso a las mejores aproximaciones en normas de Frobenius generalizadas por [75, 73]. Capítulo 3 Inversa aproximada basada en el producto escalar de Frobenius En este Capítulo se presentan unos precondicionadores de construcción en paralelo para la resolución de sistemas de ecuaciones lineales. El cálculo de estos precondicionadores se realiza mediante proyecciones ortogonales usando el producto escalar de Frobenius. Así, el problema mín||ylM — I\\F y la matriz MQ € S correspondiente a este mínimo, (siendo S cualquier subespacio de Mn{^)), se calcula explícitamente usando una fórmula acumulativa para reducir los costes computacionales cuando el subespacio S se extiende a otro que lo contenga. En cada paso los cálculos se realizan aprovechando los resultados anteriores, lo que reduce considerablemente la cantidad de tiabajo. En el Capítulo se muestran estos resultados generales para el subespacio de las matrices M tal que AM es simétrica. La principal aplicación se desarrolla para el subespacio de las matrices con un patrón de sparsidad dado el cual puede construirse iterativamente aumentando progresivamente el conjunto de entradas no nulas en cada columna. 3.1. Resultados teóricos Partiendo del sistema de ecuaciones lineales (2.1) donde A es una matriz de orden elevado, sparse y no singular. Como se ha expuesto en el Capítulo anterior (2), en general, la convergencia de los métodos iterativos basados en subespacios de Krylov no está asegurada o puede ser muy lenta. Para mejorar su comportamiento, consideramos una matriz de precondicionamiento M y transformamos el sistema (2.1) en cualquiera de los siguientes problemas equivalentes, MAx = Mb (3.1) AMy = b; X = My (3.2) 60 Inversa aproximada es decir, sistemas precondicionados por la izquierda o por la derecha, respectivamente (ver [87,97]). Usaremos aquí el precondicionamiento por la derecha y M debe seleccionarse tal que AM esté cercana a la identidad en cierto sentido. Esta cercanía puede medirse usando normas matriciales, por ejemplo, aquellas normas inducidas por las normas p (1 < p < oo) de M" o la norma Frobenius. Hemos preferido la norma Frobenius por dos razones. La primera es teórica. La norma Frobenius es prehilbertiana, lo que no ocurre con otras normas usuales como ||-||i,||-||2o||-||ooEsto nos permite usar la Teoría de Aproximación en espacios prehilbertianos. La segtmda razón es de naturaleza computacional. Las columnas de M pueden calcularse y usarse independientemente, es decir, en paralelo [38,41,61], ya que PM-/|||.= ¿||(AM-J)e,||^; e,. = (0,...,0,lo,...,Of (3.3) Así, para im subespacio vectorial S de A^„(M), el problema a resolver es rmn\\AM - I\\F = \\AN - I\\F (3.4) y se considera N un buen precondicionador del sistema (2.1) si es una inversa aproximada de A (con respecto a la norma Frobenius) en el sentido estricto, es decir, \\AN — /||F < 1. En general, esta condición no se satisface. De hecho, para una matriz dada A, siempre es posible encontrar un subespacio S que no contenga ninguna inversa aproximada de A estrictamente. En el siguiente teorema se resume la conveniencia de usar una aproximada inversa como precondicionador [8, 71]. Teorema 29 Sea A,M e Mn{R) tal que \\AM - I\\F < 1Entonces, (i) Mes no singular y M = Zil-AMY .m-O A (ii) AMes definida -positiva donde K2{AM) representa el número de condición de la matriz AM relativo a la norma 2. Como punto de partida, consideremos la norma Frobenius con pesos, \fA e Mn(M-) •• \\A\\w = tr{AWA''), {W simétrica y definida positiva) Sea K c {1,2,..., n} x {1,2,..., n} el patrón de sparsidad de M. Definimos el espacio de matrices con esa estructura sparse como, y S = {MG Mnm/rriij = 0; V(2, j) ^ K} Resultados teóricos 61^ Axelsson [8] estudia el problema precondicionado por la izquierda min||M^-/||^ = ||iVM-/||^ (3.5) y lo convierte en, {N*AWA').. = {WA')..-MhJ) e K (3.6) \\N*A-I\\l, = tr [(/ - N''A)W] (3.7) De forma similar, se resuelve el problema precondicionado por la derecha rmn||^M-J||^ = ||AiV-J||2^ (3.8) que da lugar a, (A'ANW).. = {A'W)..;V(¿, j) e K (3.9) pAT - /| |2j, = ir [(/ - AN)W] (3.10) En particular, para W = I, iA'AN)ij = aji;^(i,j)&K (3.11) \\AN-I\\l = ntr{AN) (3.12) Nuestro próximo objetivo es generalizar las ecuaciones (3.11) y (3.12) para un subespacio arbitrario «S de MnO^) (sección 3.1.1). A continuación, en la sección 3.2 obtenemos las expresiones explícitas para la matriz N y \\AN — I\\jp de una base dada del subespacio S, caracterizando aquellas que nos dan inversas aproximadas. Finalmente, en la sección 3.3 aplicamos estos resultados a algtmos precondicionadores usuales (diagonal, sparse, simétricos y otros). 3.1.1. Mejor aproximación y producto escalar de Frobenius En general, una norma matricial inducida no es prehilbertíana, ni aún cuando las normas vectoriales relacionadas lo sean. En efecto, esto puede demostrarse con el siguiente contraejemplo para la norma 2, ya que \\A\\, = \\B\\, = \\A + B\\,=^\\A-B\\,^1 62 Inversa aproximada Sin embargo, la norma de Frobenius con pesos trivialmente proviene del producto escalar {A, B)^ = tr{AWB'), VA, B G M„(M) (3.13) y, en particular, para VF = /, ||-||^ es prehilbertiana y, n n {A, B)p = tr{AB') = J]^a¿,¿)¿,-, VA, B e >Í„(M) (3.14) En lo que sigue, el término ortogonalidad se referirá al producto escalar de Frobenius (•, •)p. Entonces, el problema (3.4) definido en el espacio Prehilbertiano (A^„(R), (•, •)p) se reduce a encontrar la proyección ortogonal de la matriz / sobre el subespacio AS. Teorema 30 Sea A e Mn{R) una matriz no singular y S un subespacio vectorial de Mn{R). Entonces, la solución al problema (3.4) será, tr (A^ANM^) = tr (AM) ; VM E S (3.15) y- \\AN - I\\l = ntr{AN) (3.16) Demostración. Como dim(j4<S) = dim(íS) < oo, el teorema de la mejor aproximación en subconjtmtos cerrados y convexos asegura la existencia y vmicidad de la solución al problema (3.4) y el teorema de la proyección ortogonal caracteriza la mejor aproximación AN de / en el subespacio AS por la condición {AN - J, AM)p = 0;\/MeS (3.17) que es equivalente a, tr{A'ANM') = tr{AM) ; ^M e S (3.18) Además, \\AN - ifp = {IAN, I)p + {AN - 7, AN)p - n - tr{AN) El teorema anterior generaliza los restdtados (3.11) y (3.12), obtenidos para el subespacio de las matrices con un patrón de sparsidad dado K C {l,2,...,n}x {1,2,..., n}, en cualquier subespacio vectorial de Mni^)- S-Inversa generalizada y complemento ortogonal de Frobenius 63 3.2. <S-Inversa generalizada y complemento ortogonal de Frobenius 3.2.1. <S-Inversa generalizada En esta sección generalizamos la obtención del mejor precondicionador (en norma Frobenius) en im subespacio arbitrario <S de Mn (K) , es decir, la Sinversa generalizada para la matriz A del sistema lineal (2.1), donde A es una matriz no singular, sparse y no necesariamente simétrica. Sea S xm subespacio vectorial de A^„ (K). Sea una descomposición en valores singulares de la matriz A. Entonces, la solución al problema de optimización (3.4) viene dada por: N = VQU^ } mín ||EP - 7||j. = ||SQ -/||i. (3.19) Procedamos así: \\AN-I\\F = \\J:Q-I\\F (3.20) A = UY:V^- (3.21) A^ = yE[/^ (3.22) AA'^ = UE'^U'^ (3.23) A^A = FEV^ (3.24) Entonces, el problema de minimización (3.4) resulta. 5ngpM-/ii,={^:^f;,} = mín \\Ui:V^VXU^-I xevTsu " = ,,mm^„ \\UEXU^ - j||^ (3.25) xevTsu = mín WU'LV'^VXU'^ -UU^ xevTsu " = mín IISX — JIL xeV^su Así, mínpM-/||f= mín IIEX - JIL (3.26) Mes" " xevTsu ^ ||^iV-J||^ = ||Ey-/i|^ (3.27) 64 Inversa aproximada Luego, A = UEV^ ^ N = VYU^ con mín IISX -/|L = IISF -/|U (3.28) 3.2.2. Complemento ortogonal de Frobenius Resolver el problema de minimización (3.4) consiste en encontrar una matriz AN — I que sea ortogonal a la matriz AM en el sentido de la norma de Frobenius, es decir < AN - /, AM >F^ O, M e S, esto es, AN -le (AS)-^. Proposición 31 V5 C MniR) : {ASy = {A-^f S^ (3.29) En consecuencia, AN-I e {AS)^ = {A'^fs^ ^AN-I^ {A'^fQ; Q e S^ =^ (3.30) A^AN = A^ + Q;QeS^ (3.31) lo que constituye la Ecuación Normal Generalizada. Sea dim(5) = p y{Ei,..., Ep] tma base de S, entonces, {A^AN,Ei)^ = {Á',Ei)^ + {Q,Ei)p;\/i=\,2,...,p (3.32) Por tanto, {A^AN,Ei)p = {Á^,Ei)p,;yi = 1, 2, ...,p (3.33) constituyen las nuevas Ecuaciones Normales Generalizadas. En conclusión, N y \\AN — I\\p pueden determinarse por la ecuación normal a partir de S-^ A continuación presentamos el complemento ortogonal de algunos subespacios particulares. 3.2.2.1. Subespacio de las matrices cuadradas 5 = Mn{R) =^S^ = {0} A'^AN = A'^^N = A-'YA-'eS (3.34) S-Inversa generalizada y complemento ortogonal de Frobenius 65^ 3.2.2.2. Subespacio de las matrices sparse S = {Me MnW/mij = O-MhJ) ^K}, KC{1,2, ..,n} x {1,2, ...,n} Á^AN = Á^ + Q (3.35) {A^AN, Eij)^ = {Á^, Eij)^ ; V(z, j) e K (3.36) ((^^A)iV)..=a,,- V(¿,i)eif (3.37) 3.2.2.3. Subespacio de las matrices simétricas S = 5„(R) =^S^ = 7ín(E) A^AN = A^ + H (3.38) NA^A = A-H (3.39) Sumando ambas ecuaciones tenemos, A^AN + NA^A = A + A^ (3.40) 3.2.2.4. Subespacio de las matrices hemisimétricas S = 7Í„(M) ^S^ = 5„(E) A^AN = A^ + S (3.41) -A^^^A = A + S (3.42) Sumando ambas ecuaciones tenem.os, A^AN + ATA^A = A^ - A (3.43) 3.2.2.5. Subespacio de las matrices que simetrizan a A 5 = |M € Mn{^)/ {AMf = AM^ =^AS = Sn=^ (AS)^ = Hn AN-Ie (AS)^ (3.44) AN -I = H (3.45) AN-I = -H (3.46) Sumando ambas ecuaciones tenemos, 2 {AN -I) = Q=^AN-I = 0=^N^A-^yA-^eS (3.47) 66^ Inversa aproximada 3.1.1.6. Subespacio de las matrices que antisimetrizan a A 5 = {M 6 Mn{^)/ {AMf = -AMJ =^ AS = nn=^ (AS)^ = 5„ AN-Ie [AS)^ (3.48) AN-I = S (3.49) -AN -I = S (3.50) Restando ambas ecuaciones tenemos, 2AN = 0=^N = 0 (3.51) 3.2.2.7. Subespacio de las matrices simétricas que simetrizan a A S=¡M e Mni^)/M^ = My {AMf = AMJ => =^ AS= ¡AM G Mn{R)/M^ = M y {AMf = AM\ = ASn n Sn (3.52) AN-Ie {ASf = {ASn n Snf = {ASnf + {Snf = H^A + Hn (3.53) AN -1 = HiA + H2 (3.54) De otro modo, AN-Ie {ASf = (^<S„ n Snf = {A-^Y Un + Un (3.55) ^iV - / = (^-^)^ Hi + fÍ2 (3.56) A^AAT = Á^^Hx + Á^H2 (3.57) 3.3. Expresiones explícitas para los mejores precondicionadores Otra forma de determinar N y \\AN - I\\p es a partir de las expresiones explícitas usando una base adecuada del subespacio S de precondicionadores. Aplicación a algunos precondicionadores usuales 73^ Demostración. Sea A = UT,V^ la descomposición en valores singulares de A. Usando la ecuación (3.79), \\AN -I\\% = ||Í7E [{V\A + A')V) QP]V^- if F i||2 \U [S {V\A + A')V) 0 P - U'V] V = II (Ey't/E + E^f/V) © P - {E^U'V + U'VE^) 0 P = II [{EV'U - U'VH) S] 0 P||^ = \\[{U* {A - A') U) E] QP i2 \F = E Por lo tanto. (W{A-A)U\jaj af + a] 2 F '' yiUHA-A^)U)l i<j af + a] ¿ ^ (c/* (^ - ^*) í7)^. < ll^iv - Jlll < ¿ E (í^' (^ - ^') u)l y como ^iU^A-A')U)l = l\\A-ArF i<j obtenemos ¿^||yl - A^^ <\\ANI\\l < ^J\A - A% • En el siguiente corolario se resumen algunas consecuencias directas de este teorema Corolario 38 (/) Si 2 \\A\\2 < \\A — -<4^||^, entonces no hay ninguna inversa aproximada simétrica de A, (ii)sik2{A) ~ 1,entonces ¡ \\A - A^\\^ \\A\\-' ~ \\AN - I\\p c^^\\AA^\\^ \\A~%, (üi) i \\A - A^W^ \\A\\-' c. \\AN - I\\^ c^l\\AA^^ \\A-' ||^. El corolario 38 afirma que las cotas en (3.80) están más cercanas cuanto más se acerca el número de condición k2{A) a la unidad. Como un caso particular, si A es ortogonal, entonces \\AN-I\\^ = ^\\A-A^^yN = l{A + Á^) (3.81) Estos resultados teóricos pueden obtenerse también aplicando el teorema 32. Siguiendo el mismo procedimiento que para el caso simétrico, el mejor precondicionador hemisimétrico del sistema (2.1) es N = V[{V^{Á^-A)V)QP]V^ (3.82) 74 Inversa aproximada 3.4.4. Precondicionador M tal que AM sea simétrica Sea S ^ {M E A1„(E)/(AM)^ = AM}. Si la matriz M se factoriza como M = XÁ^{X e Mn(M)) (3.83) directamente se sigue que (AM)'^ = AM si y sólo si X'^ = X. Así, «S = {XA^/X e MnO^),X'^ = X} y, para aplicar el teorema 33 se puede obtener inmediatamente una base de S de la base canónica del subespacio de las matrices simétricas de orden n S = span l^{Ei,iA'^}l^ U {{Eij + Ej,i) A^}.^.) (3.84) donde Eij es la matriz cuadrada de orden n cuya única entrada no nula es Bij = 1. Note que las matrices de la base anterior de S sólo contienen iina fila de A^ o dos filas en posiciones permutadas. Además, partiendo de S' = span [{EÍ^ÍA^}"'^^ y aumentando la base iterativamente, siempre podemos obtener precondicionadores MQ = XQÁ^ tan cercanos a A~^ G S como se requiera, preservando la simetría de ^A^ tal que \\AN - I\\p < \\AA^ - l\\p pues A^ e S'. El mayor interés de este ejemplo radica en la posbilidad de usar el método CG para resolver las ecuaciones normales generalizadas del sistema (2.1), si Xo es simétrica definida positiva, AXoA*y = b; x = XoA^y (3.85) Capítulo 4 Valores propios y valores singulares en la proyección ortogonal de Frobenius Los resultados de la teoría de operadores completamente continuos en espacios de Hilbert determinan relaciones particulares para los valores singulares y autovalores de la mejor aproximación en AS, a la identidad. En particular se establece que el menor valor singular y el módulo del menor valor propio de AN no exceden nunca de la unidad cualquiera sea el subespacio S. A continuación se realiza un análisis de convergencia (número de condición, distribución de autovalores y de valores singulares y grado de normalidad) para el precondicionador óptimo en un subespacio arbitrario <S de A^„(R). Se propone im índice alternativo de la descomposición de Schur para el grado de normalidad 4.1. Menor valor propio y menor valor singular en la proyección ortogonal de Frobenius Recordemos que la mejor aproximación AN a la identidad en el subespacio AS se caracteriza por la condición {AN - I, AM)p = 0, "ÍM eS (4.1) de donde \\AN - I\\% = ntr{AN) (4.2) |2 ||AiV||^ = tr{AN) (4.3) Denotando por {XkYl^-^ Y {'^fc}fc=i ^ los autovalores y los valores singulares de la matriz AN (ordenados en orden decreciente de sus módulos), respectiva- 76^ Valores -propios y valores singulares mente, la ecuación (4.3) se escribe como E^^ = ¿A. (4.4) k=l k=l Es decir, para la mejor aproximación AN a la matriz identidad en el sentido de Frobenius, la suma de los cuadrados de los valores singulares coincide con la suma de los autovalores. Los siguientes resultados establecen con más precisión la relación entre los autovalores y los valores singulares de AN {Ay N reales) Lema 39 Sea S un subespacio vectorial cualquiera de Aí„(E) y mín \\AM - I\\p = \\AN - I\\p ; S C Mn{^)- (4.5) Entonces: n Yl>^l<Y:\>^k? <Y.'^l = Y.>^k<Y.\>^k\<Y.<'k (4.6) *;=! fc=l fc=l fc=l fc=l k=\ Demostración. Aplicando el Teorema de la mayorante de Weyl [106] m m ^\Xkf<J2^k (P>0 ;m=l,2,...,n) (4.7) fc=i fc=i para m = nyp= 1, p = 2yla ecuación (4.4), se concluye la demostración • En particular, si AN es simétrica, entonces la igualdad sustituye a todas las n n desigualdades en (4.6) con la única excepción de ^ A^ < ^ |Afe|. Más aún, si fc=i fc=i AN simétrica y definida positiva, entonces las seis sumas en (4.6) coinciden. Teorema 40 El menor valor singular y el módulo del menor valor propio de AN no son nunca mayores que 1 Demostración. La ecuación \\AN - I\\l = n-ír(AAr) y ellema39 nos Uevan a, n O < IIAiV - 7||^ = n-J2 \^nf < n(l - \\nf) < n(l-a^) (4.8) de donde, O < 1 - |A„|' =^ |AJ < 1 O < 1 - CT2 =^ CTn < 1 • El resultado anterior puede ser utilizado para al interpretar la proximidad a la uiúdad del módulo del menor valor propio (menor valor singular) en términos de la calidad del precondicionador. De paso, proporcionamos una demostración alternativa para el teorema anterior. Menor valor propio y menor valor singular 77^ Teorema 41 | A„| ( o an) determinan por si mismas la calidad del precondicionador N. Con más precisión, lím II^A^ - /||^ = lím \\AN - I\\p = O (4.9) |An|—1 cr„^l Demostración, De la ecuación (4.8), se obtiene directamente |A„| ^ 1 =^ ||AiV -/||^-^ O Una consecuencia inmediata del teorema 40 en relación con el precondicionador N definido por: S={ME Mn{^)/AM = (AM)'^}; mín \\AM - I\\p = \\AN - /||^ nos indica que la generalización de la ecuación normal A^Ay = b; x = A^y definida en el Capítulo 3 (apartado 3.4.4) constituye una mejora real de dicha ecuación. Con más precisión. Corolario 42 (i) Si el menor valor singular de A es mayor que la unidad, entonces, en cualquier subespacio So de Mn(^) que contenga a A^, es posible (y sencillo) encontrar una matriz NQ tal que: WANO-IWP^WAA'^-IW^ (4.10) (ii) Si el menor valor singular de A es mayor que la unidad, entonces, en cualquier subespacio So de S2 que contenga a A^, es posible (y sencillo) encontrar una matriz No tal que: ANQ simétrica y \\ANo~I\\p^\\AA^-l\\p (4.11) Demostración, (i) Sea So C Aí„(R) tal que So3 A^ y consideremos el problema: min^\\AM-I\\^=\\ANo-I\\p (4.12) Procedamos por reducción al absurdo. Si ^A^o € ^o tal que \\ANo — I\\p $ \\AA^ — /||^, entonces, el mínimo (4.10) se alcanza en A^, es decir, NQ = Á^. En virtud del teorema 40 \Xn{ANo)\ < 1 =^ \XniAA^)\ < 1 =^ a„(^) < 1 en contra de la hipótesis. • Con esto hemos probado que en todos aquellos casos en que an{A) > 1, podemos precondicionar el sistema Ax = b con ima matriz NQ ^ Á^ que simetrice a A, es decir, ANQ = (ANo)"^, y que, a la vez, mejore a la propia A^ como inversa aproximada: \\ANo - I\\p % |i^^^ - I\FFinalmente, establecemos tina condición necesaria para que la matriz A tenga una inversa aproximada en el subespacio S. Para este propósito se usan los siguientes resultados (ver [64]). 7S Valores propios y valores singulares Lema 43 (Horn) Sean A,Be Mn{R), y sea f{x) una función no decreciente de variable real en [O, oo) y convexa mediante la sustitución x = e* (—00 <t < 00). Denotamos los valores singulares ordenados de A, B,y AB por ai{A) >•• > (Tn{A) > O, (ji(5) > • • • > an{B) > O, y (Ji{AB) > > an{AB) > O,entonces m m J2fi<^k{AB))<J2fi'^k{A)ak{B)), Vm = l,2,...,n (4.13) A;=l fc=l Esto nos lleva al siguiente teorema. Teorema 44 Sea A G Mn{R) una matriz no singular y sea S cualquier subespacio vectorial de Mn{R)- Si S contiene una inversa aproximada de A entonces, E [1 - ^^ (^) ^l (^0)] < 1 (4.14) fc=i Demostración. Usando el Lema de Horn para /(x) = x^, tenemos WANWl = J24 {AN) <J2^¡ (A) al (N) (4.15) fc=i fc=i y de las ecuaciones (4.2) y (4.3), n 1 > \\AMo - I\\l = nWAMoWl > E [^ " ^^ (^) ^^ (^0)] (4.16) fc=i • 4.2. Análisis de convergencia Cuando M es una inversa aproximada de A, la desigualdad: (ver [8]) basada en el lema de Banach, muestra que el número de condición de AM puede hacerse arbitrariamente próximo a 1 con tal que \\AM — /||^ sea suficientemente pequeña. La proximidad de K2{AM) a la imidad, así como el grado de normañdad de AM y la concentración de sus autovalores y valores singulares en tomo a 1, son condiciones esenciales para la convergencia de la mayoría de los métodos iterativos [60]. A continuación, estudiamos el comportamiento de la matriz AN con respecto a estos cuatro parámetros de convergencia. La distribución de los valores singulares de ^A'^ se resume a continuación. Análisis de convergencia 79 Teorema 45 Sea S un subespacio vectorial cualquiera de Mn{R) y min \\AM - I\\^ = \\AN - /||^ (4.18) Entonces, los valores singulares de AN se encuentran en el intervalo [l-PiV-/||2,l+piV-/y (4.19) y satisfacen la relación n Y,{l-a,f<\\AN-I\\l (4.20) jfc=i Además, si \\AN — /II2 < 1 entonces ^'^^"^ ^ l-\\AN-lt '•^•^'^ Demostración. Como bien es sabido [90], los valores singulares dependen continuamente de sus argumentos y |(7fc - 1| < \\AMo - /II2 , VA: = 1, 2,..., n (4.22) Además, del lema 39 y de la ecuación (4.2) obtenemos, n n n J] (1 - Ukf = 5] (1 + Afc - 2ak) <Y.{l-\k) = \\AN - I\\l (4.23) fc=l fc=l k-l Por otra parte, la distribución de los autovalores de AN se presenta en el siguiente teorema, en el que también se propone y evalúa tm estimador del grado de normalidad de AN. Teorema 46 Sea S un subespacio vectorial cualquiera de Mn(R) y mm\\AM - ly = \\AN - I\\^ (4.24) Entonces los autovalores deANse encuentran dentro de un círculo de radio \ | AN — /11 ^ y centro 1 y satisfacen n Y,\í-h\'<\\AMo-I\\l (4.25) A;=l Además, I ¿ (|A,| - a,)' < \ \\AM,f^ (1 - Cn) (4.26) "til so Valores propios y valores singulares Demostración. De la inecuación de Weyl, el lema 39 y la ecuación (4.2) se sigue ^\1Xkf = ¿ (1 + \hf - 2Xk) < 5^ (1 - Afc) = \\AMo - I\ \l fc=l k=l k=l estimación obtenida por Grote y otros [60] usando una descomposición de Schur de la matriz AN — I. Por otro lado, de la inecuación de Weyl y las ecuaciones (4.3) y (4.4), obtenemos E (lAfcl - (7kf = ¿ (iA,|' + Afc - 2 |A,| (Tfc) < 2 ¿ (Afe - |Afc| ak) k=l fc=l fc=l < 2 ÍE Afc) (1 - ^n) = 2 II^Moll^ (1 - an) m n La expresión ^ E (I Ajt I — o"*;) se propone aquí como im estimador del grado k=l de normalidad de AN basándonos en la conocida caracterización de los operadores Hneales normales en función de sus valores propios y singulares [55]. Capítulo 5 Inversa aproximada sparse Actualmente, los métodos basados en subespacios de Krylov son las herramientas más eficientes para resolver sistemas de ecuaciones lineales (2.1). Sin embargo, estos métodos deben usarse con precondicionadores muy a menudo. Recientemente, el uso de inversas aproximadas se ha convertido en tina buena alternativa para los precondicionadores implícitos debido a su naturaleza paralelizable; para más detalles ver, [20,38,41,61], y también vm estudio comparativo completo en [24]. En este Capítulo, se construye una matriz inversa aproximada usando el producto escalar de Frobenius (sección 5.1). Aunque el estudio se ha desarrollado para inversas aproximadas precondicionando por la derecha, los resultados por la izquierda son directos. De hecho, la mayoría de los experimentos en este capítulo se han ejecutado precondicionando por la izquierda. En la sección 5.2, se propone un Algoritmo de la Inversa Aproximada mejorada (lAI) no sólo para obtener un mejor precondicionador a partir de una inversa aproximada dada y usarlo combinado con un m.étodo iterativo basado en los subespacios de Krylov, sino incluso para resolver (2.1). El apartado 5.2.1 es un resumen de algunas propiedades teóricas de la inversa aproximada sparse y la mejorada. Finalmente, en la sección 5.3 se ilustra la eficiencia de estos precondicionadores y del algoritmo lAI con algunos experimentos numéricos. 5.1. Cálculo de inversas aproximadas sparse Sea S C Mn, él subespacio de las matrices M donde se busca una inversa aproximada explícita con un patrón de sparsidad desconocido. La formulación del problema es: encontrar MQ E S tal que Mo = argmín \\AM - J||^ (5.1) Mes Además, esta matriz irúcial MQ pudiera ser posiblemente una aproximada inversa de A en un sentido estricto, es decir, según la definición dada por [71], II^Mo - /||^ = £ < 1 (5.2) 82 Inversa aproximada sparse Existen dos razones para esto. Primera, la ecuación (5.2) permite asegurar que Mo es no singular (lema de Banach), y segimda, esta será la base para construir un algoritmo explícito para mejorar MQ y resolver la ecuación (2.1). Un trabajo reciente de Grote y otros [61] nos da un algoritmo eficiente para obtener una inversa aproximada tan cercana a una matriz no singular A como se requiera. Hemos seguido dicha técnica pero variando el método de selección de las entradas en MQ y el algoritmo usado para resolver el problema (5.1). La construcción de MQ se realiza en paralelo, independizando el cálculo de cada columna. Aunque nuestro algoritmo nos permite comenzar desde cualquier entrada de la columna k, se acepta comúnmente el uso de la diagonal como primera aproximación. Además, la expresión del precondicionador diagonal óptimo es bien conocida. La siguiente entrada a considerar se selecciona dentro del conjunto de entradas candidatas, el cual se define siguiendo el criterio propuesto por Grote y otros [61]. Sea r^ el residuo correspondiente a la columna fc-ésüna, Tfc = Anik - Ck (5.3) y sea Ik el conjunto de índices de las entradas no nulas en rk, es decir, J^ — {i e {1,2, ...,n} /rik i^ 0}.Si£fc = {/ € {1,2, ...,n} /m/fc 7^ O}, entonces la nueva entrada se busca en el conjunto Jk = {3 e £% / aij ^ 0,\/i elk}- En realidad, las únicas entradas consideradas en rrik son aquellas que afectan las entradas no nulas de r^. En lo que sigue, asumimos que Ck U {j} — {i\, i\., •••,ip^} es no vacío, siendo pk el número actual de entradas no nulas de nik, y que ip^ = j, para todo j E JkPara cada j, calculamos, donde, para todo k, det (GQ) = 1 y Gf es la matriz de Gram de las columnas ¿1,2^,..., ^f de la matriz A con respecto al producto escalar euclídeo, Df es la matriz que resulta de reemplazar la última fila de la matriz Gf por akik,akik,...,akik, con 1 < I < PkSe selecciona el índice jk que minimiza el valor de \\Amk — 6^112Esta estrategia (ver [56]) define el nuevo índice seleccionado jk atendiendo solamente al conjtinto Ck, lo que nos lleva a ixn nuevo óptimo donde se actualizan todas las entradas correspondientes a los índices de CkEsto mejora el criterio de [61] donde el nuevo índice se selecciona manteniendo las entradas correspondientes a los índices de CkAsí rrik se busca en el conjunto Sk = {mk e W/rriik = 0; Vz ^ £¿ U {jk}} Vdet(Z?f) „ ._. ¿^ det {Gl,) det (Gf) Experimentos numéricos 89 ^í'^:''""!r5^=<Z'-'""~ :•), - :%-. . .: ^^i«v;^ ; ,'i ".í. ••":.: !• m^;':.^ 11! S 1 ];ír¿f m- • _ ^-'—..-^' • 1 m^. á:r '\i ••'• % ^'"' '!! ;if lí.J^l - , . .*: ': Hi: .^ fe. 5 • ít - 'S Wk •.i^jr. ^-ij^^jií." i^pin ^.:S. 1 1- ^ ..i .k .'; "">.', •" • • • ".."^"«^.i. • i '..hi . •": .,.-3^ ^ -m ^^ :^ .:H;»JS •-fifí M. 1 ^' i k ^ • ..• . n •. -*c;""! :'•" .: ' itf. •f III '^!á «- i» - % • 'rrt: .- I' ' "^ •••-; • : k.. _.. •"-,.; i^ l( ^ .í-= h '**ii'. -- % 1 rr :: •• EÍSiSy., '"."" *•*>- Jl ' • ji- .J"*"*" i « » iit " ^ i ^i^. — \ Figura 5.3: Estructura sparse de la inversa aproximada de la matriz convdifhor con £fc = 0,05 y máx nz{mok) = 50. s-:tí • _jr^. ^_.;':T'Í^ IIII (I ^- l||.»r.|^,.. !fW". íTi r l...iM ;.i ^ •••'!!!• lili ll'^l ^íí=;= ^j^-_ -f^- 5^-=^=—^^^^=^- -'^—• íi^l—; '!."??!" "^^^—•"•íS '' " ':.:.:'!'flR"iiiiii''i Figura 5.4: Estructura sparse de la inversa aproximada de la matriz convdifhor con £k = 0,05 y máx nz{mok) = 200. 90 Inversa aproximada sparse £k 0,6 0,4 0,3 0,2 k 100 100 85 58 Iter. 31 2 1 1 nz{Mo) 3087 10353 11025 27964 nz(Mo) nz(A] 0,21 0,73 0,78 1,98 \\MoA-iy 26,76 13,66 11,44 8,89 Tabla 5.2: Resultados de convergencia para oilgen con GMRES precondicionado por la izquierda. Bk = 0,6 £k = 0,5 Sk = 0,4 Sk = 0,3 £k = 0,2 Mo = 7 MQ-^ = ILU{Q) Iter. 122 73 50 37 25 408 38 nz^Mo) 1297 1880 2709 4294 10072 1000 3750 nz(Mo) nz(A) 0,35 0,50 0,72 1,14 2,69 0,26 1,00 ||Mo^-J||^ 13,33 11,53 9,39 7,38 5,03 68,01 Tabla 5.3: Resultados de convergencia para sherman con QMRCGSTAB precondicionado por la izquierda. £k = 0,5 Sk = 0,4 £k = 0,3 Sk = 0,2 MQ-' = ILUiOi) Iter. 795 656 270 142 218 nz{Mo) 19227 39956 78003 209186 86562 nz(Mo) nz(A) 0,22 0,46 0,90 2,42 1,00 ||MoA-J||^ 51,36 41,54 31,64 21,64 Tabla 5.4: Resultados de convergencia para isla con BiCGSTAB precondicionado por la izquierda. BiCGSTAB QMRCGSTAB GMRES(6) Left-IAI Right-IAI pores Iter. 5 5 2 3 5 convdifhor Iter. 7 7 2 4 6 Tabla 5.5: Comparación de los resultados de convergencia para pores con métodos de Krylov e Mí. Experimentos numéricos 91 -2 ^^• £ -4 -10 20 1 1 1 \ t\%. "' '-^ •n \\^ ; i \^ i ^ ^ \; "n ) h •r i '^ -' r 1 N A/^ ••"•• , i --rx • • ?; i ••-. •, 1 1 "Unprec." "ILUO" "Tol=0.05 nz=200" •'Tol=0.05" - "Tol=0.2" •'Tol=0.3" -•-• "Tol=0.4" •Tol=0.5" ' - 1 1 40 60 80 Number of iterations 100 120 Figura 5.5: Comportamiento de los prencondicionadores con BiCGSTAB para convdifhor. —1 "Unprec." "ILUO" "Tol=0.6" "Tol=0.5" "Tol=0.4" "Tol=0.3" "Tol=0.2" -10 li m y It 1 F 1 il i^ % H 200 400 600 Number of iterations 800 1000 1200 Figura 5.6: Comportamiento de los prencondicionadores con BiCGSTAB para isla. Capítulo 6 Efecto de la reordenación en precondicionadores del tipo inversa aproximada sparse El objetivo de este Capítulo es estudiar el efecto de la reordenación sobre nuestra versión del precondicionador SPAI (inversa aproximada sparse) propuesta por Grote y otros [61], cuyos aspectos tanto teóricos como computacionales se han analizado en [56,84]. Presentamos los resultados sobre el efecto de la reordenación no sólo en la cantidad de entiadas en los factores de las inversas, sino también en el número de pasos del método iterativo (ver [49]). Aunque la inversa A"'^ está normalmente llena, independientemente de la reordenación seleccionada, mostramos experimentalmente cómo el etecto fill-in de las inversas aproximadas sparse depende de la ordenacüon de A. Benzi y otros [25] han realizado un estudio similar para inversas aproximadas factorizadas. También es un punto de partida interesante, el artículo de Bervzi, Szyld and van Duin [21] sobre el efecto de la reordenación para la factorización incompleta en la convergencia de los métodos basados en subespacios de Krylovs no sim.étricos. En primer lugar, se resume el algoritmo de la inversa aproximada sparse en la sección 6.1. En la sección 6.2, se discuten algimas consideraciones sobre las técnicas de reordenación tales como Grado Mínimo, Mínimo Vecino, Reverse Cuthill-Mckee, [1], [54]. Finalmente, se resuelven algxmos experimentos numéricos para mostrar el efecto de los algoritmos de reordenación sobre la convergencia del método de Bi-CGSTAB [104] para la solución de sistemas de ecuaciones lineales no simétricos donde se utilizan dichas inversas aproximadas como precondicionadores. Se presentan los resultados de los sistemas cuyas matrices peretenecen a la colección Harwell-Boeing [43] y otios que provienen de la discretización por el Método de elementos finitos de diferentes problemas de contomo. 94 Efecto de la reordenación 6.1. Algoritmo de la inversa aproximada sparse Sean rl = m].A—e\. el residuo correspondiente a la fila A; de M e J^ el conjunto de índices de entradas no nulas en r|., es decir, J^ — {ie {l,2,...,n} / Tki ^ 0}. Si £fc = {Z G {1,2,..., n} / TTiik ^ 0}, entonces la nueva entrada se busca en el conjunto Jk — {j e £^ / Oj, 7^ O, V? G Xfc}. Así, las únicas entradas consideradas en mi son aquellas que afectan las entradas no nulas en r|.. Asumimos que -Cfc U {j} = {¿1, ¿2, •••, íp^} es no vacío, con pk el número actual de entradas no nulas de m^. y i^^ = j, para todo j & JkPara cada j, calculamos Pk irr.'.l c'll^-1 V [d'^t(A')]' ,m,A e,|k-l ^2^<j,^(g._^),^,(g,) (6.1) '1! '2' donde, para todo A;, det (GQ) = 1 y Gf es la matriz de Gram de las filas de la matriz A con respecto al producto escalar euclídeo, Df resulta de reemplazar la última fila de la matriz G\ por aikk,aikk,:.,aikk, con 1 < I < pkSe selecciona el índice jk que minimiza el valor de \\ml.A — e^W^. Esta estrategia define el nuevo índice seleccionado jk atendiendo solamente al conjunto £fc, lo que nos lleva a un óptimo donde todas las entradas correspondientes a los índices de Ck se actualizan. Así mi se busca en el conjunto Sk = {mi e W/mik = 0; V¿ ^ £fc U {jk}}, y ^ det (Pf) 'mk= > / ^ \ r ^x ^ (6.2) .¿^det(Gf_Odet(Gf) ' donde fh\ es el vector de entradas no nulas z| (1 < /i < /)• Cada tina de ellas se obtiene evaluando el determinante correspondiente que restdta de reemplazar la última fila en det (Gf) por e^, con I <l <pk6.2. Algunos comentarios sobre reordenación Hemos considerado diferentes técnicas de reordenación para mostrar el efecto de la reordenación sobre la solución iterativa de sistem^as de ecuaciones lineales usando los precondicionadores SPAI. La ordenación original corresponde a las matrices provenientes de la aplicación del Método de Elementos Finitos con maUas no estructuradas y refinamiento adaptativo de mallas. Los algoritmos de reordenación han sido expuestos en el capítulo 2, sección 2.5. Ahora nuestro principal objetivo es investigar si la reordenación reduce la cantidad de entradas en el precondicionador SPAI y si el uso de los precondicionadores SPAI m.ejora la convergencia del método iterativo. Sea P la matriz de permutación asociada a un algoritmo de reordenación. Como {P'^AP)~'^ = P'^A'^P, es decir, la inversa de la matriz reordenada es la Algunos comentarios sobre reordenación 95 reordenada de la matriz inversa, cuando reordenamos una matriz, su inversa aproximada debe tender a la reordenada de la inversa. Si la tolerancia de la inversa aproximada está dada por e, en el subespacio S C M„(M), mín \\MA - I\\^ = \\NA - I\\^ < e (6.3) entonces, mín WM'P^AP - /|L = mín \\PM'P^A - /IL = mín \\MA - I\\ < e (6.4) Sea S' un subespacio de M„(E) correspondiente al mismo número de entradas no nulas que S, donde la inversa aproxim^ada óptima es obtienida para dicho número de entradas no nulas. Cabe destacar, también, que el número de entradas no nulas en S es el mismo que el de P^SP. En este caso obtenemos. N'P^AP - /IL = mín WM'P'^AP - /ÍL < mín ¡¡M'P'^AP - /IL < £ (6.5) Evidentemente, el número de entradas no nulas necesarias en S' será menor o igual que el número de entradas no nulas en P'^SP y en S. Concluimos que la reordenación reduce la cantidad de entradas no nulas en la inversa aproximada para una tolerancia dada s, o, al menos, no lo aumenta. Debido a los resultados dados en (6.5), los precondicionadores tipo inversa aproximada reordenados adquieren mejores propiedades desde el punto de vista de su utUización para mejorar la convergencia de los métodos iterativos. La cercanía del número de condición de M'P'^AP a 1 se caracteriza por. , , r . 1+\\M'P'^AP-I\l MM'P-AP)<^_\^,^,^^_¡¡^^ (6.6) La desviación de la normalidad de {M'P'^AP) está acotada por, i¿(|A,| - a,f < - WM'P'^APt (1 - a„) (6.7) siendo {Xk}^^^, {crk}^^i los autovalores y valores singulares de P'^APM' (en secuencia de módulos decreciente). Finalmente, la agrupación de los autovalores y valores singulares estará dada por. ¿(l-cT,)^<||M'P^^P-/||^ (6.8) k=l n Y,{1 - \k? < WM'P'^AP - /||^ (6.9) fc=i 96 Efecto de la reordenación 6.3. Experimentos numéricos Comenzamos con el estudio de un problema cuya matriz pertenece a la colección Harwell-Boeing: orsregl. Ésta es xma matriz de simulación de im almacenamiento de petróleo para una malla de21x21x5de tamaño n = 2205 y entradas no nulas nz = 14133. Las tablas 6.1, 6.2 y 6.3 muestran los resultados obtenidos con el algoritmo BiCGSTAB precondicionado para la ordenación Original, Grado Mínimo y Cuthill-Mckee Inverso de la matriz A, respectivamente. El comportamiento de ILU(O) es comparado con varios precondicionadores SPAI correspondientes a diferentes niveles de llenado. Presentamos el número de iteraciones, el número de entradas de M y la norma de Frobenius de la matriz residuo. Los valores de nz para los precondicionadores SPAI en los casos reordenados llegan a ser menores que el caso de la ordenación original a partir de Sk = 0,2. Sin enibargo, en los otros casos al menos la norma de la matriz residuo se reduce con la reordenación, y por tanto, no contradicen los resultados teóricos anteriores. Por otro lado, el número de iteraciones del BiCGSTAB siempre disminuyó cuando usamos reordenación, excepto para el primer SPAI, que dio inestabilidad a la convergencia del algoritmo. También observamos en estos experimentos que, con requerimientos de almacenamiento similares a ILU(O), se produce una convergencia más rápida con SPAI Sk = 0,3, o incluso ek = 0,2 si se reordena con Cuthill-Mckee Inverso). La figura 6.1 representa el comportamiento del algoritmo BiCGSTAB precondicionado cuando ILU(O) y SPAI(0.3) se construyen después de la reordenación. La principal conclusión es que SPAI puede competir con ILU en máquinas paralelas. Además, si aplicamos un algoritmo de reordenación adecuado a las dos estrategias, esta competitividad se mantiene. Precondicionador Sin precondicionar ILU(O) SPAI £fc = 0,6 SPAI Sk = 0,5 SPAI Sk = 0,4 SPAI £fe = 0,3 SPAI £fc = 0,2 Iter. 548 47 299 169 86 59 37 nz{M) 2205 14133 3087 6615 10353 11025 31782 nz{M)/nz{A) 0,16 1,00 0,22 0,47 0,73 0,78 2,25 \\MA-I\\F 854454 — 26,77 22,31 13,67 11,45 8,56 Tabla 6.1: Resultados de convergencia para orsregl con el orden original y BiCGSTAB precondicionado por la izquierda. Experimentos numéricos 97 Precondicionador Sin precondicionar ILU(O) SPAI Sk = 0,6 SPAI £k = 0,5 SPAI £fc = 0,4 SPAI Ek = 0,3 SPAI Ek = 0,2 Iter. >2205 381 1396 73 36 21 14 nz{M) 2205 14133 3593 7199 10678 13850 19694 nz{M)/nz{A) 0,16 1,00 0,25 0,51 0,76 0,98 1,39 \\MA-I\\p 854454 — 25,62 20,93 13,02 8,24 5,36 Tabla 6.2: Resultados de convergencia -para orsregl con Grado Mínimo y BiCGSTAB precondicionado por la izquierda. Precondicionador Sin precondicionar ILU(O) SPAI Sk = 0,6 SPAI £k = 0,5 SPAI £k = 0,4 SPAI Ek = 0,3 SP^I Ek = 0,2 Iter. >2205 26 810 88 27 19 12 nz{M) 2205 14133 4152 7005 9304 10608 13322 nz{M)/nz{A) 0,16 1,00 0,29 0,50 0,66 0,75 0,94 \\MA-I\\F 854454 — 23,84 19,57 11,62 8,08 5,74 Tabla 6.3: Resultados de convergencia para orsregl con Reverse CutHill McKee y BiCGSTAB precondicionado por la izquierda. El segtindo ejemplo (convdifhor) es un problema de convección difusión definido en [O, i] X [0,1] por la ecuación. du ViK ox = F donde Í;I = 10^ {y - |) (a; - x^) (| - x), K = 10-^-10^,7 F = 10^-1. La matriz corresponde a una malla no estructurada de elementos finitos con n = 1960 y nz = 13412. Las tablas 6.4, 6.5, 6.6 y 6.7 son similares a las del problema anterior. La reducción del número de entradas en la SPAI es del 40-50 % para los casos de reordenación con Grado Mínimo y Cuthill-Mckee Inverso. El algoritmo de Mílúmo Vecino no afecta a nz. Además, el número de iteraciones del BiCGSTAB se redujo mediante las renumericaciones del 60-70 %. Como estamos interesados en el efecto de la reordenación de A en las características de los precondicionadores SPAI, en las figuras 6.3,6.4,6.5 y 6.6 se muestra la estructura sparse de la matriz M con Sk = 0,3 para ordenación Original, Grado Mínimo, CuthillMckee Inverso y Mínimo Vecino, respectivamente. Las entradas no nulas se representan por un ptmto. La estructuras correspondientes representan matri- 98 Efecto de la reordenación -6 -10 -12 I V'-,'. \\ ; . .''. "ii \\ "•-..•••.. tí \ \ • \ \ \ Mi-. \ ; \ i 1 Á,....A-------. -\ v^^^^ . n 1 I ILU(O) SPAÍ(0.3) -— ILU(0)-mdg SPAI(0.3)-mdg " ILU(0)-rcm ; SPAI(0.3)-rcm - - - 1 1 20 40 60 Iterations 80 100 120 Figura 6.1: Comparación del comportamiento de BiCGSTAB-ILU(O) y BiCGSTABSPAI con reordenación para orsregl. ees llenas, como cabía esperar. Sin embargo, se advierte cierto paralelismo con la estructura de A para las diferentes reordenaciones. La reducción del ancho de banda y del perfil llevada a cabo por el algoritmo de Cuthill-Mckee Inverso en la matriz A se conserva de alguna manera en la matriz M, incluso cuando existe una tendencia a explotar algunas entradas fuera del perfil. Esto queda ilustrado claramente en la figura 6.5. Los patrones de las matrices SPAI correspondientes a Grado Mínimo y Mínimo Vecino también conservan en parte las estructuras de la matriz A reordenada, respectivamente, aún cuando nuestro algoritmo SPAI no tiene por qué producir matrices M con estructura simétrica. En la figura 6.2 se compara la convergencia de BiCGSTAB-SPAI(0.2) para todos estas reordenaciones. Claramente, las reordenaciones producidas por los algoritmos de Grado Mínimo y Cuthill-Mckee Inverso tiene un efecto beneficioso en la velocidad de convergencia del algoritmo BiCGSTAB-SPAI. Experimentos numéricos 99 Precondicionador Sin precondicionar ILU(O) SPAI Ek = 0,6 SPAl£fc = 0,4 SPAI Ek = 0,3 SPAI £fc = 0,2 SPAI £fc = 0,1 Iter > 1960 74 414 302 171 83 21 nz{M) 1960 13412 3161 10693 21734 54406 167678 nz{M)/nz{A) 0,15 1,00 0,24 0,80 1,62 4,06 12,50 \\MA-I\\F 13279 — 22,42 16,99 12,78 8,70 4,36 Tabla 6.4: Resultados de convergencia para convdifhor con el orden originial y BiCGSTAB precondicionado por la izquierda. Preconditioner Unprecond. ILU(O) SPAI Ek = 0,6 SPAI Ek = 0,4 SPAI Ek = 0,3 SPAI Ek = 0,2 SPAI £fc = 0,1 Iter > 1960 57 166 99 68 40 21 nz{M) 1960 13412 2617 6255 11461 26992 92864 nz{M)/nz{A) 0,15 1,00 0,20 0,47 0,85 2,01 6,92 \\MA-I\\F 13279 — 19,82 15,18 11,42 7,78 3,85 Tabla 6.5: Resultados de convergencia convdifhor con Grado Mínimo y BiCGSTAB precondicionado por la izquierda. Preconditioner Unprecond. ILU(O) SPAI Ek = 0,6 SPAI Ek = 0,4 SPAI Ek = 0,3 SPAI £fc = 0,2 SPAI £fc = 0,1 Iter 1477 31 144 92 66 41 18 nz{M) 1960 13412 2510 6126 11355 26270 88093 nz{M)/nz{A) 0,15 1,00 0,19 0,46 0,85 1,96 6,57 \\MA-I\\F 13279 — 19,51 15,51 11,67 7,98 4,01 Tabla 6.6: Resultados de convergencia para convdifhor con Reverse CutHill McKee y BiCGSTAB precondicionado por la izquierda. 106 Efecto de la reordenación ^¿ vvO^'v*'. >NSSSS;'.- •^^^ ^S^;: :-í?ÍWS-. •"••'N^'»-» • ^'^'^'^íív-s. • • ^^^w^v . • •'»^w^;^ • •~'O0cC-. ^'>5ís>; •^*^: " **'Vtfi'.. . •..•ííss> .. "^^t •^55^. ^ ^ Figura 6.10: Patrón de sparsidad de la matriz SPAI(0.3) reordenada con Reverse Cuthill-Mckee para cuaref. :../' rX.-. y :.•. ..... - - 'X; - ... •'• ' VI . •" ' ^ .? - ^^ Figura 6.11: Patrón de sparsidad de la matriz SPA1(0.3) reordenada con Mínimo Vecino para cuaref. Experimentos numéricos 107 Concluimos, como en el problema anterior y en otros llevados a cabo no incluidos aquí, que la estructura sparse de SPAI parece partir de una estructura similar a la típica de A obtenida de la reordenación, y tiende a una matriz llena a medida que aumentamos la precisión. La figura 6.7 no aporta diferencias significativas a la figura 6.2 del segundo problema. Los algoritmos de Grado Mínimo y de Cuthill-Mckee Inverso son preferibles al Mínimo Vecino o a la ordenación Original. Sin embargo, hemos notado que si nz aumenta, las diferencias entre Grado Mínimo y Cuthill-Mckee Inverso son más apreciables a favor del primero. Capítulo 7 Conclusiones y líneas futuras El carácter prehilbertiano de la norma matricial de Frobenius permite obtener en el Capítulo 3 el mejor precondicionador N para el sistema (2.1) en cada subespacio S de A<„(E) mediante la proyección ortogonal de la identidad sobre el subespacio AS. Consecuentemente, se obtiene la ecuación (3.4), que representa la distancia mínim.a d(I, AS). Estas consideraciones hacen posible el desarrollo de expresiones explícitas tanto para A'^ como para | l^áAT—/| |i? usando una base ortogonal de AS. Posteriormente, el método de ortogonalización de Gram-Schmidt generaliza estas fórmulas para cualquier base de S. Además, el coste computacional de las últimas expresiones se reduce considerablemente usando la descomposición de AS como suma directa de espacios mutuamente ortogonales ASj. La aplicación de los resultados anteriores al caso de los precondicionadores sparse nos Ueva evidentemente a expresiones cuyos cálculos son inherentemente paralelos, puesto que las columnas Uj de N pueden obtenerse independientemente unas de otras. Asimismo, el patrón de sparsidad del precondicionador N se captura automáticamente (no es un patrón de sparsidad fijo) puesto que cada nueva entrada n¿j es equivalente a extender el subespacio precondicionador <S a S®span {Mjj}. Por otra parte, el mejor precondicionador simétrico N se obtiene, de forma alternativa, de la condición AN-I e {ASn)^ = HnA Entonces, usando la descomposición en valores singulares de la matriz A, se obtienen las cotas superiores e inferiores de H^A'^ — I\\F • Estas cotas dependen de 11^ — ^*| liT y están más cercanas cuanto más cercano esté el número de condición «2 {A) a la unidad. Las diferentes expresiones obtenidas para \\AN — I\\F directamente, caracterizan aquellas matrices A para las cuales el subespacio precondicionador S contiene inversa aproximada en el sentido estricto, es decir, \\AN — I\\F < 1El algoritmo propuesto en el Capítulo 5 permite calciolar tm precondicionador sparse MQ de una matriz sparse no simétrica A. Como las columnas (o filas) lio Conclusiones y líneas futuras de MQ se obtienen independientemente unas de otras, los cálculos pueden realizarse en paralelo. El patrón de sparsidad de estos precondicionadores se construye dinámicamente partiendo del diagonal aumentando el número de entradas no nulas. Evidentemente, este precondicionador puede ser competitivo, si se trabaja en paralelo, con los tradicionales precondicionadores implícitos. Esto dependerá en gran medida del problema y de la arquitectura del ordenador. No obstante, se ha demostrado de forma teórica y práctica la eficacia de este precondicionador en la m^ejora la convergencia de los m.étodos iterativos. Asimismo, el algoritmo lAI propuesto calcula un precondicionador mejor a partir de una inversa aproximada dada. Sólo es necesario calcular productos matriz-vector, por lo que éste puede ser usado directamente para mejorar la efectividad de la aproximada inversa en la convergencia de los métodos iterativos. Se muestran los resultados teóricos sobre esta mejora. Finalmente, el algoritmo JAI permite resolver sistemas lineales en un número de pasos conocidos a priori para una tolerancia prefijada. La cantidad necesaria de productos matriz-vector depende directamente de la torerancia y de la calidad de la inversa aproximada inicial, pero ntmca del número de incógnitas, lo cual es una propiedad bastante interesante. Se ha estudiado por muchos autores el efecto de la reordenación sobre la convergencia de los métodos basados en subespacios de Krylov con precondicionamiento para resolver sistemas de ecuaciones lineales no simétricos. Este efecto beneficioso se muestra en, [96], [21] para precondicionadores basados en factorización incompleta. En el Capítulo 6 hemos probado experimentalmente que las técnicas de reordenación tienen efectos beneficiosos en la ejecución de inversas aproximadas sparse utilizadas como precondicionadores en los métodos iterativos basados en subespacios de Krylov. La reducción del número de entradas no nulas debido a la reordenación permite obtener inversas aproxim.adas sparse con una exactitud similar a las obtenidas sin reordenación, pero con menores requerimientos de almacenamiento y coste computacional. Además, la reordenación produce precondicionadores con mejores cualidades puesto que generalmente se reduce el número de pasos para alcanzar la convergencia de los métodos iterativos. Los experimentos numéricos, pues, parecen indicar que un orden adecuado puede mejorar la eficiencia de las inversas aproximadas aquí propuestas. Sin embargo, deben realizarse más investigaciones sobre la reducción de las entradas no nulas en MQ para un número similar de iteraciones de los métodos de Krylov. También sería interesante estudiar el comportamiento de diferentes algoritmos de reordenación para un sistema de ecuaciones dado, proveniente de la discretización de problemas similares (por ejemplo, nos hemos centrado aquí en la ecuación de convección-difusión). Proponemos la investigación del efecto de otras técnicas de reordenación que tengan en cuenta las entradas numéricas de A (ver [79], [83]). Aún cuando estas técnicas son caras, en el caso en que haya que resolver muchos sistemas de ecuaciones lineales con la misma matriz, dichas técnicas pueden ser competitivas en máquinas en paralelo. Bibliografía [1] P. ALMEIDA, Solución directa de sistemas sparse mediante grafos, Departmento de Matemáticas, Universidad de Las Palmas de Gran Canaria, España, 1990. [2] F. ALVARADO Y H. DAG, Sparsified and incomplete sparse factored inverse preconditíoners, in: Proceedings ofthe 1992 Copper Mountain conference on Iterative Methods, 1, 9-14, april 1992. [3] W.E. ARNOLDI, The Principie of Minimized Iteration in the Solution of the Matrix Eingenvalue Problem, Quart. Appl. Math., 9, pp.17-29 (1951). [4] E. ASPLUND, Inverses of matrices {a¿j} which satisfy a¿j = O íor j > i + p, Mathematica Scandinavia,?,57-60,1959. [5] O. AXELSSON, Incomplete block-matrix factorization preconditioning methods. The ultímate answer?, /. Comp. Appl. Math., 12,13, 3-18,1985. [6] O. AxELSSON, A Restarted Versión of a Generalized Preconditioned Conjúgate Gradient Method, Comunications in Applied Numerical Methods, 4, 521-530,1988. [7] O. AXELSSON, A survey of preconditioned iterative methods for linear Systems of algébrale equations, BIT, 25,166-187,1995. [8] O. AXELSSON Iterative Solution Methods, Cambridge University Press, Cambridge, U.K., 1994. [9] O. AXELSSON, S. BRINKKEMPER Y V.P. IL'IN, On some versions of incomplete block-matrix factorization iterative methods, Lin. Alg. Appl, 58,3-15, 1984. [10] O. AXELSSON Y I. KAPORIN, Error norm estimation and stopping criteria in preconditioned conjugated gradient iterations, Num. Lin. Alg. Appl, 8, 265-286,2001. [11] O. AXELSSON Y L. Yu. KOLOTILINA, Diagonally compensated reduction and related preconditioning methods, Num. Lin. Alg. Appl, 1, 155-177, 1994. 112 Bibliografía [12] M. BENSON, Iterativa solution oflarge scale linear systems, Master's thesis, Lakehead University, Thunder Bay, Canadá, 1973. [13] M.W. BENSON, P.O. FREDERICKSON, Iterative solution of large sparse linear systems arising in certain multidünensional approximation problems, militas Math., 22,127-140,1982. [14] M. BENZI, A direct row-projection methodfor sparse linear systems. Department of Mathematics. Norh Carolina State University, Raleigh, NC, 1993. [15] M. BENZI, J.K. CULLUM Y M. TÚMA, Robust approximate tnverse preconditioning for the conjúgate gradient m.ethod. Los Alamos Technical Report LA-UR-99-2899,1-15,1999. [16] M. BENZI Y G.H. GOLUB, Botmds for the enfries of matrix functions with applicatíons to preconditioning, BIT, 39,3, 417-438,1999. [17] M. BENZI, J.C. HAWS Y M. TÚMA, Preconditioning highly indefinite and nonsummetric matrices. Los Alamos Technical Report LA-UR-994857, 1-25, 1999. ^ [18] M. BENZI, W. JOUBERT Y G. MATEESCU, Numerical experiments with paraUel orderings for ILU preconditioners, Elect., Transact., Num. Anal, 8,88114,1999. [19] M. BENZI, J. MARÍN, M. TÚMA, A two -level parallel preconditioner based on sparse approximate inverses, Iter. Meth. Sci. Comp. II, D.R., Kincaid et. al. editores, 1-11,1999. [20] M. BENZI, C.D. MEYER Y M. TÚMA, A Sparse Approximate Inverse Preconditioner for the Conjúgate Gradient Method, SIAM J. Sci. Comput., 17, 5,1135-1149,1996. [21] M. BENZI, D.B. SZYLD Y A. VAN DUIN, Orderings for incomplete factorization preconditioning of nonsymmetric problems, SIAM}. Sci. Comput., 20,5,1652-1670,1999. [22] M. BENZI Y M. TUMA, A sparse approximate inverse preconditioner for nonsymmetric linear systems, SIAM /. Sci. Comput., 19,3,968-994,1998. [23] M. BENZI Y M. TUMA, Sparse matrix orderings for factorized inverse preconditioners, in: Proceedings ofthe 1998 Cooper Mountain Conference on Iterative Methods, marzo 30-abril 3,1998. [24] M. BENZI Y M. TÚMA, A comparativa study of sparse approximate inverse preconditioners, Appl. Num. Math., 30,305-340,1999. Bibliografía 123 [25] M. BENZI Y M. TÚMA, Orderings for factorized sparse approximate inverse preconditioners, SIAM f. Sci. Comput., to appear. [26] L. BERGAMASCHI, G. PINI Y F.SARTORETTO, Approximate inverse preconditíoning in the parallel solutíon of sparse eigenproblems, Num. Lin. Alg. Appl, 7,99-116,2000. [27] E. BODEWIG, Matrix calculus, 2nd revised and enlarged edition, Interscience, New York, 1959. [28] R. BRIDSON Y W.P. TANG, Ordering, anisitropy and factored sparse approximate inverses, Preprint, Department of Computer Science, UniversityofWaterloo., 1998. [29] B. CARPENTIENRI, I.S. DUFF Y L. GIRAUD, Sparse pattem selecion strategies for robust Frobenius-norm minimization preconditioners in electromagnetism, Num. Lin. Alg. Appl, 7, 667-685,2000. [30] L. CESARI, Sulla risoluzione dei sistemi di equazione lineari per approssimazioni succesive, Atti. Accad. Naz. Lincei, Rend. Cl. Sci. Fis. Mat. Nat., 25, 422-428,1937. [31] T.F. CHAN, E. GALLOPOULOS, V. SIMONCINI, T. SZETO Y C.H. TONG, A Quasi-Minimal Residual Variant of the Bi-CGSTAB Algorithm for Nonsymmetric Systems, SIAM]. Sci. Statist. Comput, 15,338-247,1994. [32] E. CHOW, Robust preconditioning for sparse linear systems, Department of Computer Science, University of Minnesota, Minneapolis, MN, 1997. [33] E. CHOW, A priori sparsity pattems for parallel sparse approximate inverse preconditioners, SIAM }. Sci. Comput, 21,5,1804-1822,2000. [34] E. CHOW Y Y. SAAD, Approximate inverse techniques for blockpartitioned matrices, SLAMJ. Sci. Comput, 18,6,1657-1675,1997. [35] E.,CHOW Y Y. SAAD, Experimental study of ILU preconditioners for indefínite matrices, /. Comp. Appl. Math, 86,387-414,1997. [36] E. CHOW Y Y. SAAD, Approximate inverse preconditioners via sparsesparse iterations, SLAMJ. Sci. Comput, 19,3,995-1023,1998. [37] P. CONCUS, G.H. GOLUB Y G. MEURANT, Block predondicionint for the conjúgate gradient method,SÍAM/. Sci. Stat. Comp., 6,220-252,1985. [38] D.J.F. COSGROVE, J.C. DÍAZ Y A. GRIEWANK, Approximate inverse preconditionings for sparse linear systems, Int. f. Comput. Math., 44, 91-110, 1992. 114 Bibliografía [39] E.H. CUTHILL Y J.M. MCKEE, Reducing the Bandwidth of Sparse Symmetric Matrices, Proc. 24th National Conference ofthe Association for Computing Machinen/, Brondon Press editores, New Jersey, 157-172,1969. [40] H. DAG, Iterative methods and parallel computation for power systems. Department of Electrical Engineering, University of Wisconsrn, Madison, WI, 1996. [41] J.C. DÍAZ Y C.G. MACEDO JR., FuUy vectorizable block preconditionings with approximate inverses for non-symmetric systems of equations, Int. fourn. Num. Meth. Eng., 27,501-522,1989. [42] P.R DUBOIS, A. GREENBAUM Y G.H. RODRIGUE, Approximating the inverse of a matrix for use in iterative algorithms on vector processors, Computing, 22, 257-268,1979. [43] I.S. DUFF, R.G. GRIMESY Y J.G. LEWIS, Sparse matrix test problems, ACM Trans. Math. Software, 15,1-14,1989. [44] A.C.N. VAN DUIN Y H. WIJSHOFF, Scalable parallel preconditioning with the sparse approximate inverse triangular systems, Preprint, Computer Science Departmente, Uiúversity of Leiden, Leiden, the Netherland, 1996. [45] L.C. DUTTO, The effect of ordering on preconditioned GMRES algorithm, for solving the compressible Navier-Stokes equations, Int. J. Num. Meth. Eng., 36,457-497,1993. [46] H.C. ELMAN, A stability analysis of incomplete LU factorization, Math. Comp., 47,191-217,1986. [47] A.M. ERISMAN Y W.F. TINNEY, On computing certain elements of the inverse of a sparse matrix, Comm. ACM, 18,177-179,1975. [48] M.R. FlELD, An efficient parallel preconditioner for the conjúgate gradient algorithm, Hitachi Dublin Laboratory Technical Report HDL-TR-97-175, Cublin, Ireland, 1997. [49] E. FLÓREZ, M.D. GARCÍA, L. GONZÁLEZ Y G. MONTERO, The effect of orderings on sparse approximate inverse preconditioners for non-symmetric problems, Adv. Engng. Software, 33, 7-10,611-619,2002. [50] P.O. FREDERICKSON, Fast approximate inversión of large sparse linear systems, Math. Report, 7, Lakehead University, Thvmder Bay, Ganada, 1975. [51] M. GALÁN, G. MONTERO Y G. WINTER, Variable GMRES: an Optimizing Self-Configuring Implementation of GMRES(k) with Dynamic Memory AUocation, Tech. Rep. ofCEANI, Las Palmas, 1994. Bibliografía 215 [52] M. GALÁN, G. MONTERO Y G. WINTER, A Direct Solver for the Least Square Problem Arising From GMRES(k), Com. Num. Meth. Eng., 10, 743749,1994. [53] A. GEORGE, Computer Implementation of the Finite Element Method, Report Stan CS-7,1-208,1971. [54] A. GEORGE Y J.W. LIU, The Evolution of the Minimum Degree Ordering Algorithms, SIAM Rev., 31,1-19,1989. [55] I.C. GOHBERGH Y M.G. KREIN, Introduction to the Theory of Linear Nonselfadjoint Operators in Hilbert Space, Ed. Translations of Mathematical Monographs, 18, American Mathematical Society, Providence, Rhode Island, U.S.A., 1991. [56] L. GONZÁLEZ, G. MONTERO Y E. FLÓREZ, Aproxímate Inverse Freconditíoners Using Frobenius Inner Product I: Theorical Results, S*'* ILAS Conference, Barcelona, 1999. [57] N.I.M. GOULD Y J.A. SCOTT, Sparse approximate-inverse preconditioners using norm-minimization techniques, SIAM, }. Sci. Comput., 19, 2, 605-625,1998. [58] G.A. GRAVVANIS, Approximate inverse banded matrix techniques, Eng. Comp., 16, 3, 337-346,1999. [59] G.A. GRAVVANIS, The convergence rate and complexity of explicit preconditioned conjúgate gradient methods based on approximate inverse banded matrix techniques, Neural, Parallel & Scientific Computations, 9,355368,2001. [60] M. GROTE Y H.D. SIMÓN, Parallel Preconditioning and Approximate Inverses on the Connection Machine, Sixth SIAM Conference on Parallel Processing for Scientific Computing, 2,519-523, Philadelphia, 1992. [61] M. GROTE Y T. HUCKLE, Parallel preconditioning with sparse approximate inverses, SIAMf. Sci. Comput. 18,3,838-853,1997. [62] I. GUSTAFSSON Y D. LlNDSKOG, Completely paralleHzable preconditioning methods, Num. Un. Alg. Appl, 2,447-465,1995. [63] M.R. HESTENES Y E. STIEFEL, Methods of Conjúgate Gradients for Solving Linear Systems, ]our. Res. Nat. Bur. Sta. 49,6,409-436,1952. [64] R.A. HORN Y C.R JOHNSON, Topics in Matrix Analysis, Cambridge University Press, Cambridge, U.K., 1991. 116 Bibliografía [65] R.A. HORN Y C.R. JOHNSON, Tapies in Matrix Analysis. Ed. Cambridge University Press, Cambridge, 1999. [66] T.K. HUCKLE, Approximate sparsity pattems for the inverse of a matrix and preconditioning, in: Proceedings ofthe 15th IMACS World Congress 1997 on Scientific Computation, Modelling and Applied Mathematics, R.,Weiss and W. Schonauer editores, 2,569-574,1997. [67] T.K. HuCKLE, Efficient computation of sparse approximate inverses, Num. Lin. Alg. Appl, 5,57-71,1998. [68] O.G. JOHNSON, C.A. MICCHELLI Y G. PAUL, Polynomialpreconditioning for conjúgate gradient calculations, SIAM J. Num. Anal, 20,362-376,1986. [69] LE. KAPORIN, New convergence results and preconditioning strategies for the conjúgate gradient method, Num. Lin. Alg. Appl, 1,179-210,1994. [70] LE. KAPORIN, High quality preconditioning of a general summetric positive defínite matrix based on its U'^U + U^R + BFUdecomposition, Num. Lin. Alg. Appl, 5,483-509,1998. [71] C.T. KELLEY, , Iterative methodsfor Linear and Nonlinear Equations, Frontiers in Applied Mathematics, SIAM, Philadelphia, 1995. [72] S.A. KHARCHENKO, L.YU. KOLOTILINA, A.A. NIKISHIN Y A. YU. YEREMIN, A robust AINV-type method for constructing sparse approximate inverse preconditioners in factored form, Num. Lin. Alg. Appl, 8,165-179, 2001. [73] L.Y. KOLOTILINA, On approximate inverses of block H-matrices, Num. Anal Math. Mod., Moscú, 1989. [74] L. Yu. KOLOTILINA, A.A. NIKISHIN V A.Yu. YEREMIN, Factorized sparse approximate inverse preconditionings IV: simple approaches to rising efficiency, Num. Lin. Alg. Appl, 6,515-531,1999. [75] L. Yu. KOLOTILINA Y A. Yu. YEREMIN, On a family of two-level preconditionings of the incomplete block factorization type, Sov. J. Numer. Anal Math. Modelling, 1,292-320,1986. [76] L. Yu. KOLOTILINA Y A. Yu. YEREMIN, Factorized sparse approximate inverse preconditioning I. Theory, SIAM }. Matrix Anal Appl, 14, 45-58, 1993. [77] L. Yu. KOLOTILINA Y A. Yu. YEREMIN, Factorized sparse approximate inverse preconditioning II. Solution of 3D FE systems on massively paraUel computers, Internat. }. High Speed Comp., 7,191-215,1995.