PCA与梯度上升法
Posted 周虽旧邦其命维新
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了PCA与梯度上升法相关的知识，希望对你有一定的参考价值。
1、什么是PCA
主成分分析(Principal Component Analysis):
- 一个非监督的机器学习算法
- 主要用于数据的降维
- 通过降维，可以发现更便于人类理解的特征
- 其他应用：可视化，降噪
假设存在一根直线，将所有的点都映射在该条指直线上，这样的话点的整体分布和原来的点的分布就没有很大的差异（点和点的距离比映射到x轴或者映射到y轴都要大，区分度就更加明显），与此同时所有的点都在一个轴上（理解成一个维度），虽然这个轴是斜着的。用这种方式将二维降到了一维度
那么如何找到这个让样本间距最大的轴？
如何定义样本间间距？ 使用方差(Variance)
方差越大代表样本之间越稀疏，方差越小代表样本之间越紧密。
移动坐标轴，使得样本在每一个维度均值都为0：
demean之后的方差最大其实就是求映射后每个点到(0,0)的距离最大再求和，假设降维后轴的方向是w=(w1, w2)，Xi是映射前的向量，Xi(project)是映射后的向量
目标函数即：
对目标函数求梯度：
转化为：
转化为：
由于最终转换的结果是一个1行m列的矩阵，而我们想要得到一个n行1列的矩阵，所以还要进行一次转置
2、梯度上升法求解PCA
2.1 求数据的第一个主成分
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] + 2 + np.random.normal(0, 10., size=100)
# plt.scatter(X[:,0],X[:,1])
# plt.show()
# 均值归零函数
def demean(X):
return X - np.mean(X, axis=0)
X_demean = demean(X)
# print(X_demean)
# print(np.mean(X_demean, axis=0))
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.001):
res = np.empty(len(w))
for i in range(len(w)):
w_1 = w.copy()
w_2 = w.copy()
w_1[i] += epsilon
w_2[i] -= epsilon
# 2*epsilon必须用括号包起来，否则就不是除以2倍的epsilon而是除以2乘以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)
i_iter = 1
while i_iter < n_iters:
gradient = df(w=w, X=X)
last_w = w
# gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
# 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
# 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
w = w + eta * gradient
w = direction(w) # 注意1，每次求一个单位向量
# abs求绝对值
if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
break
i_iter += 1
return w
initial_w = np.random.random(X.shape[1]) # 注意2：不能用0向量开始
eta = 0.0001
# print(gradient_ascent(df_debug, X=X_demean, initial_w=initial_w, eta=eta))
print(gradient_ascent(df_math, X=X_demean, initial_w=initial_w, eta=eta, n_iters=100))
注意：
1、 w为0本来就是一个极值点，只不过这个极值点对应的是最小值的位置，但是我们求的是最大值的位置。在这个最小值的位置上它的梯度值也为0，所以此时把
零向量代进去是不可以的，所以开始的时候随机化一个向量就可以了
2、不能用StandardScaler标准化数据，因为PCA过程本来就是要求一个轴，使得所有的样本映射到那个轴之后样本的方差最大，一旦我们将样本的数据进行标准
化了，样本的方差就为1了，方差的最大值就不存在了，因为在标准化的过程中把样本的方差给打掉了，这样子的话就求不出来PCA真正想最大化的结果了
2.2 求数据的前n个主成分
根据推导过程可以发现，去除第一个主成分剩下的向量X’和第一主成分向量是垂直的关系，而两个垂直的向量相乘等于0，所以可以根据这个特性去验证我们求出的结果是否正确。
代码实现：
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] + 2 + np.random.normal(0, 10., size=100)
# 均值归零函数
def demean(X):
return X - np.mean(X, axis=0)
X_demean = demean(X)
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 first_component(X, initial_w, eta, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
i_iter = 1
while i_iter < n_iters:
gradient = df(w=w, X=X)
last_w = w
# gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
# 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
# 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
w = w + eta * gradient
w = direction(w) # 注意1，每次求一个单位向量
# abs求绝对值
if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
break
i_iter += 1
return w
initial_w = np.random.random(X.shape[1]) # 注意2：不能用0向量开始
eta = 0.01
w = first_component(X=X, initial_w=initial_w, eta=eta)
print(w)
# X2 = np.empty(X.shape)
# for i in range(len(X)):
# X2[i] = X[i] - X[i].dot(w) * w
X2 = X - X.dot(w).reshape(-1,1) * w
# plt.scatter(X2[:,0],X2[:,1])
# plt.show()
w2 = first_component(X2, initial_w=initial_w, eta=eta)
print(w2)
print(w.dot(w2))
# 求前n个主成分
def first_n_component(n, X, eta=0.01, n_iters=1e4, epsilon=1e-8):
X_pca = X.copy()
X_pca = demean(X_pca)
res = []
for i in range(n):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(X=X_pca, initial_w=initial_w, eta=eta)
res.append(w)
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return res
print(first_n_component(2, X))
2.3 高维数据映射为低维数据
编写PCA.py文件
import numpy as np
class PCA:
def __init__(self, n_components):
self.n_components = n_components
# 数据的前n个主成分，公式中的Wk
self.components_ = None
# 取得数据的前n个主成分
def fit(self, X, eta=0.01, n_iters=1e4):
# 均值归零函数
def demean(X):
return X - np.mean(X, axis=0)
X_demean = demean(X)
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 first_component(X, initial_w, eta, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
i_iter = 1
while i_iter < n_iters:
gradient = df(w=w, X=X)
last_w = w
# gradient是对每个维度求偏导得到的列表，如果偏导数为负则w的这个维度加上一个负值，降维后的方差趋于变大
# 如果偏导数为正，则w的这个维度加上一个正值，降维后的方差趋于变大，因此w加上导数值，降维后的方差趋于变大
# 在eta合适的情况下，随着循环进行，导数值逐渐趋近0，eta是常数，降维后的方差的变化量会越来越小
w = w + eta * gradient
w = direction(w) # 注意1，每次求一个单位向量
# abs求绝对值
if (abs(f(w=w, X=X) - f(w=last_w, X=X)) < epsilon):
break
i_iter += 1
return w
X_pca = demean(X)
# 数据的前n个主成分，公式中的Wk
self.components_ = np.empty(shape=(self.n_components, X_pca.shape[1]))
for i in range(self.n_components):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(X=X_pca, initial_w=initial_w, eta=eta, n_iters=n_iters)
self.components_[i] = w
X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w
return self
def transform(self, X):
# 将X映射到各个主成分上去，结果就是公式中的Xk
return X.dot(self.components_.T)
def inverse_transform(self,X):
# 将X反向映射回原来的特征空间，参数是公式中的Xk，结果是公式中的Xm
return X.dot(self.components_)
测试：
import numpy as np
import matplotlib.pyplot as plt
from common.PCA 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)
print(pca.components_)
pca_reduction = pca.transform(X)
# print(pca_reduction)
X_restore = pca.inverse_transform(pca_reduction)
print(X_restore)
plt.scatter(X[:,0], X[:,1], color='r')
plt.scatter(X_restore[:,0], X_restore[:,1])
plt.show()
3、scikitlearn中的PCA
scikitlearn中的PCA使用的不是梯度上升法求解出来的，而是数学方法。
3.1 scikitlearn中PCA基本使用
import numpy as np
import matplotlib.pyplot as plt
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)
X_reduction = pca.transform(X=X)
print(X_reduction.shape)
X_restore = pca.inverse_transform(X=X_reduction)
plt.scatter(X[:,0], X[:,1], color='r')
plt.scatter(X_restore[:,0], X_restore[:,1])
plt.show()
3.2 用scikitlearn中PCA分析digits数据集
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA
import time
digits = datasets.load_digits()
X = digits.data
y = digits.target
X_train,X_test,y_train,y_test = train_test_split(X, y, random_state=666)
# print(X_train.shape)
start_time = time.time()
knn_clf = KNeighborsClassifier(n_neighbors=5)
knn_clf.fit(X_train, y_train)
print(knn_clf.score(X_test,y_test))
end_time = time.time()
print(end_time - start_time)
start_time = time.time()
pca = PCA(n_components=41)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
knn_clf = KNeighborsClassifier(n_neighbors=5)
knn_clf.fit(X_train_reduction, y_train)
print(knn_clf.score(X_test_reduction, y_test))
end_time = time.time()
print(end_time - start_time)
# explained_variance_ratio_涵盖了原有数据的方差的比例
# print(pca.explained_variance_ratio_)
# 查看每个特征在原数据方差的占比
# pca = PCA(n_components=X_train.shape[1])
# pca.fit(X_train)
# print(pca.explained_variance_ratio_)
# 保留占比超过0.99的前n个主成分
pca = PCA(n_components=0.99)
pca.fit(X_train)
print(pca.n_components_)
3.3 用scikitlearn中PCA分析mnist数据集
import numpy as np
from sklearn.datasets import fetch_openml
from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA
import time
# fetch_mldata()不再能够使用是因为其所依赖的资源不再适用。将之替换为fetch_openml()
# 需要注意的是这个替换并不是一个无缝替换。例如mnist数据集需要改为mnist_784，具体数据集是在https://www.openml.org/
mnist = fetch_openml("mnist_784")
# print(mnist)
X,y = mnist['data'],mnist['target']
# print(X.shape)
X_train = np.array(X[:60000], dtype=float)
y_train = np.array(y[:60000], dtype=float)
X_test = np.array(X[60000:], dtype=float)
y_test = np.array(y[60000:], dtype=float)
print(X_train.shape, y_train.shape)
print(X_test.shape, y_test.shape)
start_time = time.time()
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
end_time = time.time()
print("模型训练时间：", end_time - start_time)
start_time = time.time()
print(knn_clf.score(X_test, y_test))
end_time = time.time()
print("模型预测时间：", end_time - start_time)
# 使用PCA降维
pca = PCA(0.9)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
start_time = time.time()
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
end_time = time.time()
print("模型训练时间：", end_time - start_time)
start_time = time.time()
print(knn_clf.score(X_test_reduction, y_test))
end_time = time.time()
print("模型预测时间：", end_time - start_time)
特征脸方法就是将PCA方法应用到人脸识别中，将人脸图像看成是原始数据集，使用PCA方法对其进行处理和降维，得到“主成分”——即特征脸，然后每个人脸都可以用特征脸的组合进行表示。这种方法的核心思路是认为同一类事物必然存在相同特性（主成分），通过将同一目标（人脸图像）的特性寻在出来，就可以用来区分不同的事物了。人脸识别嘛，就是一个分类的问题，将不同的人脸区分开来。
以上是关于PCA与梯度上升法的主要内容，如果未能解决你的问题，请参考以下文章