SVD

3 minute read

개요

이번 시간에는 차원 축소, 추천 시스템 등에 사용되는 특이값 분해 즉 SVD에 대해 다뤄보겠습니다.

선형대수학에서는 행렬의 대각화(diagonalization) 개념이 소개됩니다. non-singular 정방 행렬을 그 eigenvalues와 eigenvectors로 분해하여 표현하는 기법입니다. 행렬을 보다 basic한 성분들의 곱으로 나타낸다는 컨셉을 non-square 행렬로 확장한 것이 바로 SVD입니다.

본 포스트에서는 reduced SVD와 그 알고리즘에 쓰이는 Gram-Schmidt process를 직접 구현해보겠습니다.


필요한 라이브러리와 데이터 불러오기

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy

from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler
from sklearn.metrics.pairwise import cosine_similarity

임의의 데이터를 생성합니다. 샘플 즉 행의 개수는 1,000개입니다.

data = make_classification(n_features = 6, n_redundant=0,
                           n_samples=10**3, weights=[0.9], random_state= 42, )
data_x = data[0]
data_y = data[1]
print(data_x)
[[ 0.13556342  0.72851864  0.63760581  0.3042709   0.08448742  0.25876981]
 [-0.63002498  0.7319555  -0.19692181  1.78875994 -1.64024454 -0.46458108]
 [ 1.18673462  1.43633541 -0.23811044 -1.14142277 -1.98512688  0.04600275]
 ...
 [ 0.39579885  0.26355081  1.28409579  1.02826725 -0.97134377 -2.42639207]
 [ 1.20660815 -1.52716786 -1.05508937  0.88695219 -0.93969063  0.12626774]
 [ 0.75951489 -0.29175863 -0.05982614  1.65199039 -1.12225881  0.4935602 ]]



Gram-Schmidt process

Gram-Schmidt process는 서로 독립인 n개의 벡터를 각각에 대해 수직으로 만드는 기법입니다. 구체적인 과정은 아래와 같습니다.

$\large \displaystyle \mathrm{V} = \begin{pmatrix} \mathbf{v}_1 &\mathbf{v}_2 & \cdots & \mathbf{v}_p \end{pmatrix}$ , $i \in {1,2,\dots,p}$ 에 대해 $\mathbf{v}_i$는 i-번째 컬럼 벡터입니다. 이 $\mathbf{v}_i$를 직교 벡터인 $\mathbf{u}_i$로 바꿉니다. $\left\langle \mathbf{u}_i, \mathbf{v}_j \right\rangle$는 두 벡터의 내적입니다.

캡처1

$\large \therefore \mathrm{U} = \begin{pmatrix} \mathbf{u}_1 & \mathbf{u}_2 &\cdots & \mathbf{u}_p \end{pmatrix}$

이렇게 만들어진 $\large \mathrm{U}$ 행렬의 서로 다른 컬럼 벡터들은 직교하게 됩니다. 각 컬럼 벡터를 그 벡터의 norm으로 나누어 주면, 벡터들은 정규 직교 기저가 됩니다.

알고리즘으로 구현하면 아래와 같습니다. fit_transform 시 행렬의 컬럼을 벡터로 하며 normalizing 즉 정규화하는 것을 default로 두었습니다.

class gram_schmidt:
        
    def _proj(self, u, v):
        return (np.vdot(v, u) / np.vdot(u,u))*u
    
    def fit_transform(self, X, col_vec = True, normal = True):
        X = X.astype(float)
        if col_vec:
            mat = X.copy()
        else:
            mat = (X.T).copy()
        
        N = mat.shape[1]
        mat_orth = np.array([]).reshape(mat.shape[0], -1)
        for n in range(N):
            u = mat[:, n:n+1].copy()
            if n ==0:
                mat_orth = np.hstack((mat_orth,u))
            else:
                for i in range(n):
                    u -= self._proj(mat_orth[:, i:i+1], mat[:, n:n+1])
                mat_orth = np.hstack((mat_orth,u))
        
        if normal:
            result = mat_orth / np.linalg.norm(mat_orth, axis=0)
            if col_vec:
                return result
            else:
                return result.T
        else:
            if col_vec:
                return mat_orth
            else:
                return mat_orth.T

Gram-schmidt process를 통해 임의의 두 벡터를 정규 직교화해 봅시다.

col_vecs = np.array([[2,1], [1,2]])
gs= gram_schmidt()
col_vecs_orthonorm = gs.fit_transform(col_vecs)
x = np.arange(-1,1,.00001)
fig, ax = plt.subplots(1,1, figsize=(10,10))
ax.set_aspect(1)
ax.plot(np.vstack(([0,0], col_vecs[:,0]))[:,0], np.vstack(([0,0], col_vecs[:,0]))[:,1],  label='v1')
ax.plot(np.vstack(([0,0], col_vecs[:,1]))[:,0], np.vstack(([0,0], col_vecs[:,1]))[:,1], label='v2')
ax.scatter(x, np.sqrt(1-x**2), s=.1, color='k', label='half unit circle')

