Metodos de Monte Carlo

e de Otimizacao

Tereza Mendes

Sao Carlos, 28 de Julho de 2003

Monte Carlo Estocastico

e.g. Sistema complexo evolu-

indo segundo regra estocas-

tica:

• automato celular

• agregacao limitada por

difusao

Otimizacao Determinıstico

e.g. Terreno com 24 arvores

produz 600 macas por arvo-

re; para cada arvore a mais a

producao por arvore diminui

12 macas

f(x) = (24+x)×(600−12x)

df/dx = 0 ⇒ x = 13

Monte Carlo Estocastico

Determinıstico

usado no calculo de integrais definidas (e.g.

area)

Estocastico

Otimizacao Determinıstico

metodos probabilısticos como recozimento

simulado, algoritmos geneticos, CA

Monte Carlo Estocastico

Otimizacao Determinıstico

Referencias

•

A Guide to Monte Carlo Simulations in Sta-

tistical Physics, Landau & Binder (Cambridge,

2000)

• Monte Carlo Methods in Statistical Physics, Newman &

Barkema (Oxford, 1999)

• Monte Carlo Methods in Statistical Mechanics: Foundations

and New Algorithms, Sokal (1996),

http://citeseer.nj.nec.com/sokal96monte.html

• Notas de aula - Transicao de fase e fenomenos crıticos

(2002), http://lattice.if.sc.usp.br/

Aplicacoes

• Mecanica Estatıstica: comportamento

macroscopico (termodinamica) a partir da

descricao microscopica de sistemas como

fluidos/gases, modelos de materiais magneticos,

sistemas biologicos; Fenomenos crıticos; Sistemas

complexos.

• Materia Condensada: sistemas quanticos,

polımeros, fluidos complexos, propriedades

condutoras/magneticas.

• Cromodinamica Quantica (QCD): teoria quantica

de campos que descreve a forca nuclear como

interacao forte entre quarks e gluons; Formulacao

de Rede ↔ Mecanica Estatıstica.

Mecanica Estatıstica

A probabilidade de uma configuracao S para um

sistema fısico em equilıbrio a temperatura T (ensemble

canonico) e dada em termos de sua hamiltoniana H(S)

pela distribuicao de Boltzmann

P (S) =e−βH(S)

Z; Z =

∫

dS e−βH(S); β = 1/KT

Medias termodinamicas de observaveis (e.g. energia:

E =< H(S) >) sao dadas por < A >=∫

dS A(S)P (S).

O Modelo de Ising

“spins” de dois estados preferem estar alinhados

S Si j

-J +J

S Si j

H(S) = −J∑

<i,j>

Si Sj − H∑

i

Si

• Energia: E =< H(S) >

• Calor especıfico: CV = ∂E/∂T

• Magnetizacao: M =<∑

i Si >

• Suscetibilidade: χ = ∂M/∂H

Simulacoes Numericas

• experimentos que nao podemos/queremos realizar

(avioes, guerra nuclear, evolucao)

• problemas sem solucao analıtica; sistemas

complexos: nao-linearidade, fenomenos crıticos.

Nota: independencia de detalhes (universalidade)

• comportamento dinamico (e.g. automato celular),

regras de evolucao locais ⇒ processamento

paralelo. Nota: dinamica pode ser introduzida

• tratamento teorico, com aspectos experimentais

(dados, erros, “medidas” no tempo)

Metodos de Monte Carlo

Geradores de Numeros Aleatorios

Anyone who considers arithmetical methods of

producing random digits is, of course, in a state of sin.

John Von Neumann (1951)

gerador = prescricao algebrica que produz

sequencia de numeros ri com distribuicao desejada

(em geral uniforme em [0,1]) dada uma semente.

Nota: esta sequencia e determinıstica: operacao

repetida a partir do mesmo ponto inicial gera a

mesma sequencia ⇒ numeros pseudo-aleatorios.

Exemplo: rand()

! inicializacao

iseed = 2003

call srand(iseed)

