코딩 더 매트릭스 - 10장 직교화

필립 클라인의 저서 코딩 더 매트릭스 10장 직교화

1. 여러개의 벡터에 직교하는 투영을 찾고 이 때 여러 개의 벡터가 반드시 직교해야 하는 이유를 이해합니다.
2. 임의의 벡터 리스트를 바탕으로 서로 직교하는 생성자들을 찾는 법을 이해합니다.
3. 직교여공간을 공부하고 직합과의 관계를 파악합니다.
4. QR 인수분해를 이해하고 응용법을 공부합니다.
5. 해의 근사를 구하는 최소제곱법을 이용해 문제를 해결하는 방법을 이해합니다.

복수의 벡터들에 직교하는 투영

벡터 $v$가 벡터들의 집합 $S$내의 모든 벡터에 직교하면 $S$에 직교합니다. 예를 들어 벡터 $v = [2,0,-1]$이고 $S = Span\{[0,1,0], [1,0,2]\}$이라 하겠습니다. $\langle v, [0,1,0] \rangle = 0$, $\langle v, [1,0,2] \rangle = 0$이고 이것은 $v$가 $S$내의 모든 벡터와 직교함을 의미합니다. $S$내의 임의의 벡터를 $\alpha [0,1,0] + \beta [1,0,2]$라 할때 다음이 성립합니다.

정리하면 벡터 $v$가 벡터 $a_1, ..., a_n$에 직교할 필요충분조건은 $v$가 $Span\{a_1, ..., a_n\}$ 내의 모든 벡터에 직교하는 것입니다. $v$가 $a_1, ..., a_n$에 직교한다고 가정하고 $w$는 $Span\{a_1, ..., a_n\}$내의 임의의 벡터라 하겠습니다. $w$는 다음과 같이 쓸 수 있습니다.

$\langle v, w \rangle$에 직교의 성질을 적용하면,

따라서 $v$는 $w$에 직교합니다. 거꾸로 $v$가 $Span\{a_1, ..., a_n\}$내의 모든 벡터에 직교한다고 할 때 $Span\{a_1, ..., a_n\}$는 $a_1, ..., a_n$을 포함하므로 $v$는 $a_1, ..., a_n$에 모두 직교합니다.

서로 직교하는 벡터들의 리스트에 직교하는 b의 투영

내적을 공부하면서 벡터 $v$에 따른 $b$의 투영을 구하는 프로시저 project_orthogonal을 만들었습니다. 이번에는 벡터들의 리스트 $vlist$에 직교하는 벡터를 찾는 프로시저를 만들어보겠습니다. 한가지 주의해야 할 것은 $vlist$가 모두 직교하는 벡터들의 리스트여야 한다는 것입니다. 아래의 예시를 통해 그 이유를 살펴보겠습니다.

서로 직교하는 벡터들의 리스트를 $[v_1, v_2]$라 하겠습니다. $v_1$을 따른 $b$의 직교하는 투영은 $b^{\perp v_1} = b - \sigma_1 v_1$이라 쓸 수 있습니다. $v_2$를 따른 $b^{\perp v_1}$의 직교하는 투영은 다음과 같습니다.

$\langle b - \sigma_1 v_1 - \sigma_2 v_2, v_1 \rangle$, $\langle b - \sigma_1 v_1 - \sigma_2 v_2, v_2 \rangle$이 실제 0인지 검산해보겠습니다.

벡터들의 리스트 $[v_1, v_2]$가 서로 직교하지 않는다면 위 식은 성립하지 않았을 것입니다. 벡터들이 리스트가 서로 직교한다면 $[v_1, ..., v_n]$에 대해 직교하는 벡터는 다음 프로시저로 찾을 수 있습니다.

def project_orthogonal(b, vlist):
    for v in vlist:
        b = b - project_along(b, v)
    return b

예를 들어 $b = [1,1,1]$, $vlist = [v_1 = [0,2,2], v_2 = [0,1,-1]]$라 할때 $vlist$에 직교하는 벡터를 구하는 과정을 따라가 보겠습니다. 각 이터레이션에서 얻어지는 b를 $b_i$라 하겠습니다. 먼저 $v_1$을 따르는 $b$의 직교투영은 다음과 같습니다.

다음 이터레이션에서 $v_2$에 직교하는 투영을 찾습니다.

한 이터레이션을 $i$라 하면 다음이 성립합니다.

• $b_i$ 는 $vlist$에 속하는 $i$개 벡터들에 직교합니다.

• $b - b_i$ 는 $vlist$ 에 속하는 $i$개 벡터들의 생성에 속합니다.

