Cubic Spline

Keywords: Cubic Spline

Cubic spline은 복수개의 샘플을 통과하는 부드러운 곡선을 만들되 각 곡선이 3차함수로 이루어지도록하는 보간법의 일종이다. 별다른 논쟁의 여지 없는 심플한 설명이고, 필자 역시 이 정도로만 알아두고 적절한 library들을 통해 사용해오고는 했다. 그런데, 얼마전 우연히 별도의 library를 통하지 않고 cubic spline을 구하는 코드를 고칠 일이 생겼다. 별 것 아니겠지 하고 시작했는데, 의외로 재미있고 생각할 것이 꽤 있다 싶었기에 포스팅한다.

Goal

아래 그림에서 적색과 같이 총 N개의 샘플 \((x_{i}, y_{i})\) 이 주어졌다고 해보자. Cubic spline의 목표는 3차 함수를 사용해서 주어진 N개의 점을 모두 통과하는 부드러운 곡선을 구하는 것이다. 예를들어, 아래 그림의 청색 라인과 같이 말이다. “부드러운”이라는 표현이 다분히 비과학적인데, 모든 점에서 2번 미분이 가능하다는 의미 정도로 기억해두자.

많은 사람들이 오해하는 사실 중 하나는 cubic spline이 마치 \(f: x \rightarrow y\) 를 직접 구해준다고 생각하는 것이다. 그런데, 조금만 생각해보아도 이것은 수학적으로 불가능한 것을 알 수 있다. 예를 들어, 당장 바로 위에서 언급한 그림만 보아도 적색 곡선은 양함수 형태로 나타낼 수 없지 않은가? 그래서, cubic spline은 주로 호장 재매개화를 통해서 \((s_{i}, x_{i}), \,(s_{i}, y_{i})\) 으로 샘플을 찢고(?) 각각에 대해 spline을 활용하여 \(f_{x}: s \rightarrow x\), \(f_{y}: s \rightarrow y\)을 만드는 방법을 사용한다.

Algorithm

이제 \((x_{i}, y_{i})_{i=1...N}\)이 주어졌다고 가정하고 cubic spline을 어떻게 만들 수 있는지 차근차근 알아보도록하자.

Step 1: x, y 각각에 대한 cubic spline 문제로 분할하기

\((x_{i}, y_{i})_{i=1...N}\) 을 \(f: x \rightarrow y\) 로 모델링해버리면 양함수 형태를 구할 수 없을 수도 있으므로, 호장 재매개화를 통해 \(x,\,y\) 각각에 대한 cubic spline문제로 변환한다. 다시말해, \((x_{i}, y_{i})_{i=1...N}\) 을 \((s_{i}, x_{i})_{i=1...N}, \,(s_{i}, y_{i})_{i=1...N}\) 로 찢어두라는 이야기이다. 이제 두 샘플군에 대해 각각 cubic spline \(f_{x}: s \rightarrow x\), \(f_{y}: s \rightarrow y\) 을 구할 것이고, 마지막에 \((f_x(s_i), f_y(s_i))_{i=1...N}\)과 같이 복원하도록 할 것이다. \(f_{x}: s \rightarrow x\), \(f_{y}: s \rightarrow y\) 을 유도하는 과정은 서로 동일하므로 이후의 기술에서는 \(f_{x}: s \rightarrow x\) 를 구하는 방법에 대해서만 이야기하도록 하겠다.

Step 2: N - 1 개의 3차 함수를 설정하기

\((s_{i}, x_{i})_{i=1...N}\) 에는 \([s_{i}, s_{i+1})_{i=1...N-1}\)과 같이 N-1개의 구간이 존재할 것이다. i번째 구간 \([s_{i}, s_{i+1})\)이 3차 함수 \(f_{x}^{i} = a_{i}s^3 + b_{i}s^2 + c_{i}s + d_{i}\)를 통해 이어진다고 가정하자. 우리가 할일은 이제 모든 \(a_i\), \(b_i\), \(c_i\), \(d_i\)들을 알아내면 되는 것이다. 이들을 알아내기 위해서 잠시 초등학생으로 돌아가보자. 알아내야하는 값은 총 \(4 * (N - 1)\)개이다. 따라서, 우리는 \(4 * (N - 1)\)개의 조건식을 찾아주면 된다 (미지수가 N개일 때, 식이 N개 있으면 풀린다라고 선생님이 말씀하시지 않았던가). 이제 아래에 이어질 섹션에서 이 \(4 * (N - 1)\)개의 조건(식)을 찾아주도록 하겠다.

