[Linear Analysis] 특이값(Singular Value)

특이값분해(Singular Value Decomposition)

내용

• 비정방행렬의 고유값분해
• 특이값

비정방행렬의 고유값분해

고유 분해, 스펙트럴 분해 등은 가역적인 정방 행렬을 대상으로 합니다. 이에 반해 특이값을 이용하는 특이값 분해(Singular Value Decomposition)는 비정방 행렬을 정방행렬로 변형하여 분해할 수 있는 방법으로 선형 대수의 계산에서 가장 유용하게 사용되고 있습니다.

비정방행렬은 행렬에 그 행렬의 전치행렬과의 행렬곱으로 정방행렬로 만들 수 있습니다.

import numpy as np
import numpy.linalg as la
from sympy import *

np.random.seed(1)
X=np.random.randint(0, 11, (3, 4))
print(X)

[[5 8 9 5]
 [0 0 1 7]
 [6 9 2 4]]

XTX=np.dot(X.T, X)
print(XTX)

[[ 61  94  57  49]
 [ 94 145  90  76]
 [ 57  90  86  60]
 [ 49  76  60  90]]

위 결과와 같이 전치행렬과의 행렬곱으로 생성한 정방행렬은 대칭행렬이 됩니다. 이 대칭행렬의 고유행렬은 다음과 같이 전치행렬과 역행렬이 같으므로 정규직교(Orthonormal)행렬이 됩니다. 이 결과는 정방행렬이면서 대칭행렬이므로 행렬의 대각화에 기반을 둔 분해(decomposition)가 가능합니다.

d, P=la.eig(XTX)
np.allclose(P.T, la.inv(P))

True

정방행렬 A의 고유값과 고유벡터의 관계와 고유벡터는 단위벡터인 점을 적용하면 식 1이 성립합니다. 즉, 행렬 A와 고유벡터(v)와의 내적의 크기는 고유벡터에 대응하는 고유값의 크기와 같아집니다.

\begin{align}Av_1& =\lambda_1v_1, \Vert{v_1}\Vert=1\\\tag{식 1} \Vert{Av_1}\Vert& =\Vert{\lambda_1v_1}\Vert\\& = \vert{\lambda_1}\vert \Vert{v_1}\Vert\\ \\& = \vert{\lambda_1}\vert \end{align}

위와 같은 정방 행렬의 특성은 비정방 행렬을 정방행렬로 변환 후 적용할 수 있습니다.

예 1)

ℝ³ 벡터 x가 행렬 A에 의해 ℝ²로 선형변환이 이루어진다면 ||Ax||가 최대가 되는 단위벡터 x를 계산해 봅니다.

$$A=\begin{bmatrix}4& 11& 14\\8 & 7& -2 \end{bmatrix}$$

행렬 Ax의 최대크기를 결정하기 위해 고유값과 고유벡터를 사용할 수 있습니다. 그러나 A는 정방행렬이 아니므로 전치행렬과의 내적곱으로 정방행렬(A1)로 전환해야 합니다. 이 전환된 행렬 A1은 식 2와 같이 대칭행렬이 됩니다.

\begin{align}A1& =A^TA\\\tag{식 2}(A^TA)^T& =A^T(A^T)^T\\& = A^TA\\\therefore\, A1^T& =A1 \end{align}

비정방행렬인‖Ax‖의 최대값은 식 3과 같이 결정할 수 있습니다.

\begin{align}\Vert{Ax}\Vert & = (Ax)^T(Ax)\\\tag{식 3}& =x^TA^TAx\\& = x^T(A^TA)x\end{align}

위 식의 x^T(A^TA)x는 이차 형식으로서 ||x||² = 1(x^Tx = 1)의 제한 조건에서의 A^TA의 최대값은 이 대칭 행렬의 가장 큰 고유값이 됩니다(제한된 최대 최소 참조). 그 고유값에 대응하는 단위 고유 벡터(v₁)와의 곱인 ||Av₁||은 최대가 됩니다.

A=np.array([[4,11,14], [8, 7, -2]])
A1=A.T@A
print(A1)

[[ 80 100  40]
 [100 170 140]
 [ 40 140 200]]

d, P=la.eig(A1)
print(d.round(3))

[360.   0.  90.]

대칭행렬의 고유벡터의 L₂ norm은 1입니다. 즉, 단위벡터입니다.

for i in range(3):
    print("고유벡터 v%d의 norm: %.1f" %(i, la.norm(P[:,i])))

