[파이썬과 선형대수] 특잇값 분해의 계층화, Reduced SVD 정의와 직접 구현하기

https://aaaaaaaaaaayowooji.tistory.com/55

[파이썬과 선형대수]특잇값 분해의 정의, svd 안쓰고 직접 구현해보기

이전 글에서 특잇값 분해(SVD) 정의에 대해 알아보고 파이썬으로 구현도 해보았다.
SVD의 진정한 의미는 이번 글에서 다루어볼 계층화(layer) 가 진짜라고 생각한다.

특잇값 분해의 계층화

특잇값 분해를 하면,
$$A=\sigma_1u_1v_1^T+\sigma_2u_2v_2^T+\sigma_3u_3v_3^T+...+\sigma_mu_mv_m^T=\sum_{i=1}^{min(M,N)} \sigma_iu_iv_i^T$$
로 나타낼 수 있겠다.

특잇값 분해의 의미

공식을 대충 보면 벡터의 곱으로 간단히 나타낸 것 같지만 $u_iv_i^T$는 벡터의 외적이다. 계산 결과는 $M\times N$ 행렬이다. 여기에 $\sigma$는 크기 순으로 정렬해놓은 상태다.

다시 생각해보자. $A$라는 행렬은 $M$개(엄밀히 따지면 $min(M,N)$이지만 $M$이라고 가정)의 $M \times N$행렬들을 쌓아서 모두 더한 결과인 것이다.
추가로 $\sigma$는 내림차순으로 작아진다고 했다. 앞쪽일수록 $A$를 구성하는 요소에 영향을 더 미치게 될 것이다. 반대로 뒤쪽일수록 영향은 미미하다.

특잇값 분해의 의미를 알면 당연한 질문이 따라온다.

특잇값이 아주 작다면 무시해도 되는가?

 

이와 관련된 이야기는 아래에 등장할 것이다.

파이썬으로 확인하기

코드

#13-5_2
import numpy as np
import matplotlib.pyplot as plt

m,n=10,8

A=np.random.randn(m,n)
U,s,V=np.linalg.svd(A)
S=np.zeros(A.shape)
np.fill_diagonal(S,s)

fig,axs=plt.subplots(2,6,figsize=(15,5))

layer=np.zeros((m,n,10))
err=np.zeros(5)

#=====================plotting========================
axs[0,0].imshow(A,cmap='gray')
axs[0,0].set_title(f"A")
axs[0,0].set_axis_off()

axs[1,0].set_axis_off()

for i in range(10):

    if i<5:
        #각 layer plot
        layer[:,:,i]=s[i]*U[:,i].reshape(m,1)@V[i,:].reshape(1,n)
        axs[0,i+1].imshow(layer[:,:,i],cmap='gray')
        axs[0,i+1].set_title(f"layer{i+1}")
        axs[0,i+1].set_axis_off()
    else:
        #layer1부터 layer n까지 합을 plot
        layer[:,:,i]=np.sum(layer[:,:,0:i-4],axis=2)
        err[i-5]=np.linalg.norm(A-layer[:,:,i])
        axs[1,i-4].imshow(layer[:,:,i],cmap='gray')
        axs[1,i-4].set_title(f"layer1:{i-4}/err={np.round(err[i-5],2)}")
        axs[1,i-4].set_axis_off()



plt.show()

결과

해당 결과는 $A$행렬을 계층화하여 각 계층을 나타낸 것이다. 크기가 $10 \times 8$이기 때문에 실제로 특잇값은 8개가 나오겠지만 크기관계상 5개까지만 나타내보았다.

또한 그 아랫줄에는 레이어 1층부터 n층까지 더한 행렬을 나타내본 결과다. 단순히 $A$와 $u_{1}\sigma_{1}v_1^T$의 오차는 7.43이지만 실제로 계층들이 더해지면서 $A$와 $\sum_{i=1}^{5}u_{i}\sigma_{i}v_i^T$의 오차는 2.34로 줄었다. 실제로 layer1부터 layer5까지 더한 결과와 $A$를 비교해보면 거의 비슷함을 알 수 있다.

결과를 plotting하는 과정에서 더 간단한 알고리즘이 존재하겠음을 깨달았지만 편의상 결과가 잘 나오는 만큼 코드는 냅두었다.

---

Economic SVD(경제적 특잇값 분해)

