In [1]:

import numpy as np
import pandas as pd
from numpy import linalg as la
from copy import deepcopy

np.set_printoptions(precision=3)

PCA

PCA의 목적은 한마디로 말하면, 고차원의 원본 데이터를 저차원으로 압축하고 싶은데, 데이터의 분산(데이터의 특성)을 최대한 보존하여 낮추고자 할때 사용하는 방법이다.

좀 더 구체적으로 말하면, 데이터의 분산 특성인 공분산행렬로 데이터분포를 가장 설명할 수있는 eigen vector와 eigen value들을 찾는 것이 가능하다.

데이터 행렬 에 대한 covariance(공분산행렬)이 데이터의 분산 특성을 나타내는 행렬인데, 이 공분산행렬을 선형대수적 특성을 이용하여 데이터의 분포를 가장 잘 표현하는 $n$개(데이터의 차원 수)의 eigen vector와 eign value로 표현 할 수 있다.

그런데, 데이터의 분산 특성을 가장 잘 드러내는 (eigen value가 가장 큰) $k$개의 eigen vector들을 골라($k$ 개의 principal axis) $n$개의 bases으로 표현할 수 있는 $n$차원 안의 모든 데이터를 $k$ 개의 bases를 이루는 $k$ 차원의 데이터 공간으로 Projection 함으로써 압축한다.

eigen vector의 의미는 어떤 matrix $\mathbf{A}$ 를 통해 vector $\mathbf{x}$ 를 linear transform 했을때, 방향은 변하지 않고, 크기만 변하는 벡터를 의미한다. 이때 변하는 크기가 eigen value를 의미한다.
다음 그림은 이해를 돕기위해 $\mathbf{A}\mathbf{x} = \mathbf{b}$, where $\mathbf{A} \in \mathbb{R}^{m \times n}$에서의 space정보들이다.
eigen vector $\mathbf{e}$는 $\mathbf{A}\mathbf{e}$와 같이 linear transform 되어도, 크기만 변하는 $\mathbf{A} \mathbf{e} = \lambda \mathbf{e}$ 와 같은 특성을 지닌다. 따라서, 위의 그림에서 column space의 한 vector인 $\mathbf{e}$는 row space의 $\mathbf{A} \mathbf{e}$로 transform되어도 방향이 그대로 보존 된다.

Proof

정말로 $n$ 차원의 데이터를 $k$개의 eigen vector들을 골라 압축하는게 최선으로 데이터의 (분산) 특성을 보존하는 것일까? 수학적으로 증명해보자.

name	<div style="width:290px">notation</div>	dimension
data matrix	$\mathbf{X}$	$\mathbb{R}^{m \times n}$
eigen vectors	$\mathbf{E} = [\mathbf{e}_1, \mathbf{e}_2, …, \mathbf{e}_k]$	$\mathbb{R}^{n \times k}$
eigen values	$\boldsymbol{\lambda}_{[k \times k]} = \begin{bmatrix} \lambda_1 & 0 & … & 0 \\0 & \lambda_2 & … & 0 \\0 & 0 & … & 0 \\… & … & … & … \\0 & 0 & … & \lambda_k \\ \end{bmatrix}$	$\mathbb{R}^{k \times k}$

일단, $n$ 차원의 $m$ 개의 데이터들 들을 $k$차원으로 압축했다고 생각하자.
그러면, $n$ 차원에서 $k$차원으로 Projection 하는 $\mathbf{V} \in \mathbb{R}^{n \times k} = [\mathbf{v}_1, \mathbf{v}_1, …, \mathbf{v}_k]$를 생각할 수 있다.
(여기서, $\Vert \mathbf{v}_i \Vert = 1, \forall i$ 이다.)

편의상 $\mathbb{E}(\mathbf{X}) = 0$으로 가정한 상황에서
원본 데이터의 covariance 를 보존한다는 것은
Projection 된 $\mathbf{X} \mathbf{V} \in \mathbb{R}^{m \times k}$의 covariance 인 $cov(\mathbf{X}\mathbf{V}) \in \mathbb{R}^{k \times k}$를 최대화하는 말과 같다.
따라서, 정리하면 문제는 다음과 같다.

