[C++] EigenRand: Eigen용 Random Library 개발

Eigen는 Random 지원이 빈약하다

최근 c++로 tomotopy라는 토픽모델링 툴을 개발하면서 벡터화 가속을 위해서 Eigen이라는 라이브러리를 가져다 썼습니다. Eigen은 여러 곳에서 널리 사용되는 선형대수 연산용 C++ 라이브러리로, 사실상 이쪽 업계의 표준 아닌 표준이라고 할 수 있습니다. 오랫동안 검증되고 최적화되었기 때문에 Eigen 라이브러리만 가져다 쓰는 것으로도 충분히 속도 향상을 이룰 수 있었습니다.

다만 여러 확률 분포를 이용하는 토픽 모델링의 특성상 코드 내에서 확률 분포 내에서 임의의 숫자를 샘플링하는 작업을 굉장히 자주 반복해야하는데 불행히도 Eigen에는 랜덤 관련 함수 지원이 크게 부족했습니다. 일례로 현재 3.3.7버전에서 제공하는 Random함수는 다음 한 가지가 전부입니다.

// 100x100짜리 랜덤 행렬을 생성. 값의 범위는 [-1, 1] Eigen::MatrixXf mat = Eigen::MatrixXf::Random(100, 100);

Eigen의 Random는 여러가지 문제점을 가지고 있습니다. 먼저 rand() 및 srand()라는 오래된 c함수를 사용합니다. (이건 아마 Eigen을 c++98과 호환시키려고 그런것으로 보입니다. c++에 random 라이브러리가 추가된게 c++11부터이니 어쩔 수 없던 문제겠습니다). rand() 함수가 가지는 문제에 대해서는 이미 널리 알려졌기 때문에 여기서 자세히 다루지 않겠습니다. 간략하게 정리하자면, 함수가 생성하는 난수의 품질이 좋지 못하고, 전역변수를 사용하기 때문에 멀티스레딩에 있어 불리하며, 벡터화하기 어렵다는 것입니다. 

그리고 Eigen::Random은 실수에 대해서는 오직 [-1, 1] 범위에서, 정수에 대해서는 전체 범위에서의 균등 분포만을 생성한다는 또 다른 한계가 있습니다. 만약 정규분포를 따르는 난수를 생성하고 싶다면? 각 정수가 등장할 확률이 다른 분포에서 난수를 뽑고 싶다면? 방법이 없습니다.

그래서 Eigen에서 NullaryExpr을 이용해 C++11 버전 랜덤 생성기를 사용하는 방법을 별도로 문서화해두기는 했습니다. 다음과 같지요.

std::mt19937_64 urng{ 42 }; // 난수 생성기 // 생성된 난수를 [0, 1) 범위의 랜덤 실수로 바꿔주는 녀석 std::uniform_real_distribution<> uniform_real; // 다음 코드는 NullaryExpr를 이용해 랜덤 행렬을 생성함 Eigen::MatrixXf mat = Eigen::MatrixXf::NullaryExpr(100, 100, [&](){ return uniform_real(urng); });

NullaryExpr는 행렬의 각 요소마다 넘겨진 함수를 호출하여 [0, 1) 범위의 값을 생성하여 대입합니다. 코드가 조금 길어지긴 했지만, 이제 생성되는 난수를 쉽게 통제할 수 있고 전역변수로부터도 자유롭습니다. uniform_real_distribution 대신에 다른 난수 분포를 이용하면 정규분포도, 카이제곱 분포도, 감마 분포도 생성할 수 있습니다. (C++11 표준에 포함된 난수 분포에 대해서는 여기를 참조하세요.)

그런데 또 다른 문제가 있습니다. NullaryExpr는 각 요소를 하나씩 하나씩 처리하기 때문에 Eigen의 큰 장점 중 하나인 벡터화를 통한 속도 향상이 적용되지 않습니다. SSE나 AVX처럼 여러 개의 수를 한꺼번에 처리하는 능력이 있는 CPU에서는 4개, 8개 혹은 16개씩을 동시에 처리함으로써 빠른 연산이 가능한데, NullaryExpr를 사용하면 그게 안된다는 것이지요.

그래서 tomotopy를 개발할때 내부적으로 벡터화 가능한 Eigen용 난수 생성 custom expression을 작성하여 이를 사용했었습니다. 그런데 개발해놓고 보니 tomotopy 코드 속에만 박혀있기는 아깝고, 누구나 쉽게 쓸수 있게하는게 좋을거 같아서 아예 독립시키기로 하였습니다.

EigenRand를 소개합니다

