Day 10_PCA_MNIST

# -*- coding: utf-8 -*-
"""Day 14_PCA_MNIST.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1wu4_7BeGAp9YqSR2L1jW5SMigl2GXmmA
"""
 
# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
import gzip, sys, os
 
 
 
if sys.version_info[0] == 2:
    from urllib import urlretrieve
else:
    from urllib.request import urlretrieve
 
def download(filename, source='http://yann.lecun.com/exdb/mnist/'):
    print("Downloading %s" % filename)
    urlretrieve(source + filename, filename)
 
def load_mnist_images(filename):
    if not os.path.exists(filename):
        download(filename)
    
    # 첫 16바이트를 제외하고 60000 * 784 개의 바이트를 읽는다.
    with gzip.open(filename, 'rb') as f:
        data = np.frombuffer( f.read(), np.uint8, offset = 16 )
 
    # 60000행 784열로 리쉐잎. 한개의 숫자 = 784개의 요소(열 하나)로 표현된다
    data = data.reshape(-1,784)
    
    return data / np.float32(256)
 
## MNIST 데이터 로드
train_data = load_mnist_images('train-images-idx3-ubyte.gz')
 
# 공분산 행렬 계산, rowvar 피처간의 분산을 계산, 정규화 된다.
Sigma = np.cov(train_data, rowvar = 0, bias = 0)
 
# 공분산 행렬의 대각선에는 각 피처별 분산이 들어있다.
coordinate_variances = np.sort(Sigma.diagonal())
 
# Symetric인 공분산 행렬에 대한 고유벡터에 데이터를 투영했을때의 분산(고유값)을 저장
eigenvector_variances = np.sort( np.linalg.eigvalsh( Sigma ) )
 
# k개의 고유벡터좌표계에 투영했을때의 전체 분산의 손실비율을 계산
 
# 전체 분산의 합대비 각 피처까지의 분산의 누적합 비율을 계산.
total_coordinate_variance = np.cumsum(coordinate_variances)
total_coordinate_variance = total_coordinate_variance / total_coordinate_variance[783]
 
# 각 고유벡터들(분산)의 누적비율을 계산.
total_eigenvector_variance = np.cumsum(eigenvector_variances)
total_eigenvector_variance = total_eigenvector_variance/total_eigenvector_variance[783]
 
# 분산의 손실 계산 결과 시각화.
# 원래의 좌표계에서 K개의 컬럼만 가지고 계산한 분산의 총합보다
# K개의 고유벡터의 분산의 총합이 훨씬 높다( 잃어버린 분산이 적다 )
 
plt.plot(np.arange(1,784), total_coordinate_variance[784:0:-1], 'b-', lw=2)
plt.plot(np.arange(1,784), total_eigenvector_variance[784:0:-1], 'r-', lw=2)
plt.xlabel('projection dimension', fontsize=14)
plt.ylabel('fraction of residual variance', fontsize=14)
plt.xlim(0,784)
plt.ylim(0.0,1.0)
plt.legend(['coordinate directions', 'PCA directions'], fontsize=14)
plt.show()
 
# 구체적으로 어떤 정보가 손실 되었는지 알아보자.
 
## 고유값 오름차순으로, 고유벡터는 유닛 벡터로 정규화 된 상태. 하나의 컬럼이 하나의 고유벡터이다.
## 각 고유벡터도 고유값과 마찬가지로 왼쪽부터 오른쪽으로 분산에 대한 오름차순으로 정렬 되어있음에 주의한다.
 
eigenvalues, eigenvectors = np.linalg.eigh(Sigma)
 
# U를 784행 k열의 가장 큰 분산(고유값)을 가진 고유벡터 매트릭스라 하고
# U( U.T @ X ) 를 하게 되면 X를 K 차원상의 점으로 표현한 상태 ( U.T @ X ) 에서
# 다시 원래의 차원수를 가지도록 복원할 수 있다.
 
## 784행의 벡터를 k차원으로 줄이고 다시 784행의 벡터로 복원하는 함수.
def projection_and_reconstruction(k):
    U = eigenvectors[:, (784 -k) : 784]
    P = np.dot(U, U.T)
    return P
 
## 784열의 숫자 이미지 데이터 하나를 시각화 하는 함수.
def show_digit(x):
    # 모든 요소가 [0,255] 사이에 존재하도록 조정.
    for i in range(784):
        x[i] = max(0.0, x[i])
        x[i] = min(255.0, x[i])
    
    #시각화
    plt.axis('off')
    plt.imshow(x.reshape((28,28)), cmap=plt.cm.gray,interpolation = 'bilinear')
    plt.show()
    return
 
## 이미지( 한 행)을 받아서 K_list 에 있는 차원 수로 감소시킨 결과를 보여주는 함수
def show_effect_of_PCA(x, k_list):
    print ("원본:")
    show_digit(x)
    for k in k_list:
        if (k > 0) and (k < 784):
            print("투영된 차원의 수 ", k) 
            P = projection_and_reconstruction(k)
            show_digit(P.dot(x))
 
show_effect_of_PCA(train_data[1,], [100, 50, 25, 10])
 

출처 및 참고자료 : edx -  Machine Learning Fundamentals_week_9 Programming Assignment.1