Step 3: 통과조건에 대한 식을 만들기

당연하게도 구하는 cubic spline은 주어진 모든 점을 통과해야하고 각 샘플점에서 연속이어야 한다. 이것으로부터 우리는 아래와 같이 \(2 * (N - 1)\) 개의 조건을 얻을 수 있다.

\[\begin{align*} f_{x}^{1}(s_{1}) & = a_{1} s_{1}^3 + b_{1} s_{1}^2 + c_{1} s_{1} + d_{1} = x_{1} \\ f_{x}^{1}(s_{2}) & = a_{1} s_{2}^3 + b_{1} s_{2}^2 + c_{1} s_{2} + d_{1} = x_{2} \\ f_{x}^{2}(s_{2}) & = a_{2} s_{2}^3 + b_{2} s_{2}^2 + c_{2} s_{2} + d_{2} = x_{2} \\ f_{x}^{2}(s_{3}) & = a_{2} s_{3}^3 + b_{2} s_{3}^2 + c_{2} s_{3} + d_{2} = x_{3} \\ ... \\ f_{x}^{N-1}(s_{N-1}) & = a_{N-1} s_{N-1}^3 + b_{N-1} s_{N-1}^2 + c_{N-1} s_{N-1} + d_{N-1} = x_{N-1} \\ f_{x}^{N-1}(s_{N}) & = a_{N-1} s_{N}^3 + b_{N-1} s_{N}^2 + c_{N-1} s_{N} + d_{N-1} = x_{N} \\ \end{align*}\]

그리고, 향후 기술을 위해 위의 식들을 행렬 형태로 정리해놓도록 하자. 눈치챘겠지만 이렇게 조건식들을 행렬형태로 이어 붙여서 마지막에 역행렬을 똭! 하고 곱해서 계수들을 구할 것이다.

\[\left[ \begin{array}{ccccccccccc} s_{1}^3 & s_{1}^2 & s_{1} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{3}^3 & s_{3}^2 & s_{3} & 1 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N-1}^3 & s_{N-1}^2 & s_{N-1} & 1 \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N}^3 & s_{N}^2 & s_{N} & 1 \\ \end{array} \right] \left[ \begin{array}{c} a_{1} \\ b_{1} \\ c_{1} \\ d_{1} \\ a_{2} \\ b_{2} \\ c_{2} \\ d_{2} \\ ... \\ a_{N-1} \\ b_{N-1} \\ c_{N-1} \\ d_{N-1} \\ \end{array} \right] = \left[ \begin{array}{c} x_{1} \\ x_{2} \\ x_{2} \\ x_{3} \\ ... \\ x_{N-1} \\ x_{N} \\ \end{array} \right]\]

Step 4: 1계 미분에 대한 식을 만들기

우리는 모든 점에서 부드럽다고 말할 수 있는 (2번 미분가능한) spline을 만들고 있다. 따라서, 주어진 모든 구간에서 1계 미분 또한 연속이어야 하고, 이는 아래와 같이 \(N-2\)개의 조건식으로 나타낼 수 있다.

\[\begin{align*} f_{x}^{\prime 1}(s_{2}) = f_{x}^{\prime 2}(s_{2}) & \Rightarrow 3 a_{1} s_{2}^2 + 2 b_{1} s_{2} + c_{1} = 3 a_{2} s_{2}^2 + 2 b_{2} s_{2} + c_{2} \\ f_{x}^{\prime 2}(s_{3}) = f_{x}^{\prime 3}(s_{3}) & \Rightarrow 3 a_{2} s_{3}^2 + 2 b_{2} s_{3} + c_{2} = 3 a_{3} s_{3}^2 + 2 b_{3} s_{3} + c_{3} \\ f_{x}^{\prime 3}(s_{4}) = f_{x}^{\prime 4}(s_{4}) & \Rightarrow 3 a_{3} s_{4}^2 + 2 b_{3} s_{4} + c_{3} = 3 a_{4} s_{4}^2 + 2 b_{4} s_{4} + c_{4} \\ & ... \\ f_{x}^{\prime N-2}(s_{N-1}) = f_{x}^{\prime N-1}(s_{N-1}) & \Rightarrow 3 a_{N-2} s_{N-1}^2 + 2 b_{N-2} s_{N-1} + c_{N-2} = 3 a_{N-1} s_{N-1}^2 + 2 b_{N-1} s_{N-1} + c_{N-1} \\ \end{align*}\]

