# 體驗主成份分析

This work by Jephian Lin is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/).

$\newcommand{\trans}{^\top}
\newcommand{\adj}{^{\rm adj}}
\newcommand{\cof}{^{\rm cof}}
\newcommand{\inp}[2]{\left\langle#1,#2\right\rangle}
\newcommand{\dunion}{\mathbin{\dot\cup}}
\newcommand{\bzero}{\mathbf{0}}
\newcommand{\bone}{\mathbf{1}}
\newcommand{\ba}{\mathbf{a}}
\newcommand{\bb}{\mathbf{b}}
\newcommand{\bc}{\mathbf{c}}
\newcommand{\bd}{\mathbf{d}}
\newcommand{\be}{\mathbf{e}}
\newcommand{\bh}{\mathbf{h}}
\newcommand{\bp}{\mathbf{p}}
\newcommand{\bq}{\mathbf{q}}
\newcommand{\br}{\mathbf{r}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\bz}{\mathbf{z}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bv}{\mathbf{v}}
\newcommand{\bw}{\mathbf{w}}
\newcommand{\tr}{\operatorname{tr}}
\newcommand{\nul}{\operatorname{null}}
\newcommand{\rank}{\operatorname{rank}}
%\newcommand{\ker}{\operatorname{ker}}
\newcommand{\range}{\operatorname{range}}
\newcommand{\Col}{\operatorname{Col}}
\newcommand{\Row}{\operatorname{Row}}
\newcommand{\spec}{\operatorname{spec}}
\newcommand{\vspan}{\operatorname{span}}
\newcommand{\Vol}{\operatorname{Vol}}
\newcommand{\sgn}{\operatorname{sgn}}
\newcommand{\idmap}{\operatorname{id}}
\newcommand{\am}{\operatorname{am}}
\newcommand{\gm}{\operatorname{gm}}$

from lingeo import random_int_list

## Main idea

In statistics, data is often stored as a matrix (two-dimensional table) such that
each row represent a sample, while
each column represent a feature.
For example, the rows can be the students and the columns can be the scores in different subjects.

Let $X$ be an $N\times d$ matrix for some data and $\bx_1,\ldots,\bx_N\in\mathbb{R}^d$ its rows.
Thus,

$$
\mu = \frac{1}{N} \sum_{i = 1}^N \bx_i
$$

is the center of these data points.

One may **center** the data by creating a new matrix $X_0$ whose rows are $\bx_1 - \mu, \ldots, \bx_N - \mu$.
Let $J$ be the $N\times N$ all-ones matrix.
Then one way to describe this is $X_0 = (I - \frac{1}{N}J)X$.

The **principal component analysis** (or **PCA**) is
a method of finding the most important directions of the data.
These directions (vectors) are called the **principal components** of the data.

##### Algorithm (principal component analysis)

Input: a data represented by an $N\times d$ matrix $X$ and the desired number $k$ of principal components
Output: the principal components represented by the column vectors of a matrix $P$

1. Let $J$ be the $N\times N$ all-ones matrix.
Define $X_0 = (I - \frac{1}{n}J)X$.
2. Let $U\Sigma V\trans$ be the singular value decomposition of $X_0$.
3. Let $P$ be the $d\times k$ matrix obtained by the first $k$ columns of $V$.

Once the principal components are found.
One may project all data points onto the plane spanned by the principal components.
This can be done by the $d\times k$ matrix $Y = X_0P$,
whose rows are the new data points contains the essential information of the original data.

## Side stories

- hidden words

## Experiments

##### Exercise 1

執行以下程式碼。

### code
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

np.random.seed(0)
print_ans = False

mu = np.random.randint(-3,4, (2,))
cov = np.array([[1,1.9],
                [1.9,4]])
X = np.random.multivariate_normal(mu, cov, (20,))
mu = X.mean(axis=0)
X0 = X - mu
u,s,vh = np.linalg.svd(X0)
P = vh.T[:,:1]

plt.axis('equal')
plt.scatter(*X.T)
plt.scatter(*mu, c="red")
plt.arrow(*mu, *(P[:,0]), head_width=0.3, color="red")

xs = X0.dot(P)
ys = np.zeros_like(xs)
plt.scatter(xs, ys)

pretty_print(LatexExpr(r"V^\top ="))
print(vh)

if print_ans:
    print("red point: center")
    print("red arrow: first principal component")
    print("orange points: projection of blue points onto the red arrow")
    print("red arrow = first row of VT = first column of V")

##### Exercise 1(a)

若藍點為給定的資料。
改變不同的 `seed`，說明紅點、紅箭頭、橘點的意思。

當 `seed = 0` 時可得到以下題目資訊。

**Ans:**

:::warning

- [x] 標點
- [x] 藍點的常態分布 --> 藍點的分布趨勢
:::

紅點為藍點們的中心點，\
紅箭頭為藍點的分布趨勢，\
橘點為藍點以紅點為中心紅箭頭為基軸投影在紅箭頭上的點。

##### Exercise 1(b)

令 $X_0$ 的列向量為將藍點資料置中過後的資料點。
若 $X_0 = U\Sigma V\trans$ 為其奇異值分解。
說明紅箭頭與 $V\trans$ 的關係。

:::warning

- [x] $V\trans$轉置後的$columns$為紅箭頭之向量 --> $V\trans$ 的列為紅箭頭之向量。
:::

**Ans:**\
$V\trans$ 的第一列為紅箭頭之向量。

## Exercises

##### Exercise 2

令 $X$ 為一 $N\times d$ 矩陣，而其列向量為一筆資料。
令 $J$ 為 $N\times N$ 的全一矩陣。

##### Exercise 2(a)

說明 $\frac{1}{N}JX$ 每列的意義。

**Ans:**

因為 $J$ 為一全一矩陣，故 $\frac{1}{N}JX$ 每一列是相同的。
觀察 $\frac{1}{N}JX$ ，會發現其每一列為 $X$ 的所有列相加後除以 $N$ 。
類似於算術平均數。

##### Exercise 2(b)

說明 $(I - \frac{1}{N}J)X$ 每列的意義。

**Ans:**

令
$$
X = \begin{bmatrix}
- & \br_1 & - \\
~ & \vdots & ~ \\
- & \br_n & -
\end{bmatrix}
$$
觀察 $(I - \frac{1}{N}J)X$ ，會發現其第 $i$ 列為，以 $i=1$ 為例
$$
\frac{N-1}{N}\br_1 - \frac{1}{N}(\br_2 + \br_3 + ... + \br_n) =
\frac{N}{N}\br_1 - \frac{1}{N}(\br_1 + \br_2 + ... + \br_n).
$$
故 $(I - \frac{1}{N}J)X$ 的第 $i$ 列為 $X$ 的第 $i$ 列減掉 $X$ 每一列的平均值。
相當於一個置中的動作。

##### Exercise 3

執行以下程式碼。

### code

import numpy as np
np.random.seed(0)
X = np.random.randint(-3, 4, (5,2))
mu = X.mean(axis=0)
X0 = X - mu
print("X =\n", X)
print("mu =\n", mu)
print("X0 =\n", X0)

##### Exercise 3(a)

解釋 `mu = X.mean(axis=0)` 的意義。

**Ans:**
前方的程式碼將此 $X$ 設定為 $5\times2$ 的矩陣，其中，$5$ 個列都各為一個樣本。
那麼 `mu = X.mean(axis=0)` 可以說是，將 `mu` 設為 $X$ 中 $5$ 個樣本的算術平均數，另外 `mu` 為 $1\times2$ 的矩陣。

##### Exercise 3(b)

解釋 `X0 = X - mu` 的意義。

**Ans:**
將 $X_0$ 的值設為 $X$ 的每個列分別扣 `mu` 的矩陣。

##### Exercise 4

執行以下程式碼。

### code

import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

np.random.seed(0)

mu = np.random.randint(-3,4, (3,))
cov = np.array([[3, -0.9, -1.9],
                [-0.9, 11, -9.9], 
                [-1.9, -9.9, 12]])
X = np.random.multivariate_normal(mu, cov, (100,))
mu = X.mean(axis=0)
X0 = X - mu
u,s,vh = np.linalg.svd(X0)
P = vh.T[:,:2]

ax = plt.axes(projection='3d')
ax.set_xlim(-5,5)
ax.set_ylim(-5,5)
ax.set_zlim(-5,5)

ax.scatter(*X.T)
ax.scatter(*mu, c="red")
ax.quiver(*(mu[:,np.newaxis] + np.zeros((2,))), *P, color="red")

%matplotlib inline
Y = X0.dot(P)
plt.axis('equal')
plt.scatter(*Y.T)

##### Exercise 4(a)

說明兩段程式碼產出的圖的關係。

**Ans:**

第一張圖為在三維空間中尚未整理過的原始資料

第二張則是透過主成份分析將資料投影至兩紅色向量張出的二維空間


##### Exercise 4(b)

使用 `X.shape` 來觀察一個陣列的形狀、
並逐句解釋下方程式碼所做的事情。

mu = X.mean(axis=0)
X0 = X - mu
u,s,vh = np.linalg.svd(X0)
P = vh.T[:,:2]

**Ans:**

:::warning

- [x] shape of X = (1215, 100) --> `shape of X = (1215, 100)`，其它程式碼如 mu 也改成 `mu`；英數遇到中文前後要留空格
:::

透過執行以下代碼

print("shape of X =", X.shape)

得到 `shape of X = (1215, 100)` ,
可知形狀為在100維空間中的1215個點，

mu = X.mean(axis=0)

對陣列`X`各行求算數平均值算數並帶入`mu` ,

X0 = X - mu

將陣列`X`各行中的每一項減去當行的平均數`mu`後代入`X0`,相當於將圖形置中，

u,s,vh = np.linalg.svd(X0)

對`X0`進行奇異值分解，

P = vh.T[:,:2]

取出`vh`中的前兩項。

##### Exercise 5

**挑戰**
以下程式碼讀進一筆資料。
這筆資料實際上是貼在高維度中的一個二維平面，
且它在這個平面上排出一個一個英文字。

請找出這個英文字是什麼。

X = np.genfromtxt('hidden_text.csv', delimiter=',')
print("shape of X =", X.shape)

**Ans:**

透過執行以下代碼

import numpy as np
import matplotlib.pyplot as plt

X = np.genfromtxt('hidden_text.csv', delimiter=',')
print("shape of X =", X.shape)
mu = X.mean(axis=0)
X0 = X - mu
u,s,vh = np.linalg.svd(X0)
print(u.shape,s.shape,vh.shape)
P = vh.T[:,:2]
Y = X0.dot(P)
plt.axis('equal')
plt.scatter(*Y.T)

可得出圖形為旋轉 180° 的 "Taiwan" (?)


:::success
Nice work
:::

##### Exercise 6

令 $X_0$ 及 $P$ 為主成份分析演算法中的矩陣。
說明 $X_0P$ 與 $X_0PP\trans$ 的列向量所代量的意義。

:::warning

- [x] $V^\top$ 的前 $k$ 個向量 --> $V$ 的前 $k$ 個行向量
- [x] 則此向量即為 ${\bf x}_i$ 向量本身。--> 則此向量即為 ${\bf x}_i$ 向量在 $\vspan\{\bv_1,\ldots, \bv_k\}$ 的投影。
:::

Ans:
$X_0$ 為一 $N \times D$ 矩陣，將其寫為
$$
X_0 = \begin{bmatrix} - & {\bf x}_1 & - \\
~ & \vdots & ~ \\ - & {\bf x}_N & - \\
\end{bmatrix}
$$
$P$ 則為 $A$ 的奇異值分解 $A = U\Sigma V^\top$ 中，$V$ 的前 $k$ 個行向量所組成的 $d \times k$ 矩陣，將其寫為
$$
P = \begin{bmatrix} | & & | \\
{\bf v}_1 & \cdots & {\bf v}_k \\
| & & | \\
\end{bmatrix}
$$
則
$$
X_0P =
\begin{bmatrix} - & {\bf x}_1 & - \\
~ & \vdots & ~ \\ - & {\bf x}_N & - \\
\end{bmatrix}
\begin{bmatrix} | & & | \\
{\bf v}_1 & \cdots & {\bf v}_k \\
| & & | \\
\end{bmatrix}=
\begin{bmatrix}
\langle {\bf x}_1, {\bf v}_1 \rangle & \langle {\bf x}_1, {\bf v}_2 \rangle & \cdots & \langle {\bf x}_1, {\bf v}_k \rangle \\
\vdots & \vdots & & \vdots \\
\langle {\bf x}_i, {\bf v}_1 \rangle & \langle {\bf x}_i, {\bf v}_2 \rangle & \cdots & \langle {\bf x}_i, {\bf v}_k \rangle \\
\vdots & \vdots & & \vdots \\
\langle {\bf x}_N, {\bf v}_1 \rangle & \langle {\bf x}_N, {\bf v}_2 \rangle & \cdots & \langle {\bf x}_N, {\bf v}_k \rangle \\
\end{bmatrix}
$$
因此，$X_0P$ 的第 $i$ 列的 $k$ 個元素就紀錄了 ${\bf x}_i$ 向量與 $k$ 個主成分向量的內積，即紀錄了其在 $k$ 個主成分向量上的分量。
$X_0PP\trans$ 的第 $i$ 列則為
$$
\begin{bmatrix}
\langle {\bf x}_i, {\bf v}_1 \rangle & \langle {\bf x}_i, {\bf v}_2 \rangle & \cdots & \langle {\bf x}_i, {\bf v}_k \rangle \\
\end{bmatrix}
\begin{bmatrix} - & {\bf v}_1 & - \\
~ & \vdots & ~ \\ - & {\bf v}_k & - \\
\end{bmatrix}=
\langle {\bf x}_i, {\bf v}_1 \rangle {\bf v}_1 + \cdots + \langle {\bf x}_i, {\bf v}_k \rangle {\bf v}_k
$$
若 ${\bf v}_1, \ldots , {\bf v}_k$ 皆為單位向量，則此向量即為 ${\bf x}_i$ 向量在 $\vspan\{\bv_1,\ldots, \bv_k\}$ 的投影。

:::info
目前分數 = 6 × 檢討 = 7
:::