Simulated Annealing [algorithm] by 바죠

Simulated Annealing (SA)
http://en.wikipedia.org/wiki/Simulated_annealing
http://mathworld.wolfram.com/SimulatedAnnealing.html

 Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. (1983). "Optimization by Simulated Annealing". Science 220, 671–680.


모사된 풀림, 풀림 시늉이라고 번역되는 것이 적절하다.


global optimization(특정한 목적 함수를 최소화하는 것, 다양한 변수들에 대해서 목적 함수를 최소화 하는 것)을 수행하는 하나의 방법이다. [목적 함수를 최대화 하는 것은 최소화 하는 것과 사실상 동일한 문제이다. g(x):=-f(x) 처럼 새로운 함수를 정의해 버리면 되기 때문이다.] 아주 오래된 이론이다. 많이 회자되는 알고리듬이다. 기본적인 알고리듬이다. 이것으로 모든 것을 제대로 풀수는 없는 실정이다. 현실적으로 보다 더 우수한 알고리듬들이 제안되어 있는 상황이다. 실제 문제 풀이에서는 많이 사용된다.  하지만, 역사적으로 중요한 의미를 가지고 있다. 또한, 이것의 변형들이 많이 존재한다. 실제 온도를 바꾸지 않는다면, 주어진 온도에서 시스템이 가질 수 있는 상태들, 열역학적으로 허용된 상태들을 이론적으로는 모두 샘플링할 수 있다. 주어진 온도에서 가능한 상태를 모두 샘플링하는 Monte Carlo 방법을 그대로 사용하는 것이다. 풀림 시늉에서는 온도를 천천히 내리는 단계적인 추가 계산들이 있다.


하나의 구조(입자)가 특정 구조(위치)에서 에너지를 주고, 또 다른 구조에서 또 다른 에너지를 줄 경우,  온도에 따라서 이 새로운 구조를 수용하거나 거부한다.   상황이 바뀌면 목적 함수도 변화한다. 이 변화된 목적 함수값과 현재의 온도를 직접적으로 비교한다.  다시 말해서 주어진 온도에서 어느 정도의 목적 함수 변화는 수용할 수도 있고 수용 못할 수도 있다.  온도가 높을 수록 목적 함수가 다소 크게 변화된 것을 허용할 수 있게 된다.  온도가 굉장히 낮을 경우, 목적 함수가 증가하는 것은 매우 받아 들이기 어려운 것이다.  주어진 온도에서는 메트로폴리스(Metropolis) 알고리듬을 그대로 수용한다.  주어진 온도에서 충분히 시뮬레이션을 한다. 그리고 난 다음에 온도를 낮추어준다.  다만, 온도를 높은 온도에서 낮은 온도로 낮추어준다. 실제 금속의 풀림을 알고리듬으로 변환한 것이다. 그래서 이름이 simulated annealing이다.

0과 1사이의 랜덤 넘버(확률에 해당하는 것)를 활용한 알고리듬이다. 특정 확률 이상일 경우 변화를 받아 들이는 것이다.  이 때, 계산된 확률과 랜덤 넘버를 비교하게 된다.  확률적으로 변화를 받아 들인다는 뜻이다.  Monte Carlo 계산에 기본이 되는 알고리듬이다. 분명히 uphill search를 허용한다. 이는 확률적으로 허용된다.  반면, downhill search는 무조건 허용한다.

메트로폴리스(Metropolis) 알고리듬에서는, 주어진 온도 T1 에서 conformation이 변환될 때,  즉, x ---> x' 으로 변환 가능한지를 따지게 된다. 이 때, conformation x' 은 굉장히 일반적으로 만들어진 것일 수 있다. 랜덤으로 만들어져도 좋다.  물론, 너무 많이 거부되는 방식으로 만들어지면 안된다.  50% 정도는 거부되고 나머지는 허용될 수 있게 만들어지면 이상적이라고 할 수 있겠다.

x'을 어떻게 만드느냐에 따라서 알고리듬이 달라지고 확장될 수 있다.  유사하지만 굉장히 다른 알고리듬이 만들어 질 수 있다는데 주의해야 한다.

