Geometry & Recognition

주어진 점집합을 기술하는 직선을 얻는 방법 중에 하나로 각각의 점들이 직선에서 벗어난 거리의 제곱을 더한 값(square of residual)을 최소화시키는 기울기와 절편을 찾는 최소자승법(least square method: linear regression)이 있다. 점집합이 $\{(x_i, y_i\}$를 직선 $y=ax +b$로 피팅을 하는 경우 직선에서 벗어난 정도(residual)는 직선까지의 거리를 사용할 수도 있고 또는 $y$ 값의 차이를 이용할 수도 있다. 먼저 해를 closed form으로 쓸 수 있는 $y$ 값의 offset을 residual로 사용하자. 그러면 fitting error는 각 점에서 residual의 제곱의 합으로 주어진다. 

$$ R^2(a, b) = \sum_i  \left| y_i - (ax_i + b) \right|^2 $$

여기서 주어진 점집합의 moment를 $$ S_{xx} = \sum x_i^2 , ~S_{yy} = \sum y_i^2, ~S_{xy} = \sum x_i y_i, ~S_x = \sum x_i, ~S_y = \sum y_i$$로 놓으면 

\begin{align}   R^2(a,b) &= S_{yy} + a^2 S_{xx} + n b^2 - 2a S_{xy} + 2ab S_x - 2 b S_y\end{align}

이다. 주어진 점집합을 피팅하는 직선의 파라미터는 타원의 방정식으로 주어짐을 알 수 있고, 주워진 타원상의 임의의 파라미터에 해당하는 직선의 fittign error는 항상 일정한 값을 가짐을 알 수 있다.

좌표의 원점을 각 성분의 평균값만큼 이동하면 $S_x = S_y = 0$이 되어 더 식이 단순해진다. 분산
$$ \sigma_x^2 = \frac{1}{n} S_{xx} - \bar{x}^2 ~\to ~\sigma_{x}^2  =\frac{1}{n} \vec{x}\cdot \vec{x}$$

$$  \sigma_y^2 = \frac{1}{n} S_{yy} - \bar{y}^2 ~\to~ \sigma_y^2 = \frac{1}{n} \vec{y}\cdot \vec{y}$$ 공분산

$$\text{cov}(x,y)= \frac{1}{n} S_{xy} - \bar{x}\bar{y} \to \text{cov}(x,y) = \frac{1}{n} \vec{x}\cdot\vec{y}$$

을 써서 표현하면,

\begin{align} R^2 &= n\sigma_x^2 a^2 -2n \text{cov}(x,y) a + n b^2 + n \sigma_y^2 \\ &=n \sigma_x^2 \left( a - \frac{\text{cov}(x,y)}{\sigma_x^2 } \right)^2 + n b^2 + n \frac{\text{cov}(x,y)^2}{\sigma_x^2} + n \sigma_y^2 \\ &\ge n \frac{\text{cov}(x,y)^2}{\sigma_x^2} + n \sigma_y^2   \end{align} 

따라서 residual을 최소로 만들어 주는 직선의 기울기 $a$와 $y-$절편 $b$는

\begin{align} a &= \frac{\text{cov}(x, y) }{\sigma_x^2 } \\ b&=0\end{align}

로 주어지는데, 이 직선은 원점을 통과하는 직선이 된다. 원래의 좌표계로 돌아가면 기울기는 원점의 이동에 무관하므로 변화가 없고 직선이 $(\bar{x}, \bar{y})$을 통과해야 하므로 절편 $b$값은

$$b = \bar{y} - a \bar{x}$$

로 주어진다. 상관계수 $r(x,y)$를 이용하면 Fitting이 잘된 정도를 정량적으로 표현이 된다.

$$ r = \frac{\text{cov}(x,y) }{\sqrt{ \sigma_x^2  \sigma_y^2 } }$$

즉,

$$ R = n\sigma_y^2 ( r^2 + 1)$$

로 주어진다.

최소 자승법 문제는 추정치와 실제의 측정치와의 차이를 최소로 만드는 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 ,  \quad 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}(\underbrace{w_0, w_1,...}_{\text{nonzero}},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{gather}  \mathbf{v}_k^T . \mathbf{x} = \mathbf{u}_k^T . \mathbf{b} / w_k, \quad w_k \ne 0, \\                    \mathbf{v}_k^T . \mathbf{x} = 0, \quad w_k = 0, \\  \rightarrow ~~\mathbf{x} = \sum _{w_k \ne 0 } ( \mathbf{u}_k^T . \mathbf{b} / w_k)  \mathbf{ v} _k , \end{gather}

이어서 위의 해를 구할 수 있다. 이 해는 $|\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;
}

일반적인 이차곡선은 다음의 이차식으로 표현이 된다:

$$ F(x, y)=ax^2 + bxy + cy^2 +d x + ey + f=0$$

6개의 계수는 모두 독립적이지 않고 어떤 종류의 이차곡선인가에 따라 제약조건이 들어온다. 주어진 점 데이터 $\{ (x_i, y_i) | i=1,2,...,n\}$를 피팅하는 이차곡선을 각각의 데이터 점에서 대수적 거리 ($=|F(x_i, y_i)|$)의 제곱을 최소로 하는 조건하에서 찾도록 하자:

