본문 바로가기

Problem Solving/알고리즘

PS: 뉴비를 위한 Simulated Annealing 입문 (2)

INTRO

국내에서 가장 널리 사용되는 온라인 저지인 BOJ에는 Simulated Annealing으로 풀 수 있는 유명한 문제가 하나 있습니다. 바로 동전 뒤집기 2 (https://www.acmicpc.net/problem/2582) 입니다. 많은 사람들이 KOI에도 시뮬레이티드 어닐링 같은 특수한 알고리즘을 시용해서 푸는 문제가 나오냐고 의아해 하는데, 저 당시에는 서브태스크가 아닌 테스트 케이스 당 부분점수가 있었고, 따라서 정해가 적절한 가지치기를 동반한 백트래킹이었다고 추측해 볼 수 있을 것 같습니다. 하지만, 고생고생해서 백트래킹에 휴리스틱 함수를 적용시키는 것 보다는 시뮬레이티드 어닐링을 배워서 써먹는 쪽이 아무래도 훨씬 빠르겠죠. 시뮬레이티드 어닐링은 이러한 종류에 문제에 더 일반적으로 적용시킬 수 있기도 하고 말입니다. 그런 고로, 이 문제를 시뮬레이티드 어닐링으로 한 번 풀어 보려 합니다.

가볍게 시작해 봅시다

문제를 요약하면 가로가 n, 세로가 n개인 정사각형 형태의 동전들을 한 줄 단위로 적절히 뒤집어서, 뒷면이 나오는 동전의 수를 최소화 하는 것이라고 할 수 있습니다. n<=32 조건을 보니 누가 봐도 다항 시간에 동작하는 효율적인 해결책은 없을 것 같고, dp를 하자니 2^32 * (모종의 수)라는 시간 복잡도가 나올 것 같은데 역시 상상하기도 싫습니다. 시뮬레이티드 어닐링을 적용한다고 생각해 보면, 일단 상태를 어떻게 정의하고 변형시킬 지를 결정해야 합니다. 가장 먼저 해야 하는 관찰은 한 줄을 두 번 이상 뒤집을 필요는 없다는 것이죠. 두 번 뒤집으면 원래대로 돌아갈 테니 이는 자명합니다. 그러면 모든 세로줄과 모든 가로줄에 대해서 뒤집거나/뒤집지 않거나를 상태로 두고, 이에 대해 시뮬레이티드 어닐링을 적용할 수 있습니다. 상태가 정의 되었으니 이제 인접 상태를 정의할 순서입니다. 직관적으로, 임의의 줄을 하나 무작위로 골라서 뒤집힌 여부를 토글하는 것으로 인접 상태를 나타낼 수 있을 겁니다. 그리고 남은 것은 평가 함수를 결정하는 것이죠. 역시 뒤로 뒤집힌 동전의 개수를 단순히 세서 평가 함수로 사용할 수 있을 겁니다. 그러면 남은 것은 열심히 코딩을 하고... 볼츠만 상수를 정하고... 버그를 잡고.... 만고의 노력 끝에 내면 WA가 나옵니다. 뭐 상수를 잘 조절하면 될 수도 있을 것 같은데, 그러기 보다는 좀 더 효율적인 방법이 있는지 살피는 게 맞지 않을까요?

약간의 최적화

 조금의 관찰을 더하면, 모든 가로줄과 세로줄에 대해서 뒤집어 볼 필요가 전혀 없습니다! 가로줄의 뒤집힌 여부가 모두 결정되었다고 합시다. 그러면 세로줄을 뒤집는 시행에 대해서는, 서로 겹치는 부분이 없으므로 뒤를 향한 동전의 수 변화가 독립이 됩니다. 따라서 가로줄이 모두 정해진 순간 세로줄은 greedy하게, 더 많은 동전이 앞을 향하도록 뒤집어 줄 수 있고 이 때가 최적임이 자명합니다. 따라서 기존에 크기가 64인 상태를 고려했다면, 이제는 크기가 32인 상태를 고려할 수 있습니다. 이제 탐색해야 하는 공간의 크기가 2^32배 줄었으니, 무난하게 제출하고 맞았습니다!!를 받을 수 있습니다.

뉴비를 위한 글이 아닌 것 같은데요

 여기까지 본인이 올린 포스트의 설명을 듣고 별도의 추가적인 공부 없이 동전 뒤집기에서 AC를 받으셨다면, 뉴비가 아니므로 뒤로 가기를 누르시면 됩니다. 랜덤이 뭔지 잘 모르는 뉴비들에게는 더 친절한 설명이 필요하다는 의견을 반영하여, 코드를 보면서 천천히 설명을 하기로 했습니다. 다음은 제가 해결한 동전 뒤집기 2의 소스 코드입니다.

#include <bits/stdc++.h>
#include <random>
#define INF 987654321

using namespace std;

int scoring() {
	//TODO: score present state and return
}

double t=1, d=, k=, lim=; //TODO: set constants and seed 
int ret=INF;
std::mt19937_64 seed();
std::uniform_real_distribution<double> rng(0, 1);

void simulated_annealing() {
	double e1, e2;
	while (t > lim) {
		e1 = scoring();
     	//TODO: change state 1 to 2
		e2 = scoring();

		double p = exp((e1 - e2) / (k*t));
		if (p < rng(seed)) //TODO: return back to state 1 
		t *= d;

		ret = min(ret, scoring());
	}
}

int main() {
	scanf("%d", &n);
	//TODO: implement it
	simulated_annealing();
	printf("%d", ret);
	return 0;
}

물론 전 별로 친절한 사람이 아니기 때문에 가운데 내용을 모두 빼고 implementation is left as an exercise to the reader을 시전해 보았습니다. 안타깝게도, 문제에서 요구하는 내용을 모두 빼고 보니 시뮬레이티드 어닐링의 핵심적인 골조가 더 잘 보이게 되었습니다. 딱 봐도 중요해 보이는 simulated_anneling() 함수부터 분석해 봅시다.

void simulated_annealing() {
	double e1, e2;
	while (t > lim) {
		e1 = scoring();
     	//TODO: change state 1 to 2
		e2 = scoring();

		double p = exp((e1 - e2) / (k*t));
		if (p < rng(seed)) //TODO: return back to state 1 
		t *= d;

		ret = min(ret, scoring());
	}
}

 

 일단 t가 lim이라는 변수보다 클 때까지 계속해서 아래 과정을 반복하는군요. 느낌 상 t가 현재 온도, lim이 임계 온도라는 것을 알 수 있습니다. 반복문 내부에서는 e1과 e2에 대해서 각각 scoring() 함수를 호출하는 것을 볼 수 있습니다. 당연히도 scoring 함수는 평가 함수라는 것을 추측해 볼 수 있겠군요. 따라서 e1은 현재 상태의 평가 함수 값e2는 인접 상태의 평가 함수 값임을 알 수 있겠습니다. p에는 우리가 처음 시뮬레이티드 어닐링을 배울 때 보던 바로 그 식이 적혀있습니다. 따라서 p는 인접 상태로 전이할 확률이 되겠군요. 그 후 현재 온도에 d를 곱하는데, 아마도 d가 온도 감률인 것 같습니다. 이 과정을 한번 수행할 때마다 ret이라는 변수에 현재 해의 출력을 저장해서 그 중 최솟값을 출력하네요. 간단합니다. 여기서 rng(seed)라는 함수를 볼 수 있는데, PS만 하신 분들은 이게 뭔지 잘 모르실 겁니다.
 랜덤한 수를 출력하는 방법은 C++ 기준으로 보통 두 가지가 있습니다. 하나는 C의 srand() 함수와 rand() 함수를 사용하는 방법이 있고, 나머지 하나는 random 헤더에 정의된 메르센 트위스터를 사용하는 방법입니다. rand() 함수가 더 직관적이고 사용하기 쉬움에도 불구하고, 시스템에 따라서 발생하는 수의 범위가 달라질 수도 있고 그리 좋은 랜덤 함수도 아닌 데다가 정수만 생성할 수 있기 때문에 좋은 난수를 발생시키려면 메르센 트위스터를 사용해야 합니다. 그러나, 귀찮기 때문에 작은 범위 내에서 난수를 생성해야 할 때는 rand()함수를 이용해도 별 상관이 없습니다. 시뮬레이티드 어닐링에서 메르센 트위스터로 난수를 생성할 때는 mt19937을 이용할 수 있습니다. 자세한 내용은 구글 갓한테 물어보면 알 수 있으니, 여기서는 실전압축해서 사용법을 소개하겠습니다. 

1. std::mt19937이나, std::mt19937_64(64비트 변수용) 변수 x를 선언한다. 이때 생성자 변수로는 시드를 전달한다. 시드에는 아무 숫자를 임의로 대입할 수도 있고, 정밀 시간초를 전달해도 되지만 별 차이가 없으므로 코드포스처럼 핵을 당할 위험이 있는 게 아니라면 본인의 마음에 드는 숫자를 넣는 게 좋다.

2. std::uniform_int_distribution<정수 자료형> 이나 std::uniform_real_distribution<실수 자료형> 변수 y를 선언한다. 이때 생성자 변수로는 생성될 난수의 (최솟값, 최댓값)을 전달한다.

3.  y(x)가 난수를 반환한다. 이때 x가 y의 시드 역할을 한다.

참 쉽죠?

std::mt19937_64 seed(1119);
std::uniform_real_distribution<double> rng(0, 1);
//중략
if (p < rng(seed))

 아, 이제 seed라는 변수의 괄호 값을 채워 줄 수 있습니다. 제 생일이 11월 19일이므로 1119를 대입했고, 이를 통해 함수 중 rng(seed)가 0부터 1까지의 임의의 double형 난수를 반환한다는 것을 알 수 있겠네요. 따라서 if문 내부는 1-p의 확률로 동작함을 알 수 있겠습니다.

 개괄적인 구조를 이해했으니 이제 내부를 좀 채워 주고 싶습니다. 일단 우리가 문제의 상태를 어떻게 나타낼 것인지 고민해 봅시다. 글 처음에는 각 줄이 뒤집힌 상태로 정의한다고 했지만, 어차피 평가 함수 구현할 것을 생각하면 줄이 뒤집힌 상태를 저장하는 것과 그냥 현재 동전들의 상태를 저장하는 것이 크게 다를 바 없어 보입니다. 그래서 그냥 grid[][]라는 2차원 배열로 현재 동전이 앞인지, 뒤인지 표시하겠습니다. 아, 그러면 평가 함수는 바로 만들 수 있겠군요? 뒷면으로 표시된 동전을 1이라 하고 앞면을 0이라고 합시다. 그러면 scoring() 함수를 바로 채워 줄 수 있습니다.

int scoring() {
	int s = 0;
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			s += grid[i][j];
	return s;
}

 좋습니다. 평가 함수가 완성되었습니다. 다음으로 채워야 할 부분은 상태 전이입니다. TODO: change state 1 to state 2라고 되어 있는 부분입니다. 현재 상태인 grid 배열에서 임의의 한 줄을 뒤집어 주면 되겠습니다. 원래 상태를 따로 저장해 두고, rand() 함수를 이용해서 임의의 한 줄을 잡은 후 모두 뒤집어 줍니다. 뒤집고 나서는, 이 뒤집고 나서의 상태가 나머지 방향의 줄에서도 최적이 될 수 있도록 따로 greedy하게 처리해 주어야 합니다. 이를 위해서 함수들을 추가로 조금 정의해 두었습니다. 

void turn(int x, int y) {
	if (y == 0) {
		for (int i = 0; i < n; i++)
			grid[x][i] = 1 - grid[x][i];
	}
	else {
		for (int i = 0; i < n; i++)
			grid[i][x] = 1 - grid[i][x];
	}
}

void func() {
	for (int i = 0; i < n; i++) {
		int s = 0;
		for (int j = 0; j < n; j++)
			s += grid[j][i];
		if (s > (n / 2)) turn(i, 1);
	}
}

void simulated_annealing() {
	double e1, e2;
	int ori[40][40];
	while (t > lim) {
		e1 = scoring();
		for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) ori[i][j] = grid[i][j];
		int pos = rand() % n;
		turn(pos, 0); func();
		e2 = scoring();

		double p = exp((e1 - e2) / (k*t));
		if (p < rng(seed)) //TODO!
		t *= d;

		ret = min(ret, scoring());
}

 이제 코드가 꽤 채워졌으니 시뮬레이티드 어닐링의 작동 과정을 어느 정도 분석할 수 있습니다. 온도가 임계 온도보다 높은 동안 현재 상태의 점수 e1과 나중 상태의 e2를 비교한 다음, 온도를 감률에 따라 감소시키게 됩니다. 그리고 현재 답이 현재 답안보다 좋으면 갱신해 주네요. 딱 한 곳 빈 부분은 e1와 e2에 따라 실제로 상태를 바꿀 것인지 확인하는 부분입니다. 원래는 p의 확률로 새로운 상태로 바꿔 주어야 하지만 이미 지금 상태를 바꾸었으니 1-p의 확률로 원래 상태로 되돌려 주어도 동일한 작동 과정을 거칩니다. 따라서 다음과 같이 채워줄 수 있습니다.

if (p < rng(seed)) for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) grid[i][j] = ori[i][j];

