# 7-3 代码实现主成分分析问题

```python
import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100, 2))
X[:,0] = np.random.uniform(0., 100, size=100)
X[:,1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 10., size=100)

plt.scatter(X[:,0], X[:,1])
plt.show()
```

![](http://windmissing.github.io/images/2019/100.png)

## 第一步：demean

```python
def demean(X):
    return X - np.mean(X, axis=0)

X_demean = demean(X)
plt.scatter(X_demean[:,0], X_demean[:,1])
plt.show()
```

![](http://windmissing.github.io/images/2019/101.png)

## 第二步：梯度上升法

```python
def f(w, X):
    return np.sum((X.dot(w)**2)) / len(X)

def df_math(w, X):
    return X.T.dot(X.dot(w)) * 2. / len(X)

def df_debug(w, X, epsilon=0.0001):
    res = np.empty(len(w))
    for i in range(len(w)):
        w_1 = w.copy()
        w_1[i] += epsilon
        w_2 = w.copy()
        w_2[i] -= epsilon
        res[i] = (f(w_1, X) - f(w_2, X)) / (2 * epsilon)
    return res

# 把向量单位化
def direction(w):
    return w / np.linalg.norm(w)

def gradient_ascent(df, X, initial_w, eta, n_iters=1e4, epsilon=1e-8):
    w = direction(initial_w)
    cur_iter = 0
    while cur_iter < n_iters:
        gradient = df(w, X)
        last_w = w
        w = w + eta * gradient
        w = direction(w)
        if(abs(f(w, X)) - abs(f(last_w, X)) < epsilon):
           break
        cur_iter += 1
    return w
```

**注意1：epsilon取值比较小，因为w是方向向量，它的每个维度都很小，所以epsilon也要取很小的值**\
**注意2：每次计算出w后要对其单位化**\
如果每次计算出w后不做单位化的工作，算法也可以工作，因为w本身也是代方向的。\
但这样会导致搜索过程不顺畅。\
因为如果不做单位化，w应该是公式要求的w偏大的，这就要求eta值非常小。\
而eta值小又会导致循环次数非常多，性能就会下降。\
因此遵循公式的假设条件，每次都让w成为方向向量。

## 训练和绘制结果

```python
initial_w = np.random.random(X.shape[1])
eta = 0.001
gradient_ascent(df_debug, X_demean, initial_w, eta)
w = gradient_ascent(df_math, X_demean, initial_w, eta)
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.plot([0, w[0]*30], [0, w[1]*30], color='r')
plt.show()
```

![](http://windmissing.github.io/images/2019/102.png)\
**注意3：w不能是零向量。因为w=0本身也是在极值点上，是极小值点，此时梯度也会0**\
**注意4：不能使用StandardScaler标准化数据。**\
因为本算法的目标就是让方差最大。\
一但对数据做了标准化，样本的方差就肯定是1了，不存在方差最大值。

## 另一个更极端的例子

```python
X2 = np.empty((100, 2))
X2[:,0] = np.random.uniform(0., 100, size=100)
X2[:,1] = 0.75 * X2[:, 0] + 3.

X2_demean = demean(X2)
w2 = gradient_ascent(df_debug, X2_demean, initial_w, eta)

plt.scatter(X2_demean[:, 0], X2_demean[:, 1])
plt.plot([0, w2[0]*30], [0, w2[1]*30], color='r')
plt.show()
```

![](http://windmissing.github.io/images/2019/103.png)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://windmising.gitbook.io/liu-yu-bo-play-with-machine-learning/src/chapter7-3/7-3.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.