$$ \epsilon = \sum_{i} |F(x_i , y_i )|^2 \rightarrow \text{min}$$

타원으로 피팅하는 경우 여러 가지 제약조건을 줄 수 있지만(e.g.: $a+c=1$, $\sqrt{a^2+b^2+...+f^2 }=1$ 등등) 여기서는 이차 곡선이 타원이기 위해서는 2차 항의 계수로 만드는 행렬 $\left(\begin{array}{cc} a & b/2\\b/2 & c\end{array}\right)$의 determinant (= $ac- b^2/4>0$)가 양수이어야 한다는 조건에서 (타원을 회전+평행 이동시키면 표준형 타원으로 만들 수 있는데, 이때 두 주축의 계수 곱이 determinant다. 따라서 타원이기 위해서는 값이 양수여야 된다) 

$$ 4ac - b^2 = 1$$

로 선택하도록 하자. 이 제약조건을 넣은 타원 피팅은 다음 식을 최소화하는 계수 벡터 $\mathbf{a}=[a,b,c,d,e,f]^T$을 찾는 문제가 된다:

\begin{gather}    \epsilon = \mathbf{a}^T. \mathbf{S}. \mathbf{a}   -\lambda ( \mathbf{a}^T.\mathbf{C}.\mathbf{a}-1)\end{gather}

여기서 $6\times 6$ scattering matrix $\mathbf{S}$는

$$ \mathbf{S}= \mathbf{D}^T. \mathbf{D}.$$

$n\times 6$ design matrix $\mathbf{D}$는 

$$ \mathbf{D}=\left[ \begin{array}{cccccc} x_1^2& x_1 y_1 & y_1^2 & x_1 & y_1 & 1\\ x_2^2& x_2 y_2 & y_2^2 & x_2 & y_2 & 1\\ & & \vdots & &  \\ x_n^2& x_n y_n & y_n^2 & x_n & y_n & 1 \end{array}\right] .$$

그리고, $6\times 6$ 행렬 $\mathbf{C}$는

$$ \mathbf{C}= \left[ \begin{array} {cccccc} 0&0&2&0&0&0\\ 0&-1&0&0&0&0 \\ 2&0&0&0&0&0\\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \end{array} \right]  . $$

$\epsilon$을 $\mathbf{a}^T$에 대해 미분하면 최소 제곱 문제는 Lagrange multiplier를 고윳값으로 하는 고유 방정식의 해를 구하는 문제로 환원이 된다.

$$ \mathbf{S}.\mathbf{a} = \lambda \mathbf{C}.\mathbf{a}$$

그리고 주어진 고유값 $\lambda$에 해당하는 normalized 고유 벡터가 $\mathbf{a}$일 때 

$$\epsilon = \mathbf{a}^T. \mathbf{S}.\mathbf{a} = \lambda$$

이므로 최소의 양의 고윳값에 해당하는 고유 벡터가 해가 된다. Silverster의 law of inertia를 이용하면 위의 고유 방정식에서 양의 고윳값은 딱 1개만 존재함을 보일 수 있다. 고유값 $\lambda$에 해당하는 고유 벡터가 $\mathbf{a}$일 때 임의의 상수 $\mu$ 배를 한 $\mu\mathbf{a}$도 같은 고유값을 갖는 고유벡터다. normalized 고유벡터는 제약조건 $\mu^2 \mathbf{a}^T . \mathbf{C}. \mathbf{a} =1$을 만족하도록 크기를 조정하면 $ \mathbf{\tilde{a}} = \mathbf{a} / \sqrt{\mathbf{a}^T. \mathbf{C}. \mathbf{a}}$로 주어진다. 

이 일반화된 고유방정식의 풀이는 먼저 positive definite인 $\bf S$의 제곱근 행렬 $\bf Q=S^{1/2}$을 이용하면 쉽다. $\bf S$의 고유벡터를 구하면 eigendecomposition에 의해 $\bf  S = R\Lambda R^T$로 쓸 수 있으므로 제곱근 행렬은 $\bf  Q = R \Lambda ^{1/2} R^T$임을 쉽게 확인할 수 있다. 원래의 고유값 문제를 $\tt Q$을 이용해서 표현하면 

$$ \bf Q a = \lambda Q^{-1} C a = \lambda Q^{-1} C Q^{-1} Qa$$ 이므로 더 다루기 쉬운 $\bf Q^{-1} C Q^{-1}$의 고유값 문제로 환원이 됨을 알 수 있다.

추가로, 고유 방정식은 $\mathbf{C}$가 $3\times 3$ 크기의 block으로 나누어질 수 있으므로 $[a,b,c]^T$ 에 대한 고유 방정식으로 줄일 수 있어서 쉽게 해결할 수 있다. 물론 scattering matrix을 구성하는 moment의 개수가 많아서 matrix 연산 패키지를 사용하지 않는 경우 코드가 길어지게 된다. 

Note: 이 내용은 다음 논문을 정리한 것이다. Andrew Fitzgibbon, Maurizio Pilu, and Robert B. Fisher, Direct Least Square Fitting of Ellipses". IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 21, NO. 5, MAY 1999. 

https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/ellipse-pami.pdf

구현은  https://kipl.tistory.com/566