\[\begin{align} \text{maximize } cov(\mathbf{X}\mathbf{V})& \\ cov(\mathbf{X}\mathbf{V}) &=\frac{1}{m - 1}(\mathbf{X}\mathbf{V})^T (\mathbf{X}\mathbf{V})\\ &= \mathbf{V}^T (\frac{\mathbf{X}^T \mathbf{X}}{m - 1}) \mathbf{V} \\ &= \mathbf{V}^T \Sigma \mathbf{V} \end{align}\]

이 문제에서 $cov(\mathbf{X}\mathbf{V})$는 $\Vert \mathbf{V} \Vert$가 클수록 값이 크므로 제약사항으로 다음과 같은 조건이 추가된다. $\Vert \mathbf{V} \Vert = 1$

제약사항들(constraints)이 있는 최소/최대 문제를 풀때 사용하는 Lagrange Multiplier Method(라그랑주 승수법)을 이용하면 문제는 다음과 같이 풀린다. $\begin{align} L &= \mathbf{V}^T \Sigma \mathbf{V} - \boldsymbol{\lambda}_{[k\times k]}(\Vert \mathbf{V} \Vert - 1)\\ \frac{\partial L}{\partial V} &= 2\Sigma \mathbf{V} - 2 \boldsymbol{\lambda}{[k\times k]} \mathbf{V} = 0 \\ \Sigma \mathbf{V} &= \boldsymbol{\lambda}_{[k \times k]} \mathbf{V} \end{align}$

위의 조건(eigen vector의 정의)을 만족하는 $\mathbf{V}$는 $k$개의 eigen vectors이며, $\lambda_i \forall i$값들은 eignen value가 됨을 알 수 있다. $\therefore \mathbf{V} = \mathbf{E}$ 또한, $\lambda_i$ 가 클수록 $\mathbf{X} \mathbf{e_i}$ 가 크기 때문에, top $k$ 개의 $\lambda_{1, 2, …, k}$ 를 선택하고, 이에 대응되는 $\mathbf{e}_{1, 2, …, k}$를 선택하면 된다.

최종적으로 $\mathbf{e}_{1, 2, …, k}$ 가 span 하는 공간으로 데이터 샘플들을 Projection 하는 것이 원본 분산을 최대로 보존하는 방법이 된다.

Step 1. Centering

feature, 즉 변수(행)별로 평균을 $0$으로 centering한 행렬 $X$를 만듬
(좌표계의 원점이 평균 벡터와 일치하도록 만듬)

$m$: the number of data(sample)
$n$: the dimension of each datagen value를 의미한다.

In [2]:

m, n = 5, 3
X = np.random.randint(10, size = (m, n)).astype('float64')
print("* 원본 data")
print(X)

* 원본 data
[[5. 2. 8.]
 [8. 3. 2.]
 [1. 3. 6.]
 [7. 3. 3.]
 [4. 1. 2.]]

In [3]:

print("* 각 feature 별 평균: {}".format(np.mean(X, axis=0)))
X -= np.mean(X, axis=0)
print(X)

* 각 feature 별 평균: [5.  2.4 4.2]
[[ 0.  -0.4  3.8]
 [ 3.   0.6 -2.2]
 [-4.   0.6  1.8]
 [ 2.   0.6 -1.2]
 [-1.  -1.4 -2.2]]

In [4]:

np.matmul(X.T, X) / (m - 1)  # covarince of X

array([[ 7.5 ,  0.5 , -3.5 ],
       [ 0.5 ,  0.8 ,  0.15],
       [-3.5 ,  0.15,  7.2 ]])

Step 2. Find Covariance Matrix

공분산행렬은 다음 식으로 만들 수 있음, 여기서 $m - 1$로 나눈 이유는 $m$개의 데이터로 이루어진 샘플인 표본 분산에 대한 공분산 행렬이기 때문이다.
($m - 1$ degree of freedom, 왜냐하면 마지막 하나의 데이터는 모평균으로 맞추기 위해) $\Sigma = cov(X)= \frac{X^T X} {m − 1} \propto X^T X$