실제로 데이터를 다루게 되면 행렬의 크기가 매우 커진다고 한다. (M의 크기가 백만의 단위로 넘어가는...)
행렬의 크기를 줄일 수 있다면 메모리 효율성과 계산 속도를 향상시킬 수 있을 것이다.
그 방법으로 대두되는 것이 all zero rows or columns를 삭제하는 것이다.

정의

우선 $S$는 원래 아래와 같이 대각행렬에 all zero rows or columns들이 붙어있는 꼴이었다.
$M>N$인 경우
$$\begin{bmatrix}\sigma_1&0&0\\ 0&\sigma_{2}&0 \\ 0&0&\sigma_3\\0&0&0\\0&0&0\end{bmatrix}$$

$M<N$인 경우
$$\begin{bmatrix}\sigma_1&0&0&0&0\\ 0&\sigma_{2}&0&0&0 \\ 0&0&\sigma_3&0&0\end{bmatrix}$$

여기서 0으로만 이루어진 행이나 열을 삭제해 대각 행렬로 만든다.
기존에 $M>N$인 경우라면 $U$행렬의 $N+1$열들은 어차피 $\sigma$ 값이 0이기 때문에 의미가 없다. 그래서 아예 지워버리자는 것이 Economic SVD(경제적 특잇값 분해)이다.

기존 full SVD에서 $U,S,V$로 표기했지만 Economic SVD에서는 구별하기 위해 $\hat U, \hat S, \hat V^T$로 표기한다.

행렬식으로 나타내보면,

 

$M>N$인 경우에
$$\hat A=\begin{bmatrix}u_1&u_2&u_3&\cdots&u_n\end{bmatrix}\begin{bmatrix}\sigma_1&0&0&...&0\\ 0&\sigma_{2}&0&...&0 \\ 0&0&\sigma_3&...&0\\0&0&0&...&\sigma_n\end{bmatrix}\begin{bmatrix}v_1&v_2&v_3&\cdots &v_n\end{bmatrix}^T$$

$$\hat U:M\times N$$
$$\hat S:N\times N$$
$$\hat V^T:N\times N$$

 

$N>M$인 경우에
$$\hat A=\begin{bmatrix}u_1&u_2&u_3&\cdots&u_m\end{bmatrix}\begin{bmatrix}\sigma_1&0&0&...&0\\ 0&\sigma_{2}&0&...&0 \\ 0&0&\sigma_3&...&0\\0&0&0&...&\sigma_m\end{bmatrix}\begin{bmatrix}v_1&v_2&v_3&\cdots&v_n\end{bmatrix}^T$$
$$\hat U:M\times M$$
$$\hat S:M\times M$$
$$\hat V^T:M\times N$$

파이썬으로 확인하기

코드

import numpy as np

sizes=[(5,3),(3,5),(4,4)]

for size in sizes:
    A=np.random.randn(size[0],size[1])
    U,s,V_T=np.linalg.svd(A,full_matrices=False)
    S=np.diag(s)
    print(f'U의 shape:{U.shape},S의 shape:{S.shape},V^T의 shape:{V_T.shape}')

svd함수에서 사용한 full_matrices 인자에 대해서는 아래 글에 자세하게 작성해놓았으니 참고하면 된다.
https://aaaaaaaaaaayowooji.tistory.com/52

[파이썬 라이브러리]numpy.linalg 함수 분석(2)(qr,eig,svd)

 

결과

U의 shape:(5, 3),S의 shape:(3, 3),V^T의 shape:(3, 3)
U의 shape:(3, 3),S의 shape:(3, 3),V^T의 shape:(3, 5)
U의 shape:(4, 4),S의 shape:(4, 4),V^T의 shape:(4, 4)

---

Reduced SVD(축소형 특잇값 분해)

위에서 경제형 SVD에 대해 알아보았다. 경제형 SVD에서는 특잇값이 작다고 무시하지는 않았다. 특잇값이 0인 행이나 열에 대해서 의미가 없기 때문에 삭제했을 뿐이다.

축소형 특잇값 분해(Reduced SVD)에서는 무시해도 될 만한(negligible) 특잇값을 무시한다.
$k<min(M,N)$인 $k$번째 특잇값이 있을 때, $k+1$이후의 특잇값들이 예를 들어 $1e-5$보다 작은 값들이어서 원래 행렬 $A$의 재건(reconstructed A) 에 영향을 끼치지 않는다고 생각이 되면 무시하겠다는 것이다.
행렬식으로 보면,
$$\tilde A=\begin{bmatrix}u_1&u_2&u_3&\cdots&u_k\end{bmatrix}\begin{bmatrix}\sigma_1&0&0&\cdots&0\\ 0&\sigma_{2}&0&...&0 \\ 0&0&\sigma_3&...&0\\0&0&0&...&\sigma_k\end{bmatrix}\begin{bmatrix}v_1&v_2&v_3&\cdots&v_n\end{bmatrix}^T
$$

