오랜만에 글을 써 봅니다. 근래에 학업 공부를 하고 그러다 보니까 블로그에 공부한 걸 잘 안 올렸네요.(나름 자료를 찾아보고 틈틈히 보긴 했는데 쩝)


  2018년 겨울방학에 들어서고, 꽤 전부터 생각해오던 CUDA를 이용한 프렉탈 그리기를 직접 프로그래밍했습니다. 중학교 때 프렉탈을 처음 알게 되고 나서, 그 신비함과 아름다움에 빠져서 언젠가는 꼭 프렉탈을 그려내는 프로그래밍을 하겠다고 다짐을 했는데(그 옛날에 julia set의 상수마다 이미지를 만들어서 1TB 하드를 꽉 채우겠다는 생각도 해본 거 같습니다. 뭐 정작 여러 가지 그려보다 보니 이쁜 모양이 나오는 julia set 상수는 그렇게 많진 않았지만요)겨울방학 들어서고 하게 되었네요.


  요즘엔 안드로이드 프로그래밍을 거의 안 하고(앱 유지보수 정도...) 딥 러닝 등을 Python으로 하고 있습니다. 그래서 이번 프렉탈 그리기는 제가 옛날부터 주로 쓰던 Java AWT라던가 그런 걸로 만들지 않았고, Python 3에다가 PyQt5와 PyCUDA를 활용해서 만들어 봤습니다. 아래 이미지들이 제가 지금 제작한(그리고 제작하고 있는) 앱의 스크린샷입니다.

 

