[Linear Algebra] 특이값분해(Singular Value Decomposition)

특이값분해(Singular Value Decomposition)

관련된 내용

• 특이값

m×n 차원의 행렬 A의 분해는 식 1과 같은 대각 행렬(Σ)와 관련됩니다.

$$\tag{식 1}\Sigma=\begin{bmatrix}D& 0\\0 & 0 \end{bmatrix}$$

식 1에서 D는 A^TA의 고유값(≠ 0)을 대각원소로 구성한 대각행렬이고 차원은 r×r(r ≤ Min(m, n))입니다. 또한 Σ의 차원은 행렬 A의 차원과 같습니다.

행렬 A에 대해 A^TA의 0이 아닌 고유벡터와 고유값을 각각 v_i와 λ_i라 하고 각 고유벡터가 정규직교관계이면 식 2와 같이 나타낼 수 있습니다. 그러므로 고유벡터에 행렬 A와의 내적 결과인 Av₁, Av₂, …, Av_r 벡터들은 서로 직교합니다. 이들을 정규화(normalization)한 벡터들을 u₁, u₂, ..., u_r은 라고 하면 이 벡터들은 모두 정규 직교 벡터가 됩니다.

\begin{align}u_i & =\frac{Av_i}{\Vert{Av_i}\Vert}\\\tag{식 2}& = \frac{Av_i}{\sigma_i}\\ Av_i&=\sigma_i u_i\\ & 1\le i \le r \end{align}

식 2의 AV는 고유값이 0이 되는 모든 부분까지 확장하여 식 3과 같이 나타낼 수 있습니다.

\begin{align}AV&=\begin{bmatrix}Av_1&\cdots& Av_r& 0& \cdots& 0 \end{bmatrix}\\&= \begin{bmatrix} \sigma_1 u_1&\cdots& \sigma_r u_r& 0& \cdots& 0 \end{bmatrix}\\ \tag{식 3}& = \begin{bmatrix} u_1&\cdots& u_r& 0& \cdots& 0 \end{bmatrix} \begin{bmatrix} \sigma_1&\cdots& 0 &0& \cdots & 0\\ \vdots& \ddots & \vdots & \vdots & \ddots & \vdots\\ 0& \cdots & \sigma_r & 0 & \cdots & 0\\ \vdots& \ddots & \vdots & \vdots & \ddots & \vdots\\ 0& \ddots & 0 & 0 & \ddots & 0\end{bmatrix}\\ & = U\Sigma\end{align}

고유행렬인 V는 정규직교행렬(orthogonal matrix)이므로 식 3은 식 4와 같이 나타낼 수 있습니다. 식 3에서 0인 고유값을 포함하였으므로 V 역시 고유값 0에 대응하는 고유벡터를 포함합니다.

\begin{align}AV&=U\Sigma\\\tag{식 4}A&=U\Sigma V^{-1}\\&=U\Sigma V^T \end{align}

식 4를 특이값 분해(singular value decomposition, SVD)라고 하며 U의 각 열벡터를 좌특이벡터(left singular vectors), V를 우특이 벡터(right singular vectors)라 합니다.

예 1)

행렬 A의 특이값 분해를 실행해 봅니다.

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

A의 특이값 분해를 위해 A^TA의 고유값과 고유벡터를 계산합니다. 이 정방행렬은 대칭행렬로서 고유행렬은 정규직교기저입니다.

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.]

0이 아닌 고유값에 대응하는 특이값(σ)과 그 행렬을 생성합니다. 특이값 행렬은 A과 같은 차원입니다.

r=d[[0,2]]
V=P[:,[0,2]]
sigma=np.sqrt(r)
print(sigma.round(3))

[18.974  9.487]

sigmaMat=np.c_[np.diag(sigma), np.zeros((2,1))]
print(sigmaMat.round(3))

[[18.974  0.     0.   ]
 [ 0.     9.487  0.   ]]

0이 아닌 고유값에 대응하는 고유벡터들로 구성된 행렬 U는 다음과 같습니다.

u1=(A@V[:,0])/sigma[0]
print(u1.round(3))

[-0.949 -0.316]

u2=(A@V[:,1])/sigma[1]
print(u2.round(3))

[ 0.316 -0.949]

U=np.c_[u1, u2]
print(U.round(3))

[[-0.949  0.316]
 [-0.316 -0.949]]

행렬 A 분해에 대한 위 식을 검정합니다. 위 과정에서 0이 아닌 고유값들만을 고려했기 때문에 A에 대한 분해는 고유값 0에 대응하는 열 벡터를 V에 첨가하여야 합니다(객체 V1).

V1=np.c_[V, P[:, 1]]
print(V1.round(3))