마찬가지로 Economic SVD에서는 구별하기 위해 $\hat U, \hat S, \hat V^T$로 표기한 것처럼 Reduced SVD에서는 구별하기 위해 $\tilde U, \tilde S, \tilde V^T$로 표기한다.

$$\tilde{U}:M\times k$$
$$\tilde S:k\times k$$
$$\tilde V^T:k\times N$$

더이상 정규직교행렬이 아니다?

$$\tilde U \tilde U^{T}\neq I \quad \tilde U^{T} \tilde U= I$$
$$\tilde V \tilde V^{T}\neq I \quad \tilde V^{T} \tilde V= I$$
$\tilde U$는 $M \times k$이다. 정규직교 행렬은 열끼리 행끼리 크기가 1이고 직교한다.
$M \times k$이라는 것을 생각해보면 $k$개의 열끼리는 직교가 유지되고 각 열의 크기는 1이지만 행끼리는 그렇지 않다. 따라서 $\tilde U^{T}\tilde U$의 계산은 크기가 모두 1인 행과 열의 내적으로 계산되므로 $I$가 되는 반면에 $\tilde U \tilde U^{T}\neq I$는 $M \times M$행렬이 될 뿐이다.

$$\tilde U \tilde U^{T}=\begin{bmatrix}u_1&u_2&u_3&\cdots& u_k\end{bmatrix}\begin{bmatrix}-u_1^T-\\-u_2^T-\\-u_3^T- \\ \vdots \\-u_k^T-\end{bmatrix}$$
$$\tilde U^T \tilde U=\begin{bmatrix}-u_1^T-\\-u_2^T-\\-u_3^T- \\ \vdots \\-u_k^T-\end{bmatrix}\begin{bmatrix}u_1&u_2&u_3&\cdots& u_k\end{bmatrix}$$

$\tilde V \tilde V^T$ 역시 $\tilde U \tilde U^T$와 유사하게 계산해보면 그 이유를 알 수 있다.

파이썬으로 확인하기

위 사례를 파이썬으로 구현해보았다.

코드

import numpy as np
import matplotlib.pyplot as plt

A=np.random.randint(0,20,(6,8))

U,s,V_T=np.linalg.svd(A)
k=4
U_til=U[:,:k]
print(f'shape of U_til:{U_til.shape}')
print(f'U_til@U_til.T:\n{np.round(U_til@U_til.T,2)}')
print(f'U_til.T@U_til:\n{np.round(U_til.T@U_til,2)}\n')

V_T_til=V_T[:k,:]
print(f'shape of V_T_til:{V_T_til.shape}')
print(f'V_til_T@V_til:\n{np.round(V_T_til@V_T_til.T,2)}')
print(f'V_til_T.T@V_til:\n{np.round(V_T_til.T@V_T_til,2)}')

결과

shape of U_til:(6, 4)
U_til@U_til.T:
[[ 0.49  0.26 -0.27  0.02  0.32  0.09]
 [ 0.26  0.84  0.07  0.08 -0.23  0.06]
 [-0.27  0.07  0.72  0.21  0.02  0.28]
 [ 0.02  0.08  0.21  0.72  0.19 -0.34]
 [ 0.32 -0.23  0.02  0.19  0.65  0.19]
 [ 0.09  0.06  0.28 -0.34  0.19  0.58]]
U_til.T@U_til:
[[ 1.  0. -0. -0.]
 [ 0.  1.  0.  0.]
 [-0.  0.  1. -0.]
 [-0.  0. -0.  1.]]

shape of V_T_til:(4, 8)
V_til_T@V_til:
[[ 1. -0.  0.  0.]
 [-0.  1.  0. -0.]
 [ 0.  0.  1.  0.]
 [ 0. -0.  0.  1.]]
