Linear Least Square Fitting: perpendicular offsets

주어진 데이터를 가장 잘 피팅하는 직선을 찾기 위해서는 일반적으로 각 데이터 점의 \(y\)값과 직선이 예측하는 \(y\)값의 차이, 즉 residual의 제곱합을 최소화하는 최소자승법을 사용한다. 이 방법은 계산이 간단하고 대부분의 경우 좋은 결과를 제공하지만, 오차를 \(y\)축 방향으로만 측정한다는 한계가 있다. 따라서 직선의 기울기가 매우 큰 경우에는 실제 점과 직선 사이의 최단거리는 작더라도 \(y\,\)방향 residual은 크게 나타날 수 있으며, 데이터에 포함된 outlier의 영향도 더욱 커져 올바른 직선 피팅을 방해할 수 있다.

이러한 문제를 완화하는 간단한 방법 중 하나는 \(y\,\)방향 residual 대신 데이터 점과 직선 사이의 수직거리(최단거리)를 오차로 사용하는 것이다. 이 경우 오차는 직선의 기울기에 의존하지 않는 기하학적 거리로 정의되므로, 기울기가 큰 직선에서도 보다 안정적인 피팅 결과를 얻을 수 있다. 이러한 방법을 Orthogonal Least Squares 또는 Total Least Squares(TLS)라 하며, 특히 \(x\)축과 \(y\)축 모두에 측정 오차가 존재하는 경우에 적합하다.

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

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

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

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

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $(\theta, s)$을 구해야 한다. $$ L = \sum_i d_i^2 = \sum_i \big( x_i \sin \theta - y_i \cos \theta +s \big)^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_i x_i y_i + s \cos\theta \sum_i x_i + s \sin \theta \sum_i y_i = 0$$

$$ \frac{1}{2}\frac{\partial L}{\partial s}=\sin \theta \sum_i x_i -\cos \theta \sum_i y_i  + N s=0$$

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

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

두 번째 식은 원점에서 직선까지 거리 $s$가 0임을 의미하므로 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$ \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$$이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$  \boxed{\left(-A+ \sqrt{A^2+ (2B)^2} \right)  (x-\bar{x}) = 2B (y - \bar{y})   }  $$

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

$$ \bar{x} = \frac{1}{N} \sum_i x_i, \quad \bar{y} = \frac{1}{N} \sum_i y_i$$

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

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

단, $B=0$이면 위 식을 그대로 사용하지 말고 별도로 처리해야 한다. $A>0$이면 최적 직선은 $y=\bar y$인 수평선이고, $A<0$이면 $x=\bar x$인 수직선이다. 특히 $A=B=0$이면 중심화한 데이터의 $x$, $y$ 방향 분산이 같고 공분산이 0이다. 즉 데이터가 2차 모멘트 수준에서 모든 방향으로 동일하게 퍼져 있어 특별한 주축이 존재하지 않는다. 이 경우 질량중심을 지나는 모든 직선이 동일한 수직거리 제곱합을 가지므로, 유일한 최적 직선을 정할 수 없다. 원 위에 균등하게 놓인 점들, 정사각형 또는 정삼각형의 꼭짓점들이 대표적인 예다. 이 결과는 PCA를 적용한 결과와 동일하다. 중심화한 데이터의 공분산 행렬에서 가장 큰 고윳값에 대응하는 고유벡터가 최적 직선의 방향을 결정하며, 가장 작은 고윳값의 고유벡터는 그 직선의 법선 방향을 결정한다. 따라서 PCA 기반 피팅은 수직선을 포함한 모든 방향의 직선을 자연스럽게 처리한다.  https://kipl.tistory.com/211.

 

다만 TLS 역시 거리의 제곱합을 최소화하므로 outlier에 본질적으로 강인한 방법은 아니다. outlier의 영향을 줄이려면 Huber 또는 Tukey 가중치를 이용한 반복 가중 최소자승법, 혹은 RANSAC과 같은 robust fitting 방법을 사용해야 한다. 또한 통계에서 Deming regression은 일반적으로 $x$, $y$ 측정오차의 분산비를 반영하는 직교회귀를 가리킨다. 여기서 다룬 TLS는 두 축의 오차 분산이 같다고 보는 Deming regression의 특수한 경우로 해석할 수 있다.

PCA Line Fitting