짠! tomotopy에서 사용된 난수 생성 함수들을 긁어모아서 EigenRand라는 이름을 붙이고 라이브러리화시켰습니다.

https://github.com/bab2min/EigenRand

처음에는 난수 생성기랑 균등분포, 정규분포 생성기능 밖에 없었는데, 라이브러리화시키려고 하니깐 욕심이 생겨서 c++11 표준에 들어간 다른 난수 분포들도 하나 둘씩 추가했습니다. 아직 c++11 표준에 있는 모든 분포가 지원되는 건 아니지만, 이제 어느 정도 모양은 갖춰진거 같아서 이렇게 공개합니다~!

EigenRand 역시 Eigen처럼 header-only 라이브러리이기 때문에 그저 헤더 파일을 include하는 것만으로 사용이 가능합니다. 사용법도 간단하지요.

std::mt19937_64 urng{ 42 }; // 난수 생성기 // 표준 정규분포를 따르는 100x100짜리 랜덤 행렬을 생성. Eigen::MatrixXf mat = Eigen::Rand::normal<Eigen::MatrixXf>(100, 100, urng); // 또 numpy와 유사하게 -Like 함수들이 있습니다. // mat과 타입/크기가 같고, 감마 분포를 따르는 랜덤 행렬을 생성합니다. Eigen::MatrixXf mat2 = Eigen::Rand::gammaLike(mat, urng);

참 쉽죠?

C++11의 난수 생성기 외에도 Eigen::Rand::Vmt19937_64라는 벡터화된 난수 생성기도 지원합니다. 이 난수 생성기는 벡터 명령어 지원여부에 따라 64bit 난수를 동시에 2개 혹은 4개까지 생성할 수 있습니다. 현재 지원하는 기능들에 대해서는 여기를 확인하시면 되겠습니다.

속도는 얼마나 빠를까

아무리 기능이 많아도 속도가 느리면 소용이 없겠지요. 여러 시스템에서 속도 평가한 결과를 공유드립니다.

다음은 Intel(R) Xeon(R) Platinum 8171M CPU @ 2.60GHz (Ubuntu 16.04, gcc7.5)에서 테스트한 결과입니다. 표 안의 수치들은 백만개의 난수를 생성할때 소요된 시간을 초 단위로 표현한 것입니다.

	Eigen	C++ std	EigenRand
			(No Vect.)	(SSE2)	(SSSE3)	(AVX)	(AVX2)
balanced	9.0	-	5.9	1.5	1.4	1.3	0.9
balanced(double)	8.7	-	6.4	3.3	2.9	1.7	1.7
chiSquared	-	80.5	249.5	64.6	58.0	29.4	28.8
discrete(int32)	-	-	14.0	2.9	2.6	2.4	1.7
discrete(fp32)	-	-	21.9	4.3	4.0	3.6	3.0
discrete(fp64)	-	72.4	21.4	6.9	6.5	4.9	3.7
exponential	-	31.0	25.3	5.5	5.3	3.3	2.9
extremeValue	-	66.0	60.1	11.9	10.7	6.5	5.8
gamma(0.2, 1)	-	207.8	211.4	54.6	51.2	26.9	27.0
gamma(5, 3)	-	80.9	60.0	14.3	13.3	11.4	8.0
gamma(10.5, 1)	-	81.1	248.6	63.3	58.5	29.2	28.4
lognormal	-	66.3	55.4	12.8	11.8	6.2	6.2
normal(0, 1)	-	38.1	28.5	6.8	6.2	3.8	3.7
normal(2, 3)	-	37.6	29.0	7.3	6.6	4.0	3.9
randBits	-	5.2	5.4	1.4	1.3	1.1	1.0
uniformReal	-	12.9	5.7	1.4	1.2	1.4	0.7
weibull	-	41.0	35.8	17.7	15.5	8.5	8.5

	C++ std	EigenRand
		(No Vect.)	(SSE2)	(SSSE3)	(AVX)	(AVX2)
Mersenne Twister(int32)	4.7	5.6	4.0	3.7	3.5	3.6
Mersenne Twister(int64)	5.4	5.3	4.0	3.9	3.4	2.6

이건 AMD Ryzen 7 3700x CPU @ 3.60GHz (Windows 10, MSVC2017)에서 테스트한 결과입니다.

	Eigen	C++ std	EigenRand
			(SSE2)	(AVX)	(AVX2)
