담금질 기법¶
담금질 기법(Simulated Annealing, SA)은 함수의 전역 최적해를 근사하는 무작위 알고리즘입니다. 탐색 과정에서 어느 정도의 무작위성을 사용하기 때문에 무작위 알고리즘이라고 부르며, 따라서 같은 입력에 대해서도 출력이 달라질 수 있습니다.
문제¶
상태 $s$의 에너지를 계산하는 함수 $E(s)$가 주어집니다. 우리가 할 일은 $E(s)$가 최소가 되는 상태 $s_{best}$를 찾는 것입니다. SA는 상태가 이산적이고 $E(s)$에 지역 최솟값이 여러 개 있는 문제에 적합합니다. 외판원 문제(Travelling Salesman Problem, TSP)를 예로 들어 보겠습니다.
외판원 문제(TSP)¶
2차원 공간에 있는 노드들의 집합이 주어집니다. 각 노드는 $x$ 좌표와 $y$ 좌표로 나타냅니다. 노드들을 방문할 순서를 정할 때, 그 순서대로 방문하면 이동 거리가 최소가 되는 순서를 찾는 것이 과제입니다.
동기¶
담금질(annealing)은 재료를 가열한 뒤 식혀서 내부 원자들이 내부 에너지가 최소인 배열로 재배치되게 하는 금속 공정이며, 그 결과 재료의 성질이 달라집니다. 상태는 원자들의 배열이고, 내부 에너지는 최소화할 함수입니다. 원자들의 원래 상태는 내부 에너지의 지역 최솟값이라고 생각할 수 있습니다. 재료의 원자들을 재배치하려면 전역 최솟값에 도달할 수 있도록 내부 에너지가 최소가 아닌 영역을 지나가게끔 재료에 동기를 부여해야 합니다. 이 동기는 재료를 더 높은 온도로 가열하여 부여합니다.
담금질 기법은 말 그대로 이 과정을 시뮬레이션합니다. 임의의 상태(재료)에서 시작하여 높은 온도를 설정합니다(가열합니다). 이제 알고리즘은 높은 온도라는 동기를 받았으므로 현재 상태보다 에너지가 더 높은 상태도 받아들일 준비가 되어 있습니다. 이 덕분에 알고리즘이 지역 최솟값에 갇히지 않고 전역 최솟값을 향해 나아갈 수 있습니다. 시간이 흐르면 알고리즘이 식으면서 에너지가 더 높은 상태를 거부하고, 지금까지 찾은 가장 가까운 최솟값으로 이동합니다.
에너지 함수 E(s)¶
$E(s)$는 최소화(또는 최대화)해야 하는 함수입니다. 이 함수는 모든 상태를 실수에 대응시킵니다. TSP의 경우, $E(s)$는 상태에 있는 노드의 순서대로 한 바퀴 원을 그리며 이동한 거리를 반환합니다.
상태¶
상태 공간은 에너지 함수 $E(s)$의 정의역이며, 상태는 상태 공간에 속하는 임의의 원소입니다. TSP의 경우 모든 노드를 방문하기 위해 택할 수 있는 모든 경로가 상태 공간이며, 이 경로들 중 어느 하나를 상태로 볼 수 있습니다.
이웃 상태¶
이전 상태와 가까운, 상태 공간 안의 상태입니다. 이는 보통 원래 상태에 간단한 변환을 적용하여 이웃 상태를 얻을 수 있다는 뜻입니다. 외판원 문제의 경우, 무작위로 노드 2개를 골라 현재 상태에서 두 노드의 위치를 맞바꾸어 이웃 상태를 얻습니다.
알고리즘¶
무작위 상태 $s$에서 시작합니다. 매 단계에서 현재 상태 $s$의 이웃 상태 $s_{next}$를 고릅니다. $E(s_{next}) < E(s)$이면 $s = s_{next}$로 갱신합니다. 그렇지 않으면 $s_{next}$로 이동할지 $s$에 머무를지를 결정하는 확률 수락 함수 $P(E(s),E(s_{next}),T)$를 사용합니다. 여기서 T는 온도이며, 처음에는 높은 값으로 설정되고 단계마다 천천히 감소합니다. 온도가 높을수록 $s_{next}$로 이동할 가능성이 커집니다. 동시에 모든 반복에 걸쳐 가장 좋은 상태 $s_{best}$도 기록해 둡니다. 수렴하거나 시간이 다 될 때까지 진행합니다.
극댓값이 여러 개인 이 함수의 최댓값을 찾는 담금질 기법(simulated annealing)을 시각적으로 나타낸 것입니다.
온도(T)와 감쇠(u)¶
시스템의 온도는 알고리즘이 더 높은 에너지를 가진 상태를 얼마나 기꺼이 받아들이는지를 수치로 나타냅니다. 감쇠는 알고리즘의 "냉각 속도"를 수치로 나타내는 상수입니다. 냉각 속도가 느릴수록(더 큰 $u$) 더 좋은 결과를 낸다고 알려져 있습니다.
확률 수락 함수(PAF)¶
$P(E,E_{next},T) = \begin{cases} \text{True} &\quad\text{if } \mathcal{U}_{[0,1]} \le \exp(-\frac{E_{next}-E}{T}) \\ \text{False} &\quad\text{otherwise}\\ \end{cases}$
여기서 $\mathcal{U}_{[0,1]}$는 $[0,1]$ 위의 연속 균등 난수입니다. 이 함수는 현재 상태, 다음 상태, 온도를 입력받아 불리언 값을 반환하며, 이 값은 탐색이 $s_{next}$로 이동할지 $s$에 머무를지를 알려 줍니다. $E_{next} < E$ 이면 이 함수는 항상 True를 반환하며, 그렇지 않더라도 $\exp(-\frac{E_{next}-E}{T})$의 확률로 이동할 수 있다는 점에 유의하십시오. 이 확률은 깁스 측도에 해당합니다.
bool P(double E,double E_next,double T,mt19937 rng){
double prob = exp(-(E_next-E)/T);
if(prob > 1) return true;
else{
bernoulli_distribution d(prob);
return d(rng);
}
}
코드 템플릿¶
class state {
public:
state() {
// Generate the initial state
}
state next() {
state s_next;
// Modify s_next to a random neighboring state
return s_next;
}
double E() {
// Implement the energy function here
};
};
pair<double, state> simAnneal() {
state s = state();
state best = s;
double T = 10000; // Initial temperature
double u = 0.995; // decay rate
double E = s.E();
double E_next;
double E_best = E;
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
while (T > 1) {
state next = s.next();
E_next = next.E();
if (P(E, E_next, T, rng)) {
s = next;
if (E_next < E_best) {
best = s;
E_best = E_next;
}
E = E_next;
}
T *= u;
}
return {E_best, best};
}
사용 방법:¶
상태 클래스의 함수들을 상황에 맞게 채우십시오. 최솟값이 아니라 전역 최댓값을 찾으려는 경우, $E()$ 함수가 최대화하려는 함수의 음수를 반환하도록 하고 마지막에 $-E_{best}$를 출력하십시오. 아래 매개변수들을 필요에 맞게 설정하십시오.
매개변수¶
- $T$ : 초기 온도입니다. 탐색을 더 오래 실행하려면 더 높은 값으로 설정하십시오.
- $u$ : 감쇠율입니다. 냉각 속도를 결정합니다. 냉각 속도가 느릴수록(u 값이 클수록) 실행 시간이 길어지는 대신 대체로 더 좋은 결과를 얻을 수 있습니다. 반드시 $u < 1$이어야 합니다.
루프의 반복 횟수는 다음 식으로 주어집니다.
$N = \lceil -\log_{u}{T} \rceil$
$T$와 $u$를 고르는 팁: 국소 최솟값이 많고 상태 공간이 넓다면 냉각 속도가 느려지도록 $u = 0.999$로 설정하십시오. 그러면 알고리즘이 더 많은 가능성을 탐색할 수 있습니다. 반면 상태 공간이 더 좁다면 $u = 0.99$로도 충분합니다. 확신이 서지 않는다면 안전하게 $u = 0.998$ 이상으로 설정하십시오. 알고리즘 한 번의 반복에 드는 시간 복잡도를 계산하고, 이를 이용해 TLE를 피할 수 있는 $N$ 값을 어림한 다음, 아래 공식을 사용하여 $T$를 구하십시오.
$T = u^{-N}$
TSP 구현 예시¶
class state {
public:
vector<pair<int, int>> points;
std::mt19937 mt{ static_cast<std::mt19937::result_type>(
std::chrono::steady_clock::now().time_since_epoch().count()
) };
state() {
points = {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}} ;
}
state next() {
state s_next;
s_next.points = points;
uniform_int_distribution<> choose(0, points.size()-1);
int a = choose(mt);
int b = choose(mt);
s_next.points[a].swap(s_next.points[b]);
return s_next;
}
double euclidean(pair<int, int> a, pair<int, int> b) {
return hypot(a.first - b.first, a.second - b.second);
}
double E() {
double dist = 0;
int n = points.size();
for (int i = 0;i < n; i++)
dist += euclidean(points[i], points[(i+1)%n]);
return dist;
};
};
int main() {
pair<double, state> res;
res = simAnneal();
double E_best = res.first;
state best = res.second;
cout << "Length of shortest path found : " << E_best << "\n";
cout << "Order of points in shortest path : \n";
for(auto x: best.points) {
cout << x.first << " " << x.second << "\n";
}
}
알고리즘의 추가 변형:¶
- TLE를 방지하기 위해 while 루프에 시간 기반 종료 조건을 추가하십시오
- 위에서 구현한 감쇠는 지수 감쇠입니다. 필요에 따라 언제든지 다른 감쇠 함수로 바꿀 수 있습니다.
- 위에서 제시한 확률 수용 함수는 지수의 분자에 있는 $E_{next} - E$ 항 때문에 에너지가 더 낮은 상태를 수용하는 쪽을 선호합니다. 이 항을 제거하기만 하면 PAF가 에너지 차이와 무관해집니다.
- 에너지 차이 $E_{next} - E$가 PAF에 미치는 영향은 아래와 같이 지수의 밑을 키우거나 줄여서 증가시키거나 감소시킬 수 있습니다:
bool P(double E, double E_next, double T, mt19937 rng) { double e = 2; // set e to any real number greater than 1 double prob = pow(e,-(E_next-E)/T); if (prob > 1) return true; else { bernoulli_distribution d(prob); return d(rng); } }