$b - b_i = \sigma_1 v_1 + ... + \sigma_i v_i$는 $b = b_i + \sigma_1 v_1 + ... + \sigma_i v_i$로 바꿔쓸 수 있고 선형방정식으로 표현된 것이므로 아래와 같이 행렬로 표현 가능합니다. 이때 $i$는 파이썬의 인덱스에 맞춰 0부터 시작하도록 변경하겠습니다.

행렬-벡터 곱셈에서 선형방정식의 계수들인 $\sigma_i$를 얻는 새로운 프로시저 aug_project_orthogonal를 작성해 보겠습니다.

def aug_project_orthogonal(b, vlist):
    sigmadict = { len(vlist) : 1 }
    for i, v in enumerate(vlist):
        sigma = (b * v) / (v * v) if v * v > 1e-20 else 0
        sigmadict[i] = sigma
        b = b - sigma * v
    return (b, sigmadict)

생성자들의 직교집합 만들기

프로시저 project_orthogonal는 인수로 서로 직교하는 vlist를 받아야 했습니다. 임의의 벡터들 $v_1, ..., v_n$을 바탕으로 서로 직교하는 $Span\{v_1, ..., v_n\}$의 생성자들은 어떻게 찾을 수 있을까요?

프로시저 orthogonalize는 벡터 리스트 $v_1, ..., v_n$을 받아 서로 직교하고 $Span\{v_1, ..., v_n\} = Span\{v_1^*, ..., v_n^* \}$을 만족하는 ${v_1}^*, ..., {v_n}^*$을 반환합니다.

def orthogonalize(vlist):
    vstarlist = []
    for v in vlist:
        vstarlist.append(project_orthogonal(v, vstarlist))
    return vstarlist

orthogonalize로부터 얻은 $v_1^*, ..., v_n^*$이 $Span\{v_1, ..., v_n\} = Span\{v_1^*, ..., v_n^* \}$을 만족하는지 증명해 보겠습니다. 프로시저의 이터레이션 과정을 따라가면 우선 $Span\{v_1\} = Span\{v_1^* \}$입니다. $v_2^* = v_2 - \sigma_1 v_1^*$ 이며 $Span\{v_1^*, v_2^* \} = Span\{v_1^*, v_2 - \sigma_1 v_1^*\} = Span\{v_1, v_2\}$이 선형결합에 의해 성립합니다. $Span\{v_1^*, v_2^*, v_3^* \} = Span\{v_1^*, v_2 - \sigma_1 v_1^*, v_3 - \sigma_{22} v_2^* - \sigma_{21} v_1^*\} = Span\{v_1, v_2, v_3\}$ 역시 선형결합에 의해 성립합니다.

orthogonalize를 행렬방정식으로 아래와 같이 표현할 수 있습니다.

orthogonalize 의 사용

서로 직교하는 영이 아닌 벡터들은 일차독립입니다. 일차독립임을 증명하기 위해 다음과 같은 식을 살펴보겠습니다.

일차독립이라면 위 식을 만족하는 계수 $\alpha_1, ..., \alpha_n$이 모두 0이어야 합니다. 식의 양변에 $v_1^*$과의 내적을 취해보겠습니다.

$v_1^*$은 영벡터가 아니므로 $\parallel v_1^* \parallel$은 0이 아닙니다. 따라서 $\alpha_1$는 0입니다. 마찬가지로 $v_2^*, ..., v_n^*$에 대해 내적을 취하면 계수 $\alpha_1, ..., \alpha_n$은 모두 0이어야 식을 만족합니다. 따라서 서로 직교하는 벡터들 $v_1^*, ..., v_n^*$은 일차독립입니다.

orthogonalize를 이용해 벡터 리스트들이 생성하는 공간의 기저를 찾을 수 있습니다. orthogonalize의 인수인 벡터리스트 $[v_1, ..., v_n]$과 반환값인 $[v_1^*, ..., v_n^*]$에 대해 $Span\{[v_1, ..., v_n]\} = Span\{[v_1^*, ..., v_n^*]\}$ 이므로 아래와 같이 한번의 조건을 추가하면 기저를 찾을 수 있습니다.

from vecutil import zero_vec

def find_basis(vlist):
    vstarlist = orthogonalize(vlist)
    return [ v for v in vstarlist if v != zero_vec(v.D) ]

기저를 찾을 수 있으므로 기저의 개수인 랭크도 구할 수 있습니다.

앞서 살펴봤던 행렬방정식,

