사원수와 회전(Quaternion and Rotation)

[경고] 아래 글을 읽지 않고 "사원수와 회전"을 보면 바보로 느껴질 수 있습니다.

1. 사원수

2. 좌표계 기반 벡터

3. 회전 행렬

4. 오일러 각

사원수(四元數, quaternion)는 물리학에 벡터(vector)란 개념을 선물해준 고마운 존재이다. 사원수가 아름답기 때문에 그동안 수많은 찬사를 받았지만, 사원수는 치명적인 약점이 존재한다. 우리 직관을 너무 벗어난 사원수의 대수적 특성은 사실 멀리하고 싶은 그리움이다. 우리만 이런 딜레마를 느낄까? 당연히 아니다. 사원수가 널리 퍼진 19세기말부터 대수 기반의 사원수가 아닌 직관적인 벡터 개념을 만들기 위한 경쟁이 시작되었다. 사원수에 벡터가 이미 포함되어 있었지만, 완벽한 사원수 대수를 버리고 어딘지 부실하게 벡터만 강조한 벡터 해석학(vector analysis)이 1881년기브스 42세, 조선 고종 시절에 출현했다[1], [2]. 벡터 해석학은 미국 최초의 공학 박사이자 예일대학교(Yale University) 교수인 기브스Josiah Willard Gibbs(1839–1903)가 만들었다. 사원수라는 엄밀한 수학 체계를 어려워하는 예일대 학생들을 위해 벡터 개념을 간단히 사용할 수 있도록 기브스 교수는 좌표계 기반 벡터에 대한 자체 교재를 만들어서 학생들을 가르쳤다. 좌표계 기반 벡터 교재는 영국에 있는 헤비사이드Oliver Heaviside(1850–1925)에게 1888년 무렵에 전달되었다. 깐깐한 헤비사이드가 기브스의 벡터 개념을 칭찬했지만, 헤비사이드는 이미 1884년헤비사이드 34세, 조선 고종 시절에 기브스와는 독립적으로 사원수로 기술된 맥스웰 방정식(Maxwell's equations)을 자신만의 벡터 기반 맥스웰 방정식으로 간략화했다. 그뒤 기브스는 너무 바빠서 새로운 벡터 개념을 다듬을 시간이 없었지만, 기브스의 제자인 윌슨이 벡터 해석학[3]이란 멋진 책을 써서 1901년에 출판했다. 이 교재로 인해 사원수라는 개념은 물리학자의 손을 떠나 원래 있어야 할 수학자에게 돌아갔다. 요즘 물리학자는 기브스와 헤비사이드가 제안한 좌표계 기반 벡터를 사용해서 사물의 움직임을 계산한다.

[그림 1] 회전축 $\hat e$에 대한 3차원 공간의 회전(출처: wikipedia.org)

[그림 2] 공간 회전의 사원수 표현식을 위한 3차원 좌표계(출처: wikipedia.org)

사원수는 수학자가 발견한 교환 법칙이 성립하지 않는 최초의 대수 체계라서 존재 가치가 분명히 있다. 하지만 사원수와 경쟁하는 벡터 개념이 너무 직관적이라서 사원수는 다수의 사랑을 다시 받기는 어렵다. 그럼에도 불구하고 3차원 회전 연산(rotation operation)만 보면, 사원수의 회전 표현식이 벡터나 행렬 공식보다 확실히 아름답다.

[3차원 공간 회전을 위한 사원수 표현식]

단위 벡터 $\hat k$를 회전축으로 놓고, 3차원 벡터 $\bar v$를 $\theta$만큼 회전시킨 벡터 $\bar v_\text{rot}$는 다음처럼 표현된다.

                  (1)

여기서 단위 벡터(unit vector)는 스칼라가 $0$인 단위 사원수(unit quaternion)이다.

[증명]

사원수 표현식을 증명하기 위한 사원수의 벡터 항등식은 다음과 같다.

                              (2)

                              (3)

회전을 표현하기 위해 사용한 사원수 $\bf q$의 크기는 $\theta$에 관계없이 항상 $1$이다.

                              (4)

그래서 식 (1)처럼 ${\bf q}^{-1}$ = ${\bf q}^*$이 성립한다. 벡터 $\bar v_\text{rot}$의 크기도 $|{\bf q} \bar v{\bf q}^*|$ = $|{\bf q}| |\bar v| |{\bf q}^*|$ = $|\bar v|$처럼 보존된다. 사원수 $\bf q$를 식 (1)에 대입한 후, 식 (2)와 (3)을 이용해 정리한다.

                              (5)

여기서 삼각 함수의 합차 공식에 의해 $\cos^2 (\theta/2) - \sin^2 (\theta/2)$ = $\cos \theta$, $2 \sin^2 (\theta/2)$ = $1 - \cos \theta$이다. 식 (5)는 로드리그의 회전 공식(Rodrigues' rotation formula)과 동일하므로, 3차원 공간의 회전 표현식이 증명된다.

______________________________

공간 회전에 대한 사원수 표현식을 증명할 때, 로드리그의 회전 공식과 비교한 부분이 약간 어색해 보일 수 있다. 하지만 역사적으로 보면, 로드리그의 회전 공식이 나온 직후에 사원수가 제안되었으므로 우리의 접근 방식은 타당하다. 

[참고문헌]

[1] M. J. Crowe, "A history of vector analysis," University of Louisville, 2002. (방문일 2020-09-18)

[2] E. B. Wilson, "Reminiscences of Gibbs by a student and colleague," Bull. Amer. Math. Soc., vol. 37, no. 6, pp. 401–416, 1931.

[3] E. B. Wilson and J. W. Gibbs, Vector Analysis: A Text-book for the Use of Students of Mathematics & Physics: Founded Upon the Lectures of J. W. Gibbs, Scribner, 1901.

블록 행렬(Block Matrix)

[경고] 아래 글을 읽지 않고 "블록 행렬"을 보면 바보로 느껴질 수 있습니다.

1. 행렬

[그림 1] 행렬 원소의 정의(출처: wikipedia.org)

행렬(matrix)의 원소는 [그림 1]처럼 스칼라(scalar)로 정의함이 원칙이다. 하지만 행렬의 원소가 매우 많은 경우는 행렬 원소의 배치에 따라 구획(block)을 나누어서 부분 행렬(submatrix)로 정의하면 편리하다. 행렬의 원소를 구획으로 구분해 부분 행렬로 만든 후, 부분 행렬의 배치로 전체 행렬을 구성하는 행렬을 블록 행렬(block matrix)이라 한다. 블록 행렬의 원소는 행렬, 벡터(vector), 스칼라 등이 될 수 있다. 예를 들어, 행렬 $\bf M$을 2행 2열로 구획한 블록 행렬로 표현할 수 있다.

                  (1)

여기서 ${\bf A}, {\bf B},{\bf C},{\bf D}$는 $\bf M$의 부분 행렬이다. 혹은 $m \times n$ 형태로 구획한 블록 행렬도 다음처럼 정의할 수 있다.

                  (2)

여기서 ${\bf M}_{ij}$는 $\bf M$의 부분 행렬이다. 블록 행렬의 연산은 일반 행렬의 연산과 거의 유사하다. 다만 블록 행렬의 원소는 일반적으로 행렬이기 때문에, 각 원소를 계산할 때 교환 법칙을 쓰면 안 된다. 예를 들면, 가우스–요르단 소거법(Gauss–Jordan elimination)을 이용해 식 (1)의 역행렬을 블록 행렬 형태로 유도할 수 있다.

             (3)

따라서 식 (1)의 역행렬을 다음처럼 공식화할 수 있다.

             (4)

부분 행렬이 대각선 위치에만 있으면 블록 대각 행렬(block diagonal matrix)이라 한다.

                  (5)

역행렬 정의 ${\bf A A}^{-1}$ = $\bf I$를 이용하면 블록 대각 행렬의 역행렬은 다음과 같다.

                  (6)

블록 대각 행렬의 행렬식은 각 부분 행렬에 대한 행렬식의 곱이다.

                  (7)

식 (7)은 라플라스 전개(Laplace expansion)를 이용해 쉽게 증명할 수 있다. 블록 대각 행렬의 대각합(trace)도 각 부분 행렬에 대한 대각합의 합으로 표현할 수 있다. 

                  (8)

부분 행렬을 삼각 행렬 형태로 배치하면 블록 삼각 행렬(block triangular matrix)이 된다. 블록 삼각 행렬은 블록 상삼각 행렬(block upper triangular matrix)과 블록 하삼각 행렬(block lower triangular matrix)로 나눈다.

                  (9)

                  (10)

행렬 $\bf A$가 순열 행렬(permutation matrix) $\bf P$에 의해 다음처럼 블록 상삼각 행렬이 된다면, $\bf A$를 가약 행렬(可約行列, reducible matrix)이라 부른다.

                  (11)

여기서 ${\bf B},{\bf C},{\bf D}$는 $\widetilde{\bf A}$의 부분 행렬, ${\bf B}$와 ${\bf D}$는 정방 행렬(square matrix)이다. 행렬 $\bf A$의 양쪽에 순열 행렬을 곱하면, 순열(permutation)에 따라 행과 열이 동시에 바꾸어진다. 이로 인해 행이 바뀌면 열도 같은 규칙으로 바뀌어서 $\widetilde{a}_{ij}$ = $a_{\sigma_i \sigma_j}$가 성립한다. 여기서 $\sigma$는 $\bf P$를 구성하는 순열이다. 가약 행렬이 될 수 없으면 기약 행렬(旣約行列, irreducible matrix)이라 한다. 기약 행렬 개념은 프로베니우스Ferdinand Georg Frobenius(1849–1917)가 1912년프로베니우스 63세, 일제 식민지 시절에 제안했다. 가약 행렬의 개념을 이해하기 위해 변형된 연립 방정식 $\widetilde{\bf A}{\bf x}$ = $\bf k$를 부분 행렬 관점으로 써본다[1].

                  (11)

여기서 $\bf x$ = $[{\bf x}_1~{\bf x}_2]^T$, $\bf k$ = $[{\bf k}_1~{\bf k}_2]^T$이다. 가약 행렬의 조건에 의해 ${\bf B}$와 ${\bf D}$는 정방 행렬이므로, 식 (11)의 오른쪽 식은 잘 정의된다. 행렬 $\bf A$가 축약 가능[혹은 가약]이면, 처음에는 연립된 ${\bf x}_1$과 ${\bf x}_2$를 순열 행렬로 분리해서 더 풀기 쉬운 두 개의 행렬 방정식으로 간략화할 수 있다. 연립 방정식을 축약 혹은 간략화하지 못해서 ${\bf x}_1$과 ${\bf x}_2$를 있는 그대로 풀어야 한다면, ${\bf A}$는 기약이라 한다. 식 (11)을 보면 가약 행렬의 의미를 충분히 이해할 수 있다. 하지만 왜 굳이 프로베니우스는 큰 특징이 보이지 않는 기약 행렬을 정의했을까? 프로베니우스가 제안한 기약 행렬은 역행렬의 존재성을 판정하는 좋은 기준이다. 물론 행렬식을 계산하면 역행렬이 있을지 판정할 수 있지만 노력이 너무 많이 들어간다. 그래서 행렬이 대각 우세(diagonal dominance)인지 혹은 기약인지를 관찰해서 역행렬의 존재성을 판단한다.

[참고문헌]

[1] R. S. Varga, Matrix Iterative Analysis, 2nd ed., Berlin: Springer-Verlag, 1991.

행렬의 반복법(Iterative Method of Matrix)

[경고] 아래 글을 읽지 않고 "행렬의 반복법"을 보면 바보로 느껴질 수 있습니다.

1. 행렬의 고유치와 고유 벡터

2. 행렬 노름과 조건수

[그림 1] 고유 벡터를 구하기 위한 반복법(출처: wikipedia.org)

행렬(matrix) 이론은 잘 개발되어 있기 때문에 원하는 답을 언제나 정확하게 얻을 수 있다. 하지만 행렬의 차원이 커지면 계산 시간이 매우 빠르게 증가해서, 행렬 알고리즘(matrix algorithm or algorism)은 있으나 없는 것과 마찬가지인 상황이 빈번하게 생긴다. 그래서 매우 큰 행렬을 연산할 때는 정확함보다는 주어진 시간에 근사해를 얻는 방법이 필요하다. 행렬은 현실 문제를 푸는 합리적인 방법론이라서, 큰 차원의 행렬을 연산하기 위한 효율적인 반복법(iterative method)이 잘 발달되어 있다.

반복법으로 역행렬(inverse matrix)을 간단히 계산할 수 있. 원칙적으로 역행렬은 가우스–요르단 소거법(Gauss–Jordan elimination)으로 정확히 계산하지만, 행렬의 차원이 커지면 어마어마한 계산 시간으로 인해 반복법을 쓸 수밖에 없다. 역행렬의 반복법을 위해 행렬 $\bf T$에 대한 부분 합 ${\bf S}_m$을 구한다.

                  (1)

만약 $m$이 커질 때 ${\bf T}^m$이 $0$으로 간다면, 식 (1)을 이용해 역행렬을 다음처럼 표현할 수 있다.

                  (2)

따라서 어떤 행렬 $\bf A$의 역행렬은 식 (2)에 ${\bf T}$ = ${\bf I} - {\bf A}$를 대입해 계산할 수 있다. 식 (2)와 같은 방법을 쓰면, 다양한 역연산(inverse operator)을 무한 급수(infinite series)로 정의할 수 있다. 이 경우, 식 (2)에 나타나는 무한 급수를 노이만 급수(Neumann series)라 한다. 노이만 급수는 1877년노이만 45세, 조선 고종 시절에 노이만Carl Neumann(1832–1925)이 제안하였다. 이외에도 노이만은 노이만 경계 조건(Neumann boundary condition)으로도 유명하다. 노이만의 아버지Franz Ernst Neumann(1798–1895)는 자기장(magnetic field)에 나오는 자기 벡터 포텐셜(magnetic vector potential)과 노이만 공식(Neumann formula)을 제안하였다.

식 (2)에 등장한 노이만 급수의 수렴성을 판정한다. 임의의 열 벡터 $\bf x$를 곱해도 식 (2)의 우변에 있는 무한 급수는 수렴한다고 가정한다. 이 조건을 이용해 열 벡터 $\bf x$를 $n$개의 고유 벡터(eigenvector) ${\bf x}_i$에 대한 선형 결합으로 바꾸어서 노이만 급수에 곱한다.

                  (3)

여기서 $\lambda_i$는 고유 벡터 ${\bf x}_i$에 대한 고유치(eigenvalue), $\bf T$의 차원은 $n \times n$이다. 식 (3)의 마지막식이 모든 경우에 수렴하려면, 최대 고유치의 크기가 $1$보다 작아야 한다. 따라서 노이만 급수의 수렴 조건은 행렬 $\bf T$에 대해 $\Lambda_\text{max} < 1$이다. 여기서 $\Lambda_\text{max}$는 크기가 가장 큰 고유치[= $\max(|\lambda_i|)$]를 의미한다.

보통 가우스 소거법(Gaussian elimination)으로 푸는 연립 방정식 ${\bf Ax}$ = $\bf b$도 반복법으로 해결할 수 있다. 행렬 $\bf A$를 나누어서 $\bf A$ = ${\bf S} - {\bf T}$라고 둔다. 그러면 연립 방정식의 해인 열 벡터 $\bf x$를 다음과 같은 반복법으로 표현할 수 있다[1].

                  (4)

여기서 $k$가 커지면 ${\bf x}_k$는 $\bf x$에 수렴한다고 가정한다. 식 (4)의 왼쪽과 오른쪽 등식을 서로 빼서 반복법의 오차 벡터 ${\bf e}_k$에 대한 관계식을 만든다.

                  (5)

여기서 ${\bf e}_k$ = ${\bf x} - {\bf x}_k$이다. 다음 단계로 식 (3)에 있는 위쪽 등식처럼 오차 벡터를 고유 벡터의 선형 결합으로 바꾼다.

                  (6)

여기서 $\lambda_i$는 ${\bf S}^{-1}{\bf T}$의 제$i$번째 고유치, ${\bf S}^{-1}{\bf T}$의 차원은 $n \times n$이다. 따라서 모든 고유치가 $|\lambda_i| < 1$을 만족하면, 식 (4)에서 얻은 반복법의 결과는 해 $\bf x$에 수렴한다. 혹은 더 쉽게 크기가 가장 큰 $\Lambda_\text{max}$가 $1$보다 작으면 반복법은 수렴한다. 또한 반복법의 수렴 속도는 최대 크기의 고유치 $\Lambda_\text{max}$에 달려있다. 그래서 로랑 급수(Laurent series)의 수렴 반경(radius of convergence)처럼 $\Lambda_\text{max}$를 스펙트럼 반경(spectral radius)이라 부른다.

                  (7)

수렴 조건을 표현할 때, 주파수 범위를 뜻하는 스펙트럼이 쓰여서 약간 이상하다. 하지만 주파수에 대한 푸리에 급수(Fourier series)처럼 고유 벡터를 이용해 임의의 열 벡터를 분해한 후 반복법의 수렴 특성을 분석하므로, 스펙트럼 반경은 나름 타당한 용어이다. 또한, 스펙트럼 반경은 노이만 급수의 항이 무한대로 갈 때의 점근적 수렴성을 나타내는 특성으로 인해 점근 수렴 인자(asymptotic convergence factor)라고 부르기도 한다. 

식 (4)에 기반한 반복법을 사용하려면, $\bf S$의 역행렬을 구해야 한다. 큰 행렬의 역행렬은 또 다른 문제점을 만들기 때문에, $\bf S$의 선택은 반복법의 효율에 결정적인 영향을 끼친다. 그래서 효율적으로 $\bf S$의 역행렬을 구하기 위해, $\bf S$는 보통 대각 행렬(diagonal matrix)이나 삼각 행렬(triangular matrix)로 택한다. 이로 인해 $\bf S$를 반복법을 위한 전조건자(preconditioner)라 부른다. 행렬 $\bf S$를 대각 행렬로 간편하게 선택하는 반복법은 야코비 방법(Jacobi method)이라 한다. 야코비 방법을 쓰기 위해 ${\bf A}, {\bf S}, {\bf T}$를 다음과 같이 정한다.

                  (8a)

                  (8b)

식 (8)을 식 (4)에 대입해서 정리하면, 야코비 방법의 반복 공식을 얻는다.

             (9a)

                  (9b)

여기서 $\bf S$ = $\bf D$, $\bf T$ = $-({\bf L} + {\bf U})$, $\bf A$ = ${\bf S} - {\bf T}$ = ${\bf D}+{\bf L}+{\bf U}$, 그리고 ${\bf D}, {\bf L}, {\bf U}$는 각각 대각, 하삼각, 상삼각 행렬(diagonal, and lower and upper triangular matrices)이다. 다음처럼 $\bf S$를 하삼각 행렬(lower triangular matrix)로 선택한 반복법은 가우스–자이델 방법(Gauss–Seidel method)이 된다.

                  (10)

식 (10)을 식 (4)에 넣어서 $x_i^{(k+1)}$을 중심으로 모으면 가우스–자이델 방법을 위한 공식이 만들어진다.

                  (11a)

                  (11b)

여기서 $\bf S$ = ${\bf D} + {\bf L}$, $\bf T$ = $-{\bf U}$, $\bf A$ = ${\bf S} - {\bf T}$ = ${\bf D}+{\bf L}+{\bf U}$이다. 야코비 방법과 가우스–자이델 방법을 가중치 $\omega$로 합쳐서 수렴 속도가 빠른 새로운 반복법을 만들 수 있다. 이 반복법은 축차 가속 완화(逐次加速緩和, successive over-relaxation, SOR) 방법 혹은 SOR 방법이라 부른다[5]. 축차는 다음 순서라는 뜻이며, 축차 가속은 다음 순서에도 동일한 완화값을 써서 반복법이 가속된다는 의미이다. 수학에서 완화(緩和, relaxation)는 엄격하지 않은 느슨한(relax) 방법으로 답을 찾는 전략, 혹은 풀기 쉬운 문제로 바꾸어서 답을 대충 근사하는 방식을 뜻한다. 완화를 계속 반복하면 대부분 답에 수렴하지만, 발산할 수도 있어서 완화를 적용할 때는 수렴성을 꼭 확인해야 한다. 최적화(optimization) 기법 중에서 공간 사상(space mapping)은 완화에 기반을 두고 있다.

[그림 2] 완화 인자 $\omega$ 변화에 대한 스펙트럼 반경 $\rho$의 특성(출처: wikipedia.org)

먼저 SOR 방법을 적용하기 위해 행렬 $\bf A$를 다음처럼 $\bf A$ = ${\bf D} + {\bf L} + {\bf U}$로 나눈다.

                  (12)

여기서 $\bf D$ = $\operatorname{diag}(a_{11}, a_{22}, \cdots, a_{nn})$, ${\bf L}$과 ${\bf U}$는 LU 분해와 전혀 관계없이 $\bf A$의 좌하단 및 우상단 원소로 각각 구성한다. 식 (12)에 나온 ${\bf L}$과 ${\bf U}$는 각각 엄격 하삼각 및 상삼각 행렬(strictly lower and upper triangular matrices)이 된다. SOR 방법은 하삼각 행렬을 쓰는 가우스–자이델 방법이 중심이며, 가중치 $\omega$를 가우스–자이델 방법에 곱한다. 추가적으로 $x_i^{(k+1)}$을 구할 때, 야코비 방법적 요소를 추가하기 위해 대각 행렬에 가중치 $1-\omega$를 곱해 $x_i^{(k+1)} - x_i^{(k)}$를 수정 항으로 사용한다.

                  (13)

여기서 가중치 $\omega$는 완화 인자(relaxation factor)라 부른다. 완화 인자 $\omega$가 $1$이면, SOR 방법은 가우스–자이델 방법과 동일하다. 완화 인자가 $\omega < 1$ 및 $\omega > 1$인 값을 가지는 경우는 각각 감속 완화(減速緩和, under-relaxation)와 가속 완화(加速緩和, over-relaxation)라 이름 붙인다[6]. 완화 인자 $\omega$는 보통 $0 < \omega < 2$ 사이에서 있으며,[과감하게 $\omega$를 복소수로 확장하기도 한다.] 통상적으로 [그림 2]처럼 $1 < \omega < 2$의 범위에서 스펙트럼 반경이 최소가 되도록 $\omega$ = $\omega_\text{opt}$를 선택한다. 여기서 $\omega_\text{opt}$는 최적 완화 인자(optimal relaxation factor)이다. 식 (13)을 $x_i^{(k+1)}$에 대한 공식으로 표현하면 다음과 같다.

                  (14)

식 (13), (14)에서 $\omega$가 $0$ 가까이 갈 때의 성질은 무엇일까? 우리가 원하는 해답을 찾기 위해 노이만 급수인 식 (2)를 식 (13)에 적용해서 정리한다.

                  (15)

우세 항 $\omega$만 남기고 고차 항을 모두 없애서 더욱 간략화한다.

                  (16)

반복을 통해서 식 (16)이 수렴하게 되는 해는 야코비 방법에 의한 식 (9)가 된다. 

                  (17)

따라서 $\omega \to 0$으로 간 경우에 SOR 방법은 야코비 방법에 수렴한다.

[참고문헌]

[1] G. Strang, Linear Algebra and its Applications, 4th ed., Brooks/Cole, 2006.

[2] D. M. Young, Iterative Solution of Large Linear Systems, New York: Dover Publications, 1971.

[3] R. S. Varga, Matrix Iterative Analysis, 2nd ed., Berlin: Springer-Verlag, 1991.

[4] W. Auzinger and J. M. Melenk, Iterative Solution of Large Linear Systems, 2017. (방문일 2020-09-04)

[5] D. M. Young, Iterative Methods for Solving Partial Difference Equations of Elliptic Type, Ph.D. Thesis, Harvard University, 1950. (방문일 2020-09-04)

[6] P. Wild and W. Niethammer, "Over and underrelaxation for linear systems with weakly cyclic Jacobi matrices of index $p$," Linear Algebra Appl., vol. 91, pp. 29–52, Jun. 1987.

[7] A. Hadjidimos, "Successive overrelaxation (SOR) and related methods," J. Comput. Appl. Math., vol. 123, no. 1–2, pp. 177–199, Nov. 2000.