V_til_T.T@V_til:
[[ 0.72 -0.01  0.06 -0.14 -0.03  0.13  0.4  -0.06]
 [-0.01  0.64 -0.11  0.11  0.21  0.4  -0.04  0.07]
 [ 0.06 -0.11  0.45  0.24 -0.15 -0.01  0.1   0.38]
 [-0.14  0.11  0.24  0.49  0.3  -0.15  0.14  0.17]
 [-0.03  0.21 -0.15  0.3   0.5  -0.15  0.17 -0.2 ]
 [ 0.13  0.4  -0.01 -0.15 -0.15  0.44 -0.08  0.15]
 [ 0.4  -0.04  0.1   0.14  0.17 -0.08  0.35 -0.04]
 [-0.06  0.07  0.38  0.17 -0.2   0.15 -0.04  0.4 ]]

A의 재건(Reconstruction of A)

$k$개의 특잇값 만으로 다시 $\tilde U \tilde S \tilde V^T$를 진행해주면 $A$와 동일하지는 않지만 유사하게 재건할 수 있다.
말로만 하는 것보다 그래프를 보는 것이 더 이해가 빠르다.

파이썬으로 확인하기

코드

import numpy as np
import matplotlib.pyplot as plt

A=np.random.randint(0,20,(6,8))

U,s,V_T=np.linalg.svd(A)

err_list=np.zeros(A.shape[0])
for k in range(1,A.shape[0]):

    U_til=U[:,:k]
    V_T_til=V_T[:k,:]
    s_til=s[:k]
    S=np.diag(s_til)

    err_list[k-1]=np.linalg.norm(A-U_til@S@V_T_til)

plt.plot(range(1,s.shape[0]+1),s/np.sum(s)*100,'ks-')
plt.plot(range(1,s.shape[0]+1),err_list,'ro-')
plt.legend(['variation(%)','reconstruction error'])
plt.xlabel("number of singular values")

plt.show()

결과

행렬 $A$는 $6 \times 8$이다. 검정색 그래프는 각 특잇값을 전체 특잇값의 합으로 나눈 값이다. 즉, 특잇값의 중요도라고 생각할 수 있다. 중요도가 커질수록 $A$를 재건하는 데에 영향을 많이 미친다.
빨간색 그래프는 $A$를 재건했을 때 실제 $A$와의 프로베니우스 오차다.

이를 통해 더 큰 $k$값을 지정해 재건할수록 오차가 더 작아짐을 알 수 있다. 다시 말해서 실제 $A$와 유사해진다는 것이다.
이런 경우에 $k$은 오차 그래프가 급 완만해지는 5정도로 잡으면 적당하다.

Eckart-Young Theorem

위에서처럼 $k$을 몇으로 잡느냐는 중요한 문제다. 적당한 수를 잡지 않으면 별의미 없는 특잇값까지 포함되어 메모리를 잡아먹거나 계산이 느려지고 반대로는 제대로 된 재건을 할 수도 없다.

이에 대한 이론으로 Eckart-Young Theorem이다.
$$||A-A_k||_F\leq||A-B||_F$$
주어진 행렬 $A$에 대해, rank가 $r$이하인 $A_k$ 중에서, 프로베니우스 노름으로 $A$와 가장 가까운 행렬 $A_k$를 구하는 것이 이 이론의 중점이다.

그 중에서 최적의 $k$를 구하는 방법에는 여러 가지가 있다.
대표적인 방법으로 위에서 보인 그래프가 있다. Scree Plot이라고 하는 이 그래프는 특잇값의 크기를 순서대로 시각화한 것으로 특잇값이 급격히 감소하는 엘보우 포인트
를 $k$로 잡을 수 있겠다.

다른 방법으로 $\sum_{i=1}^{r} \sigma^2$을 총 에너지라고 했을 때, 에너지 비율은 ${\sum_{i=1}^{k} \sigma^2}/{\sum_{i=1}^{r} \sigma^2}$이 특정값(90% 혹은 95%) 이상이 되는 $k$을 선택하는 방법이 있다.

---

 

특잇값 분해의 계층화에 대해 다양한 예시를 살펴보았다.
하다보니 많은 예를 들려는 욕심에 글이 많이 길어졌다. 하지만 더욱 더 SVD의 본질에 다가갔음을 느끼는 계기가 되었다.

---

출처

서적

개발자를 위한 실전 선형대수학/한빛미디어/마이크,코헨 지음
Linear Algebra and its Application 5th Ed/PEARSON/David C.Lay, Steven R.Lay

WebSite

https://python.quantecon.org/svd_intro.html#eckart-young-theorem

Singular Value Decomposition (SVD)

 

https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=3