을 만족하는 프로시저 aug_orthogonalize을 작성해보겠습니다. aug_orthogonalize은 계수 $\sigma_i$ 벡터를 같이 반환하는 프로시저 aug_project_orthogonal를 이용합니다.

def aug_orthogonalize(vlist):
    vstarlist = []
    sigvecs = []
    for v in vlist:
        vstar, sigdict = aug_project_orthogonal(v, vstarlist)
        vstarlist.append(vstar)
        sigvecs.append(sigdict)
    return (vstarlist, sigvecs)

직교여공간

벡터공간 $W$와 그의 부분공간 $U$에 대해 다음을 $U$의 $W$에 대한 직교여공간이라 합니다.

$V$가 $U$의 $W$에 대한 직교여공간이라고 할때 $U \cap V$는 영벡터입니다. $V$와 $U$에 모두 속하는 임의의 벡터 $u$가 있다고 할때 $\langle u, u \rangle = 0$이고 이는 $u$의 norm도 0임을 뜻하므로 $u$는 영벡터입니다. 이는 직합의 정의와 같습니다. 따라서 $U$의 $W$에 대한 직교여공간이 $V$이면 아래와 같이 쓸 수 있습니다.

평면의 법선

법선(normal)은 수직이라는 의미입니다. 평면이 두 개의 3-벡터 $u_1$과 $u_2$의 생성으로 명시된다고 해보겠습니다. 이 경우 평면은 2차원의 벡터공간 $U$입니다. $n$은 $U$에 직교하는 영이 아닌 벡터라고 하면 $Span\{n\}$은 $\Bbb{R}^3$에 속하는 $U$의 직교여공간의 부분공간입니다. $dim\text{ }U + dim\text{ }Span\{n\} = dim\text{ }\Bbb{R}^3$이므로 $dim\text{ }Span\{n\} = 1$이고 이는 $Span\{n\}$이 일차원 직선임을 의미합니다. 이 직선이 $U$의 법선입니다.

한편 평면을 아핀 hull로서 아래와 같이 표현할 수 있습니다.

아핀 hull은 평면의 $u_1$만큼의 평행이동이므로 $Span\{n\}$은 그대로 법선입니다.

평면을 선형방정식에 대한 해집합으로 표현하면 아래와 같습니다.

위 식에 대응하는 동차선형방정식을 아래와 같이 쓰면,

벡터 $[a,b,c]$는 해집합의 소멸자입니다. 동시에 내적이 0이므로 $[a,b,c]$를 평면의 법선의 한 후보라 할 수 있습니다.

직교여공간 계산하기

직교여공간을 어떻게 찾을 수 있을까요? 벡터공간 $W$의 기저가 $[w_1, ..., w_n]$, 부분공간 $U$의 기저가 $[u_1, ..., u_m]$일때 $[u_1, ..., u_m, w_1, ..., w_n]$을 프로시저 orthogonalize의 인수로 대입함으로서 직교여공간의 기저를 찾을 수 있습니다. orthogonalize의 이터레이션 $m$까지 얻은 서로 직교하는 벡터리스트는 $[u_1^*, ..., u_m^*]$이고 $Span\{[u_1, ..., u_m]\} = Span\{[u_1^*, ..., u_m^*]\}$이 성립합니다. 이터레이션 $m$이후 얻게되는 $[w_1^*, ..., w_n^*]$은 $[u_1^*, ..., u_m^*]$에 모두 직교하고 $Span\{[u_1, ..., u_m]\} = Span\{[u_1^*, ..., u_m^*]\}$이므로 $[w_1^*, ..., w_n^*]$는 $U$에 직교합니다. 따라서 이는 $U$의 $W$에 대한 직교여공간의 기저라 할 수 있습니다. 단, $[w_1^*, ..., w_n^*]$에서 영벡터를 제외해야 합니다. 직교여공간을 찾는 프로시저 find_orthogonal_complement를 작성해 보겠습니다.

def find_orthogonal_complement(U_basis, W_basis):
    vstarlist = orthogonalize(U_basis + W_basis)
    return [ v for v in vstarlist[len(U_basis):] if v != zero_vec(v.D) ]

예를 들어 $\Bbb{R}^3$에 속하는 $Span\{[8,-2,2], [0,3,3]\}$의 직교여공간을 찾아보겠습니다.

from vecutil import list2vec

U_basis = [
    list2vec([ 8,-2,2 ]),
    list2vec([ 0,3,3 ]),
]
W_basis = [
    list2vec([ 1,0,0]),
    list2vec([ 0,1,0 ]),
    list2vec([ 0,0,1 ]),
]
for v in find_orthogonal_complement(U_basis, W_basis):
    print(v)

