Lab.3: Solución de sistema de ec. lineales con Montecarlo / Metropolis

Algoritmos Paralelos

Glen Rodríguez

Aplicaciones

Redes eléctricas Redes de tránsito Cadenas de Markov Solución de algunos tipos de

ecuaciones diferenciales parciales

Problema

Obtener vector x (tamaño nx1) que cumple: Ax=b.

Donde b es un vector conocido del mismo tamaño de x, y A es una matriz de tamaño nxn que cumple: A es NO SINGULAR Se puede escoger una matriz M no

singular tal que MA=I-L, y que |λ(L)| < 1 para todos los eigenvalues λ de L (no singular también)

Algoritmo

λ es un número real llamado eigenvalue si hay solución de Lu= λu

El problema se puede ver como: x=Lx + Mb = Lx + f

Si M=I, se puede demostrar que la sumatoria ΣLi converge a A-1 si i∞

En los ‘50s se demostró que un “random walk” puede aproximar esta sumatoria.

Algoritmo x0=0; xm=b+Lxm-1 donde xm=m-esima

aproximación a x

l ij

Pseudocódigo1. Calcular pij, vij. Se puede asumir que pi=constante

(ej.:0.01). Entonces se puede hacer pij iguales en una fila i =(1-pi)/numero de pij ≠0 en fila i=ti; vij=lij/ti. O se puede hacer pij= (1-pi)* | lij | /Σ |lij| en fila i; vij= Σ |pij| en fila i / (1-pi)

2. For elem=1 to n3. x_sum=0;4. for k=1 to m5. v_prd=1;; ult=entero entre 1 y n tal que lult,elem ≠0;

x_est=bult/pult

6. Escoger r= real aleatorio unif. entre 0 y 17. While r<=1-pult

8. Asociar r a un valor “escog” usando las pult,j 9. v_prd= v_prd* vult,escog

Pseudocódigo

9. x_est=v_prd*bescog/pescog

10. ult=escog11. end while12. x_sum= x_sum+x_est;13. end for14. x[elem]= x_sum/m;15. end for

Que se puede paralelizar: cualquiera de los for.Mejoras: una vez conocido x[1],…x[q], usar esos

valores para x[q+1],…,x[n]

Redes eléctricas

Sólo necesitamos 3 leyes: Ohm: ΔV=IR entre los extremos de un

conductor de resistencia R 1ra de Kirchhoff: Σ Iin = Σ Iout en cada

unión entre 2 o más conductores 2da de Kirchhoff: Σ ΔV = ΣFuentes en

un loop ó lazo cerrado.

Ejemplo

Redes de tráfico Se usa solo la 1ra Kirchhoff.

Resolviendo Ec. de Poisson con Montecarlo

Sea la ecuación de Poisson:

Con fuente sinusoidal : -2π2 sen(π x) sen(π y) Y cero en todas las fronteras En diferencias

finitas:

Solución con “random walks” De un punto (i,j) comenzar el “random

walk”.La probabilidad de moverse a cualquiera de los 4 puntos vecinos es ¼.

Generar un número aleatorio para escoger el vecino.

Añadir g(x,y) en la nueva posición Repetir hasta llegar a una frontera Esto fue solo UNA “random walk” Después de N “random walks” el estimado de

u(i,j) es:

Paralelización de esta solución Hacer el siguiente proceso para todos los

puntos Comenzar en una frontera

De afuera hacia adentro Fila por fila Actualizar (intercambiar) data en el límite entre

procesadores Actualizar el walk dentro de los puntos de un

procesador Si el walk sale de un procesador, informarle al otro

procesador

Metrópolis Se usa para integrales con funciones de peso. Genera un random walk pero guiado por una función de

peso w(x). Ejemplo en 2-D comenzando de un punto (xi,yi):

1. Escoger δ, el step size

2. Generar dos aleatorios R1, R2 uniformes en el rango [-1,1]

3. El nuevo punto será1. xT

i+1 = xi + δR1

2. yTi+1 = yi + δR2

4. Evaluar el ratio de w en el punto actual vs. el punto anterior1. r= w(xT

i+1 ,yTi+1) / w(xi,yi)

Metrópolis5. Si r>1 aceptar el nuevo punto y regresar a (2)

1. xi+1 = xTi+1 ; yi+1 = yT

i+1

6. Si r<1 aceptar el nuevo punto con probabilidad r. O sea, generar un aleatorio uniforme η en [0, 1] y aceptar el punto solo si r>η (hacer lo mismo que en el paso 5)

7. Caso contrario, rechazar el nuevo punto y volver a (2)

El step size no debe ser ni muy chico ni muy grande.

Los puntos obtenidos se usan luego para integrar, como en MonteCarlo (Constante*n/N).

La paralelización es similar a la discutida para las integrales con Montecarlo