ax.plot(np.vstack(([0,0], col_vecs_orthonorm[:,0]))[:,0], np.vstack(([0,0], col_vecs_orthonorm[:,0]))[:,1], label='v1_orthonorm', linewidth = 4)
ax.plot(np.vstack(([0,0], col_vecs_orthonorm[:,1]))[:,0], np.vstack(([0,0], col_vecs_orthonorm[:,1]))[:,1], label='v2_orthonorm', linewidth = 4)
plt.legend(fontsize=12)
plt.show()

png

검정 실선은 반지름이 1인 반원입니다. 임의의 벡터 $\mathbf{v}_1$과 $\mathbf{v}_2$가 각각 녹색과 빨간색 벡터로 정규 직교화된 것을 확인할 수 있습니다.


reduced SVD

SVD는 $m \times n$ 행렬을 세 가지 성분으로 분해하는 기법입니다. Full SVD와 다르게 reduced SVD는 좌측 특이벡터 $U$를 non-square 행렬로 계산하여 계산량을 단축합니다.

  • Reduced(thin) SVD

$\large X_{n\times p} = U\; S \; V^* $

$\large \displaystyle U_{n\times p} ={ \begin{pmatrix} \frac{X V_1}{s_1} &\frac{X V_2}{s_2} &\cdots & \frac{X V_p}{s_p} \end{pmatrix} } = X V S^{-1}$

$\large S_{p\times p}$는 $X^* X$ 행렬의 0보다 큰 고윳값을 대각 성분으로 하는 대각 행렬입니다. 대각 성분은 크기의 내림차순으로 정렬됩니다.

$\large {V}_{p\times p}$는 $X^*X$ 행렬의 정규 직교 고유벡터를 컬럼으로 하는 행렬입니다.

순서상으로 $S$와 $V$ 행렬을 구한 후 $U$ 행렬을 계산합니다.

위 알고리즘을 numpy로 구현하면 아래와 같습니다.

def reduced_svd(data, ):
    dim = data.shape[1]
    rank = np.linalg.matrix_rank(data)
    relu = np.vectorize(lambda x: np.real(x) if np.real(x)>=0 else .0)
    
    #eval_u, evec_u = np.linalg.eig(data.dot(data.T), )
    eval_v, evec_v = np.linalg.eig(data.T.dot(data), )
    
    gs = gram_schmidt()
    #evec_u_gs = gs.fit_transform(evec_u)
    evec_v_gs = gs.fit_transform(evec_v)
        
    s = eval_v.copy()
    s = np.sqrt(relu(s))
    s1 = np.sort(s)[::-1]
    if dim > rank:
        s1[-(dim-rank):] = 0
    
    S = np.eye(dim)*s1
    
    v_idx = np.sqrt(relu(np.real(eval_v))).argsort()[-dim:][::-1]
    v = evec_v_gs[:, v_idx ]
    
    u = data.dot(v)/s1
    
    return u, S,  v

앞서 불러온 데이터를 SVD로 분해해보겠습니다. 먼저 numpy 내장 함수를 사용하고 이후 직접 구현한 함수를 적용하여 그 결과를 비교해봅시다.

  • Sample data

캡처2 먼저 샘플 행렬을 Standardization합니다.

$\large \displaystyle \mathbf{X}_{\normalsize standardized} = \begin{pmatrix} \frac{X_1 - \mu_1}{\sigma_1} & \frac{X_2 - \mu_2}{\sigma_2} & \cdots &\frac{X_p - \mu_p}{\sigma_p} \end{pmatrix} $ , $\mu_i$와 $\sigma_i$는 각각 $X_i$($i \in {1,2,\cdots,p}$)의 평균과 표준편차.

scaler = StandardScaler()
data_x_std = scaler.fit_transform(data_x)

numpy로 아래와 같이 계산할 수 있습니다. s는 1차원 array이므로 각 elements를 대각 성분으로 갖는 대각 행렬을 생성합니다. 세 번째로 return되는 행렬은 $V$가 아니라 $V^*$임에 주의합니다.

u, s, vh = np.linalg.svd(data_x_std, full_matrices=False, )
S = np.zeros((data_x_std.shape[1], data_x_std.shape[1]), float)
np.fill_diagonal(S, s)

계산 결과 $\large X \simeq U S V$임을 얻습니다.

np.allclose(np.dot(u, np.dot(S, vh)), data_x_std )
True

이번에는 제가 구현한 코드를 검증해봅시다. 아래와 같이 제대로 계산되는 것을 볼 수 있습니다.

uu, ss, vv = reduced_svd(data_x_std)
np.allclose(np.dot(uu, np.dot(ss, vv.T)), data_x_std )
True

$\large \therefore X \simeq U S V$

이번 시간에는 Gram-schmidt process와 SVD에 대해 알아보았습니다. 다음 포스트에서는 SVD를 이용한 머신러닝 기법에 대해 다루겠습니다.

source of teaser


Scroll to Top

Leave a comment