>>>

     0     1      2
-------------------
 0.111 0.222 -0.222

         0        1        2
----------------------------
 -8.33E-17 1.67E-16 5.55E-17

        0        1        2
---------------------------
 8.33E-17 5.55E-17 1.67E-16

결과값의 두번째와 세번째는 부동소수점 오차를 고려한 영벡터입니다.

QR 인수분해

서로 직교하는 벡터들은 만약 그들의 norm이 모두 1이면 정규직교(orthonormal)한다고 합니다. 행렬은 만약 그 열들이 정규직교하면 열-직교(column-orthogonal)라고 하고 정방 열-직교 행렬은 직교행렬이라고 합니다.

정규직교하는 열들을 가지는 행렬 $Q$가 있다고 할때 $Q^TQ$를 살펴보겠습니다.

$i = j$일 때 $q_i^* \cdot q_j^* = 1 (\because 정규직교)$이고 $i \neq j$일 때 $q_i^* \cdot q_j^* = 0$입니다. 따라서 $Q^TQ$는 $n \times n$ 단위행렬입니다. 따라서 $Q^T$와 $Q$는 서로 역행렬입니다.

QR인수분해는 행렬을 두 행렬 $Q$와 $R$의 곱 $QR$로 표현하는 방법입니다. 여기서 행렬 $Q$와 $R$을 찾는 것은 프로시저 aug_orthogonalize로부터 얻는 두 벡터 리스트에서 시작합니다. aug_orthogonalize의 첫번째 반환값인 vstarlist는 서로 직교하는 열벡터들의 리스트 $[v_1^*, ..., v_n^*]$이고 두번째 반환값 sigvecs는 계수 $\sigma_i$와 1로 이루어진 열벡터들의 리스트입니다. 이 때 sigvecs 열벡터 리스트는 상삼각행렬을 이룹니다. vstarlist를 $Q$에, sigvecs를 $R$에 대응시키면 남은 작업은 $Q^TQ$가 단위행렬이 되도록 $Q$를 정규화하는 것입니다. 정규화는 다음 과정으로 이루어집니다.

• 열벡터 $v_1^*, ..., v_n^*$을 가지는 행렬의 열들 각각을 각각의 norm으로 나누어 줍니다.
• 삼각행렬의 각 행에 대응하는 $v_1^*, ..., v_n^*$의 norm을 곱해 줍니다.

$QR$이 행렬의 곱셈이므로 삼각행렬 $R$에서는 도트곱에 대응하는 위치에 norm을 곱해 준다고 생각하면 이해가 쉽습니다.

단 norm이 0이 되는 상황을 피하기 위해(영벡터가 없어야 함을 뜻함) $A$의 열들이 일차독립이어야 QR인수분해가 가능합니다. 일차독립 조건이 충족된다면 다음이 성립합니다.

이는 프로시저 orthogonalize의 인수 벡터 리스트와 결과 벡터 리스트의 생성에 대한 식 $Span\{v_1, ..., v_n\} = Span\{v_1^*, ..., v_n^*\}$과 같은 원리입니다.

QR 인수분해로 $Ax = b$ 풀기

QR 인수분해를 이용해 행렬 방정식 $Ax = b$를 풀어보겠습니다.

위 식의 양변에 $Q^T$를 곱하면 다음을 얻습니다.

$Q^TQ$는 단위행렬이므로 다음을 얻습니다.

이 때 $R$은 상삼각행렬이므로 후진대입법을 이용하면 $x$값을 찾을 수 있습니다.

만약 행렬 $A$의 행의 개수가 열의 개수보다 많다면 $Ax = b$의 해가 없을 수 있습니다. 이 경우 $A$의 열들의 선형결합 중 $b$에 가장 가까운 벡터를 찾음으로서 문제를 해결할 수 있습니다.. 직교화는 이 문제를 풀 수 있는 솔루션을 제공합니다. $b$에 가장 가까운 벡터는 $b^{\parallel v}$라 할 수 있고 이것은 $A$의 열공간상에 대한 $b$의 투영입니다.

한편 가장 가까운 벡터 $b^{\parallel v}$를 찾기 위해 다음 정리가 유용합니다. $Q$가 열-직교 기저라 하고 $V = Col\text{ }Q$라 할 때 정의역이 $Q$의 행라벨 집합과 동일한 임의의 벡터 $b$에 대해 $Q^Tb$는 $b^{\parallel v}$의 좌표표현이고 $QQ^Tb = b^{\parallel v}$입니다. $b^{\parallel v}$는 $V$내에 있으므로 다음과 같이 쓸 수 있습니다.