[[-0.333 -0.667 -0.667]
 [-0.667 -0.333  0.667]
 [-0.667  0.667 -0.333]]

np.allclose(A,  U@sigmaMat@V1)

True

위의 과정은 numpy.linalg.svd() 함수로 대체할 수 있습니다.

U, s, VT=la.svd(A)
print(U. round(3))

[[ 0.949 -0.316]
 [ 0.316  0.949]]

print(s.round(3))

[18.974  9.487]

print(VT. round(3))

[[ 0.333  0.667  0.667]
 [ 0.667  0.333 -0.667]
 [-0.667  0.667 -0.333]]

예 2)

행렬 B의 특이값 분해를 실행합니다.

$$B=\begin{bmatrix}1& -1\\2& 2\\2& -2 \end{bmatrix}$$

linalg.svd() 함수를 적용합니다.

B=np.array([[1,-1],[-2,2],[2, -2]])
U,s,VT=la.svd(B)
print(U.round(3))

[[-0.333  0.667 -0.667]
 [ 0.667  0.667  0.333]
 [-0.667  0.333  0.667]]

print(s.round(3))

[4.243 0.   ]

sMat=np.r_[s.reshape(1,-1), np.zeros((2,2))]
print(sMat.round(3))

[[4.243 0.   ]
 [0.    0.   ]
 [0.    0.   ]]

print(VT.round(3))

[[-0.707  0.707]
 [ 0.707  0.707]]

print(np.isclose(B,U@sMat@VT))

[[ True  True]
 [ True  True]
 [ True  True]]

위 결과는 특이값 분해에 대한 좌특이행렬, 우특이행렬을 반환합니다. 그러나 특이값은 1개인데 반해 좌특이 벡터는 3개로서 불일치를 보여줍니다. 이를 확인하기 위해 위 과정을 다시 복기합니다.

B^TB(=B1)의 고유값과 고유행렬을 계산합니다.

B1=B.T@B 
print(B1)

[[ 9 -9]
 [-9  9]]

d, P=la.eig(B1)
print(d)

[18.  0.]

0이 아닌 고유값을 특이값으로 하므로 행렬 B의 특이값은 1개입니다. 그러므로 그에 대응하는 고유벡터로부터 계산되는 U는 1개의 벡터로 구성됩니다.

s=np.sqrt(d[0])
print(s.round(3))

4.243

sMat1=np.array([[s, 0],[0,0],[0,0]])
print(sMat.round(3))

[[4.243 0.   ]
 [0.    0.   ]
 [0.    0.   ]]

u1=B@(P[:,0].reshape(-1,1))/s
print(u1.round(3))

[[ 0.333]
 [-0.667]
 [ 0.667]]

U는 정규직교행렬로서 정방행렬이어야 합니다. 즉, 행렬의 차원이 3×3이어야 합니다. 그러므로 u2, u3를 생성하기 위해 u1의 선형결합으로부터 영공간을 사용합니다(식 5).

\begin{align}\tag{식 5}& u1_{\text{linsys}}=\begin{bmatrix}0.333 & -0.667& 0.667\\0. & 0. & 0. \\0. & 0. & 0. \end{bmatrix}\\& \begin{bmatrix}0.333 & -0.667& 0.667\\0. & 0. & 0. \\0. & 0. & 0. \end{bmatrix}\begin{bmatrix} x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix} 0\\0\\0\end{bmatrix} \end{align}

위 선형결합의 표준행렬(u1_linsys)의 영공간은 다음과 같습니다.

u1_linsys=np.r_[u1.reshape(1,-1), np.zeros((2,3))]
print(u1_linsys.round(3))

[[ 0.333 -0.667  0.667]
 [ 0.     0.     0.   ]
 [ 0.     0.     0.   ]]

u1_base=Matrix(u1_linsys).nullspace()
u1_base=u1_base[0].row_join(u1_base[1])
u1_base

$\begin{bmatrix}2.0& -2.0\\1 & 0\\ 0 & 1 \end{bmatrix}$

위 영공간은 u1과의 선형결합의 해집합의 기저벡터들입니다. 이 벡터들의 정규직교기저벡터는 QR 분해의 Q가 됩니다.

Q,R=la.qr(np.array(u1_base, dtype=float))
print(Q.round(3))

[[-0.894  0.298]
 [-0.447 -0.596]
 [-0.    -0.745]]

u1과 위 결과의 결합으로 좌특이행렬이 됩니다.

U=np.c_[u1, Q]
print(U.round(3))

[[ 0.333 -0.894  0.298]
 [-0.667 -0.447 -0.596]
 [ 0.667 -0.    -0.745]]

위 결과들로부터 B=UΣV^T를 확인합니다.

print(np.isclose(B, U@sMat@P.T))

[[ True  True]
 [ True  True]
 [ True  True]]