Linear Algebra
///
///
🤿

Some Math Basics often used in Photogrammetry

생성일
2022/02/10 02:54
태그
SVD
Least-squares
Skew-symmetric
3d geometry를 공부하다보면 수학적인 지식을 통해 해결하는 경우가 많다, 미적분,선형대수 등등.. 오늘은 그 중에서도 가장 필요한 부분 행렬 분해, 해 구하기, 회전 등등에서 필요한 지식들을 정리해 보도록 하겠습니다.

Solving Ax=bAx = b

이러한 형태의 문제를 풀 때 크게 3가지의 경우로 나뉘어 진다.

case1. 행렬 AA가 sqaure matrix이고 full rank를 가지는 경우

가장 최고의 상황이며, 이러한 경우 unique solution을 가지게 된다.
사용할수 있는 방법은
1.
Gaussian Elimination
2.
Inversion of AA : x=A−1b\boldsymbol{x}=A^{-1} \boldsymbol{b}
3.
Cholesky decomposition chol⁡(A)=LL⊤\operatorname{chol}(A)=L L^{\top} 이 때, 행렬 LL은 lower triangular matrix이고, Ly=b and L⊤x=yL \boldsymbol{y}=\boldsymbol{b} \text { and } L^{\top} \boldsymbol{x}=\boldsymbol{y}를 푼다. 일반적으로 행렬 AA가 sparse한 경우 유용하다.
4.
QR decomposition
5.
Conjugate gradients

case2. 행렬 AA가 over-determined된 경우

가장 흔하게 볼 수 있는 상황이다. unique한 solution이 존재하지 않는다. 이러한 경우 곧바로 Ax=bA x=b를 풀기 보다는 ∥Ax−b∥\|A x-b\|를 최소화 하는 xx를 찾으면 된다. 식으로 표현하면 아래와 같다
x∗=arg⁡min⁡x∥Ax−b∥\boldsymbol{x}^{*}=\arg \min _{\boldsymbol{x}}\|A \boldsymbol{x}-\boldsymbol{b}\|
가장 일반적인 least squares approach이며, 해를 구하면 아래와 같이 나온다.
x=(A⊤A)−1A⊤b\boldsymbol{x}=\left(A^{\top} A\right)^{-1} A^{\top} \boldsymbol{b}

case3. 행렬 AA가 under-determined된 경우

무수히 많은 해가 존재한다. (당연하다, 미지수의 개수가 식의 개수보다 많으니까) 이러한 경우에는 Ax=bA x=b를 푸는 xx를 찾되, ∥x∥\|\boldsymbol{x}\|를 최소화하면 된다, 해는 아래와 같이 나온다.
x=A⊤(AA⊤)−1b\boldsymbol{x}=A^{\top}\left(A A^{\top}\right)^{-1} \boldsymbol{b}

Solving Ax=0Ax = \mathbf{0}

Ax = 0 인 방정식을 풀때가 많다 (DLT, Camera Calibration, fundamental/Essential M)
이러한 Homogenous System Ax=0Ax = \mathbf{0} 을 푸는 것은 x≠0x \neq 0인 non-trivial solution을 찾는 것이다. 일반적으로 해가 여러개이기 때문에 underdetermined된 경우라고 생각할 수 있으며, null⁡(A)\operatorname{null}(A)의 원소를 구하는 것과 같다.
일반적으로, sqaured matrix AA에 대해 아래 식을 만족한다.
dim⁡(A)=dim⁡(null⁡(A))+rank⁡(A)\operatorname{dim}(A)=\operatorname{dim}(\operatorname{null}(A))+\operatorname{rank}(A)
•
rank⁡(A)\operatorname{rank}(A) 개의 non-zero eigenvalue를 가지게 된다.
•
dim⁡(null⁡(A))\operatorname{dim}(\operatorname{null}(A)) 개의 zero eigenvalue를 가지게 된다.
eigenvalue의 기본 개념 Aν=λνA \nu=\lambda \nu을 다시 떠올려 본다면, eigenvalue가 0이라면 아래 식을 만족하게 된다.
Aν=0ν=0A \nu=0 \nu=0

그렇다! 결국 우리가 해야 할 것은 eigenvalue 0 에 해당하는 eigenvector를 찾는 것이다!

따라서 임의의 행렬 A∈Rm×nA \in \mathbb{R}^{m \times n}에 대해서 ATAA^{T}A 하여 symmetric하게 만든 후 EVD 하여 구한 eigenvalue 값이 AA에대한 svd의 singular value 와 연관된다는 사실을 생각하면 해를 구할 수 있게 된다. 따라서, 우린 이제부터 svd하여서 singular value, vectors를 구하면 된다.