In [5]:

C = np.cov(X, rowvar=False)
print(C)

[[ 7.5   0.5  -3.5 ]
 [ 0.5   0.8   0.15]
 [-3.5   0.15  7.2 ]]

Step 3. Find Eigen Values and Vectors

공분산행렬(C)을 기반으로 고유값(l)과 고유벡터(principal_axes) 구하기

In [6]:

l, principal_axes = la.eig(C)
print(l)
print(principal_axes)

[10.86   3.913  0.727]
[[ 0.723  0.681 -0.118]
 [ 0.026  0.144  0.989]
 [-0.69   0.718 -0.087]]

고유값을 높은 순으로 정렬하고, 이에 대응하는 고유벡터와 순서를 맞춤

In [7]:

idx = l.argsort()[::-1]
idx

array([0, 1, 2])

In [8]:

l, principal_axes = l[idx], principal_axes[:, idx]
print(l)
print(principal_axes)

[10.86   3.913  0.727]
[[ 0.723  0.681 -0.118]
 [ 0.026  0.144  0.989]
 [-0.69   0.718 -0.087]]

In [9]:

# np.matmul(principal_axes.T, principal_axes) 

Step 4. Dimensionality Reduction

차원축소 예 (고유값을 기준으로 가장 큰 $k$개의 고유 벡터 선택)

principal axis(principal_compoents) 를 구함

In [10]:

# k = 2 
print(principal_axes[:, :2])

[[ 0.723  0.681]
 [ 0.026  0.144]
 [-0.69   0.718]]

principal axis 를 기반으로 원본데이터 $X$에 대한 principal components 를 구함
그리고 top 2 개의 components 들에 대해서 dimensionality reduction을 수행

In [11]:

mapped_data = X.dot(principal_axes)
print(mapped_data)

[[-2.634  2.672 -0.725]
 [ 3.703  0.548  0.431]
 [-4.119 -1.344  0.909]
 [ 2.29   0.586  0.462]
 [ 0.76  -2.462 -1.077]]

In [12]:

mapped_data_reduced = X.dot(principal_axes[:, :2])
print(mapped_data_reduced)

[[-2.634  2.672]
 [ 3.703  0.548]
 [-4.119 -1.344]
 [ 2.29   0.586]
 [ 0.76  -2.462]]

In [13]:

# 새로운 변수 z1 과 z2 (5개의 데이터에대한 각각의 dimension component를 갖고있는 vector)
z1 = mapped_data_reduced[:,0]
z2 = mapped_data_reduced[:,1]
print(z1)
print(z2)

[-2.634  3.703 -4.119  2.29   0.76 ]
[ 2.672  0.548 -1.344  0.586 -2.462]

Example 2

In [14]:

n, d = 4, 2
X = np.array([[0,-2],[3,-3],[2,0],[-1,1]]).astype('float64')
print(X)

[[ 0. -2.]
 [ 3. -3.]
 [ 2.  0.]
 [-1.  1.]]

In [15]:

# 각 row마다 빼진다. 알아서 broadcasting 되어 계산됨
X -= np.mean(X, axis=0)
print(X)

[[-1. -1.]
 [ 2. -2.]
 [ 1.  1.]
 [-2.  2.]]

\[\Sigma =cov(X)= \frac{X^T X}{m−1}\]

Covariance Matrix C(or $\Sigma$) can be computed by $\frac{X^TX}{m - 1}$

In [16]:

print(np.matmul(X.T,X) / (n - 1))

[[ 3.333 -2.   ]
 [-2.     3.333]]

In [17]:

C = np.cov(X, rowvar=False)
print(C)

[[ 3.333 -2.   ]
 [-2.     3.333]]

In [18]:

l, principal_axes = la.eig(C)
print(l)
print(principal_axes)

[5.333 1.333]
[[ 0.707  0.707]
 [-0.707  0.707]]

In [19]:

idx = l.argsort()[::-1]
l, principal_axes = l[idx], principal_axes[:, idx]
print(l)
print(principal_axes)

[5.333 1.333]
[[ 0.707  0.707]
 [-0.707  0.707]]

In [20]:

mapped_data = X.dot(principal_axes)
print(mapped_data)

[[-2.220e-16 -1.414e+00]
 [ 2.828e+00 -4.441e-16]
 [ 2.220e-16  1.414e+00]
 [-2.828e+00  4.441e-16]]

SVD

SVD(Sigular Value Decomposition) means that
$A$ can be decomposed by $U \in \mathbb{R}^{m \times k}, \Sigma \in \mathbb{R}^{k \times k}, V \in \mathbb{R}^{n \times k}$ mathmatically, where $k \le d $.

$\begin{align} A &= U \Sigma V^T \\ &= \begin{pmatrix} \mid & \mid & {} & \mid \\\\ \vec u_1 & \vec u_2 &\cdots &\vec u_k \\\\ \mid & \mid & {} & \mid \end{pmatrix} \begin{pmatrix} \sigma_1 & & & \\ & \sigma_2 & & \\ & & \ddots & \\ & & & \sigma_k \end{pmatrix} \begin{pmatrix} ㅡ & \vec v^T_1 & ㅡ \\ ㅡ & \vec v^T_2 & ㅡ \\ &\vdots& \\ ㅡ & \vec v^T_k & ㅡ \end{pmatrix} \\ &= \sigma_1 \vec u_1 \vec v_1^T + \sigma_2 \vec u_2 \vec v_2^T +\cdots+ \sigma_k \vec u_k \vec v_k^T \text{, where } -1 \le \vec{u}_i \vec{v}_i^T \le 1 \end{align}$ 이로써 알 수 있는 점은 $\mathbf{A}$의 원소 값들은 $\vec{u}_i \vec{v}_i^T$ 의 linear combination 에 의해 결정되며 $ \alpha_i$ 에 의해 크기가 결정된다는 점이다.

$AA^T$ 또는 $A^TA$ 은 real symmetric matrix이며 다음과 같은 특성들을 지닌다. 한글 정리된 블로그 1 블로그 2

symmetric 이므로 모든 eigen value 모두 실수(특히 양수).
- real symmetrix 이므로 positive definite. $\therefore \forall \lambda >0 $
모든 eigen vector가 1차 독립(full rank)이므로 대각화가 가능.
- 특히 eigen vector들이 서로 orthogonal 하므로 orthogonally diagonalizable(직교대각화가 가능)하다.

정리하면 다음과 같다.

직교행렬 $\mathbf{Z}$ 는 $\mathbf{Z} \mathbf{Z}^T = \mathbf{I}$ 를 만족 $\because \mathbf{Z}^T = \mathbf{Z}^{-1}$

$U$는 $AA^T$ (즉, $A^T$의 공분산 행렬) 의 고유벡터들 이며 orthogonal matrix(직교행렬)이다.

$V$는 $A^TA$ (즉, $A$의 공분산 행렬) 의 고유벡터들 orthogonal matrix(직교행렬)이다.

$\Sigma \Sigma^T$ 또는 $\Sigma^T \Sigma= \lambda$ 따라서, singular value (특이값) $\sigma = \sqrt{\lambda}$

의미 분석

$\mathbf{A}$ 는 linear trasnformd을 하는 행렬이고, $U, \Sigma, V$로 3개의 linear transform을 하는 행렬로 나눠질 수 있다. 그런데, 기하학적으로 $V^T, U$는 크기는 변하지 않으며 방향 변환을 $\Sigma$는 방향은 변하지 않으며 크기 변환만 수행하는 연산을 의미한다.

또한, 의미적으로 $\mathbf{A}$ 를 $n$개의 feature를 지닌 $m$개의 샘플로 구성된 데이터 행렬이라고 생각했을때, 데이터의 특성(데이터의 분산표현)은 다음과 같다.

행공간에 대한 데이터 특성 $\mathbf{A}^T\mathbf{A}$
열공간에 대한 데이터 특성 $\mathbf{A}\mathbf{A}^T$