! numero aleatorio

r = rand()

• semente: numero inteiro

• a cada passo um novo inteiro e produzido e usado

como semente para o passo sucessivo

• inteiros sao convertidos em reais em [0,1]

• numero aleatorio em [a, b]: a + rand() * (b-a)

Caracterısticas de um bom gerador

• distribuicao dos ri e uniforme ⇒ testes

• perıodo muito maior do que o comprimento da

sequencia usada na simulacao

• possıvel armazenar a cada momento a semente

associada a um numero da sequencia

• sequencia produzida a partir de uma semente e

a mesma em computadores diferentes

• tempo para geracao dos ri e o menor possıvel

Metodo Congruencial Linear

in+1 = (a in + c) mod m

rn+1 =dble(in+1)

dble(m)

onde a, c e m sao numeros inteiros fixos. Para um

gerador com perıodo longo e necessario m muito

grande. Nota: m pode ser maior do que o maior

inteiro que pode ser armazenado no computador.

Nestes casos resolve-se o problema fatorizando m.

Area do Cırculo de Raio 1

lancando N pontos aleatorios uniformemente no

quadrado: x ∈ [−1, 1], y ∈ [−1, 1]

y

−1

−1

1

1

x.... .... .

..

.. ..

..

.

.

.

. .A◦

A2

=π

4=

n

N

onde n < N e o numero

de pontos contidos no

cırculo.

Integral Unidimensional

I =

∫ 1

0

f(x) dx =

∑

i f(xi)

N

com xi uniformemente distribuıdos em [0,1]. Na

verdade, para N finito

I ≡∑

i f(xi)

N

e uma variavel aleatoria, que converge para seu valor

medio I com erro proporcional a 1/√

N .

σ2I

=σ2

f

N=

<f2 > − <f >2

N

Metodos Determinısticos

f(x)

x x + h

Regra do trapezoide: estima a area

compreendida entre x e x + h como

a area do trapezoide definido pela

aproximacao linear da funcao entre

estes dois pontos.

∫ x+h

x

f(x′) dx′ =h

2[f(x) + f(x + h)] + O(h3f ′′)

erro para integral em [a, b] e O(h2) ∼ N−2

Regra de Simpson: aproximacao de 3 pontos para f(x)

⇒ erro e ∼ N−4

Comparacao

• d = 1: com 2N pontos o erro diminui por um fator

4 (trapezoide), 16 (Simpson) ou√

2 (Monte Carlo).

• d > 1: N ∼ 1/hd ⇒ erro N−2/d (trapezoide) ou

N−4/d (Simpson). Monte Carlo comeca a ser

vantagem a partir de dimensao d = 8.

• tipicamente, em mecanica estatıstica d ∼ 103.

• tempo para somar 21000 termos (Ising 3d com 10

pontos por direcao) em computador de 1 Tflops:

t = 10288s = 10270 × idade do universo

Amostragem por Importancia

I =

∫ 1

0f(x) w(x) dx =

∑

i f(xi) w(xi)

N

com xi uniformemente distribuıdos em [0,1].

Considere w(x) normalizada:∫ 10 w(x)dx. Se w(x)

for muito concentrada ha grande perda de

eficiencia. Portanto

I =

∫ 1

0f(x)w(x)dx =

∑

i f(xi)

N

onde os xi sao gerados com a distribuicao w(xi).

Amostragem de distribuicoes

Temos necessidadde de metodos para produzir xi

com distribuicao w(xi) (amostrar w). Por exemplo

a distribuicao exponencial

w(x) = e−x em [0,∞]

e amostrada gerando r uniforme em [0,1] e

tomando x = −log(1− r), ja que isso corresponde a

P (x) dx = P (r) dr ⇒ P (x) =d

dx(1−e−x) = w(x)

Distribuicao discreta

Por exemplo:

• x = 1 com probabilidade p

• x = −1 com probabilidade 1 − p

gero r uniforme em [0,1] e tomo

• x = 1 se r ≤ p