이 벡터가 $Q^Tb$임을 증명하면 정리를 만족합니다. $Q^Tb$의 한 엔트리는 다음과 같이 쓸 수 있습니다.

이는 $b^{\parallel v}$의 한 엔트리의 계수 $\alpha_j$가 $\langle q_j, b \rangle$임을 뜻합니다. 따라서 $Q^Tb$를 다음과 같이 쓸 수 있습니다.

벡터 $[\alpha_1, ..., \alpha_j]$는 $b^{\parallel v}$의 좌표표현이고 $QQ^Tb$는 $b^{\parallel v}$입니다.

$Rx = Q^Tb$에서 우변은 $b^{\parallel v}$의 좌표표현이고 양변에 $Q$를 곱하면 다음을 얻습니다.

이는 $b$에 가장 가까운 벡터 $b^{\parallel v}$를 만족하는 해 $\hat{x}$를 찾을 수 있음을 증명합니다. 이와 같이 해의 근사를 구하는 방법을 최소제곱이라 합니다.

최소제곱의 응용

최소제곱의 한 응용 예는 2차원 데이터에 가장 잘 일치하는 직선을 찾는 것입니다. 예시로 다음과 같은 표가 있습니다.

위 데이터가 선형에 가까운 모델일 때 나이($x$)에 대한 일차 함수 $f(x) = a + cx$를 뇌 중량을 가장 잘 예측하는 함수라 하겠습니다. 여기서 예측 에러의 제곱의 합을 최소화하는 $a$, $c$를 찾는 것이 목적입니다. 아래와 같이 모델을 행렬 방정식으로 정리할 수 있습니다.

실제 결과 벡터를 $b$라 할 때 위 식의 우변은 $b$에 가장 가까운 $b^{\parallel v}$라 할 수 있습니다. 최소제곱을 이용해 가설 벡터 $[a,c]$를 찾을 수 있습니다.

이외에도 결과 데이터에 가장 근접한 이차함수를 찾는 이차함수 피팅에도 최소제곱을 이용할 수 있습니다. 예를 들면 다음과 같습니다.

마찬가지로 예측 에러의 제곱의 합을 최소화하는 벡터 $[u_0, u_1, u_2]$를 찾는 것이 목적입니다.

QR 인수분해에 필요한 프로시저들

• 일차독립인 벡터 리스트를 직교화, 정규화하는 프로시저 orthonormalize

import math

def orthonormalize(L):
    Lstar = orthogonalize(L)
    norms = [ math.sqrt(v * v) for v in Lstar ]
    Lnormalized = [ v / norms[i] if norms[i] > 1e-20 else zero_vec(Lstar[i].D) for i, v in enumerate(Lstar) ]
    return Lnormalized

• 일차독립인 벡터 리스트를 QR 인수분해하는 프로시저 aug_orthonormalize

def dictlist2veclist(D, dictlist):
    veclist = []
    for d in dictlist:
        veclist.append(Vec(D, d))
    return veclist

def adjust_multipliers(v, multipliers):
    for k in v:
        v[k] *= multipliers[k]

def aug_orthonormalize(L):
    Lstar, sigvecs = aug_orthogonalize(L)
    norms = [ math.sqrt(v * v) for v in Lstar ]
    Qlist = [ v / norms[i] if norms[i] > 1e-20 else zero_vec(Lstar[i].D) for i, v in enumerate(Lstar) ]
    for v in sigvecs:
        adjust_multipliers(v, norms)
    D = set(range(len(Qlist)))
    Rlist = dictlist2veclist(D, sigvecs)
    return (Qlist, Rlist)

• QR 인수분해를 이용해 해 또는 근사해를 찾는 프로시저 QR_solve

from matutil import mat2rowdict, mat2coldict, coldict2mat

def triangular_solve(rowlist, v):
    D = rowlist[0].D
    label_list = sorted(list(D))
    u = zero_vec(D)
    for i in reversed(range(len(rowlist))):
        c = label_list[i]
        row = rowlist[i]
        u[c] = (v[i] - row * u) / row[c]
    return u

def QR_solve(A, b):
    Alist = [ v for k,v in mat2coldict(A).items() ]
    Qlist, Rlist = aug_orthonormalize(Alist)
    Q = coldict2mat(Qlist)
    R = coldict2mat(Rlist)
    x = triangular_solve(mat2rowdict(R), Q.transpose() * b)
    return x