또한 처음에 입력을 받는 부분을 구현해 주었습니다. 이 문제에서는 초기 상태를 제공해 주기 때문에 그대로 쓰면 되지만, 초기 상태를 제공하지 않는 다른 문제에서는 랜덤으로/그리디하게 찾은 근사해로부터 시작할 수 있습니다.

int main() {
	scanf("%d", &n);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			char ch;
			scanf(" %c", &ch);
			grid[i][j] = (ch == 'H' ? 0 : 1);
		}
	}
	simulated_annealing();
	printf("%d", ret);
	return 0;
}

 마지막으로는 상수를 결정해 주어야 합니다. 설명했듯이 일반적으로 온도 감률은 0.9999에서 0.9995 정도로 잡는 것이 일반적입니다.(문제에 따라 달라질 수 있음) 또한 임계 온도를 결정해 주어야 하는데, 사실 PS에서는 반복문의 횟수를 결정해 주는 것이 수행 시간을 예측하기에 더 적절하겠죠. 하지만 여기에서는 그냥 이런 식으로 했습니다(...) 마지막으로 볼츠만 상수를 정해야 하는데, 이 문제에서는 저는 10으로 설정해서 모든 데이터에 대해 올바른 정답을 출력할 수 있었습니다. 이 볼츠만 상수를 정하는 방법은 정말 많은 방법이 있는 것으로 알고 있는데 일반적인 문제에 대해 꽤 괜찮은 볼츠만 상수를 찾는 방법은 저도 모르고, 어림짐작하여 꽤 괜찮을 것 같은 값을 넣고 확인하는 것을 반복하고 있습니다. 혹시 노하우가 있으신 분이 있다면 제보 부탁드립니다.

 이상과 같은 내용을 종합하면 아래와 같은 코드가 나옵니다. 200ms 정도의 시간 안에 여유롭게 모든 데이터를 통과하는 것을 확인할 수 있습니다. 이 문제에 대해서는 비트 연산 등의 추가적인 최적화를 적용할 수 있지만, 코드를 굳이 읽기 힘들게 만들 필요가 없다고 생각하여 넣지는 않았습니다.