이때, 행공간과 열공간은 $k$개의 직교하는 bases vector로 span하며 그 벡터들은 다음과 같다.

$m$차원의 행공간의 데이터 특성에 대한 $k$개의 서로 직교하는 기저 표현을 $U$
$n$차원의 열공간의 데이터 특성에 대한 $k$개의 서로 직교하는 기저 표현을 $V$

또한, 열공간과 행공산 사이에 $k$공간의 연결고리가 있으며, 두 공간의 방향이 일치하는 상황이다.

$k$차원의 행공간과 열공간 사이의 각 차원 scale 차이가 $\Sigma$의 각 원소를 의미

결론적으로 특이값 분해(SVD)가 말하는 것은 다음과 같다.
직교하는 벡터 집합 $V$ 에 대하여, $\Sigma$ 를 이용한 선형 변환 후에 그 크기는 변하지만 여전히 직교할 수 있게 되는 그 직교 집합은 $U$라는 것이다.

In [21]:

user, item = 5, 3
A = np.random.randint(10, size = (user, item)).astype('float')
# user, item = 4, 2
# A = np.array([[2,3],[1,4],[0,0],[0,0]]).astype('float')
# A = np.array([[1,0,0,0,0],[0,0,2,0,3],[0,0,0,0,0],[0,2,0,0,0]])
print("* 원본 data")
print(A)
# print("* 각 dimension 별 평균")
# print(np.mean(A, axis=0))
# A -= np.mean(A, axis=0)
# print(A)

* 원본 data
[[6. 9. 8.]
 [1. 0. 0.]
 [1. 1. 8.]
 [0. 9. 4.]
 [6. 3. 5.]]

Use Libriary

In [22]:

# we now perform singular value decomposition of X
# "economy size" (or "thin") SVD
# U, s, Vt = la.svd(A, full_matrices=True)
U, s, Vt = la.svd(A, full_matrices=False)
V = Vt.T
print(s)
S = np.diag(s)
print(U)
print(V)
print(S)

[18.647  6.666  4.781]
[[-0.718 -0.035  0.26 ]
 [-0.02   0.066  0.17 ]
 [-0.338  0.539 -0.766]
 [-0.455 -0.724 -0.281]
 [-0.403  0.424  0.487]]
[[-0.38   0.441  0.813]
 [-0.649 -0.753  0.105]
 [-0.659  0.488 -0.573]]
[[18.647  0.     0.   ]
 [ 0.     6.666  0.   ]
 [ 0.     0.     4.781]]

In [23]:

print(np.matmul(np.matmul(U, S), V.T))

[[ 6.000e+00  9.000e+00  8.000e+00]
 [ 1.000e+00  1.275e-15  1.495e-15]
 [ 1.000e+00  1.000e+00  8.000e+00]
 [-1.608e-15  9.000e+00  4.000e+00]
 [ 6.000e+00  3.000e+00  5.000e+00]]

Trucated SVD(즉 dimensionality reduction을 함):

In [24]:

k = 2 
np.matmul(np.matmul(U[:,:k], S[:k,:k]), V.T[:k,:])

array([[ 4.99 ,  8.87 ,  8.711],
       [ 0.339, -0.085,  0.466],
       [ 3.978,  1.384,  5.903],
       [ 1.094,  9.141,  3.229],
       [ 4.105,  2.755,  6.334]])

Compute Numerically

recall that …

Note that $AA^T$ 또는 $A^TA$ 은 symmetric matrix이므로, eigen decomposition이 가능
$U$는 $AA^T$ 즉, $A^T$의 공분산 행렬의 (k개의)고유벡터들
$V$는 $A^TA$ 즉, $A$의 공분산 행렬의 (k개의)고유벡터들

\[cov(A)= \frac{A^T A}{m−1} \propto A^T A\]

Step 1. Find Covarance of $A$ and $A^T$

In [25]:

print(np.matmul(A.T, A))
print(np.matmul(A, A.T))
C1 = np.matmul(A.T, A)
C2 = np.matmul(A, A.T)

