PCA的推导与求解（二）

生信编程日常

发布于 2020-10-29 11:17:04

4490

发布于 2020-10-29 11:17:04

文章被收录于专栏：生物信息学、python、R、linux

我们知道在PCA中是最终要求w，使得以下目标函数最大。：

它其实存在数学解的，在sklearn中也是求得的数学解。不过我们也可以用梯度上升法求解。f(X)的梯度为：

也就是:

最后可化为最简形式(n * 1维)：

根据以上的推导，我们写出梯度上升的求解：

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

# 梯度
def df(w, X):
    return X.T.dot(X.dot(w)) * 2. / len(X)

def direction(w):
    return w/np.linalg.norm(w)

# 梯度上升

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

假如我们有如下数据：

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (10,8)

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

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

可以用自己实现的方法求解：

# 求解
initial_w = np.random.random(X.shape[1]) # 不能初始为0向量
eta = 0.001
print(gradient_ascent(df, X_demean, initial_w, eta))

这样可输出第一主成分。

不过我们仍可用sklearn中的PCA方法将其降维：

from sklearn.decomposition import PCA
X = np.empty((100, 2))
X[:,0] = np.random.uniform(0., 100., size = 100)
X[:,1] = 0.75 * X[:,0] + 2. * np.random.normal(0., 10., size = 100)

pca = PCA(n_components=1)
pca.fit(X)

X_reduction = pca.transform(X)
X_restore = pca.inverse_transform(X_reduction)

plt.scatter(X[:,0], X[:,1], color = 'b', alpha = 0.5)
plt.scatter(X_restore[:,0], X_restore[:,1], color = 'r', alpha = 0.5)
plt.show()