balanced	20.8	-	1.9	2.0	1.4
balanced(double)	21.7	-	4.1	2.7	3.0
chiSquared	-	153.8	33.5	21.2	17.0
discrete(int32)	-	-	2.4	2.3	2.5
discrete(fp32)	-	-	2.6	2.3	3.5
discrete(fp64)	-	55.8	5.1	4.7	4.3
exponential	-	33.4	6.4	2.8	2.2
extremeValue	-	39.4	7.8	4.6	4.0
gamma(0.2, 1)	-	128.8	31.9	18.3	15.8
gamma(5, 3)	-	156.1	9.7	8.0	5.0
gamma(10.5, 1)	-	148.5	33.1	21.1	17.2
lognormal	-	104.0	6.6	4.7	3.5
normal(0, 1)	-	48.8	4.2	3.7	2.3
normal(2, 3)	-	48.8	4.5	3.8	2.4
randBits	-	4.2	1.7	1.5	1.8
uniformReal	-	30.7	1.1	1.0	0.6
weibull	-	46.5	10.6	6.4	5.2

	C++ std	EigenRand
		(SSE2)	(AVX)	(AVX2)
Mersenne Twister(int32)	5.0	3.4	3.4	3.3
Mersenne Twister(int64)	5.1	3.9	3.9	3.3

AVX2 환경에서는 최대 10배까지 빨라지는걸 확인할 수 있습니다. Eigen을 사용중이고 난수 생성을 해야하는 경우라면 굳이 C++11 표준함수를 쓰는것보다는 EigenRand를 쓰는걸 추천드립니다~

과연 확률 분포를 잘 근사할까

특정 확률 분포에서 난수를 추출하기 위해서는 복잡한 실수 연산(exp, log, tan, tgamma 등)이 필요한 경우가 많습니다. 표준 수학 함수에서는 이 함수들을 최대한 정밀하게 계산하기 위해서 노력하기 때문에 속도가 느려지는 경우가 많죠. 반면 SIMD를 이용한 벡터화된 수학 함수에서는 빠른 속도를 위해 정밀도를 다소 포기하는 경향이 있습니다. EigenRand도 난수 추출에 벡터화된 수학 함수를 사용하기 때문에 근사에서 발생하는 오차가 생각보다 클 수 있습니다. 따라서 실제 추출된 결과가 원 확률 분포와 동떨어져 있을 수도 있겠지요. 그래서 과연 난수 생성 결과가 실제 분포를 잘 근사하는지 측정해보았습니다.

특정 분포에서 생성된 난수들이 과연 원 분포를 잘 따르는지 판단하는 것은 어려운 일이기 때문에 둘을 직접 비교하기 보다, 표집된 난수의 누적확률분포와 원 분포의 누적확률분포를 비교하는것으로 대신했습니다. (나중에 알아보니 이 수치가 두 분포 간의 Earth Mover's Distance와 일치한다고 하는군요.) 둘 사이의 오차가 작을 수록 표집된 난수가 원 분포를 잘 따른다고 볼 수 있습니다.

각각의 분포로부터 총 32768개의 샘플을 표집하여 샘플의 분포를 구하고, 이 분포와 실제 확률 분포 간의 EMD를 계산했습니다. 아래 표는 서로 다른 seed값에 대해 50번씩 수행하여 평균 및 표준편차를 구한것입니다. (괄호 안이 표준편차) (balanced는 c++표준에 없기 때문에 Eigen의 Random함수로 대신했습니다.)

	C++ std	EigenRand
balanced*	.0034(.0015)	.0034(.0015)
chiSquared(7)	.0260(.0091)	.0242(.0079)
exponential(1)	.0065(.0025)	.0072(.0022)
extremeValue(1, 1)	.0097(.0029)	.0088(.0025)
gamma(0.2, 1)	.0380(.0021)	.0377(.0025)
gamma(1, 1)	.0070(.0020)	.0065(.0023)
gamma(5, 1)	.0169(.0065)	.0170(.0051)
lognormal(0, 1)	.0072(.0029)	.0067(.0022)
normal(0, 1)	.0070(.0024)	.0073(.0020)
uniformReal	.0018(.0008)	.0017(.0007)
weibull(2, 1)	.0032(.0013)	.0031(.0010)

값이 작을수록 샘플의 분포가 원 확률분포를 잘 따른다는 뜻인데, 보시면 결과가 거의 엎치락 뒤치락합니다. 표준편차를 보면 둘 간의 차이가 오차범위 내라는 것을 알 수 있습니다. 따라서 실용적으로 쓰기에는 정밀도 문제가 거의 없다고 봐도 되겠습니다.

마치며

사실 이렇게 다양한 종류의 확률 분포가 있는지 몰랐었습니다. 개발하면서 이것저것 찾아보니 참 다양한 확률분포가 여러 곳에서 쓰이고 있더라구요. 해당 분포들을 계산할 때 유용하게 잘 쓰이면 좋겠습니다.