cubic spline interpolation

   cubic spline은 주어진 점을 매끄럽게 연결하는 알고리즘이다. 왜 cubic 이냐 하면, 두 점을 잇는 곡선을 3차 다항식(a₀ + a₁x + a₂x² + a₃x³)으로 사용하기 때문이다. 즉, 다시 말해, 서로 떨어져 있는 두 점 사이를 연결해야 하는데(그래서 interpolation^[각주:1]), 그 연결하는 선을 3차 다항식으로 만들고자 하는 것이다. 이와 같은 목적에 부합하는 알고리즘 중 널리 사용되는 것의 원리는 다음과 같다. 

우리가 갖고 있는 점은 위의 녹색점 뿐이라고 하자. 그 중간을 곡선으로 매끄럽게 연결하고 싶다. 즉, 각 데이터 x_i, x_i+1은 곡선 S_i로 연결을 시키는데, 각 S_i가 3차 다항식이라고 가정하는 것이다. 이 상황에서, 각 x_i+1에서 두 곡선 S_i와 S_i+1이 부드럽게 연결되기 위해서는 어떠해야 할까? 상식 선에서 크게 벗어나지 않는 다음과 같은 조건들을 생각해 볼 수 있다.

1. 각 점에서 원래의 데이터의 값 x_i을 곡선이 지나간다.
2. 각 점에서 두 곡선이 만나야 한다. 즉 각 점 x_i 에서 두 곡선의 함수값이 같아야 한다. 
3. 각 점에서 두 곡선의 매끄럽기가 같게 한다. 즉, 각 점 x_i에서 두 곡선의 미분값이 같게 한다.
4. 두 곡선이 '더욱' 매끄럽게 만나게 하자. 즉, 각 점 x_i에서 두 곡선의 이차미분값을 같게 한다.

위 사안 중 1,2 번은 당연히 필수적인 것이고, 나머지 사안들은 '부드럽게 연결'시키기 위해 사용할 수 있는 조건들이라 할 수 있겠다. 4번은 3차 다항식의 계수값을 계산할 수 있게 하기 위해 사용하는 목적도 있다. 또한 양 끝 쪽에서의 조건도 필요한데, 가장 양 끝에서 두 번 미분한 값이 0 이거나, 한 번 미분한 값이 원래의 데이터의 미분값과 같다는 조건 둘 중 하나를 이용하곤 한다. 위 조건은 수학적으로 다음과 같이 표현할 수 있다.

위의 조건에 대하여, 각 Si를 3차 다항식으로 한 뒤 각 계수를 풀면 풀린다.

각 계수를 계산해서 3차 다항식을 구한 뒤 두 점 사이에 있는 값들을 채워 나가는 C++ 코드는 다음과 같다. (코드는 Richard L. Burden; J. Douglas Faires, Numerical Analysis, 8th edition에 있는 알고리즘을 코드로 나타낸 것이다).

bool cubic_spline(std::vector<double>* x_series, std::vector<double>* y_series, std::vector<double> *destX, std::vector<double>* destY)

{   

    int n = __min((int)x_series->size()-1, (int)y_series->size()-1);

    // Step 1.

    double *h = new double[n+1];

    double *alpha = new double[n+1];

    int i = 0;

    for(i = 0; i<=n-1; i++){

        h[i] = (*x_series)[i+1] - (*x_series)[i];

    }

 

    // Step 2.

    for(i = 1; i<=n-1;i++){

        alpha[i]= 3*((*y_series)[i+1]-(*y_series)[i])/h[i]-3*((*y_series)[i]-(*y_series)[i-1])/h[i-1];

    }

 

    // Step 3.

    double *l = new double[n+1];

    double *u = new double[n+1];

    double *z = new double[n+1];

    double *c = new double[n+1];

    double *b = new double[n+1];

    double *d = new double[n+1];

 

    l[0] = 1; u[0] = 0; z[0] = 0;

 

    // Step 4.

    for(i = 1; i<=n-1; i++){

        l[i] = 2*((*x_series)[i+1] - (*x_series)[i-1]) - h[i-1]*u[i-1];

        u[i] = h[i]/l[i];

        z[i] = (alpha[i] - h[i-1]*z[i-1]) / l[i];

    }

 

    // Step 5.

    l[n] = 1;     z[n] = 0;     c[n] = 0;

 

    // Step 6.

    for(i = n-1; i>=0; i--){

        c[i] = z[i] - u[i]*c[i+1];

        b[i] = ((*y_series)[i+1] - (*y_series)[i])/h[i] - h[i]*(c[i+1] + 2*c[i])/3;

        d[i] = (c[i+1] - c[i]) / (3*h[i]);

    }

 

    for(i = 0; i<=n-2;i++){

        double x = (*x_series)[i];

        double inc = ((*x_series)[i+1] - (*x_series)[i])*0.1;

        for(; x < (*x_series)[i+1]; x+=inc){

            double x_offset = x - (*x_series)[i];

            double Sx = (*y_series)[i] + b[i]*x_offset + c[i]*x_offset*x_offset + d[i]*x_offset*x_offset*x_offset;

            if(destX != NULL){

                destX->push_back(x);

            }

            if(destY != NULL){

                destY->push_back(Sx);

            }

        }           

    }

    delete [] h;

    delete [] alpha;

    delete [] l;

    delete [] u;

    delete [] z;

    delete [] c;

    delete [] b;

    delete [] d;

 

 

    return true;

}

위 코드에서 녹색으로 칠한 부분에 0.1 을 곱한 것은, 원래 데이터 사이를 10개로 나누기 위해서이다. 위 알고리즘을 이용한 것과 matlab에 있는 것을 비교한 그림은 다음과 같다.

위의 행이 위 코드를 이용한 것이고, 아래 행이 matlab으로 테스트 해 본 결과이다. 열린 원이 원래의 데이터 였고, 곡선은 그 데이터들을 부드럽게 연결하는 것이다.

여기까지가 cubic spline interpolation 에 대한 설명이다.



위 그림은 사실 내가 cubic spline 을 구현하면서 알게 된 문제를 보여 주기 위해 찾은 예제인데, 그것은 위의 알고리즘의 경우 monotonicity가 보장되지 않는다는 것이다. 즉, 원래 데이터는 계속 증가만 하는데 cubic spline 으로 연결한 곡선은 내려간다거나 혹은 그 반대. 이것을 해결하기 위한 것이 바로 monotone cubic Hermite interpolation 이다. 

1. inter-, 즉 '사이'값을 추정한다는 것 [본문으로]