[DAY 65] Numpy를 이용한 SVD(특이값 분해)

데이터 압축, 토픽모델링 등에서 사용되는 SVD를 구현했다.

scikit-learn에서 추상화되어 API로 기능을 제공하고 있지만 원리를 이해하기 위해 Numpy를 사용했다.

실습 요약

선형대수 모듈을 사용한 SVD(Singular Value Decomposition, 특이값 분해)
행렬 곱셈, numpy의 N-dimensional array 다루기 - indexing, slicing, transpose, reshape
대각행렬 만들기, allclose 등을 통한 행렬 비교
numpy의 clip 함수 활용

실습

import numpy as np
import matplotlib.pyplot as plt
from scipy import datasets

실습에 필요한 라이브러리 로드

img = datasets.face()
plt.imshow(img)

scipy.datasets의 face 메소드를 사용하면 가로x세로 1024x768인 라쿤 얼굴 컬러 이미지를 얻을 수 있다. plt.imshow를 통해 확인한다.

cute

img.shape  # (768, 1024, 3)
img.size  # 2359296
img.ndim  # 3

해당 이미지는 768행 1024열 3채널의 3차원 array이며, 사이즈는 768*1024*3이다. 3차원이기 때문에 ndim 값은 3이다.

img

array의 값을 확인하면 아래와 같다.

컬러 이미지는 3개의 채널의 R, G, B 값으로 색상이 결정된다.

img[:, :, 0].shape, img[:, :, 1].shape, img[:, :, 2].shape

# 실행 결과
((768, 1024), (768, 1024), (768, 1024))

각 채널은 모두 768행 1024열의 행렬 형태(2D array)다.

R, G, B 채널

R, G, B 채널의 값은 위와 같으며, 예를 들어 사진의 0행 3열 부분은 R 153, G 144, B 165의 RGB 가산혼합으로 색상이 만들어진다.

img[0, 3]  # 각 좌표마다 3채널 존재

# 실행 결과
array([153, 144, 165], dtype=uint8)

0행 0열부터 767행 1023열까지 이런 방식으로 색상이 결정되며 그것이 합쳐져 가로 1024 세로 768의 라쿤 이미지가 되는 것이다.

# 0번째 채널 : Red
plt.imshow(img[:, :, 0], cmap='Reds')

0번째 채널 값으로 그린 라쿤

 

# 1번째 채널 : Green
plt.imshow(img[:, :, 1], cmap='Greens')

1번째 채널 값으로 그린 라쿤

 

# 2번째 채널 : Blue
plt.imshow(img[:, :, 2], cmap='Blues')

2번째 채널 값으로 그린 라쿤

이런 RGB 값이 합쳐져 컬러풀한 이미지가 된다.

img_array = img / 255
img_array.min(), img_array.max()

# 실행 결과
(0.0, 1.0)

0~255 값을 255로 나누어 normalize했다. 이미지가 달라지는 것은 아니다.

실습 중에 numpy reshape과 3D array에 대한 부분이 있었는데 이해가 잘 되지 않았다.
그래도 수강생분들의 많은 도움으로 이해할 수 있었다.
numpy 3D array에 대해 정리하는 글을 하나 써야겠다.


행렬 곱셈(matrix multiplication)
두 개의 행렬에서 한 개의 행렬을 만드는 이항 연산
A행렬과 B행렬을 곱한다면, 첫째 행렬(A)의 열 개수와 둘째 행렬(B)의 행 개수가 동일해야 성립한다. 그 결과 행렬은 A행렬의 행 개수와 B행렬의 열 개수를 가진다.

https://blog.finxter.com/numpy-matmul-operator/

위 사진은 2행2열인 a와 2행 3열인 b의 곱을 보여준다.
결과 행렬은 2행 3열이 되었다.

이 세 표현은 똑같이 a와 b의 행렬 곱셈을 의미한다
1. a@b
2. np.dot(a, b)
3. np.matmul(a, b)

흑백 이미지로 변환

# gray 이미지로 변경 : 2차원
img_gray = img_array @ [0.2126, 0.7152, 0.0722]
img_gray.shape

# 실행 결과
(768, 1024)

3차원인 img_array에 np.array([0.2126, 0.7152, 0.0722])를 곱해 2차원으로 변경했다. (RGB -> gray)

plt.imshow(img_gray, cmap='gray')
plt.show()

2차원으로 변경된 gray scale의 라쿤 확인

SVD(Singular Value Decomposition)

from numpy import linalg