[[ 74.  73.  86.]
 [ 73. 172. 131.]
 [ 86. 131. 169.]]
[[181.   6.  79. 113. 103.]
 [  6.   1.   1.   0.   6.]
 [ 79.   1.  66.  41.  49.]
 [113.   0.  41.  97.  47.]
 [103.   6.  49.  47.  70.]]

Step 2. Find Eigen Values and Eigen Vectors of $Cov(A)$ or $Cov(A^T)$

$Cov(A)$ 의 고유벡터들을 구하면 $V$가 된다.

In [43]:

# C1 is a symmetric matrix and so it can be diagonalized:
l, principal_axes = la.eig(C1)
# sort results wrt. eigenvalues
idx = l.argsort()[::-1]
l, principal_axes = l[idx], principal_axes[:, idx]
# the eigenvalues in decreasing order
print ("eigenvalues = \n", np.sqrt(l))
# a matrix of eigenvectors (each column is an eigenvector)
V = principal_axes
print ("eigenvectors, which is same with V = \n", V)
# principal_components = A.dot(principal_axes)
# print ("principal_components = \n", principal_components)

eigenvalues = 
 [18.647  6.666  4.781]
eigenvectors, which is same with V = 
 [[-0.38   0.441 -0.813]
 [-0.649 -0.753 -0.105]
 [-0.659  0.488  0.573]]

In [27]:

# orthogonal for each other. very small value
print(sum(V[:, 0] * V[:, 1]))
print(sum(V[:, 0] * V[:, 2]))
print(sum(V[:, 1] * V[:, 2]))

220446049250313e-16
551115123125783e-17
43689570931383e-16

Step 3 Find $U, \Sigma$

$A = U\sum V^T$ 에서 양쪽 수식의 오른쪽에 $V$ 를 곱하면,
$AV = U\sum = [s_1 U_1, s_2 U_2, … ,s_r U_r] $ 인데
각 column마다 대응되는 singular value $s_i\vert_{i=1,2, …, r}$를 나누어주면 $U$ 를 구할 수 있다.

이때, singular value들은 V를 구할때 찾은 eigen values들이 된다.

In [28]:

singulars = np.sqrt(l)
print(singulars)
U = np.matmul(A, V)
# print(U)
for i, sing in enumerate(singulars):
    U[:,i] = U[:,i]/sing
print(U)

[18.647  6.666  4.781]
[[-0.718 -0.035 -0.26 ]
 [-0.02   0.066 -0.17 ]
 [-0.338  0.539  0.766]
 [-0.455 -0.724  0.281]
 [-0.403  0.424 -0.487]]

In [29]:

S = np.diag(singulars)
S

array([[18.647,  0.   ,  0.   ],
       [ 0.   ,  6.666,  0.   ],
       [ 0.   ,  0.   ,  4.781]])

In [30]:

A

array([[6., 9., 8.],
       [1., 0., 0.],
       [1., 1., 8.],
       [0., 9., 4.],
       [6., 3., 5.]])

Sanity Check:Reconstruct $A$ by multiplication of $U, S, V^T$

In [31]:

np.matmul(np.matmul(U, S), V.T)

array([[ 6.000e+00,  9.000e+00,  8.000e+00],
       [ 1.000e+00,  5.963e-16, -3.135e-16],
       [ 1.000e+00,  1.000e+00,  8.000e+00],
       [ 4.682e-15,  9.000e+00,  4.000e+00],
       [ 6.000e+00,  3.000e+00,  5.000e+00]])

일반적으로, 어떤 rating matrix $\mathbf{A}$을 가지고, 추정한 SVD값과의 RMSE는 최소가 됨이 증명 되어있다.
하지만, $\mathbf{A}$에 missing value 가 있을 때에는 SVD를 구할수 없어 prediction 된 $U$, $\Sigma$, $V$를 찾아야한다. 따라서, 이 predication된 SVD 는 latent factor 모델을 이용한 추천시스템에 사용된다.

