[openlecture] Simulated Annealing(모의 담금질?)에 대해 알아봅시다.

  • 일루
    일루

    icpc의 세상에는 최적해를 구하는 문제들이 많지만, 실제 세상에는 최적해를 구하기에는 상태 공간이 너무 커서 시간이 감당할 수 없을 정도로 오래 걸리는 문제들이 많습니다.

    예제로는 TSP, vertex cover나 여러 가지 실생활의 문제들이 있습니다. (어떤 것들이 있을까요? ^^;;)

    이렇게 상태 공간이 아주 넓은 문제들을 푸는 방법 중 가장 만병통치약스러운 Simulated Annealing을 알아보도록 하겠습니다.

    일단 Local Search부터 알아보도록 하겠습니다.

    초기 상태를 정한다;
    while (시간이 될때까지) {
      현재 상태에서 어느 정도 변화를 시킨 상태를 만들어내고 다음 상태라고 이름을 붙인다.
      다음 상태를 평가한다;
      if (다음 상태 평가치 >= 현재 상태 평가치) 현재 상태 = 다음 상태;
    }

    이렇게 하면 됩니다.

    이 방법의 문제점이라면 부분적인 최적해에서 벗어나지 못할 수 있다는 것인데요, 10차원 이상으로 생긴 상태공간을 생각해 볼 때, 부분적으로 점수가 좋은 지역이 있을 수 있지만 최적해와는 거리가 멀고 점수 차이도 꽤 날 수 있겠지요.

    이러한 문제를 해결하기 위해서 Simulated Annealing이라는 기법이 나왔습니다.

    이 방법은 실제 담금질을 모티브로 하여 만들어진 것으로, 철강 등을 뜨겁게 가열했다가 천천히 식어가면서 원자들이 재배열되며 강도가 좋아지는 과정을 따라서 해 보는 것이 되겠습니다. 여기서 온도가 뜨거울 수록 원자들이 그 에너지를 사용하여 부분적으로 에너지가 더 높은 상태로 이행하는 확률이 높아집니다.

    while (시간이 될때까지) {
      현재 상태에서 어느 정도 변화를 시킨 상태를 만들어내고 다음 상태라고 이름을 붙인다.
      다음 상태를 평가한다;
      차이 = 다음 상태 평가치 - 현재 상태 평가치
      if (차이 >= 0 이거나 차이 < 0 이고 차이가 용인 가능한 정도일 때) {
        현재 상태 = 다음 상태;
        현재 상태가 가장 좋은 상태이면 최고 상태에 기억해 놓는다; (옵션)
      }
    }

    이렇게 하면 됩니다.

    SA를 구성하는 여러 가지 요소에 대해서 알아봅시다.

    (1) 다음 상태를 어떻게 만들 것인가?

    현재 상태에서 어느 정도 변화를 시키면서, 평가하기는 쉽게 만들어야 합니다. 예를 들어 TSP 문제를 풀 때 SA를 사용한다고 하면, 다음 상태는 인접한 두 점을 바꾼다 정도로 할 수 있겠죠. 평가는 참 쉽습니다 ^^; 인접한 두 점이 아니라 아무 두 점을 바꾼다고 해도 됩니다. 이래도 평가가 상수 시간에 됩니다.

    약간 고급 기술로는 아무렇게나 바꿔보는 것이 아니라 좋아질 확률이 높은 방향으로 바꿔보는 것이 있습니다. (TSP의 경우에는 어떻게 하면 될까요?)

    (2) 상태를 평가한다.

    이것은 SA의 요소라기보다는 평가 함수를 빠르게 동작하게 짤 수 있으면 SA의 전체 iteration 수가 늘어나기 때문에 넣어보았습니다. 중복된 계산은 철저히 피해야 하며 계산이 복잡하게 되는 상태 변화도 피해야 하겠습니다.

    (3) 온도

    가장 쉽게는 T가 시간에 연동되게 하는 것입니다. 시간 0에서는 (절대 온도)1000도였다가 마지막에는 0도가 되면 되겠죠.

    구현상으로는 T를 구하려면 gettime을 불러야 하는데, Topcoder MM에서는 gettime을 0.01~0.1초 정도에 한 번씩 부르도록 해야 합니다. 0.01초보다 자주 부르면 system call이 과다해져서 Topcoder server에서 극도의 비효율을 보일 수 있으며, 0.1초보다 덜 자주 부르면 TLE를 불러올 수 있습니다.

    (4) 용인 가능한 차이는?

    온도가 0도일 때는 로컬 서치가 됩니다. 즉 용인 가능한 차이는 0입니다.

    보통은 exp(차이/상수/T) > rand() 면 상태 업데이트를 하도록 하는 방법을 많이 사용합니다.

    여기에서 상수를 정하는 것은 문제의 특성에 따라서 보고 임의로 정할 수도 있고, 인접한 상태 사이의 차이 정도를 보고 적절히 조절되도록 할 수도 있습니다. (Adaptive SA)

    경우에 따라서 차이가 0 이상일 때도 꼭 용인하도록 하지는 않을 수도 있습니다.

    아주 간단하게 썼는데요 특히 초보자 분들께 많은 도움이 되었으면 좋겠습니다 ^^;;


    13년 전