#include <bits/stdc++.h>
#include <random>

using namespace std;
int n;
int grid[40][40];

void turn(int x, int y) { 
	if (y == 0) {
		for (int i = 0; i < n; i++)
			grid[x][i] = 1 - grid[x][i];
	}
	else {
		for (int i = 0; i < n; i++)
			grid[i][x] = 1 - grid[i][x];
	}
}

void func() {
	for (int i = 0; i < n; i++) {
		int s = 0;
		for (int j = 0; j < n; j++)
			s += grid[j][i];
		if (s > (n / 2)) turn(i, 1);
	}
}

int scoring() {
	int s = 0;
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			s += grid[i][j];
	return s;
}

double t=1, d=0.9999, k=10, lim=0.0005;
std::mt19937_64 seed(1919);
std::uniform_real_distribution<double> rng(0, 1);

int ret = 999;

void simulated_annealing() {
	double e1, e2;
	int ori[40][40];
	while (t > lim) {
		e1 = scoring();
		for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) ori[i][j] = grid[i][j];
		int pos = rand() % n;
		turn(pos, 0);
		func();
		e2 = scoring();

		double p = exp((e1 - e2) / (k*t));
		if (p < rng(seed)) for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) grid[i][j] = ori[i][j];
		t *= d;

		ret = min(ret, scoring());
	}
}