Solution to Ax=0Ax = 0 via SVD

행렬 U,VU,V는 orthogonal matrix이고, 행렬 DD는 diagonal entries가 내림차순으로 정렬된 diagonal matrix입니다.
행렬 VTV^{T}는 앞서 다룬 행렬 DD의 각 singular values에 해당하는 singular vector를 row형태로 나열한 것이며, 이를 행렬 VV로 transpose해서 살펴보면, 행렬 VV의 마지막 column은 가장 작은 singular value에 해당하는 singular vector임을 알 수 있습니다.
따라서 행렬 VV의 마지막 column이 Ax=0A x=0을 풀기 위해 우리가 찾는 non-trivial solution이 됩니다. ( Avi=σiui=0A \mathbf{v}_{i}=\sigma_{i} \mathbf{u}_{i} = \mathbf{0} )
그렇지만, 마지막 column에 해당하는 singular value가 0이면 참 좋겠지만, 아닐 수도 있습니다. 그러나, 이 값이 ∥x∥=1\| x\| = 1인 조건 하에서 ∥Ax∥\| Ax\|를 가장 최소화 하는 값을 나타낸다는 것은 변하지 않기 때문에 이 것을 해로 여깁니다.

Least Squares Problem

least squares 방법은 observation에 근거해 파라미터나 state를 추정하는 테크닉입니다. state xx를 predicted observation에 매핑시키기 위해 state xx를 업데이트 해 나가는 과정입니다. Predicted Observation이 real observation과 가깝도록 하는게 목표입니다. 많은 분야에서 사용됩니다. (Mapping, Localization, SLAM, Calibration, Bundle Adjustment...)
Predicted measurement와 real measurement 사이의 차를 최소화하는 방향으로 식을 전개하며 Sum of the Squared Error 형태로 나타나게 됩니다.
min⁡x∑(z−f(x))2\min _{\boldsymbol{x}} \sum(\boldsymbol{z}-f(\boldsymbol{x}))^{2}
Error는 predicted 와 Actual measurement의 차이로 표현됩니다.
ei(x)=zi−fi(x)\mathbf{e}_{i}(\mathbf{x})=\mathbf{z}_{i}-f_{i}(\mathbf{x})
이 때, 이 Error는 평균이 0인 normal distribution의 분포를 가진다고 가정합니다. 이러한 분포를 표현하기 위해 information matrix Λi\boldsymbol{\Lambda}_{i}를 사용하며, 앞서 제시한 squared error of a measurement를 표현하면 아래와 같습니다.
ei(x)=ei(x)TΛiei(x)e_{i}(\mathbf{x})=\mathbf{e}_{i}(\mathbf{x})^{T} \boldsymbol{\Lambda}_{i} \mathbf{e}_{i}(\mathbf{x})
참고로, squared error of a measurement의 결과는 scalar의 형태를 띄게 됩니다.
이러한, Error Function은 linear한 형태가 아니기 때문에, 선형화 하는 과정이 필요합니다. 이 때, Taylor Expansion을 이용합니다.
ei(x+Δx)≃ei(x)⏟ei+Ji(x)Δx\mathbf{e}_{i}(\mathbf{x}+\Delta \mathbf{x}) \simeq \underbrace{\mathbf{e}_{i}(\mathbf{x})}_{\mathbf{e}_{i}}+\mathbf{J}_{i}(\mathbf{x}) \Delta \mathbf{x}
행렬 J\mathbf{J}는 Jacobian을 의미합니다.
Jf(x)=(∂f1(x)∂x1∂f1(x)∂x2…∂f1(x)∂xn∂f2(x)∂x1∂f2(x)∂x2…∂f2(x)∂xn…………∂fm(x)∂x1∂fm(x)∂x2…∂fm(x)∂xn)\mathbf{J}_{f}(x)=\left(\begin{array}{cccc}\frac{\partial f_{1}(x)}{\partial x_{1}} & \frac{\partial f_{1}(x)}{\partial x_{2}} & \ldots & \frac{\partial f_{1}(x)}{\partial x_{n}} \\\frac{\partial f_{2}(x)}{\partial x_{1}} & \frac{\partial f_{2}(x)}{\partial x_{2}} & \ldots & \frac{\partial f_{2}(x)}{\partial x_{n}} \\\ldots & \ldots & \ldots & \ldots \\\frac{\partial f_{m}(x)}{\partial x_{1}} & \frac{\partial f_{m}(x)}{\partial x_{2}} & \ldots & \frac{\partial f_{m}(x)}{\partial x_{n}}\end{array}\right)
이후 Gauss-Newton 방식에 의해 다음과 같이 반복하여 해를 구합니다.
(1) state x\mathbf{x} 근처에서 선형화 하기
ei(x+Δx)≃ei(x)+JiΔx\mathbf{e}_{i}(\mathbf{x}+\Delta \mathbf{x}) \simeq \mathbf{e}_{i}(\mathbf{x})+\mathbf{J}_{i} \Delta \mathbf{x}
(2) Linear System에 해당하는 term들 구하기
b⊤=∑iei⊤ΛiJiH=∑iJi⊤ΛiJi\mathbf{b}^{\top}=\sum_{i} \mathbf{e}_{i}^{\top} \boldsymbol{\Lambda}_{i} \mathbf{J}_{i} \quad \mathbf{H}=\sum_{i} \mathbf{J}_{i}^{\top} \boldsymbol{\Lambda}_{i} \mathbf{J}_{i}
(3) Linear System을 풀기
Δx∗=−H−1b\Delta \mathbf{x}^{*}=-\mathbf{H}^{-1} \mathbf{b}
이후 계속해서 state x\mathbf{x}를 x←x+Δx∗\mathrm{x} \leftarrow \mathrm{x}+\Delta \mathrm{x}^{*}로 업데이트 합니다.
가우스 뉴턴법 말고도 Levenberg-Marquadt 방법도 존재합니다. 이러한 non-linear least squares solving 하는 방법에 대해서는 따로 지면을 할애해 서술해 보겠습니다!