고유벡터 v0의 norm: 1.0
고유벡터 v1의 norm: 1.0
고유벡터 v2의 norm: 1.0

A^TA의 고유값 중에 최대값에 대응하는 벡터는 첫번째 고유 벡터 입니다. 즉, 다음의 고유행렬 P의 0열에 해당합니다.

Av1=A@P[:,0]
print(Av1)

[-18.  -6.]

la.norm(Av1).round(3)

18.974

특이값

예 1로부터 m×n 형태의 행렬 A에 대한 A^TA는 대칭 행렬이므로 정규직교적으로 대각화가 가능합니다. 즉, 그 대칭행렬의 고유벡터(v_i)는 정규직교 벡터이므로 ℝⁿ 차원의 정규 직교 기저가 됩니다. λ₁, λ₂, …, λ_n을 A^TA의 고유값이라고 하면 식 4가 성립합니다. (1 ≤ i ≤ n)

\begin{align}A^TAv_i& =\lambda_iv_i\\\Vert{Av_i}\Vert & = (Av_i)^T(Av_i)\\\tag{식 4} & =v_i^TA^TAv_i\\& =v_I^T(\lambda_i v_i)\\& = \lambda_i v_i^Tv_i\\& =\lambda_i \\ v_i:&\;\text{정규직교벡터} \rightarrow v_i^T=v_i^{-1} \end{align}

식 4로부터 A^TA의 고유값은 음수가 아님을 나타냅니다 (식 5).

$$\tag{식 5} λ_1 \gt λ_2 \gt … λ_n \gt 0

A^TA의 고유값(λ_i)의 제곱근을 행렬 A의 특이값(singular value)이라고 하며 식 6과 같이 σ₁, σ₂, …, σ_n으로 나타냅니다.

\begin{align}\sigma_i & =\sqrt{\lambda_i}\\\tag{식 6}& =\Vert{Av_i}\Vert\\ 1\le& i\lt n \end{align}

예 1에서 A^TA의 최대 고유값은 360이며 이 고유값의 제곱근이 행렬 A의 최대 특이값이 됩니다.

np.sqrt(d[0]).round(3)

18.974

이 값은 ||Ax||의 최대값과 같습니다. 벡터 x는 A^TA의 최대 고유값에 대응하는 고유벡터입니다.

예 2)

다음 행렬 A의 특이값을 결정합니다.

$$A=\begin{bmatrix}3& 2& 2\\2& 3& -2 \end{bmatrix}$$

A=np.array([[3,2,2],[2,3,-2]])
A1=A.T@A
print(A1)

[[13 12  2]
 [12 13 -2]
 [ 2 -2  8]]

A^TA는 대칭행렬이지만 비가역행렬입니다. 그러나 다음 결과와 같이 정규직교행렬입니다.

la.det(A1).round(3)

0

비가역 행렬임은 피벗열이 아닌 벡터의 존재를 나타냅니다. 이 피벗열의 수는 기저벡터의 수와 같으며 행렬의 급수(rank)로 확인할 수 있습니다. 다음 결과와 같이 기저벡터의 수는 2개입니다.

la.matrix_rank(A1)

2

d,P=la.eig(A1)
np.allclose(P.T, la.inv(P))

True

print(d.round(3))

[25.  0.  9.]

고유값 0일 경우 식 7이 성립되지 않습니다. 즉, 0은 고유값이 아니므로 특이값은 2개입니다. 이는 행렬 A1의 급수와 같습니다.

\begin{align}\tag{식 7}A^TAv &=\lambda v\\\lambda, v:&\,\text{고유값, 고유벡터}\end{align}

sigma=np.sqrt(d)
print(sigma.round(3))

[5. 0. 3.]

그러므로 0을 제외한 고유값 25, 9의 제곱근 5와 3이 행렬 A의 특이값이 됩니다.

이 특이값들의 수는 원시행렬의 열공간(기저공간)의 수와 같습니다. 열공간은 sympy객체.columnspace() 함수로 확인할 수 있습니다.

Matrix(A).columnspace()

[Matrix([
 [3],
 [2]]),
 Matrix([
 [2],
 [3]])]

특이값

비정방행렬 A의 A^TA는 대칭행렬입니다. 그러므로 0이 아닌 고유값들에 대응하는 고유벡터들은 정규직교기저입니다.

0이 아닌 고유값의 수가 행렬 A의 특이값의 수가 되며 이는 A의 열공간(Col A) 수 즉, 급수와 같습니다(식 8).

Rank A = 특이값의 수

(식 8)