Step 3에서 만들기 시작했던 행렬에 이 조건식들을 마저 이어붙이면 아래와 같다. 중간정산을 한번 하자면 Step 3-4를 통해 \(2 * (N-1) + N-2 = 3N - 4\)개 조건식을 확보했다.

\[\left[ \begin{array}{ccccccccccc} s_{1}^3 & s_{1}^2 & s_{1} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{3}^3 & s_{3}^2 & s_{3} & 1 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N-1}^3 & s_{N-1}^2 & s_{N-1} & 1 \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N}^3 & s_{N}^2 & s_{N} & 1 \\ \hline 3 s_{2}^{2} & 2 s_{2} & 1 & 0 & -3 s_{2}^{2} & -2 s_{2} & -1 & 0 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & 0 & 3 s_{N-1}^{2} & 2 s_{N-1} & 1 & 0 & -3 s_{N-1}^{2} & -2 s_{N-1} & -1 & 0 \\ \end{array} \right] \left[ \begin{array}{c} a_{1} \\ b_{1} \\ c_{1} \\ d_{1} \\ a_{2} \\ b_{2} \\ c_{2} \\ d_{2} \\ ... \\ a_{N-1} \\ b_{N-1} \\ c_{N-1} \\ d_{N-1} \\ \end{array} \right] = \left[ \begin{array}{c} x_{1} \\ x_{2} \\ x_{2} \\ x_{3} \\ ... \\ x_{N-1} \\ x_{N} \\ \hline 0 \\ ... \\ 0 \\ \end{array} \right]\]

Step 5: 2계 미분에 대한 식을 만들기

다음 스텝도 뭐 별것 없다. 모든 점에서 2번 미분이 가능한 spline을 만들고 싶으므로 2계 미분이 연속임을 보장하도록 아래와 같이 \(N-2\)개의 조건식을 추가하자. 이제는 적기도 귀찮다. 쿨하게 바로 행렬에 이어 붙이도록 하겠다. 지금까지 \(2 * (N-1) + N-2 + N-2= 4N - 6\)개 조건식을 확보했다.

\[\begin{align*} f_{x}^{\prime \prime 1}(s_{2}) = f_{x}^{\prime \prime 2}(s_{2}) & \Rightarrow 6 a_{1} s_{2} + 2 b_{1} = 6 a_{2} s_{2} + 2 b_{2} \\ f_{x}^{\prime \prime 2}(s_{3}) = f_{x}^{\prime \prime 3}(s_{3}) & \Rightarrow 6 a_{2} s_{3} + 2 b_{2} = 6 a_{3} s_{3} + 2 b_{3} \\ f_{x}^{\prime \prime 3}(s_{4}) = f_{x}^{\prime \prime 4}(s_{4}) & \Rightarrow 6 a_{3} s_{4} + 2 b_{3} = 6 a_{4} s_{4} + 2 b_{4} \\ & ... \\ f_{x}^{\prime \prime N-2}(s_{N-1}) = f_{x}^{\prime \prime N-1}(s_{N-1}) & \Rightarrow 6 a_{N-2} s_{N-1} + 2 b_{N-2} = 6 a_{N-1} s_{N-1} + 2 b_{N-1} \\ \end{align*}\] \[\left[ \begin{array}{ccccccccccc} s_{1}^3 & s_{1}^2 & s_{1} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{3}^3 & s_{3}^2 & s_{3} & 1 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N-1}^3 & s_{N-1}^2 & s_{N-1} & 1 \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N}^3 & s_{N}^2 & s_{N} & 1 \\ \hline 3 s_{2}^{2} & 2 s_{2} & 1 & 0 & -3 s_{2}^{2} & -2 s_{2} & -1 & 0 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & 0 & 3 s_{N-1}^{2} & 2 s_{N-1} & 1 & 0 & -3 s_{N-1}^{2} & -2 s_{N-1} & -1 & 0 \\ \hline 6 s_{2} & 2 & 0 & 0 & -6 s_{2} & -2 & 0 & 0 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & 0 & 6 s_{N-1} & 2 & 0 & 0 & -6 s_{N-1} & -2 & 0 & 0 \\ \end{array} \right] \left[ \begin{array}{c} a_{1} \\ b_{1} \\ c_{1} \\ d_{1} \\ a_{2} \\ b_{2} \\ c_{2} \\ d_{2} \\ ... \\ a_{N-1} \\ b_{N-1} \\ c_{N-1} \\ d_{N-1} \\ \end{array} \right] = \left[ \begin{array}{c} x_{1} \\ x_{2} \\ x_{2} \\ x_{3} \\ ... \\ x_{N-1} \\ x_{N} \\ \hline 0 \\ ... \\ 0 \\ \hline 0 \\ ... \\ 0 \\ \end{array} \right]\]