x ---> x' 로 변환 할 때, 에너지가 E, E'으로 바뀌게 된다. 이 때, 온도 T, beta=1./(kBT).
min[1, exp(-beta E')/exp(-beta E)] > r, 여기서, r은 랜덤 넘버이다. 0.과 1사이에 존재하는 숫자이다.

T ---> 0으로 가는 영역에서는 사실상 quench를 실행할 수 있다.  목적 함수값이 낮아지는 경우만 업데이트하는 경우. 이 경우는 사실상, 목적 함수 값이 낮아지지 않으면, x ---> x' 변환이 허용되지 않는다. 확률적으로 허용되지 않는다. 의도적으로 이러한 목적으로 활용될 수 있다. 이것은 국소 최적화를 의미하는 것이다.  물론, gradient를 활용하지는 않았다.  

simulated annealing 알고리듬에는 반드시 온도를 내려주는 방법이 존재해야만 한다. 특별하게 온도를 고정하고 지속적으로 높은 온도에서 만들어진 x를 낮은 온도 쪽으로 공급을 할 수 있다면,  결국 매우 낮은 온도 쪽에서는 사실상 simulated annealing 의 마지막 단계를 수행하는 꼴이 된다.  이렇게 여러개의 온도 분포에 대해서 각각의 Monte Carlo 시뮬레이션을 수행함과 동시에 지속적인 x의 교환을 허용한다면,  새로운 방식의 simulated annealing이라고 볼 수 있다.  여러 개의 온도 분포를 두고 온도 분포 사이의 해들을 교환할 수 있다면 이것은 새로운 방식의 simulated annealing이라고 볼 수 있다.


온도가 사실상 0 으로 수렴한 경우는 무조건 에너지가 낮아질때에만 x ---> x' 변환이 허용되는 경우에 해당한다.  특별한  에너지 최소화 과정으로 볼 수 있다.

http://en.wikipedia.org/wiki/Metropolis_algorithm



그림 출처:
http://pubs.rsc.org/en/content/articlelanding/2005/cp/b509983h

아래의 프로그램은 피팅 프로그램이다.  그런데 simulated annealing 방법을 활용한 것이다. 즉, 목적 함수를 특정 다항 함수 형식으로 고정하고 주어진 데이터와의 차이가 최소화 되도록 프로그램을 만들었다. 데이터 포인트 당 다항식과 데이터의 차이가 최소화 되도록 만드는 것이 목표이다. 다시 말해서, 아주 좋은 다항식 계수를 찾는 것이 목표이다. 이 다항식 계수는 결국, 다항식을 줄 것이고, 그 다항식은 실험 데이터를 가장 잘 표현하는 다항식이 될 것이다.

       module sa_mod
!      Written by In-Ho Lee, KRISS, December, 2012.
       implicit none
       private
       save
       integer ndata,ndeg, nseed,nsim,ntem, nsize
       integer, allocatable :: iseed(:)
       real*8 temp,obj_lowest, act,rjt
       integer, parameter :: ntry=2
       real*8 act0(0:ntry-1),rjt0(0:ntry-1)
       real*8, allocatable :: xpt(:),ypt(:)
       real*8, allocatable :: qq(:),pp(:),qq_lowest(:)
       real*8 poly

       public :: initial,final,do_sa

       contains

       subroutine initial(nsim0,ntem0,nseed0,temp0)
       implicit none
       integer nsim0,ntem0,nseed0
       real*8 temp0
       real*8 tmp
       integer k
       logical lexist

       nsim=nsim0
       ntem=ntem0
       temp=temp0
       nseed=nseed0

       call random_seed(size=nsize)
       allocate(iseed(nsize))
       do k=1,nsize
       iseed(k)=nseed+k
       enddo
       call random_seed(put=iseed)
       deallocate(iseed)

       ndata=100

       allocate(xpt(ndata))
       allocate(ypt(ndata))
       call preparation()

       ndeg=2
       allocate(pp(0:ndeg))
       allocate(qq(0:ndeg))
       allocate(qq_lowest(0:ndeg))
       qq=0.d0
       do k=0,ndeg
       call random_number(tmp)
       tmp=tmp-0.5d0
       qq(k)=tmp
       enddo

       inquire(file='fort.11',exist=lexist)
       if(lexist)then
       open(11,file='fort.11',form='formatted')
       read(11,*) k
       read(11,*) qq_lowest
       read(11,*) tmp
       close(11)
       qq=qq_lowest
                 endif

       end subroutine initial

       subroutine final()
       implicit none
       integer j

       open(7,file='fort.7',form='formatted')
       do j=1,ndata
       write(7,*) xpt(j),ypt(j)
       enddo
       write(7,*) '&'
       do j=1,ndata
       write(7,*) xpt(j),poly(xpt(j),qq_lowest,ndeg)
       enddo
       close(7)
       write(6,*) 'obj_lowest,qq_lowest'
       write(6,*) obj_lowest
       write(6,*) qq_lowest

       open(11,file='fort.11',form='formatted')
       write(11,*) ndeg
       write(11,*) qq_lowest
       write(11,*) obj_lowest
       close(11)

       deallocate(qq,pp,qq_lowest)
       deallocate(xpt,ypt)
       end subroutine final

       subroutine preparation()
       implicit none
       real*8 tmp,xi,xf,dx
       integer j

       xi=0.d0
       xf=10.d0
       dx=(xf-xi)/float(ndata-1)
       do j=1,ndata
       xpt(j)=xi+dx*float(j-1)
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*2.0d0
       ypt(j)= 0.1d0 +0.2d0*xpt(j) +0.3d0*xpt(j)**2    +tmp
       enddo
       end subroutine preparation


       subroutine do_sa()
       implicit none
       real*8 objqq,objpp,dell, tmp
       integer item,isim,itry,k

       call cal_obj(qq,objqq)
       obj_lowest=objqq ; qq_lowest=qq
       do item=1,ntem
       print*, temp,'temp'
       act=0.d0 ; rjt=0.d0 ; act0=0.d0 ; rjt0=0.d0
       do isim=1,nsim
       itry=mod(isim,ntry)

       if(itry == 0)   then
       pp=qq
       call random_number(tmp)
       k=tmp*(ndeg+1)
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*sqrt(temp)
       pp(k)=pp(k)+tmp
       elseif(itry ==1)then
       pp=qq
       do k=0,ndeg
       call random_number(tmp)
       tmp=tmp-0.5d0
       tmp=tmp*sqrt(temp)
       pp(k)=pp(k)+tmp
       enddo
                       endif

       call cal_obj(pp,objpp)
       dell=(objpp-objqq)/temp
       if(dell > 50.d0) dell=50.d0
       if(dell < -50.d0) dell=-50.d0
       call random_number(tmp)
       if(tmp < min(1.d0,exp(-dell)))then
       qq=pp ; objqq=objpp
       act=act+1.d0
       act0(itry)=act0(itry)+1.d0
                                     else
       rjt=rjt+1.d0
       rjt0(itry)=rjt0(itry)+1.d0
                                     endif

       if(objqq < obj_lowest)then
       obj_lowest=objqq ; qq_lowest=qq
                             endif

       enddo
       do itry=0,ntry-1
       write(6,'(f18.8,2x,a3,2x,2f20.3)') act0(itry)/(act0(itry)+rjt0(itry)),'(%)',act0(itry),rjt0(itry)
       enddo
       write(6,'(f18.8,2x,a3,2x,2f20.3)') act/(act+rjt),'(%)',act,rjt
       print*, temp,objqq,obj_lowest
       temp=0.9d0*temp
       if(temp <= 1.d-6) temp=1.d-6
       enddo

       end subroutine do_sa

       subroutine cal_obj(q,objq)
       implicit none
       real*8 q(0:ndeg),objq
       integer i

       objq=0.d0
       do i=1,ndata
       objq=objq+(ypt(i)-poly(xpt(i),q,ndeg))**2
       enddo
       objq=objq/float(ndata)
       end subroutine cal_obj
      real*8 function poly(x,q,ndeg)
       implicit none
       integer ndeg
       real*8 x,q(0:ndeg)
       real*8 tmp
       integer i

       tmp=0.d0
       do i=0,ndeg
       if(i == 0)then
       tmp=tmp+q(i)
                 else
       tmp=tmp+q(i)*x**float(i)
                 endif
       enddo
       poly=tmp
       end function poly

       end module sa_mod
!234567890
       program test
       USE sa_mod, ONLY : initial,final,do_sa
       implicit none
       integer nsim0,ntem0,nseed0
       real*8 temp0
       integer itemp,irate,itemq

       call system_clock(itemp,irate)


       nsim0=100000
       ntem0=10
       nseed0=1234
       temp0=1.d0

       call initial(nsim0,ntem0,nseed0,temp0)

       call do_sa()

       call final()

       call system_clock(itemq)
       write(6,'(2e15.8,2x,a6)') float(itemq-itemp)/float(irate)/60., float(itemq-itemp)/float(irate)/3600.,' min,h'
       end program test

-------------------------------------------------------------------------------------
복제품 맞바꿈 알고리듬 (super simulated annealing):
http://incredible.egloos.com/4587643


유전 알고리듬:
http://incredible.egloos.com/4856511


differential evolution:
http://incredible.egloos.com/5265158


Conformational Space Annealing (CSA):

http://incredible.egloos.com/6208916
-------------------------------------------------------------------------------------

특히, 복제품 맞바꿈 알고리듬은 super simulated annealing이라고 부를 수도 있다.  온도가 높은곳에서 다양한 구조들을 지속적으로 생성하고, 온도가 낮은곳으로 새로운 구조들을 지속적으로 공급한다. 온도가 낮은곳에서는 사실상 매우 낮은 온도에서 Simulated Annealing를 하는 것으로 볼 수 있다. 온도 분포가 있는 것을 계속 유지하고 있다.  하지만, 가장 온도가 낮은 곳에서 항상 Simulated Annealing (풀림 시늉)의 결과를 받아 볼 수 있다.

 
--------------------------------------------------------------------------------------
Monte Carlo 방법: 메트로폴리스 방법
http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm
https://en.wikipedia.org/wiki/Monte_Carlo_method

몬테칼로는 특정 지역 이름이다. http://en.wikipedia.org/wiki/Monte_Carlo
모나코라는 나라, 프랑스에 붙어 있는 조그만 왕국이다. 
모나코는 AS 모나코라는 축구팀으로 우리에게 친숙한 곳이다.
이 나라의 행정부가 있는곳이다.

몬테칼로는 주사위를 던져서 시뮬레이션을 수행하는 방법을 지칭한다. 
실제 아래의 예에서 알 수 있듯이 진짜로 주사위를 던진다.
계산된 확률이 주사위를 던져서 나온확률보다 큰지 작은지를 판단한다.
실로 놀라운일이다. 이렇게 계산을 할 수 있다는 것이 놀라운일이다.
http://en.wikipedia.org/wiki/Dice

계산 물리학, 통계 물리학 영역에서는 가장 위대한 알고리듬 중의 하나로 평가되고 사용되고 있다.
예를 들어, 주어진 온도에서 가능한 상태들을 계산할 수 있고, 온도의 함수로 내부에너지, 비열등을 직접계산할 수 있다.
주어진 해밀토니안에서 가능한 상태들을 골고루 다룰 수 있고,
각 상태들마다 에너지를 계산할 수 있다면 Monte Carlo 시뮬레이션은 가능하다.
온도에 의존하는 물리량의 계산이 직접적으로 가능하다.   

{x} --> {x'} 형태로의 구조 변환, 상태 변환을 만들 수 있어야 한다. 매우 일반적인 변환을 만들 수 있으면 더욱더 좋다.
E({x})  --> E({x'}) 위의 변환에 대응하는 에너지의 변환이 계산 가능해야 한다.

온도가 적당히 높은 경우에는, 구조 변환에 따라서 에너지가 다소 높아질 경우에도 확률적으로 그러한 변환이 받아들여진다.
물론, 에너지가 낮아지는 그러한 변환은 모두 받아들여진다. 100% 받아지게 된다.

온도가 0으로 수렴하는 경우에는, 특별히, 더 높은 에너지를 가지는 구조로의 변환은 사실상 받아들여지지 않는다.
이 경우 모든 변환이 사실상 거부된다. 이 경우는 사실상 quenching이라고 부른다. 계속해서 낮은 에너지만 찾게 된다. 
본질적으로 가장 에너지가 낮은 구조를 찾는 방법으로 활용될 수 있다. 
얼마나 효율적으로 {x} --> {x'} 형태로의 구조 변환, 상태 변환을 만들 수 있느냐가 관건이다. 

주어진 온도에서 충분히 수렴된 계산을 진행할 수 있다.
관심이 있는 모든 온도에서 동일한 계산을 수행할 수 있다.
온도의 함수로 에너지를 그려 볼 수 있다. 
E=E(T)

https://en.wikipedia.org/wiki/Grace_Kelly


핑백


최근 포토로그