선형대수 모듈 linalg를 사용한다.

# linalg.svd 로 img_gray 분해
U, s, Vt = linalg.svd(img_gray)
U.shape, s.shape, Vt.shape

# 실행 결과
((768, 768), (768,), (1024, 1024))

linalg.svd를 통해 특이값 분해되었다.

그렇다면 다시 위 3개를 행렬 곱셈하면 이미지가 복원될까?
하지만 지금은 s와 Vt가 행렬 곱셈하기 위한 shape가 성립되지 않는다.

Sigma = np.zeros((U.shape[1], Vt.shape[0]))
np.fill_diagonal(Sigma, s)
Sigma

np.zeros를 통해 0으로 채워진 768행 1024열의 행렬을 생성하고, fill_diagonal로 사선으로 s값을 채운 후 확인했다.

이제 U, Sigma, Vt를 행렬 곱셈해보자.

# 매우 근소한 차이를 보이는 것을 알 수 있다
dif = img_gray - (U @ Sigma @ Vt)
dif

특이값 분해 전 img_gray와, 특이값 분해 후 다시 행렬 곱셈한 결과의 차이가 있다.

아주 작은 차이를 확인할 수 있다. == 연산자를 통해 확인해 봐도 True가 나오지 않는다.

# linalg.norm 차이값 비교
# e-12 정도의 오차가 있는 것으로 확인
linalg.norm(dif)

# 실행 결과
1.3245284105876096e-12

linalg.norm을 통해 확인해 보면 약 e-12의 차이를 볼 수 있다.
이렇게 매우 작은 차이를 보일 때 흡사한지 확인하는 함수가 있다.

# allclose : Returns True if two arrays are element-wise equal within a tolerance
# 거의 같은지 확인하는 함수. 두 array의 차이가 특정 tolerance보다 작으면 True
np.allclose(img_gray, U @ Sigma @ Vt)

# 실행 결과
True

np.allclose 함수는 거의 같은 경우 True를 반환하는 함수다. 특정 한계보다 차이가 작으면 같다고 판단하고 True를 반환한다.

Sigma 시각화

x축은 행을, y축은 특이값 분해 전의 중요한 성분을 나타낸다. 앞쪽 행에 중요한 정보들이 많은 것을 알 수 있다.

Sigma의 대각 원소 중 상위 k개의 특이값만 사용

import os
from sys import getsizeof

base = f'data/numpy-svd/'

k = 10
approx10 = U @ Sigma[:, :k] @ Vt[:k, :]
# 객체 크기 확인
print(f'객체 크기 - 원본 : {getsizeof(img_gray)} bytes, approx10 : {getsizeof(approx10)} bytes')

fig, ax = plt.subplots(1, 2, figsize=(8, 4))
# 원본 이미지
ax[0].imshow(img_gray, cmap='gray')
ax[0].set_title('origin')
# approx10 이미지
ax[1].imshow(approx_10, cmap='gray')
ax[1].set_title('approx10')
# 사진 각각 저장
plt.imsave(f'{base}approx_ori.png', img_gray, cmap='gray')
plt.imsave(f'{base}approx_10.png', approx10, cmap='gray')
# 사이즈 확인 : 파일로 저장한 결과 압축되었다
origin_size = os.stat(f'{base}approx_ori.png').st_size
approx10_size = os.stat(f'{base}approx_10.png').st_size
print(f"파일 크기 - 원본 : {origin_size} bytes, approx10 : {approx10_size} bytes")

상위 10개의 특이값만 사용하여 복원해 보고, 사진과 사진의 사이즈를 확인한다.

객체 크기는 같지만 파일 크기가 원본에 비해 절반 이하로 절약되었다. 하지만 복원이 너무 약하게 된 것 같다.

이번엔 상위 100개의 특이값을 사용해서 확인한 결과, 원본과 비슷하게 복원된 것을 볼 수 있다. 하지만 압축률이 크지 않다.

3차원 이미지 SVD
3차원의 경우엔 2차원에 채널이 추가된 형태인 점에 유의해야 한다.
diagonal을 채널별로 한 점, clip 함수를 통해 가장 작은 값을 0, 가장 큰 값을 1이 되게 변경한 점 등 2차원과 다른 점이 있었다.

3차원의 경우 사용할 특이값 k=10으로 지정했을 경우 파일의 압축률이 2차원의 경우보다 더 작았다.

특이값 분해와 행렬 곱셈, 행, 열, 채널의 개념, numpy의 여러 함수에 대해 익혔다.