Skew-symmetric Matrix

skew-symmetric matrix는 S⊤=−SS^{\top}=-S 를 만족하는 행렬 SS 를 말합니다. 이러한 행렬은 main diagonal이 0의 값을 가지고 det⁡(S)=0\operatorname{det}(S)=0 을 만족합니다.
주로 이러한 skew-symmetric matrix를 통해 cross-product를 표현할 때 많이 사용합니다.
a×b=[a]×b=Sab\boldsymbol{a} \times \boldsymbol{b}=[\boldsymbol{a}]_{\times} \boldsymbol{b}=S_{a} \boldsymbol{b}
[a]×=Sa=[0−a3a2a30−a1−a2a10]{[\boldsymbol{a}]_{\times}}={S}_{a}=\left[\begin{array}{ccc}0 & -a_{3} & a_{2} \\a_{3} & 0 & -a_{1} \\-a_{2} & a_{1} & 0\end{array}\right]
[a1a2a3]⏟a×[b1b2b3]⏟b=[−a3b2+a2b3a3b1−a1b3−a2b1+a1b2]=[0−a3a2a30−a1−a2a10]⏟Sa[b1b2b3]⏟b\underbrace{\left[\begin{array}{l}a_{1} \\a_{2} \\a_{3}\end{array}\right]}_{\boldsymbol{a}} \times \underbrace{\left[\begin{array}{l}b_{1} \\b_{2} \\b_{3}\end{array}\right]}_{\boldsymbol{b}}=\left[\begin{array}{c}-a_{3} b_{2}+a_{2} b_{3} \\a_{3} b_{1}-a_{1} b_{3} \\-a_{2} b_{1}+a_{1} b_{2}\end{array}\right]=\underbrace{\left[\begin{array}{ccc}0 & -a_{3} & a_{2} \\a_{3} & 0 & -a_{1} \\-a_{2} & a_{1} & 0\end{array}\right]}_{S_{a}} \underbrace{\left[\begin{array}{l}b_{1} \\b_{2} \\b_{3}\end{array}\right]}_{\boldsymbol{b}}

Derivative of a Rotation Matrix