추천시스템 활용 예제는 프로도의 머릿속 블로그 에 잘 정리되어 있어 링크를 남긴다.

Reference

[1] korean blog - summary
[2] Eigendecomposition of a Matrix - wiki
[3] korean blog - Symmetric Matrix and Spectrum Theory
[4] SVD - wiki
[5] Applying SVD on Item-based Filtering - paper, IEEE 2005
[6] Matrix Factorization 정리
[7] 공돌이의 수학정리노트

Appendix

eigen values, eigen vectors들은 power iteration 을 반복하여 찾는 것이 가능하다.

power iteration 방법을 사용 $v_{t} = \frac{M v_{t-1}} { \lVert Mv_{t-1} \rVert }$

power iteration을 통해 첫번째로 수렴하는 v값이 가장 큰 eigen value와 대응되는 eigen vector이다.
이때, largest eigen value는 $v^T M v$ 이다.
Gram schmidt process에 따라 $M$ 에서 largest igen vector 성분을 다음과 같이 제거하고,
$M_{new} = M - \lambda v v^T$ power iteration을 반복하면 두번째로 큰 eigen value와 대응하는 eigen vector를 찾을 수 있다.

Theory: Gram schmidt process - wiki
Korean blog: eigen values, vectors - darkprogrammer
Another korean blog: eigen values, vectors

In [38]:

M = np.array([[1,1,1],[1,2,3],[1,3,5]]).astype("float")
print(M)

l, principal_axes = la.eig(M)
idx = l.argsort()[::-1]
l, principal_axes = l[idx], principal_axes[:, idx]
print(l) # eigen values
print(principal_axes) # eigen vectors

[[1. 1. 1.]
 [1. 2. 3.]
 [1. 3. 5.]]
[ 7.162e+00  8.377e-01 -2.958e-16]
[[-0.218 -0.886  0.408]
 [-0.522 -0.248 -0.816]
 [-0.825  0.391  0.408]]

power iteration in order to find eigen vector corresponding to the largest eigen value. $v_{t} = \frac{M v_{t-1}} { \lVert Mv_{t-1} \rVert }$

In [39]:

elipslion = 1e-8
v = np.ones((np.size(M, axis=0),1))
iteration = 0
distance = 100
while distance > elipslion:
    print("iteration", iteration, "...")
    prev_v = v
    Mv = np.matmul(M,v)
    v = Mv/la.norm(Mv, 'fro')
    iteration = iteration + 1
    distance = np.linalg.norm(v-prev_v)
    print("distance = ",np.linalg.norm(v-prev_v))
print(v)

iteration 0 ...
distance =  0.890429726059787
iteration 1 ...
distance =  0.04894577225762312
iteration 2 ...
distance =  0.005731964656124627
iteration 3 ...
distance =  0.0006704398749804303
iteration 4 ...
distance =  7.84167576513936e-05
iteration 5 ...
distance =  9.171868646813307e-06
iteration 6 ...
distance =  1.0727703710604327e-06
iteration 7 ...
distance =  1.2547456944453433e-07
iteration 8 ...
distance =  1.4675897110185787e-08
iteration 9 ...
distance =  1.716538713429754e-09
[[0.218]
 [0.522]
 [0.825]]

Find largest eigen value by solving $Mv = \lambda v = v \lambda $ $\lambda = v^T M v$

In [40]:

l = np.matmul(np.matmul(v.T, M),v)
print(l)

[[7.162]]

eliminate principal eigen vector $M_{new} = M - \lambda v v^T$

In [41]:

M = M - np.matmul(l*v,v.T)
print(M)

[[ 0.658  0.184 -0.291]
 [ 0.184  0.051 -0.081]
 [-0.291 -0.081  0.128]]

Find second largest eigen value

power iteration again
then, we can get eigen vector corresponding to 2nd largeset eigen value

In [42]:

elipslion = 1e-8
v = np.ones((np.size(M, axis=0),1))
iteration = 0
distance = 100
while distance > elipslion:
    print("iteration", iteration, "...")
    prev_v = v
    Mv = np.matmul(M,v)
    v = Mv/la.norm(Mv, 'fro')
    iteration = iteration + 1
    distance = np.linalg.norm(v-prev_v)
