Geometry & Recognition

$$\begin{gathered} \{{\mathbf{C}}_k\}=\text{argmin}(L) \\ L=\sum _{k=0}^{N-1}{|{\mathbf{P}}_k-\sum _{i=0}^n{b}_{i,n}(t_k){\mathbf{C}}_i|}^2 \end{gathered}$$

$$\begin{array}{c}{t}_k=d_k/d_{N-1},d_k=d_{k-1}+|{\mathbf{P}}_k-{\mathbf{P}}_{k-1}|\end{array}$$

$$\text{Bernstein basis polynomial: }{b}_{i,n}(t)=\left(\begin{array}{c}n\\ i\end{array}\right)t^i(1-t{)}^{n-i}=\sum _{j=0}^n{t}^j{M}_{ji}$$

$$M_{ji}=\frac{n!}{(n-j)!}\frac{(-1{)}^{j+i}}{i!(j-i)!}=(-1{)}^{j+i}{M}_{jj}\left(\begin{array}{c}j\\ i\end{array}\right)(j\ge i)$$

$$\text{note: }{M}_{jj}=\left(\begin{array}{c}n\\ j\end{array}\right),M_{ji}=0(j<i)$$

// calculate binomial coefficients using Pascal's triangle
void binomialCoeffs(int n, double** C) {
    // binomial coefficients
    for (int i = 0; i <= n; i++)
        for (int j = 0; j <= i; j++)
            if (j == 0 || j == i) C[i][j] = 1;
            else                  C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
std::vector<double> basisMatrix(int degree) {
    const int n = degree + 1;
    std::vector<double* > C(n);
    std::vector<double> Cbuff(n*n);
    for (int i = n; i-->0;) 
        C[i] = &Cbuff[i * n];
    // C[degree, k]; k=0,..,degree)
    binomialCoeffs(degree, &C[0]);
    // seting the diagonal;
    std::vector<double> M(n * n, 0); // lower triangle; 
    for (int j = 0; j <= degree; j++)
        M[j * n + j] = C[degree][j];
    // compute the remainings;
    for (int i = 0; i <= degree; i++) 
        for (int j = i + 1; j <= degree; j++)
            M[j * n + i] = ((j + i) & 1 ? -1 : 1) * C[j][i] * M[j * n + j];
    return M;
}

주어진 데이터를 잘 피팅하는 직선을 찾기 위해서는 데이터를 이루는 각 점의 $y$ 값과 같은 위치에서 구하려는 직선과의 차이 (residual)의 제곱을 최소화시키는 조건을 사용한다. 그러나 직선의 기울기가 매우 커지는 경우에는  데이터에 들어있는 outlier에 대해서는 그 차이가 매우 커지는 경우가 발생할 수도 있어 올바른 피팅을 방해할 수 있다. 이를 수정하는 간단한 방법 중 하나는 $y$값의 차이를 이용하지 않고 데이터와 직선 간의 최단거리를 이용하는 것이다. 

수평에서 기울어진 각도가 $\theta$이고 원점에서 거리가 $s$인 직선의 방정식은 

$$x\sin \theta -y\cos \theta +s=0$$

이고, 한 점 $(x_i,y_i)$에서 이 직선까지 거리는

$$d_i=|\sin \theta y_i-\cos \theta x_i+s|$$

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $\theta ,s$을 구해야 한다. $$L=\sum _i(\sin \theta x_i-\cos \theta y_i+s{)}^2$$$$(\theta ,s)=\text{argmin}(L)$$

$\theta$와 $s$에 대한 극값 조건에서 

$$\frac{1}{2}\frac{\partial L}{\partial \theta }=\frac{1}{2}\sin 2\theta \sum _i(x_i^2-y_i^2)-\cos 2\theta \sum x_i{y}_i+s\sin \theta \sum _i{x}_i+s\cos \theta \sum _i{y}_i=0$$

$$\frac{1}{2}\frac{\partial L}{\partial s}=\cos \theta \sum y_i-\sin \theta \sum _i{x}_i-Ns=0$$

주어진 데이터의 질량중심계에서 계산을 수행하면 $\sum _i{x}_i=\sum _i{y}_i=0$ 이므로 데이터의 2차 모멘트를 $$A=\sum _i(x_i^2-y_i^2),B=\sum _i{x}_i{y}_i$$로 놓으면 직선의 파라미터를 결정하는 식은

$$\begin{gathered} \frac{1}{2}A\sin 2\theta -B\cos 2\theta =0\to \tan 2\theta =\frac{2B}{A} \\ s=0 \end{gathered}$$

두 번째 식은 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$\tan \theta =\frac{-A\pm \sqrt{A^2+(2B{)}^2}}{2B}$$

두 해 중에서 극소값 조건을 만족시키는 해가 직선을 결정한다. 그런데

$$\frac{1}{2}\frac{{\partial }^{2}L}{\partial {\theta }^2}=A\cos 2\theta +2B\sin 2\theta =\pm \sqrt{A^2+(2B{)}^2}>0$$

이므로 위쪽 부호로 직선($x\sin \theta =y\cos \theta$)이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$(-A+\sqrt{A^2+(2B{)}^2})(x-\overline{x})=2B(y-\overline{y})$$

여기서 주어진 데이터의 질량중심은 원좌표계에서

$$\overline{x}=\frac{1}{N}\sum _i{x}_i,\overline{y}=\frac{1}{N}\sum _i{y}_i$$

이다. 또한 원좌표계에서 $A$와 $B$의 계산은 

$$A=\sum _i[(x_i-\overline{x}{)}^2-(y_i-\overline{y}{)}^2],B=\sum (x_i-\overline{x})(y_i-\overline{y})$$

이 결과는 데이터 분포에 PCA를 적용해서 얻은 결과와 동일하다. PCA에서는 공분산 행렬의 고유값이 큰 쪽에 해당하는 고유벡터가 직선의 방향을 결정했다. https://kipl.tistory.com/211.  또한 통계에서는 Deming regression이라고 불린다.

최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 parameter vector를 구하는 것이다.
$$\underset{\mathbf{x}}{\text{argmin}}|\mathbf{A}.\mathbf{x}-\mathbf{b}{|}^2.$$

여기서 design matrix $\mathbf{A}$는 추정에 사용이 된 basis 함수를 각각의 독립변수에서 계산한 값이고, $\mathbf{x}$는 구하고자 하는 parameter vector이며, $\mathbf{b}$는 측정값 vector이다. 예를 들면 주어진 측정값을 $n-1$차의 다항식을 이용하여서 피팅하려고 하는 경우에는 parameter는 다항식의 계수가 되며 $n$-차원의 vector space을 형성하게 되며, $\mathbf{A}$는

$$A_{ij}=(x_i{)}^j,j=0,...,n-1$$

가 일 것이다. 일반적으로 $A_{ij}$는 $x_i$에서 계산된 basis-함수의 값이 된다. 위의 식을 $\mathbf{x}$에 대해서 미분을 하면 극값을 취하면 최소자승 문제는 아래의 행렬식을 푸는 문제가 된다

$$({\mathbf{A}}^T.\mathbf{A}).\mathbf{x}={\mathbf{A}}^T.\mathbf{b}.$$

${\mathbf{A}}^T.\mathbf{A}$ 은 $n\times n$ matrix다. 이 행렬이 역행렬을 가지게 되면

$$\mathbf{x}=({\mathbf{A}}^T.\mathbf{A}{)}^{-1}.(\mathbf{A}.\mathbf{b}),$$

를 하여서 원하는 parameter vector를 얻을 수 있다. 그러나 피팅 문제에 따라 행렬 ${\mathbf{A}}^T.\mathbf{A}$가 매우 singular 해져 역행렬을 구할 수 없게 되는 경우에 종종 생긴다. 예를 들면, 저주파의 신호를 고주파 기저 함수를 이용하여서 최소자승법을 사용하고자 하는 경우 등에 이러한 문제에 부딪히게 된다. 이런 경우에는 직접적으로 ${\mathbf{A}}^T.\mathbf{A}$의 역행렬을 구하는 방법을 이용할 수 없고

$$\mathbf{A}.\mathbf{x}=\mathbf{b}$$

의 식을 $\mathbf{A}$의 SVD(Singular Value Decomposition)를 이용하여서 풀 수가 있다. $\mathbf{A}$를 SVD 하면 ${\mathbf{A}}_{m\times n}={\mathbf{U}}_{m\times n}.{\mathbf{w}}_{n\times n}.{\mathbf{V}}_{n\times n}^T$의 형태로 분해할 수 있다. 여기서 $\mathbf{w}=\text{diag}(\underset{\text{nonzero}}{\underbrace{w_0,w_1,...}},0,..,0)$로 쓰여지는 대각행렬이다. matrix $\mathbf{U}$와 $\mathbf{V}$의 column vector를 사용하면
$$\mathbf{A}=\sum _{w_k\ne 0}{w}_k{\mathbf{u}}_k\otimes {\mathbf{v}}_k^T$$

의 형태로 쓰인다. ${\mathbf{u}}_k$는 $\mathbf{U}$의 $k$-번째 열벡터이고, ${\mathbf{v}}_k$는 $\mathbf{V}$의 $k$-번째 열벡터로 각각 orthonormal basis를 형성한다. parameter 벡터를 $\{{\mathbf{v}}_k\}$ basis로 전개를 하면 영이 아닌 singularvalue에 해당하는 성분만 가지게 된다. 구체적으로 위의 $\mathbf{A}$ 분해와 ${\mathbf{u}}_j^T.{\mathbf{u}}_k={\delta }_{jk}$, 그리고 $\sum _k{\mathbf{v}}_k\otimes {\mathbf{v}}_k^T={\mathbf{I}}_{n\times n}$임을 이용하면,

$$\begin{array}{c}{\mathbf{v}}_k^T.\mathbf{x}={\mathbf{u}}_k^T.\mathbf{b}/w_k,w_k\ne 0,\\ {\mathbf{v}}_k^T.\mathbf{x}=0,w_k=0,\\ \to \mathbf{x}=\sum _{w_k\ne 0}({\mathbf{u}}_k^T.\mathbf{b}/w_k){\mathbf{v}}_k,\end{array}$$

이어서 위의 해를 구할 수 있다. 이 해는 $|\mathbf{A}.\mathbf{x}-\mathbf{b}{|}^2$를 최소화한다.

int svd(double *A, int m, int n, double* w, double *V); // from cmath libary.
void fit_func(double x, double val[], int n) {          // polynomial fitting sample;
    val[0] = 1;
    for(int i = 1; i < n; ++i)
        val[i] = x * val[i - 1];
}
#define EPSILON 1.E-8
int svd_fit(const double x[], const double y[], const int m, const int n,
            void (*fit_func)(double , double [], int ),
            double params[],
            double *error)
{
    double *A = new double [m * n];
    double *w = new double [n];
    double *V = new double [n * n];
    // evaluate design matrix;
    for (int i = 0; i < m; ++i)
        fit_func(x[i], &A[i * n + 0], n) ;

    svd(A, m, n, w, V);
    // now A becomes U matrix;
    // truncate small singular values;
    double wmax = 0;
    for (int i = 0; i < n; ++i)
        if (w[i] > wmax) wmax = w[i];
    double thresh = wmax * EPSILON;
    for (int i = 0; i < n; ++i)
        if (w[i] < thresh) w[i] = 0;
    
    // back substitution;
    double *tmp = new double [n];
    for (int j = 0; j < n; ++j) {
        double s = 0;
        if (w[j]) {
            for (int i = 0; i < m; ++i)
                s += A[i * n + j] * y[i];
            s /= w[j];
        }
        tmp[j] = s;
    }
    for (int j = 0; j < n; ++j) {
        double s = 0;
        for (int jj = 0; jj < n; ++jj)
            s += V[j * n + jj] * tmp[jj];
        params[j] = s;
    };

    //estimate error;
    *error = 0;
    for (int i = 0; i < m; ++i) {
        fit_func(x[i], &A[i * n + 0], n); //use A as a tmp buffer;
        double sum = 0;
        for (int j = 0; j < n; ++j) sum += params[j] * A[i * n + j] ;
        double err = (y[i] - sum);
        *error += err * err ;
    }
    delete[] A; delete[] w; delete[] V;
    delete[] tmp;
    return 1;
}