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

1. Lab.3: Solución de sistema de ec. lineales con Montecarlo / Metropolis Algoritmos Paralelos Glen Rodríguez

2. Aplicaciones • Redes eléctricas • Redes de tránsito • Cadenas de Markov • Solución de algunos tipos de ecuaciones diferenciales parciales

3. 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)

4. 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.

5. Algoritmo • x0=0; xm=b+Lxm-1 donde xm=m-esima aproximación a x

6. l ij

8. Pseudocódigo • 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) • For elem=1 to n • x_sum=0; • for k=1 to m • v_prd=1;; ult=entero entre 1 y n tal que lult,elem≠0; x_est=bult/pult • Escoger r= real aleatorio unif. entre 0 y 1 • While r<=1-pult • Asociar r a un valor “escog” usando las pult,j • v_prd= v_prd* vult,escog

9. Pseudocódigo • x_est=v_prd*bescog/pescog • ult=escog • end while • x_sum= x_sum+x_est; • end for • x[elem]= x_sum/m; • 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]

10. 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.

11. Ejemplo

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

13. 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:

14. 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:

15. 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

16. 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): • Escoger δ, el step size • Generar dos aleatorios R1, R2 uniformes en el rango [-1,1] • El nuevo punto será • xTi+1 = xi + δR1 • yTi+1 = yi + δR2 • Evaluar el ratio de w en el punto actual vs. el punto anterior • r= w(xTi+1 ,yTi+1) / w(xi,yi)

17. Metrópolis • Si r>1 aceptar el nuevo punto y regresar a (2) • xi+1 = xTi+1 ; yi+1 = yTi+1 • 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) • 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