#     print("distance = ",np.linalg.norm(v-prev_v))
print(v)

iteration 0 ...
iteration 1 ...
[[ 0.886]
 [ 0.248]
 [-0.391]]

Futhermore

MDS(multi-dimensional scaling)
t-SNE(t-stochastic neighbor embeeding)

Share on

Twitter Facebook LinkedIn

	Starwars	Avata	Alien	Casablanca	Titanic
Starwars	1.000000	1.000000	1.000000	-0.748783	-0.748783
Avata	1.000000	1.000000	1.000000	-0.748783	-0.748783
Alien	1.000000	1.000000	1.000000	-0.748783	-0.748783
Casablanca	-0.748783	-0.748783	-0.748783	1.000000	1.000000
Titanic	-0.748783	-0.748783	-0.748783	1.000000	1.000000

	Superman	Batman	Ironman	Amelie
u1	1.0	2.0	8.0	10.0
u2	10.0	NaN	8.0	3.0
u3	8.0	9.0	9.0	2.0
u4	4.0	5.0	9.0	7.0

Dimensionality Reduction PCA, SVD Implementation

sw Yoo

PCA

Proof

Step 1. Centering

Step 2. Find Covariance Matrix

Step 3. Find Eigen Values and Vectors

Step 4. Dimensionality Reduction

Example 2

SVD

의미 분석

Use Libriary

Compute Numerically

Step 1. Find Covarance of $A$ and $A^T$

Step 2. Find Eigen Values and Eigen Vectors of $Cov(A)$ or $Cov(A^T)$

Step 3 Find $U, \Sigma$

추천시스템 활용 예제

Reference

Appendix

Futhermore

Share on

Leave a comment

You may also enjoy

Pybullet Tutorial을 해보자.

Reeds-Shepp (RS) Path

Model Predictive Control (MPC)

RRT star

	Starwars	Avata	Alien	Casablanca	Titanic
u1	1.0	1.0	1.0	0.0	0.0
u2	3.0	3.0	3.0	0.0	0.0
u3	4.0	4.0	4.0	0.0	0.0
u4	5.0	5.0	5.0	0.0	0.0
u5	0.0	0.0	0.0	4.0	4.0
u6	0.0	0.0	0.0	5.0	5.0
u7	0.0	0.0	0.0	2.0	2.0

	u1	u2	u3	u4	u5	u6	u7
u1	1.0	1.0	1.0	1.0	-1.0	-1.0	-1.0
u2	1.0	1.0	1.0	1.0	-1.0	-1.0	-1.0
u3	1.0	1.0	1.0	1.0	-1.0	-1.0	-1.0
u4	1.0	1.0	1.0	1.0	-1.0	-1.0	-1.0
u5	-1.0	-1.0	-1.0	-1.0	1.0	1.0	1.0
u6	-1.0	-1.0	-1.0	-1.0	1.0	1.0	1.0
u7	-1.0	-1.0	-1.0	-1.0	1.0	1.0	1.0

	Starwars	Avata	Alien	Casablanca	Titanic
u1	1.0	1.0	1.0	0.0	0.0
u2	3.0	3.0	3.0	0.0	0.0
u3	4.0	4.0	4.0	0.0	0.0
u4	5.0	5.0	5.0	0.0	0.0
u5	0.0	0.0	0.0	4.0	4.0
u6	0.0	0.0	0.0	5.0	5.0
u7	0.0	0.0	0.0	2.0	2.0

	Starwars	Avata	Alien	Casablanca	Titanic
u1	1.0	1.0	1.0	0.0	0.0
u2	3.0	3.0	3.0	0.0	0.0
u3	4.0	4.0	4.0	0.0	0.0
u4	5.0	5.0	5.0	0.0	0.0
u5	0.0	0.0	0.0	4.0	4.0
u6	0.0	0.0	0.0	5.0	5.0
u7	0.0	0.0	0.0	2.0	2.0