• x = −1 se r > p

Metodo da Rejeicao

Podemos amostrar f(x) se soubermos amostrar g(x)

(nao-normalizada) tal que g(x) ≥ f(x) para todo x.

1. gera x com dist. prop. a g(x)

2. aceita com prob. f(x)/g(x)

Os x aceitos terao distribuicao f(x). Aceitacao media:

A =

∫ ∞

−∞f(x) dx

∫ ∞

−∞g(x) dx

Para a distribuicao de Boltzmann

w(S) =e−βH(S)

Z; < A > =

∫

A(S) w(S) dS

nao ha esperanca de uma amostragem direta!

Solucao: Monte Carlo dinamico. “Inventamos”

uma evolucao temporal — uma dinamica

(markoviana) — de forma que as amostras geradas

no tempo tenham distribuicao w(S).

Cadeias de Markov

Processo estocastico X0, X1, . . . , Xt tal que

P (Xt+1|x0 . . . xt) = P (Xt+1|xt)

Historia determinada pela distribuicao inicial

P (X0) e pela matriz de transicao pxy.

Note: p(2)xy probabilidade de x → y em dois passos.

Claramente∑

y pxy = 1 para todo x.

Distribuicao estacionaria w(x)∑

x

w(x) pxy = w(y) para todo y

Cadeia aperiodica:

limt→∞

p(t)xy = w(y)

i.e. a cadeia converge para w independentemente

da distribuicao inicial.

Problema: dada uma cadeia (i.e. uma matriz pxy)

determinar se existe e qual e a distribuicao w(x).

Monte Carlo: qual deve ser pxy para que a

distribuicao estacionaria seja a distribuicao

w(x) que queremos amostrar, e.g. a distribuicao

de Boltzmann.

Daı medias temporais na cadeia convergem

(quando t → ∞) para medias na distribuicao de

equilıbrio w(x) do sistema considerado. Note:

desprezo o transiente inicial.

Condicoes para pxy

(A) Irredutibilidade (ergodicidade):

para x, y existe n tal que p(n)xy

(B) A distribuicao desejada e estacionaria:∑

x w(x) pxy = w(y)

Condicao suficiente:

(B’) Balanco detalhado: w(x) pxy = w(y) pyx

∑

x

w(x) pxy =∑

x

w(y) pyx = w(y)

Programa: seguir X(t) = xi e calcular medias

< A > =

∫

A(x) w(x) dx =

∑

i A(xi)

N

Problema: amostras nao sao independentes.

Correlacao normalizada (C(0) = 1):

C(k) =< Ai Ai+k > − < Ai >2

< A2i > − < Ai >2

Para k tal que C(k) → 0 as amostras sao

efetivamente independentes.

Metodo de Metropolis

Baseado em propor um passo x → y e aceitar ou

nao.

• se w(y)/w(x) ≥ 1 aceita

• em caso contrario aceita com probabilidade

w(y)/w(x)

Note: probabilidade de aceitacao:

Axy = min

{

1,w(y)

w(x)

}

Cadeia: pxy = Txy Axy (propor e aceitar).

Em geral considero Txy = Tyx.

Satisfaz balanco detalhado:

w(x) pxy = w(x) Txy min

{

1,w(y)

w(x)

}

= w(y) Tyx min

{

w(x)

w(y), 1

}

= w(y) pyx

Distribuicao de Boltzmann:

w(x) = e−βE(x)/Z ⇒ w(y)

w(x)= e−β∆E

• aceita se ∆E = E(y) − E(x) ≤ 0

• caso contrario aceita com prob. e−β∆E

Nota: se o passo for rejeitado mantem-se o

estado anterior.

Aplicacao ao Modelo de Ising

O Modelo de Ising

“spins” de dois estados preferem estar alinhados

S Si j

-J +J

S Si j

H(S) = −J∑

<i,j>

Si Sj − H∑

i

Si

• Energia: E =< H(S) >

• Calor especıfico: CV = ∂E/∂T

