Post on 06-Mar-2016
Precondicionado con inversas aproximadasMSc. Miguel Vargas Flixmiguelvargas@cimat.mx
http://www.cimat.mx/~miguelvargas
11/10/11 1/27
Contenido
Contenido
Nmero de condicin
Precondicionadores
Algunos precondicionadores para sistemas ralos
Gradiente conjugado precondicionado
Precondicionadores con inversa aproximada rala
Inversa aproximada rala con mtodo del residuo minimo
Gradiente biconjugado precondicionado con SAI
Algunos resultados numricos
Inversa aproximada rala factorizada
Algunos resultados numricos
Otro precondicionador
Preguntas?
Bibliografa
11/10/11 2/27
Nmero de condicin
Nmero de condicin
El nmero de condicin de una matriz Ann no singular, para una norma est dado por A =AA1.
Para la norma 2,
2 A=A2A12=max Amin A
,
donde son los valores singulares de la matriz.
Para una matriz A simtrica positiva definida,
A =max A min A
,
donde son los eigenvalores de A.
Un sistema de ecuaciones A x=b es considerado bien condicionado si un pequeo cambio en los valores de A o un pequeo cambio en b resulta en un pequeo cambio en x.
Un sistema de ecuaciones A x=b es considerado mal condicionado si un pequeo cambio en los valores de A o un pequeo cambio en b resulta en un cambio grande en x.
As, matrices con un nmero de condicin cercano a 1 se dicen que estn bien condicionadas.
11/10/11 3/27
Precondicionadores
Precondicionadores
Un precondicionador de una matriz A es otra matriz M tal que M A tiene un nmero de condicin menor
M A A.
En los mtodos iterativos estacionarios (ej. Gauss-Seidel) y los ms generales mtodos del subespacio de Krylov (ej. gradiente conjugado), el precondicionador, al reducir el nmero de condicin, buscar reducir el nmero de pasos necesarios para converger a la solucin de un sistema lineal de ecuaciones.
En vez de resolver el problemaA xb=0,
se resuelve el problemaM A xb =0.
En otras palabras, un precondicionador ideal sera una matriz M tal queM A1,
es decir MA1.
11/10/11 4/27
Algunos precondicionadores para sistemas ralos
Algunos precondicionadores para sistemas ralos
JacobiM1=diag A.
Para este precondicionador se toma la solo diagonal principal de A.
Factorizacin Cholesky incompleta. A=L LT sera la factoizacin completa, tomemosM1=G lG l
T,con
G lL.Para formar G l
G l i j0 si A i j0,con
G i j=1
G j j Ai jl=1j1
G i l G j l, para i jG j j=A j jl=1
j1
G j l2 .
11/10/11 5/27
Gradiente conjugado precondicionado
Gradiente conjugado precondicionado
x0, coordenada inicialr0 bA x0, gradiente inicialp0 r0, descenso inicialk 0mientras rk
k r k
T r kp k
T A pkxk1 x kk pkr k1 rkk A pk
k r k1
T rk1r k
T r kpk1 r k1k pkk k1
x0, coordenada inicialr0 bA x0, gradiente inicialp0 M r0, descenso inicialk 0mientras rk
k r k
T M r kp k
T A pkxk1 xkk pkr k1 rkk A pk
k r k1
T M r k1r k
T M r kpk1 M r k1 k pkk k1
Con los precondicionadores de factorizacin incompleta no definimos M, sino M1, entonces hacemos qk M rk, lo que quiere decir que en cada paso hay que resolver un sistema de ecuaciones M1qk r k.
11/10/11 6/27
Precondicionadores con inversa aproximada rala
Precondicionadores con inversa aproximada rala
Vamos a describir algorimos que permitan generar una matriz M se sea una aproximacin a la inversa de A, que adems tenga la propiedad de ser rala.
El problema es que a partir de una matriz rala su inversa no es necesariamente rala.
Ejemplo de las estructuras de una matriz rala y de su inversa
11/10/11 7/27
Precondicionadores con inversa aproximada rala
Una forma de buscar la inversa aproximada se basa en minimizar la norma de Frobenius del residuo IAM
F M =IAMF2 (1)
La norma de Frobenius se define como
AF=i=1m
j=1
n
ai j2= tr AT A.
Podemos descomponer (1) en sumas desacolpladas de las 2-normas de las columnas individuales del residuo IAM, es decir
F M =IAMF2=
j=1
n
e jAm j22,
donde e j es la j-sima columna de I y m j es la j-sima columna de M .
Podemos as paralelizar la construccin del precondicionador.
El problema ahora es: como elegir las entradas de M que son significativas para construir la inversa aproximada?
11/10/11 8/27
Inversa aproximada rala con mtodo del residuo minimo
Inversa aproximada rala con mtodo del residuo minimo
E. Chow, Y. Saad. Approximate Inverse Preconditioners via Sparse-Sparse Iterations. SIAM Journal on Scientific Computing. Vol. 19-3, pp. 995-1023. 1998
El algoritmo para construir M es el siguiente:
Sea M=M 0para j 1n
Sea m j M e jpara i 1sr j e jAm j
j r j
T A r jAr j
T A r jm j m j j r jTruncar las entradas menos significativas de m j
fin parafin para
11/10/11 9/27
Inversa aproximada rala con mtodo del residuo minimo
Es un mtodo iterativo donde s es usualmente pequeo. Si no se utilizan vectores ralos, este algorimo es de orden O n2 .Para inicializar el algoritmo puede se utiliza:
M 0=G,con
=tr AG
tr AG AG T ,G puede ser elegido como G=I o como G=AT.
M 0= I es la eleccin ms econmica de calcular. Sin embargo M 0= AT produce una mejor
aproximacin inicial en algunos casos.
Vamos ahora a ver como utilizar vectores ralos para reducir la cantidad de operaciones.
Dado un vector ralo b, la notacin b indica el nmero de entradas no cero de b.
11/10/11 10/27
Inversa aproximada rala con mtodo del residuo minimo
Una tcnica para elegir las entradas m j ms significativas es truncar entradas en la direccin de descenso.
Sea M=M 0para j 1n
Sea m j M e j un vector raloSea d j un vector ralo con la misma estructura que m jpara i 1sr j e jAm jd j r j se toman solo las entradas de r j con estructura en d j
j r j
T Ad jAd j
T Ad jm j m j j d jsi m j lfil
agregar maxr jk tal que k no exista en la estructura de m jfin para
fin paraApproximate inverse via MR iteration with dropping in the search direction
Problemas:
11/10/11 11/27
Inversa aproximada rala con mtodo del residuo minimo
La estructura final de M no ser necesariamente simtrica.
No se garantiza que M no sea singular si s no es grande.
No funciona como precondicionador para el gradiente conjugado precondicionado.
Funciona para el gradiente biconjugado precondicionado.
11/10/11 12/27
Gradiente biconjugado precondicionado con SAI
Gradiente biconjugado precondicionado con SAI
El mtodo de gradiente biconjugado se basa en el mtodo de gradiente conjugado y al igual que ste sirve para resolver systemas de ecuaciones lineales
A x=b,con la diferencia que Ann no tiene que ser simtrica.
Este mtodo requiere calcular un pseudo-gradiente g k y una pseudo-direccin de descenso pk.
El mtodo se construye de tal forma que que los pseudo-gradientes g k sean ortogonales a los gradientes g k y que las pseudo-direccines de descenso pk sean A-ortogonales al las direcciones de descenso pk [Meie94 pp6-7].
Si la matriz A es simtrica, entonces el mtodo es equivalente al gradiente conjugado.
No garantiza convergencia en n pasos como lo hace el gradiente conjugado.
11/10/11 13/27
Gradiente biconjugado precondicionado con SAI
Hay que hacer cuatro multiplicaciones matrix-vector: pk A, A pk, r k M, M r k. Para hacerlo eficiente, hay que almacenar A y M dos veces cada una usando CSR y CSC.
x0, coordenada inicialr0 bA x0, gradiente inicialr0 r 0
T, pseudo-gradiente inicialp0 M r0, direccin de descenso inicialp0 r0 M, pseudo-direccin de descenso, toleranciak 0mientras rk
k r k M r kpk A pk
xk1 xkk pkr k1 rk A pkr k1 rk pk A
k r k1 M r k1rk M r k
pk1 M r k1 k1 pkpk1 r k1 M k1 pkk k1
11/10/11 14/27
Algunos resultados numricos
Algunos resultados numricos
Ecuacin de calor 2D, elementos triangulares, 62,216 elementos, 31,615 variables.Solver 1 Thread 2 Thread 4 Thread 6 Thread 8 Thread Steps MemoryBG 1.40 1.02 0.76 0.68 0.69 721 14,879,486BG + Jacobi 1.68 1.09 0.81 0.73 0.72 682 15,132,430BG + ILU(0) 1.31 1.25 1.21 1.20 1.22 190 23,313,662BG + SAI(1) 26.93 14.05 7.40 5.36 4.55 465 21,323,518
1 Thread 2 Thread 4 Thread 6 Thread 8 Thread0
5
10
15
20
25
BG
BG + Jacobi
BG + ILU(0)
BG + SAI(1)
11/10/11 15/27
Inversa aproximada rala factorizada
Inversa aproximada rala factorizada
Los precondicionadores para el gradiente conjugado requiere que al igual que la matriz A, el precondicionador M sea simtrico positivo definido.
El siguiente mtodo genera un precondicionador que cumple que M sea positivo definido y adems es posible calcularlo en paralelo.
Sea A que una matriz binaria construida a partir de A
Ai j={1 si i= j o D1/2 AD1/2 i jumbral0 otro caso ,el umbral es un valor no negativo y la matriz diagonal D es
D i i={Ai i si Ai i>01 otro caso .La matriz diagonal D funciona como un escalamiento.
Notese que la estructura de A es ms rala que la estructura deA.
11/10/11 16/27
Inversa aproximada rala factorizada
Buscamos construir un precondicionadorM=GTG ,
donde G es una matriz triangular inferior tal queGL1,
donde L es el factor Cholesky de A.
Sea S L la estructura de G , sta se construir de tal forma que tenga la misma estructura que la parte triangular inferior de A.
Buscamos entonces minimizar IG LF2, lo que puede lograrse sin conocer el factor L,
resolviendo las ecuaciones(G L LT)i j=(LT)i j, (i , j )S L,
donde S L contiene la estructura de G , la ecuacin anterior es equivalente a(G A)i j=( I )i j, i , j S L,
con G=D1G , y D es la diagonal de L, dado que {LT}i j para (i , j )S L es una matriz diagonal.Las entradas de G se calculan por renglones (tiene la misma estructura que G ). stas se pueden calcular en paralelo.
11/10/11 17/27
Inversa aproximada rala factorizada
El algoritmo para calcular las entradas de G es el siguiente:
Sea S L la estructura de Gpara i 1n
para (i , j )S Lresolver (AG i) j=i j
fin_parafin_para
Resolver (AG i) j=i j significa que, si m=(G i) el nmero de entradas no cero del rengln i de G , entonces hay que resolver un sistema SPD local pequeo de tamao mm.
(AG i) j=i j, ( i , j)S L,
Ejemplo. Si la estructura del rengln 9 de G contiene 4, 7, 8, 9, entonces, para calcular los valores de G 9 hay que resolver el sistema:
(a4 4 a4 7 a4 8 a4 9a7 4 a7 7 a7 8 a79a8 4 a8 7 a8 8 a89a9 4 a9 7 a9 8 a99)(g 9 4g 9 7g 98g 99)=(0001).
11/10/11 18/27
Inversa aproximada rala factorizada
Finalmente, para obtener GG=G D.
Como es un sistema simtrico positivo definido, podemos resolver estos sistemas utilizando factorizacin Cholesky.
Es posible eficientar la solucin de estos sistemas haciendo casos particulares para m=1, m=2, hasta m=3, y el caso general para m4.
Este precondicionador M=GTG se puede aplicar entonces al gradiente conjugado precondicionado:
Este precondicionador es simtrico positivo definido si no hay ceros en la diagonal de G l.
El algoritmo para construir el precondicionador es paralelizable.
El algoritmo es estable si A es simtrica positiva definida.
E. Chow. Parallel implementation and practical use of sparse approximate inverse preconditioners with a priori sparsity patterns. International Journal of High Performance Computing, Vol 15. pp 56-74, 2001.
11/10/11 19/27
Inversa aproximada rala factorizada
El algoritmo es el siguiente
x0, coordenada inicialr0 bA x0, gradiente inicialp0 M r0, descenso inicialk 0mientras rk
k r k
T M r kp k
T A pkxk1 xkk pkr k1 rkk A pk
k r k1
T M r k1r k
T M r kpk1 M r k1 k pkk k1
M r k se calcula con dos multiplicaciones matriz-triangular*vector,M r k=G l
T(G l r k ).No hay que resolver un sistema de ecuaciones en cada paso.
11/10/11 20/27
Inversa aproximada rala factorizada
Es posible aumentar la cantidad de entradas del precondicionador y mejorar el nmero del condicin del problema. Con el costo de aumentar la cantidad de tiempo utilizada en construir el precondicionador.
El mtodo se conoce como potencias de matrices ralas.
Se puede construir S L es a travs de las estructuras formadas porA, A2, ..., Al ,
Las potencias de A pueden ser calculadas combinando los renglones de A. Denotemos el k-simo rengln de Al por Ak ,:
l , asAk ,:
l = Ak ,:l1 A.
La estructura S L ser entonces la parte triangular inferior de Al .
11/10/11 21/27
Algunos resultados numricos
Algunos resultados numricos
Ecuacin de calor en 2D, 1813,582 elementos, con 1,005,362 variables, A =6 ' 355,782.Solver 1 thread [s] 2 threads [s] 4 threads [s] 6 threads [s] 8 threads [s] StepsCholesky 108.4 85.1 74.4 71.0 69.276CG 120.8 82.2 81.8 84.3 85.769 3,718CG-Jabobi 132.0 92.9 84.6 87.9 89.251 3,504CG-ICholesky 81.1 76.2 74.2 75.9 77.112 978CG-FSAI 116.9 74.2 66.4 57.3 55.707 1,288
Cholesky CG CG-Jabobi CG-ICholesky CG-FSAI0.000
20.000
40.000
60.000
80.000
100.000
120.000
140.000
1 thread 2 threads 4 threads 6 threads 8 threads
Tim
e [s
]
11/10/11 22/27
Algunos resultados numricos
Deformacin de slidos en 2D, 501,264 elementos, con 909,540 variables, A =18 ' 062,500.Solver 1 thread [s] 2 threads [s] 4 threads [s] 8 threads [s] Steps MemoryCholesky 227.2 131.4 81.9 65.4 3,051,144,550CG 457.0 305.8 257.5 260.4 9,251 317,929,450CG-Jacobi 369.0 244.7 212.4 213.7 6,895 325,972,366CG-IChol 153.9 121.9 112.8 117.6 1,384 586,380,322CG-FSAI 319.8 186.5 155.7 152.1 3,953 430,291,930
Cholesky CG CG-Jacobi CG-IChol CG-FSAI0
50
100
150
200
250
300
350
400
450 1 thread2 threads4 threads8 threads
Tim
e [s
]
11/10/11 23/27
Algunos resultados numricos
Deformacin de slidos en 3D, 264,250 elementos, con 978,684 variables, A =69 ' 255,522.
Solver 1 thread [s] 2 threads [s] 4 threads [s] 6 threads [s] 8 threads [s] MemoryCholesky 8,562 4,418 2,635 2,055 1,911 19,864,132,056CG-FSAI 4,410 2,703 1,631 1,484 1,410 1,440,239,572CG-Jacobi 9,568 5,591 3,444 3,207 3,278 923,360,936
Cholesky CG-FSAI CG-Jacobi0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
1 2 4 6 8
Tim
e [s
]
Cholesky CG-FSAI CG-Jacobi0
2,000,000,000
4,000,000,000
6,000,000,000
8,000,000,000
10,000,000,000
12,000,000,000
14,000,000,000
16,000,000,000
18,000,000,000
20,000,000,000
Mem
ory [
byte
s]
11/10/11 24/27
Otro precondicionador
Otro precondicionador
M. J. Grote, T. Huckle. Parallel preconditioning with sparse approximate inverses. SIAM Journal of Scientific Computing, Vol. 18-3, pp. 838853, 1997.
Tambin parte de minimizar por submatrices de Aminm je j A m j2, (2)
factorizando la matriz AA=QR0 ,
donde R es una matriz triangular superior.
Sea c=QT e j, la solucin de (2) esm j=R
1 c.
11/10/11 25/27
Preguntas?
Preguntas?
11/10/11 26/27
Bibliografa
Bibliografa
[Saad03] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003.
[Benz02] M. Benzi. Preconditioning Techniques for Large Linear Systems: A Survey. Journal of Computational Physics 182, pp418477. Elsevier Science, 2002.
[Chow98] E. Chow, Y. Saad. Approximate Inverse Preconditioners via Sparse-Sparse Iterations. SIAM Journal on Scientific Computing. Vol. 19-3, pp. 995-1023. 1998
[Chow01] E. Chow. Parallel implementation and practical use of sparse approximate inverse preconditioners with a priori sparsity patterns. International Journal of High Performance Computing, Vol 15. pp 56-74, 2001.
11/10/11 27/27