int main() {
	scanf("%d", &n);
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			char ch;
			scanf(" %c", &ch);
			grid[i][j] = (ch == 'H' ? 0 : 1);
		}
	}
	simulated_annealing();
	printf("%d", ret);
	return 0;
}

 이와 같이 시뮬레이티드 어닐링으로 문제를 해결해 보았습니다. 개인적으로 시뮬레이티드 어닐링의 강점은 평가 함수와 상태를 적절히 정의할 수만 있다면 어떤 형태의 최적화 문제에 대해서도 괜찮은 해를 제공해 주는 점이라고 생각합니다. 이 포스팅에서 제공한 스켈레톤 코드와, 문제를 해결하는 과정을 적절히 자기 방식대로 활용하여 휴리스틱 문제가 나와도 당황하지 않고 침착하게 해결해 낼 수 있는 뉴비가 되기를 기원합니다.

 

P.S. 사실 이 글을 쓰게 된 동기는 NYPC입니다. 포스팅 작성 시점에서 NYPC 본선이 얼마 남지 않은 것으로 알고 있는데, NYPC 예선에서 휴리스틱 문제가 많이 출제된 것으로 보아 본선에서도 휴리스틱 문제가 한 문제 나올 것 같다고 예상하고 있습니다. (대부분의 휴리스틱 문제를 시뮬레이티드 어닐링으로 해결할 수 있었습니다 :blobaww: ) 대부분의 문제는 적절한 결정론적 휴리스틱과 그리디를 섞어 고득점을 받을 수 있겠지만, 정 힘들다면 시뮬레이티드 어닐링은 좋은 대안이 될 수 있습니다. 잘 연습해서 좋은 결과 있기를 기원합니다.