• Magnetizacao: M =<∑

i Si >

• Suscetibilidade: χ = ∂M/∂H

Metropolis para o Modelo de Ising

Percorre a rede, para cada sıtio propoe inversao do

spin. Modificacao proposta para o spin Si e

Si → −Si. Probabilidade de aceitacao

e−β E(y)

e−β E(x)=

e+βJ Si j n.n. i Sj

e−βJ Si j n.n. i Sj= e2 βJ Si hi

uma iteracao consiste em uma “varrida” completa

da rede. Ao final da iteracao calcula-se A(S) para

a configuracao gerada, e retoma-se o processo de

geracao de novas configuracoes.

Metodo do Banho Termico

Baseado na amostragem exata da distribuicao de

probabilidade condicional obtida para o sıtio i

mantendo-se todos os outros sıtios fixos. Varre-se a

rede, escolhendo um novo valor para Si independente

do antigo com a probabilidade

P (Si) = eβJSi � j n.n. iSj × const

apos normalizacao:

• P (Si = +1) = eβJhi/(eβJhi + e−βJhi) ≡ p

• P (Si = −1) = 1 − p

Calculo de Erros

• considerar apenas amostras efetivamente

independentes (atraves da analise da

correlacao temporal); erro por desvio

padrao, jack-knife, bootstrap.

• considerar todas as amostras e estimar o

erro levando em conta as correlacoes:

binning, metodo da janela auto-consistente.

Correlacoes

Media (de Monte Carlo) de A: A =1

N

N∑

i=1

Ai

Variancia: σ2A

=σ2

A

N

[

1 + 2

N−1∑

k=1

C(k)

]

=σ2

A

N(2 τ)

onde a correlacao temporal e

C(k) =< Ai Ai+k > − < Ai >2

< A2i > − < Ai >2

e τ e o tempo de auto-correlacao para o observavel A.

Tempo de Auto-correlacao

Considere C(k) = e−k/τ , τ grande (τ << N)

1 + 2N−1∑

k=1

C(k) ≈ 2∞∑

k=0

e−k/τ − 1

≈ 2 τ

∫ ∞

0

e−udu − 1 ≈ 2τ

Portanto defino

τ ≡ 1

2+

N−1∑

k=1

C(k)

Resumo∫

f(x) dµ → 1

N

N∑

i=1

f(xi); dµ =e−βH(x)

Zdx

xi com distribuicao de probabilidade µ

• MC estatico: xi independentes (erro ∼ 1/√

N)

• MC dinamico: simulacao de uma cadeia de Markov

com a distribuicao de equilıbrio µ (erro ∼√

τ/N)

Tempo de autocorrelacao τ relacionado ao critical

slowing-down: diverge na regiao crıtica para os

algoritmos locais. Solucao: algoritmos globais ou

metodos de sobre-relaxacao.

Metodos de Aglomerados

Modelo de Potts de q estados: Si = 1, 2, ..., q

H = −J∑

<i,j>

(δSi,Sj− 1)

Algoritmo de Swendsen-Wang: coloca elos com

probabilidade p entre sıtios com Si = Sj . O

modelo fica agora escrito em termos de

aglomerados interagentes de spins conectados.

Para p = 1 − e−βJ a energia de interacao e zero →aglomerados podem ser atualizados

independentemente. Note: baseado no modelo de

Fortuin-Kasteleyn, pode ser aplicado tambem com

campo magnetico.

Algoritmo de Wolff: escolhe-se um sıtio ao acaso e

ao redor dele forma-se um unico aglomerado.

Favorece aglomerados maiores → tempos de

auto-correlacao menores. Wolff tambem

generalizou o algoritmo para modelos de spins

contınuos.

Fenomenos Crıticos

Fluidos

��

T

LIQUIDO

GAS

P ponto

critico’

’

’

Magnetos��� ��� T

T

H = 0

T

MH

+

−

c

Expoentes Crıticos

t ≡ (T − Tc)/Tc

Quando t → 0:

