Transcript
Precondicionamiento con factorizaciones incompletas Miguel Vargas-Félix
[email protected] http://www.cimat.mx/~miguelvargas
CIMAT, September 29, 2015
1/24
Número de condición El número de condición κ de una matriz A no singular, para una norma ∥⋅∥ está dado por κ(A)=∥A∥⋅∥A−1∥.
Para la norma ∥⋅∥2,
κ(A)=∥A∥⋅∥A−1∥=
σ max ( A) , σ min (A)
donde σ son los valores singulares de la matriz. Para una matriz A simétrica positiva definida, λ max (A) , κ(A)= λ min (A ) donde λ son los eigenvalores de A.
http://www.cimat.mx/~miguelvargas
2/24
El número de condición indica que tan sensible es una función a pequeños cambios en la entrada. Así, matrices con un número de condición cercano a 1 se dicen que están bien condicionadas. Un sistema de ecuaciones A x=b es considerado bien condicionado si un pequeño cambio en los valores de A o un pequeño cambio en b resulta en un pequeño cambio en x. Un sistema de ecuaciones A x=b es considerado mal condicionado si un pequeño cambio en los valores de A o un pequeño cambio en b resulta en un cambio grande en x. Al reducir el número de condición de un sistema, se puede incrementar la velocidad de convergencia del gradiente conjugado.
http://www.cimat.mx/~miguelvargas
3/24
Gradiente conjugado precondicionado Entonces, en vez de resolver el problema A x−b=0, se resuelve el problema −1
M ( A x−b )=0, −1
con M
una matriz cuadrada, la cual recibe el nombre de precondicionador. −1
−1
−1
El mejor precondicionador sería claro M =A , así x=M b, y el gradiente conjugado convergiría en un paso. −1
Al igual que la matriz A, el precondicionador M
tiene que ser simétrico positivo definido. −1
Hay dos tipos de precondicionadores, implícitos M y explícitos M . −1
En algunos casos es costoso calcular M , en general se utilizan precondicionadores con inversas fáciles de calcular o precondicionadores implicitos que se puedan factorizar. http://www.cimat.mx/~miguelvargas
4/24
El algoritmo es el siguiente: entrada: A , x 0, b,ε
r 0 ← A x0 −b −1 q0 ← M r0 p0 ← −q 0 k ← 0 mientras ∥r k∥>ε
w ← A pk T rk qk αk ← T pk w x k +1 ← x k +α k p k r k +1 ← r k +α w −1 q k +1 ← M r k +1, o resolver M q k + 1 ← r k +1 T r k +1 q k +1 βk ← r Tk q k p k +1 ← −q k +1 +βk +1 p k k ← k +1 Nótese que ahora el algoritmo require aplicar el precondicionador en cada paso. http://www.cimat.mx/~miguelvargas
5/24
Precondicionador Cholesky incompleto Éste consiste en tener un precondicionador T
M=H t D H t , dónde H t es una matriz triangular inferior rala que tiene una estructura similar a la factorización Cholesky de A. El entero t indicará el nivel de truncamiento. La estructura de H 0 será igual a la estructura de la parte triangular inferior de A. La estructura de H n será igual a la estructura de L (factorización Cholesky de A). Los valores de H t se llenarán utilizando la fórmulas de la factorización Cholesky, así H i i =1,
(
)
j−1
1 H i j= Ai j−∑ H i k H j k D k , para i> j Dj k =1 j−1
D j = A j j − ∑ H 2j k D k . k =1
Veamos un ejemplo, para una matriz A ∈ℝ http://www.cimat.mx/~miguelvargas
424×424
, la estructura de H t sería... 6/24
A=
η ( A )=2456
http://www.cimat.mx/~miguelvargas
7/24
H 0=
η ( H0 )=1436
http://www.cimat.mx/~miguelvargas
8/24
H5 =
η ( H5 )=1767
http://www.cimat.mx/~miguelvargas
9/24
H10 =
η ( H10 )=2196
http://www.cimat.mx/~miguelvargas
10/24
H 20 =
η ( H 20 )=5262
http://www.cimat.mx/~miguelvargas
11/24
H 30=
η ( H30 )=6779
http://www.cimat.mx/~miguelvargas
12/24
H 40 =
η ( H 40 )=7012
http://www.cimat.mx/~miguelvargas
13/24
L=
η ( L )=7012
http://www.cimat.mx/~miguelvargas
14/24
Factorización Cholesky simbólica truncada Para una matriz rala A, definamos a i ≝ { k >i
∣
A i k ≠0 }, i=1…n,
como el conjunto de los índices de los elementos no nulos del renglón i de la parte estrictamente triangular superior de A. T
De forma análoga definimos para la matriz H t , los conjuntos hi ≝ { k >i
∣
H i k ≠0 }, i=1…n.
Para la matriz ejemplo siguiente, el conjuto a2 contendrá { 3, 4 }; el conjunto h2 contendrá { 3, 4, 6 }.
(
A1 1 A2 1
A=
A1 2 A2 2 A3 2 A4 2
A1 6 A2 3 A3 3
A2 4 A3 5 A4 4
A5 3 A6 1 a 2= {3, 4 }
http://www.cimat.mx/~miguelvargas
A5 5 A5 6 A6 5 A6 6
) ( T
Ht =
H 11
H 12 H 22 H 23 H 24 H 33 H34 H 44
H 35 H 45 H 55
H 16 H 26 H 36 H 46 H 56 H 66
)
h2={ 3, 4, 6 }
15/24
T
Requeriremos de conjuntos r i que serán usados para registrar las columnas de H t cuyas estructuras T afectarán al renglón i de H t .
for i ← 1, 2,…, n · ri ← ∅ for i ← 1, 2,…, n · hi ← a i · for each j∈r i · · hi ← hi ∪{x∈h j∣x−i
ε
w ← A pk T rk qk αk ← T pk w x k +1 ← x k +α k p k r k +1 ← r k +α w T resolver H t D Ht q k + 1=r k +1 T r k +1 q k +1 βk ← r Tk q k p k +1 ← −q k +1 +βk +1 p k k ← k +1 http://www.cimat.mx/~miguelvargas
18/24
Pero hay un problema... Puede suceder que el precondicionador no sea SPD. Para evitar esto utilizaremos el algoritmo de Munksgaard [Munk80], que consiste en dos estrategias: 1. Perturbar la diagonal de A con un factor α, j−1
D j j =α A j j −∑ H 2j k D k , k =1
esto asegurará que el precondicionador sea SPD, el valor de α se puede encontrar por prueba y error. 2. Perturbar los pivotes para mejorar la estabilidad si éstos son negativos o cecanos a cero, si D j j ≤u
{
a j k∣ , entonces D j j = ∣ ∑ ( ) k≠ j
∑ ∣a j k∣ k≠ j
1
∑ ∣a j k∣≠0 k≠ j . si ∑ ∣a j k∣=0 si
k≠ j
Un valor de u adecuado es 0.01. http://www.cimat.mx/~miguelvargas
19/24
Algunos resultados Problem Dimension Nodes Elements Element type Variables nnz(A) Tolerance
Arc 2 50,805 49,972 Quadrilateral 101,610 1’808,980 1x10-5
45.0
39.1
40.0 35.0 30.0
23.5
25.0
18.6
Time [s]
20.0
18.6
17.7
18.9 14.4
15.0 10.0
19.6
5.6
5.5
5.0 0.0
http://www.cimat.mx/~miguelvargas
20/24
Problem Dimension Nodes Elements Element type Variables nnz(A) Tolerance
Solver Cholesky CG CG + Jacobi CG + Cholesky [0] CG + Cholesky [1] CG + Cholesky [10] CG + Cholesky [100] CG + Cholesky [1000] CG + Cholesky [10000] CG + Cholesky [100000]
http://www.cimat.mx/~miguelvargas
Arc 2 50,805 49,972 Quadrilateral 101,610 1’808,980 1x10-5
Time [s] 5.6 39.1 23.5 18.6 19.6 18.6 17.7 18.9 14.4 5.5
Steps 7723 4070 1519 1519 1482 1162 708 308 1
nnz(H) 6911289 954741 954741 990125 1572311 3434219 5715237 6911289
Memory 215100454 40768170 40768170 71330422 71330422 72179638 86152102 130837894 185582326 214287574
Memory 100.0% 19.0% 19.0% 33.2% 33.2% 33.6% 40.1% 60.8% 86.3% 99.6%
21/24
¿Paralelizable? En cada paso tenemos que resolver el sistema T
H t D Ht q k + 1=r k +1. T
Hagamos d=H t q k +1 y c=D d, entonces tendremos tres sistemas de ecuaciones T
H t c=r k +1, D d=c, H t q k +1=d. Primero tenemos que resolver H t c=r k +1 para c haciendo una sustitución hacia adelante con
(
i−1
)
1 ci = r i −∑ H i k c k . Hii k =1 Después resolver D d=c, para obtener d. T
Finalmente, resolver H t q k +1=d para q k +1 sustituyendo hacia atrás con
(
n
)
1 T qi= T d i − ∑ H i k q k . H ii k =i+1 Es muy difícil paralelizar esta sustitución [Heat91 pp454-455] debido a que existe una dependencia secuencial entre las variables. http://www.cimat.mx/~miguelvargas
22/24
¿Preguntas?
[email protected]
http://www.cimat.mx/~miguelvargas
23/24
Referencias [Benz02] M. Benzi. Preconditioning Techniques for Large Linear Systems: A Survey. Journal of Computational Physics, Vol. 182, pp. 418-477. 2002. [Golu96] G. H. Golub, C. F. Van Loan. Matrix Computations. Third edidion. The Johns Hopkins University Press. 1996. [Heat91] M T. Heath, E. Ng, B. W. Peyton. Parallel Algorithms for Sparse Linear Systems. SIAM Review, Vol. 33, No. 3, pp. 420-460. 1991. [Meie94] U. Meier-Yang. Preconditioned conjugate gradient-like methods for nonsymmetric linear systems. University of Illinois. 1992 [Munk80] N. Munksgaard. Solving sparse symmetric sets of linear equations by preconditioned conjugate gradients. ACM Transactions on Mathematical Software, Vol 6-2, pp. 206219. 1980 [Noce06] J. Nocedal, S. J. Wright. Numerical Optimization, Springer. 2006. [Saad03] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003.
http://www.cimat.mx/~miguelvargas
24/24