Step 6: 어라? 식이 2개가 모자르다

모두에서 밝혔듯이 우리는 총 \(4 * (N - 1)\)개의 식이 필요했다. 그런데, 뽑아낼 수 있는 조건은 얼추 다 뽑아낸 것 같은데, 아직도 식이 2개가 모자르다. 이 2개의 식은 당신의 입맛에 따라 만들어주면 되는데, 이 2개의 식을 어떻게 설계하냐에 따라 spline의 형태가 달라지게 된다. 유명한 것들은 나름 이름도 붙어있다 (natural spline, periodic spline, not a knot spline 등등). 필자는 그 중 가장 간단한, quadratic spline, 즉, \(f_{x}^{1}(s)\)과 \(f_{x}^{N-1}(s)\)은 2차 함수라고 가정하는 조건식을 추가해보도록 하겠다. 이 조건은 간단하게 \(a_{1} = a_{N_1} = 0\)으로 표현이 가능하다.

Step 7: 역행렬을 구하자

Step 6까지 구한 조건들을 마저 행렬로 표현하면 아래와 같다. 이제 아래 식에서 정사각행렬의 역행렬만 잘 구하면 계수들을 쉽게 구해줄 수 있다. 아마도 독자들 중에는 “저게 역행렬이 늘 있나?” 싶은 분들도 있으리라 (물론, 당연히 데이터에 같은 샘플은 없어야 하겠다). 여기서 바로 Step 6의 진가가 드러나는데, 이 2개의 조건을 잘 설정해서 정사각행렬이 singular하지 않도록 만들어줄 수 있다. 물론, 필자가 고른 quadratic spline도 그 중 하나인데, 다른 방법들도 궁금하다면 구글링해보길 바란다.

\[\left[ \begin{array}{ccccccccccc} s_{1}^3 & s_{1}^2 & s_{1} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{2}^3 & s_{2}^2 & s_{2} & 1 & 0 & ... & 0 \\ 0 & 0 & 0 & 0 & s_{3}^3 & s_{3}^2 & s_{3} & 1 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N-1}^3 & s_{N-1}^2 & s_{N-1} & 1 \\ 0 & 0 & 0 & 0 & 0 & ... & 0 & s_{N}^3 & s_{N}^2 & s_{N} & 1 \\ \hline 3 s_{2}^{2} & 2 s_{2} & 1 & 0 & -3 s_{2}^{2} & -2 s_{2} & -1 & 0 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & 0 & 3 s_{N-1}^{2} & 2 s_{N-1} & 1 & 0 & -3 s_{N-1}^{2} & -2 s_{N-1} & -1 & 0 \\ \hline 6 s_{2} & 2 & 0 & 0 & -6 s_{2} & -2 & 0 & 0 & 0 & ... & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & ... & 0 & 6 s_{N-1} & 2 & 0 & 0 & -6 s_{N-1} & -2 & 0 & 0 \\ \hline 1 & 0 & ... & ... & ... & ... & ... & ... & ... & ... & 0 \\ 0 & ... & ... & ... & ... & ... & 0 & 1 & 0 & 0 & 0 \\ \end{array} \right] \left[ \begin{array}{c} a_{1} \\ b_{1} \\ c_{1} \\ d_{1} \\ a_{2} \\ b_{2} \\ c_{2} \\ d_{2} \\ ... \\ a_{N-1} \\ b_{N-1} \\ c_{N-1} \\ d_{N-1} \\ \end{array} \right] = \left[ \begin{array}{c} x_{1} \\ x_{2} \\ x_{2} \\ x_{3} \\ ... \\ x_{N-1} \\ x_{N} \\ \hline 0 \\ ... \\ 0 \\ \hline 0 \\ ... \\ 0 \\ \hline 0 \\ 0 \\ \end{array} \right]\]

이게 뭐 구현이 대단할 것이 있겠냐만은, 혹시나 과제를 위해 밤새는 대학(원)생들을 위해 여기에 예제 코드를 남겨둔다.