앞서 다룬 skew-symmetric matrix를 이용해 Rotation Matrix의 derivative를 구할 수 있습니다. 아래 과정은 유도하는 과정을 나타낸 것입니다.
임의의 rotation matrix RR 에 대해서, RR⊤=IR R^{\top}=I 가 성립된다는 것을 알고 있습니다.
x축을 기준으로 θ\theta 만큼 회전하는 rotation matrix Rx(θ)R_{x}(\theta)를 생각해봅니다. Rx(θ)Rx⊤(θ)=IR_{x}(\theta) R_{x}^{\top}(\theta)=I를 만족합니다.
위의 식에서 양변을 θ\theta 에 대하여 미분해보면,
Rx(θ)Rx⊤(θ)=IR_{x}(\theta) R_{x}^{\top}(\theta)=I
ddθ(Rx(θ)Rx⊤(θ))=ddθI\frac{d}{d \theta}\left(R_{x}(\theta) R_{x}^{\top}(\theta)\right)=\frac{d}{d \theta} I
곱미분 성질을 이용하여 아래와 같이 서술할 수 있습니다.
ddθRx(θ)Rx⊤(θ)+Rx(θ)ddθRx⊤(θ)=0\frac{d}{d \theta} R_{x}(\theta) R_{x}^{\top}(\theta)+R_{x}(\theta) \frac{d}{d \theta} R_{x}^{\top}(\theta)=0
(AB)⊤=B⊤A⊤(A B)^{\top}=B^{\top} A^{\top} 성질을 이용하면,
ddθRx(θ)Rx⊤(θ)+(ddθRx(θ)Rx⊤(θ))⊤=0\frac{d}{d \theta} R_{x}(\theta) R_{x}^{\top}(\theta)+\left(\frac{d}{d \theta} R_{x}(\theta) R_{x}^{\top}(\theta)\right)^{\top}=0
와 같이 표현할 수 있고, ddθRx(θ)Rx⊤(θ)=S\frac{d}{d \theta} R_{x}(\theta) R_{x}^{\top}(\theta) = S 로 치환하게 되면,
행렬 SS는 S+S⊤=0S + S^{\top}=0 를 만족하는 skew-symmetric matrix임을 알 수 있습니다.
주어진 행렬 SS를 실제 rotation matrix RxR_{x}를 이용해 구해보면,
Rx=[1000cos⁡θsin⁡θ0−sin⁡θcos⁡θ]R_{x}=\left[\begin{array}{ccc}1 & 0 & 0 \\0 & \cos \theta & \sin \theta \\0 & -\sin \theta & \cos \theta\end{array}\right]
S=[0000−sin⁡θ−cos⁡θ0cos⁡θ−sin⁡θ]⏟ddθRx(θ)[1000cos⁡θ−sin⁡θ0sin⁡θcos⁡θ]⏟Rx⊤(θ)=[00000−1010]S=\underbrace{\left[\begin{array}{ccc}0 & 0 & 0 \\0 & -\sin \theta & -\cos \theta \\0 & \cos \theta & -\sin \theta\end{array}\right]}_{\frac{d}{d \theta} R_{x}(\theta)} \underbrace{\left[\begin{array}{ccc}1 & 0 & 0 \\0 & \cos \theta & -\sin \theta \\0 & \sin \theta & \cos \theta\end{array}\right]}_{R_{x}^{\top}(\theta)} = \left[\begin{array}{ccc} 0 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{array}\right]
이는 x축에 대한 unit vector ex=[100]e_{x}=\left[\begin{array}{l}1 \\0 \\0\end{array}\right]를 이용한 skew-symmetric matrix SexS_{\boldsymbol{e}_{x}}로 표현 가능합니다.
따라서 아래와 같이 서술 할 수 있으며,
ddθRx(θ)Rx⊤(θ)=Sex\frac{d}{d \theta} R_{x}(\theta) R_{x}^{\top}(\theta)=S_{\boldsymbol{e}_{x}}
양변에 Rx(θ)R_{x}(\theta)를 곱하게 되면,
ddθRx(θ)=ddθRx(θ)Rx⊤(θ)Rx(θ)⏟I=SexRx(θ)\frac{d}{d \theta} R_{x}(\theta)=\frac{d}{d \theta} R_{x}(\theta) \underbrace{R_{x}^{\top}(\theta) R_{x}(\theta)}_{I}=S_{\boldsymbol{e}_{x}} R_{x}(\theta)
정리하면, rotation matrix Rx(θ)R_{x}(\theta)의 derivative는 skew-symmetric matrix SexS_{\boldsymbol{e}_{x}}와 rotation matrix Rx(θ)R_{x}(\theta)의 곱의 형태로 표현할 수 있습니다.
ddθRx(θ)=SexRx(θ)\frac{d}{d \theta} R_{x}(\theta)=S_{\boldsymbol{e}_{x}} R_{x}(\theta)
x축 말고도 임의의 다른 축 rr에 대해서도 적용됩니다.

Reference

Some Math Basics often used in Photogrammetry (Cyrill Stachniss, 2021)
An brief, informal, and incomplete collection of math basics and tools that are often used in Photogrammetry Slides: http://www.ipb.uni-bonn.de/html/teaching...
https://www.youtube.com/watch?v=Q042jupFMbU&t=1643s
Made with