7개의 댓글이 있습니다.
  • JongMan
    JongMan

    우와 예제 소스 코드도 보여주세요 ㅋ


    13년 전 link
  • hyunhwan
    hyunhwan

    위에서 일루님이 언급하신 내용들이 반영된 최근에 제가 파고 있던 MM26을 위한 SA를 이용한 코드입니다 :)

    문제의 링크는 다음과 같습니다.
    http://www.topcoder.com/longcontest/?module=ViewProblemStatement&rd=11148&pm=8426

    문제를 요약하자면 다음과 같습니다.

    N행 N열의 실수로 이뤄진 행렬 W가 주어졌을 때, 다음의 수식의 값을 최대화 하는 0이상 N-1이하의 숫자로 이뤄진 길이 N의 순열 P를 찾는 프로그램을 작성하라.

    아래의 소스코드는 순열의 임의의 두 원소를 잡아 swap을 한 결과를 다음 상태로 간주하고 SA를 이용하여 구현된 코드입니다.

    #include <sys/time.h>
    #include <cassert>
    #include <vector>
    #include <cstring>
    #include <algorithm>
    #include <cstdio>
    #include <cmath>
    #include <ctime>
    #include <cstdlib>
    using namespace std;
    // WELL512 random number generator
    #define W 32
    #define R 16
    #define P 0
    #define M1 13
    #define M2 9
    #define M3 5
    
    #define MAT0POS(t,v) (v^(v>>t))
    #define MAT0NEG(t,v) (v^(v<<(-(t))))
    #define MAT3NEG(t,v) (v<<(-(t)))
    #define MAT4NEG(t,b,v) (v ^ ((v<<(-(t))) & b))
    
    #define V0            STATE[state_i                   ]
    #define VM1           STATE[(state_i+M1) & 0x0000000fU]
    #define VM2           STATE[(state_i+M2) & 0x0000000fU]
    #define VM3           STATE[(state_i+M3) & 0x0000000fU]
    #define VRm1          STATE[(state_i+15) & 0x0000000fU]
    #define VRm2          STATE[(state_i+14) & 0x0000000fU]
    #define newV0         STATE[(state_i+15) & 0x0000000fU]
    #define newV1         STATE[state_i                 ]
    #define newVRm1       STATE[(state_i+14) & 0x0000000fU]
    
    #define FACT 2.32830643653869628906e-10
    
    static unsigned int state_i = 0;
    static unsigned int STATE[R];
    static unsigned int z0, z1, z2;
    
    void InitWELLRNG512a () {
      int j;
      state_i = 0;
      for (j = 0; j < R; j++)
        STATE[j] = time(NULL);
    }
    
    double Random() {
      z0    = VRm1;
      z1    = MAT0NEG (-16,V0)    ^ MAT0NEG (-15, VM1);
      z2    = MAT0POS (11, VM2)  ;
      newV1 = z1                  ^ z2; 
      newV0 = MAT0NEG (-2,z0)     ^ MAT0NEG(-18,z1)    ^ MAT3NEG(-28,z2) ^ MAT4NEG(-5,0xda442d24U,newV1) ;
      state_i = (state_i + 15) & 0x0000000fU;
      return ((double) STATE[state_i]) * FACT;
    }
    
    int RandomI( int mini, int maxi ) {
      // Output random integer in the interval min <= x <= max
      // Relative error on frequencies < 2^-32
      if (maxi <= mini) {
        if (maxi == mini) return mini; else return 0x80000000;
      }
      // Multiply interval with random and truncate
      int r = int((double)(unsigned int)(maxi - mini + 1) * Random() + mini); 
      if (r > maxi) r = maxi;
      return r;
    }
    
    inline double get_time() {
        timeval tv;
        gettimeofday(&tv, 0);
        return tv.tv_sec+1e-6*tv.tv_usec;
    }
    
    int N;
    double T;
    int WW[1024][1024];
    int Q[1024][1024];
    float score_norm;
    double start_time;
    
    
    struct Permute {
      int score( int PP[] ) {
          int ret = 0;
          for(int i=0;i<N;++i) for(int j=i+1;j<N;++j) ret += WW[PP[i]][PP[j]];
          return ret;
      }
      vector<int> findOrder( vector<double> w ) {
        start_time = get_time();
        InitWELLRNG512a();
        N = (int)sqrt( (double)w.size() );
        score_norm = (N*sqrt(N*log(N)-N) / 100.0);
    
        for(int i=0;i<N;++i) 
          for(int j=0;j<N;++j) 
            WW[i][j] = (int)(w[i*N+j] / score_norm * 10000000);
        for(int i=0;i<N;++i) 
          for(int j=0;j<N;++j)
            Q[i][j] = WW[i][j] - WW[j][i];
    
        int PP[1000];
        int BP[1000];
    
        for(int i=0;i<N;++i) PP[i] = i;
        memcpy( BP, PP, sizeof(int)*N );
        int best = score( BP );
        int l, r; 
        int nIter = 0;
        int E1 = best;
    
        double start_t = 5e5;
        if( N >= 500 ) start_t = 1e5;
        const double time_limit = 29.0;
    
        while( 1 ) {
          if (!(nIter&((1<<16)-1))) {
            double elapsed = get_time() - start_time;
            if ( elapsed > time_limit ) break;
            T = start_t * ( 1 - elapsed / time_limit );
          }
          l = RandomI( 0, N-1 );
          r = RandomI( 0, N-1 );
          if( l == r ) continue;
          if( l > r ) swap(l,r);
          int E2 = Q[PP[r]][PP[l]] + E1;
          for(int i=l+1;i<r;++i) 
            E2 += Q[PP[r]][PP[i]] - Q[PP[l]][PP[i]];
          if( exp( (E2-E1) / T ) > Random() ) {
            E1 = E2;
            swap( PP[l], PP[r] );
          }
          if( best < E2 ) {
            // fprintf( stderr, "updated : %.3fsec(%d-th iteration) %f->%f(%f)\n", get_time() - start_time, nIter, best / 10000000.0 , E2 / 10000000.0, T );
            best = E2;
            memcpy( BP, PP, sizeof(int)*N );
          }
          ++nIter;
        }
    
    
        // fprintf( stderr, "nIter : %d\n", nIter );
        vector<int> ret(N);
        for(int i=0;i<N;++i) ret[i] = BP[i];
        return ret;
      }
    };
    

    13년 전 link
  • Megalusion
    Megalusion

    우와 리!베!


    13년 전 link
  • JAYS
    JAYS

    일루님이 쓰신 글과 LIBe님께서 작성하신 코드 감사합니다. 코드를 보면 시간이 지날수록 용인 가능한 차이가 커지는 것 같습니다. 이를 프로그램의 초기단계에서는 비교적 좁은 범위의 공간을 탐색하다가 시간이 지날수록 넓은 범위의 공간을 시도하는 것으로 해석해도 되는지요? 반대방향으로 (넓은 범위에서 좁은 범위) 탐색하는 것에 비해 어떤 장점이 있는지요?


    9년 전 link
  • Being
    Being

    반대로 해석하신 것 같습니다. 시간이 지날수록 용인 가능한 차이가 줄어듭니다.


    9년 전 link
  • JAYS
    JAYS

    아...... 잘못생각했네요... exp(차이/상수/T) > rand() 에서 차이가 0 이상이라고 생각했네요;; 답변 감사합니다!


    9년 전 link
  • JAYS
    JAYS

    정말 많은 도움이 되는 글인 것 같아요.


    9년 전 link
  • 정회원 권한이 있어야 커멘트를 다실 수 있습니다. 정회원이 되시려면 온라인 저지에서 5문제 이상을 푸시고, 가입 후 7일 이상이 지나셔야 합니다. 현재 문제를 푸셨습니다.