Calor especıfico C ∼ |t|−α

Parametro de ordem M ∼ |t|β

Suscetibilidade χ ∼ |t|−γ

Comprimento de correlacao ξ ∼ |t|−ν

onde: < S0 Sx > ∼ e−x/ξ

Universalidade

Transicao de Desconfinamento

• transicao de fase no inıcio do universo

• natureza da transicao

• propriedades da fase de altas temperaturas

Percolacao

• n: densidade de sıtios (ou

elos) ocupados

• np: menor n tal que a origem

pertence ao aglomerado per-

colante

P ∼ (1 − n/np)βp n → np+

S ∼ (np − n)−γp n → np

Percolacao no Modelo de Ising

simple site (or bond) percolation fails to reproduce

Ising-model exponents

geometric cluster ≥ droplet

Solution: site percolation + T-dependent bond

probability

nb = 1 − e−2 J/KT

then get agreement for exponents, and percolation

point corresponds to Tc.

Percolacao de aglomerados de spins

Modelo O(4) 3d, L = 16, 24, 32

Tamanho medio dos clusters

No ponto crıtico ξ → ∞ e ha invariancia de escala

Leis de Escala

Fs(t, H) = b−d Fs(bytt, byHH)

ξ(t, H) = b ξ(bytt, byHH)

• Fs(t,H) = td/ytΦ(H/tyH/yt) → α,β,γ,ν in terms of

yt, yH (yt = 1/ν)

• Fs(t,H) = Hd/yH Ψ(t/Hyt/yH ) get δ

(Mt=0 ∼ H1/δ) and universal scaling function

M/H1/δ = f(t/H1/∆) where ∆ = βδ = ν yH

Finite-size scaling

Fs(t, H, uj , L) = L−d Q(L1/νt, L∆/νH, Lyjuj)

Q is universal, uj are irrelevant fields: yj < 0

M = L−β/νg(L1/νt, L∆/νH, L−ω)

Tomando H = 0 e desprezando correcoes mais altas

M = L−β/νg(L1/νt)

determino o expoente crıtico por um ajuste dos

dados para diversos L a t = 0.

• Binder Cumulant: at H = 0

UL ≡ −1

3

∂4Fs

∂H4 χ−2L−d = 1 − 1

3

<M4 >

<M >2

has constant (universal) value at Tc (neglecting ω)

• χ2 Method: expand FSS form of χ around t = 0

χ = Lγ/ν [c0 + (c1 + c2L−ω)tL1/ν + c3L

−ω]

• Finite-size-scaling extrapolation: use FSS form

A(β, 2L)

A(β,L)= fA(ξ(β,L)/L) + O(L−ω)

for obsvervable A together with ξ

Problemas de Otimizacao

Problemas de Otimizacao

e.g. minimizacao do funcional H(φ), ou

equivalentemente solucao do sistema linear

Aφ = f ⇒ H(φ) =1

2φ Aφ − fφ + c

atraves de um metodo iterativo.

Metodos Locais

• Jacobi (Gauss-Seidel): minimizacao local

exata

• sobre-relaxacao: combinacao com

atualizacao “micro-canonica”

Metodos Globais

• multi-grid

• gradiente conjugado

• aceleracao de Fourier

Minimizacao Global

Problemas com muitos mınimos/maximos

• sistemas com vınculos (e.g. φ · φ = 1)

• frustracao (e.g. vidros de spins)

Metodos Estocasticos

• recozimento simulado

• algoritmo genetico

Simulacoes Numericas

• experimentos que nao podemos/queremos realizar

(avioes, guerra nuclear, evolucao)

• problemas sem solucao analıtica; sistemas

complexos: nao-linearidade, fenomenos crıticos.

Nota: independencia de detalhes (universalidade)

• comportamento dinamico (e.g. automato celular),

regras de evolucao locais ⇒ processamento

paralelo. Nota: dinamica pode ser introduzida

• tratamento teorico, com aspectos experimentais

(dados, erros, “medidas” no tempo)