(GitHub 링크: https://github.com/AMakeApp/FractalViewer)


  이러한 프렉탈을 그리는데는 사실 CPU를 써도 크게 상관은 없습니다. 단지 GPU를 쓰는 이유는 프렉탈을 그리는데 들어가는 알고리즘이 상당한 연산을 필요로 하기 때문입니다. 이 연산 과정을 알려면 위의 프렉탈들이 어떻게 만들어지는지 알아볼 필요가 있습니다.


  제가 현재 앱에 구현해 둔 프렉탈은 위 두 이미지인 Mandelbrot Set과 Julia Set입니다. 사실 이 둘은 모양이 상이함에도 상당한 상관관계가 있습니다. 왜냐하면 이 둘의 집합을 계산하는 과정이 상당히 유사하기 때문이죠. 그렇기에 좀 더 많이 알려져 있는 Mandelbrot Set에 대해 알아봅시다. Mandelbrot Set은 이름에서 보이듯이 Mandelbrot가 고안한 집합입니다. 사실 많이 알려져 있듯이 저렇게 알록달록 해야지 Mandelbrot Set은 아니고 수학적으로는 아래 이미지처럼 복소평면 상에서 검정색 점들의 집합이 Mandelbrot Set입니다.

Mandelset hires


이 집합은 일반적으로 우리가 집합을 정의할 때와는 조금 다르게 특이하게 정의되어 있는데, 다음 점화식에서 수열이 발산하지 않게 하는 복소수 $c$의 집합으로 정의되어 있습니다.

$$z_0 = 0, z_{n + 1} = z_n^2 + c$$


즉, 이 집합은 무한 집합일테고 복소평면 상의 어느 한 점에 대해 이 점이 Mandelbrot Set에 속하는지 계산하려면 각 복소수에 대해 수많은 반복 연산을 통해 수열의 발산 여부를 확인해야 합니다. 다행히도 $| z_n | > 4$라면 발산한다고 볼 수 있기 때문에 일정한 반복 연산 후 $z_n$의 절대값을 비교하면 좀 더 수월하게 계산 할 수 있습니다.(아마 이 부분이 수학적으로 증명되어 있는 거 같습니다. 이 문서를 참고하세요) 하지만 그렇다고 해도 수많은 반복 연산이 필요한 것은 사실입니다. 초중고등학교에서 배우는 것처럼 단순히 함수꼴로 구할 수 있는 것도 아니고, 각 점에 대해 일일히 노가다를 해야하는 것과 마찬가지니까요. 각 점에 대해서는 제가 대충 해본 결과 200번 정도 이상하면 눈으로 보기에는 프렉탈이 충분히 이쁘게 나옵니다. 커다란 프렉탈 이미지를 생성하려고 2000 x 2000의 프렉탈 이미지를 생성한다고 합시다. 대략 이정도면 200 * 2000 * 2000 = 8억 번 정도의 연산이 필요로 합니다. 일반적으로 한 알고리즘에 대해서 1억 번 연산에 대한 소요시간을 1초로 잡을텐데(요즘 CPU는 더 빨라서 더 많이 연산을 할 수는 있습니다) 프렉탈 이미지를 계산하기 위해서는 8초나 걸린다는 겁니다. 요즘 CPU는 중사양 CPU도 4코어 8쓰레드 정도 가질텐데, 멀티 코어 프로그래밍을 한다고 해도 최소 1초는 걸린다는 이야기가 됩니다.


  이렇게 CPU는 프렉탈 이미지를 계산하는 데는 상당히 오래 걸리기 때문에 GPU의 수많은 코어를 사용해서 더 빠르게 계산을 할 수 있습니다. 비유를 하자면 CPU는 고급 인력 4~8명이 일하는 것이라면 GPU는 초보자 수 천 명이 일을 하는 것과 비슷합니다. 프렉탈 연산이 많은 반복 연산이 필요로 하지만 정작 그 연산 자체가 그렇게 복잡하지 않기 때문에(특히 분기 구문 같은 것이 많지 않기 때문에) 이러한 경우에는 GPU에서 빠르게 계산할 수 있습니다. 이러한 이유로 프렉탈을 빠르게 그려내기 위해서 GPU를 활용했습니다.


  GPU로 그려내기로 결정하고 나서, 제 컴퓨터에 장착되어 있는 GTX970에서 GPGPU 프로그래밍을 하기 위해서 Nvidia에서 지원하는 CUDA를 활용했습니다. 여기서 CUDA는 일반적으로 CUDA C를 사용해서 C/C++과 함께 사용되는데, GUI로 이미지를 보여주면서 편하게 프로그래밍하고 싶었기 때문에 Python 3를 기반으로 CUDA를 사용할 수 있게 만들어 둔 PyCUDA를 사용했습니다. 그리고 유려한 GUI를 보여주는 Qt를 Python에 접목시킨 PyQt5를 사용했고요.


다음 글부터는 이 Fractal을 그리는데 어떻게 CUDA 프로그래밍을 했는지 제 시행착오와 함께 CUDA에 대해 이야기해보겠습니다.

 음...딥러닝쪽 게시글을 오래 전부터 안 올리고 있었는데, 사실 좀 귀찮았습니다(ㅋㅋㅋㅋㅋㅋ)

  뭐 정리는 해야할 거 같지만, 데이터만 쌓아두고 여러 모델 구축해서 돌리면서 재미있게 놀았죠. 최근에 Mathematica를 구매해서 내장 딥러닝 모듈을 통해 CIFAR-10을 학습하여 76%인가?(30분 학습했는데, 이정도면 대단한 거겠죠? ㅋㅋㅋ) 나오기도 했고, 방학 때는 소논문으로 날씨 예측을 주제로 삼아 새롭게 Keras를 공부하여 패턴 예측을 시도해보기도 하였습니다. 뭐, 나중에 간단하게 모델들 업로드하면서 글로 정리하기로 하고... 오늘은 제가 학수고대한 NCS가 어제 DHL로 도착해서 간단히 리뷰를 해보려고 합니다!


단촐하지만 위엄있는 USB 스틱!


  NCS(Neural Compute Stick)은 아마 인텔의 자회사?로 추정되는 Movidius에서 개발한 딥러닝 inference용(즉, 학습은 불가능합니다) 컴퓨팅 스틱으로 저전력으로 딥러닝을 할 수 있도록 만들어진 것입니다. 79달러이며 저는 Mouser에서 9만원 정도로 구매하였습니다. 다른 리셀러는 모르겠는데, Mouser는 60달러 이상 무료 배송이고 어차피 79달러에서 할인을 하진 않을테니 그냥 여기서 구매했습니다. 아마 여기가 가장 낫지 않을까 싶네요. 지금은 주문이 많이 밀려서 그런지 꽤 걸리는 거 같고, 제가 주문하고 오기까지 한 1달 정도(배송은 1주일도 안 걸렸을 거예요) 걸렸네요.

  Myriad 2를 기반으로 하며 성능은 대략 100Gflops 정도이며, 전력을 겨우 1W만 먹는다고 합니다(!) 그래서인지 Movidius에서 지금 RPI에서 돌아가도록 SDK 배포판을 만들어서 제공하고 있으며, Movidius에선 RPI3를 예시로 보여주지만 제가 가지고 있는 RPI B(1버전)에서도 잘 돌아가더라고요. 아쉽게도 현재 지원 중인 딥러닝 프레임워크는 Caffe 밖에 없습니다. 예제로 들어있는 AlexNet, GoogLeNet, Gender, SqueezeNet을 돌려보면 각각 283ms, 569ms, 237ms, 304ms 정도 걸립니다(Inference Time 기준). 이미지 크기가 크지만 않다면 대역폭 문제가 크게 작용하지 않아서 그런지, 제 컴퓨터에서 USB 3.0으로 할 때나, RPI에서 USB 2.0(정확히는 2.0보다 안 좋습니다. RPI는 젠더로 USB 대역폭을 제공해서...)으로 할 때나 걸리는 시간은 거의 같습니다. 생각보다 딜레이 시간이 커도 각 딜레이 시간에 맞게 카메라에서 frame을 받아서 inference하면 뭐 그래도 쓸만한 거 같네요(기본 제공 예제에 있습니다). 그리고 생각보다(?) 79달러라는 싼 가격 때문인지 많이 구매해서 쓸 수 있게(?) 여러 개의 USB 스틱을 사용할 수 있도록 지원하는 거 같네요.


  지금 고민 중인 RPI Zero W를 구매하게 된다면, 카메라도 같이 구매하고, Movidius랑 같이 엮어서 Multiple Object Detection을 해보려고 생각하고 있습니다. 잘만 응용하면 시각 장애인에게 음성으로 자신 앞에 어떤 물체가 있는지 알려줄 수 있기 때문에(실제품이 나온다고 들었긴 했는데, 한 번 제 손으로 만들어보고 싶네요 ㅎㅎ) 한 번 만들어볼 생각을 하고 있습니다. 크기도 작고, 전력도 작게 먹는데다가, 무려 RPI를 지원하기 때문에 응용할 곳은 무궁무진할 거 같네요. 나중에 TensorFlow랑 Keras만 지원되면 정말 좋을 거 같습니다 ㅎㅎ

  중3때 삼각비라는 부분에서 sin, cos, tan을 배웁니다. 특수각인 0°, 30°, 45°, 60°, 90°에서의 sin, cos, tan값을 배우죠. 고2에 오면 이과 과정 미적분 II에서 삼각함수 단원에서 좀 더 자세히 배웁니다. 60분법 대신에 호도법을 사용하며 정의역도 실수 전체로 확장하죠. 삼각함수 덧셈, 뺄셈정리를 통해 15°, 75° 등의 각을 삼각함수의 특수각 값을 이용하여 구하죠. 하지만 현실에서 각도는 특수각의 합, 차, 1/2만 존재하지 않습니다. 6°, 34.5° 등의 일반적인 각에서 우리는 삼각함수의 값을 구할 필요가 있을 때가 있죠. 그러나 우리가 배운 내용으로는 이러한 값을 구할 수 없습니다. 그러나 컴퓨터는 이러한 값들을 아주 손쉽게 계산해냅니다. 어떻게 계산해내는 걸까요?


  그 과정을 알아보기 위하여 C언어의 math.h에서 사용하는 삼각함수 라이브러리(glibc)를 뜯어보았습니다. C언어에서는 math.h를 include하여 sin, cos, tan 등의 삼각함수를 사용할 수 있습니다. 이 링크에서 아래 코드의 원본을 보실 수 있습니다.

/* k_sinf.c -- float version of k_sin.c
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */

#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
#endif

#include <math.h>
#include <math_private.h>

static const float 
half =  5.0000000000e-01,/* 0x3f000000 */
S1  = -1.6666667163e-01, /* 0xbe2aaaab */
S2  =  8.3333337680e-03, /* 0x3c088889 */
S3  = -1.9841270114e-04, /* 0xb9500d01 */
S4  =  2.7557314297e-06, /* 0x3638ef1b */
S5  = -2.5050759689e-08, /* 0xb2d72f34 */
S6  =  1.5896910177e-10; /* 0x2f2ec9d3 */

float __kernel_sinf(float x, float y, int iy)
{
    float z,r,v;
    int32_t ix;
    GET_FLOAT_WORD(ix,x);
    ix &= 0x7fffffff;           /* high word of x */
    if(ix<0x32000000)           /* |x| < 2**-27 */
       {if((int)x==0) return x;}        /* generate inexact */
    z   =  x*x;
    v   =  z*x;
    r   =  S2+z*(S3+z*(S4+z*(S5+z*S6)));
    if(iy==0) return x+v*(S1+z*r);
    else      return x-((z*(half*y-v*r)-y)-v*S1);
}

  이 kernel_sin을 사용하여 실제의 sin함수를 구현하는 부분인 s_sinf.c 코드를 보면 $-{\pi \over 4} < x < {\pi \over 4}$ 에서는 간단히 iy = 0으로 정의되어 있을 것을 볼 수 있으니 iy = 0일 때 알고리즘을 살펴봅시다. 크기가 크지 않으니 z, v, r을 계산하여 x + v * (S1 + z * r)을 계산하여 넘기는 것으로 볼 수 있습니다. 이를 식으로 정리해보면 아래와 같습니다.

$$ \sin x = x + S1 x^3 + S2 x^5 + S3 x^7 + S4 x^9 + S5 x^{11} + S6 x^{13} $$


  이제 위에 하드코딩된 값을 대입하여 이 그래프를 그려보겠습니다. Legend에 나와있는 것처럼 f(x)가 위의 함수고, sin(x)는 Mathmatica 기본 내장 함수입니다.


사실상 일치해서 sin(x) 값이 다 덮어씌운 것만 보입니다. 이제 정의역을 f(x)에 정의된 것보다 좀 더 늘려볼까요?

정의역을 $-2\pi < x < 2\pi$ 까지 늘려봤습니다. 이제는 어느정도 차이가 두드러지게 나타나네요! 실제 sin 함수를 구현하는 s_sinf.c 코드를 보시면 $-{\pi \over 4} < x < {\pi \over 4}$ 이외에는 오류처리거나 reduction을 통해 표현하는 것을 볼 수 있습니다. 왜냐면 sin은 $0 \le x \le {\pi \over 4}$ 만 알면 정의역이 모든 실수일 때 정의할 수 있기 때문입니다.


  위와 같이 다항식으로 해석함수(고등학교 과정에서는 무한히 미분 가능한 함수 정도? 자세한 건 맨 아래에 있습니다)를 표현하는 것을 테일러 급수(Taylor series)라고 합니다. 이제 테일러 급수가 어떻게 미분 가능한 함수를 표현하는지 알아봅시다. 일단 미분 가능한 함수에 대한 n차 근사식을 $p_n(x)$ 라고 합시다. 그리고 우리는 한 점 $\left(x_0,\ f(x_0)\right)$ 과 그 점에 대한 미분계수( $f'(x)$ )가 주어졌을 때 그 점에서의 접선을 나타낼 수 있습니다 공식은 아래와 같죠.

$$ p_1(x) = f(x_0) + f'(x_0)(x - x_0)$$


이 공식에서 우리는 아래와 같은 성질을 찾을 수 있습니다.

$$ p_1(x_0) = f(x_0),\ p_1'(x_0) = f'(x_0) $$


따라서 이와 같은 성질을 이용하여 이번엔 2차 다항식에서 근사식을 만들어봅시다. 위의 성질을 사용하면 아래와 같은 조건을 만족하는 $p_2(x)$ 를 구해야겠죠.

$$ p_2(x_0) = f(x_0),\ p_2'(x_0) = f'(x_0),\ p_2''(x_0) = f''(x_0) $$


이에 대해 만족하는 식을 구하면 아래와 같습니다.

$$ p_2(x) = f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 $$


이와 같은 방식을 계속 반복하면 n차 미분한 식을 $f^n(x)$ 라고 했을 때 아래와 같은 n차 근사식을 구할 수 있습니다.

$$ p_n(x) = \sum^n_{k=0}{f^k(x - x_0)^k \over k!} $$


이때 아래와 같은 표현으로 $f(x)$ 를 나타내는데, 이 표현을 테일러 정리라고 하며 $R_n$ 을 나머지 항이라고 부릅니다.

$$ \begin{align*} f(x) &= p_n(x) + R_n \\&= f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 + \cdots + {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1} \end{align*} \\ R_n(x) = {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1}$$


나머지항은 일종의 테일러 정리의 식과 실제 $f(x)$ 의 식 간의 차이(오류)를 보여주므로, $R_n(x)$ 가 $n = \infty$ 일 때 0으로 수렴하면 테일러 정리의 식과 실제 식과 같아지므로 이때 우리는 위 무한급수를 테일러 급수라고 하며, 특별히 $x_0 = 0$ 인 경우에는 매클로린 급수라고 합니다.


이러한 방식으로 테일러 급수를 유도해낼 수 있는데, 수학적으로 엄밀하게 증명해보려면 미분법과 적분법을 사용하는 방법이 있습니다. 굳이 쓸 필요는 없어 보여서 이 링크를 참고하시기 바랍니다. 검색해도 잘 나와요.


자, 이제 대략적으로 테일러 급수에 대한 내용을 정리해보았습니다. 이제 그러면 적용을 해봅시다. 위에서 예시를 든 것처럼 $\sin x$ 에 대해서 테일러 급수를 써봅시다. 아래 그래프는 $\sin x$ 의 그래프(검정)와, $x_0 = 0$ 에서 테일러 급수의 1차(빨강), 3차(주황), 5차(노랑), 7차(초록), 9차(파랑), 11차(보라) 항까지 합한 것의 그래프입니다.

보이는 것처럼 합하는 항수가 커질수록 실제 $\sin x$ 에 다가가는 것을 볼 수 있습니다.


이제 다음으로 파데 근사(Padé approximant)에 대해 써보려고 합니다. 어떤 함수를 다항식으로 표현한 것이 테일러 급수라면, 어떤 함수를 유리식으로 표현한 것이 파데 근사입니다. 테일러 급수의 일반화이며 상위호환이라고 하는 거 같네요. 음이 아닌 정수 $m,\ n$에 대하여 파데 근사의 꼴은 아래와 같습니다. 이때 $n = 0$ 이 되면 매클로린 급수와 같게 됩니다.

$$ R(x)_{m/n} = {\sum_{j=0}^m{a_jx^j} \over 1 + \sum_{i=1}^n{b_kx^k}} $$


...파데 근사 계산을 한 번 해서 정리해보고 싶었는데 내용 이해가 잘 안 되는 관계로(쿨럭) 링크는 남겨둘테니 궁금하신 분들은 여기서 살펴보세요. 계산법을 제대로 이해하긴 힘들거 같고, 대신에 Mathematica에서 PadeApproximant를 함수로 제공하기 때문에 이것으로 $x_0=0$ 일 때 계산한 $\sin x$ 를 아래 표로 정리해보았습니다.

 m\n

1(빨강)

3(주황)

5(초록)

7(파랑)

 0(0)

$x$

$x-\frac{x^3}{6}$

 $x-\frac{x^3}{6}+\frac{x^5}{120}$

$x-\frac{x^3}{6}+\frac{x^5}{120}-\frac{x^7}{5040}$

 2(0.3)

 $\dfrac{x}{1 + \frac{x^2}{6}}$

 $\dfrac{x - \frac{7x^3}{60}}{1 + \frac{x^2}{20}}$

$\dfrac{x - \frac{x^3}{7} + \frac{11x^5}{2520}}{1 + \frac{x^2}{42}}$

 $\dfrac{x - \frac{11x^3}{72} + \frac{13x^5}{2160} - \frac{x^7}{12096}}{1 + \frac{x^2}{72}}$

 4(0.7)

  $\dfrac{x}{1 + \frac{x^2}{6} + \frac{7x^4}{360}}$

$\dfrac{x - \frac{31x^3}{294}}{1 + \frac{3x^2}{49}+\frac{11x^4}{5880}}$

 $\dfrac{x - \frac{53x^3}{396} + \frac{551x^5}{166320}}{1 + \frac{13x^2}{396} + \frac{5x^4}{11088}}$

 $\dfrac{x - \frac{241x^3}{1650} + \frac{601x^5}{118800} - \frac{121x^7}{2268000}}{1 + \frac{17x^2}{825} + \frac{19x^4}{118800}}$


그래프로 나타내면 아래와 같습니다. 조금 굵은 검정 선은 $\sin x$ 그래프이고, n 옆에 써 있는 색상에 m 옆에 써 있는 어두움을 가진 선이 그 (m, n)의 파데 근사 그래프입니다.

(무슨 그래픽 기교하는 줄)

이렇게 컴퓨터는 함수를 근사하여 계산합니다. 나중에 기회가 되면 파데 근사를 좀 더 조사해서 올려보도록 하겠습니다.


p.s 참고 자료: Introduction to Numerical Analysis - Alastair Wood, 위키백과-테일러 급수, 네이버 지식백과, 네이버캐스트, 위키백과-파데 근사, WolframMathWorld-Padé Approximant

+ Recent posts