主要思想：移动坐标轴，将n维特征映射到k维上（k<n），这k维是全新的正交特征。
这k维特征称为主元，是又一次构造出来的k维特征。而不是简单地从n维特征中去除其余n-k维特征。
假设有如下数据：
上面的数据这么解读，表示有两个点：
这两个点在初始坐标系下（也就是自然基e1 = [1,0] e2 = [0,1] ）下坐标值为：
随着坐标系的不同，X1,X2 的值会不断变化，都是在坐标轴上的投影长度，要想值尽可能得分配给X1X2，借鉴最小二乘法的思想，就是让
假设：
根据点积的几何意义，可以表示：
那么：
上式其实是一个二次型：
其中U为正交矩阵，即UUT = I，
对角矩阵，奇异值1>奇异值2
把P带回去：
因为U是正交矩阵，所以令
n也是单位矩阵
继续会带：
最初的问题就转化为：
用拉格朗日乘子法计算上述条件极值，结果是当n1 = 1，n2 = 0时取到极值。
因此可以推出要寻找的主元1，即：
总结下：
同理，
协方差矩阵
按行来解读，得到
按列解读，得到
矩阵就变成了
样本方差为：
样本协方差为：
两相比较可以得到一个新的矩阵，也就是协方差矩阵：
协方差矩阵
P、Q都可以进行奇异值分解：
可见，协方差矩阵Q的奇异值分解和P相差无几，只是奇异值缩小了n倍，但是不妨碍奇异值之间的大小关系，所以在实际问题中，往往都是直接分解协方差矩阵。
python实现
基本步骤：
- 对数据进行归一化处理（代码中并不是这么做的，而是直接减去均值）
- 计算归一化后的数据集的协方差矩阵
- 计算协方差矩阵的特征值和特征向量
- 保留最重要的k个特征（通常k要小于n）。也能够自己制定。也能够选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k
- 找出k个特征值相应的特征向量
- 将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据。
1.引入包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
2.定义一个均值函数
# 计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征
def meanX(dataX):
return np.mean(dataX,axis=0)
#axis=0表示依照列来求均值。假设输入list,则axis=1
2.编写pca方法
"""
參数：
- XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征
- k：表示取前k个特征值相应的特征向量
返回值：
- finalData：參数一指的是返回的低维矩阵，相应于输入參数二
- reconData：參数二相应的是移动坐标轴后的矩阵
"""
def pca(XMat, k):
average = meanX(XMat)
m, n = np.shape(XMat)
data_adjust = []
avgs = np.tile(average, (m, 1))
data_adjust = XMat - avgs
covX = np.cov(data_adjust.T) #计算协方差矩阵
featValue, featVec= np.linalg.eig(covX) #求解协方差矩阵的特征值和特征向量
index = np.argsort(-featValue) #依照featValue进行从大到小排序
finalData = []
if k > n:
print "k must lower than feature number"
return
else:
#注意特征向量时列向量。而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值
selectVec = np.matrix(featVec.T[index[:k]]) #所以这里须要进行转置
finalData = data_adjust * selectVec.T
reconData = (finalData * selectVec) + average
return finalData, reconData
3.载入数据集的函数
# 把txt文件读取为矩阵
def read_txt_to_array(path, delimiter='\t'): # delimiter是数据分隔符
fp = open(path, 'r', encoding='utf-8')
string = fp.read() # string是一行字符串，该字符串包含文件所有内容
fp.close()
row_list = string.splitlines() # splitlines默认参数是‘\n’
data_list = [[float(i) for i in row.strip().split(delimiter)] for row in row_list]
return np.array(data_list)
4.可视化结果（k=2的情况）
def plotBestFit(data1, data2):
dataArr1 = np.array(data1)
dataArr2 = np.array(data2)
m = np.shape(dataArr1)[0]
axis_x1 = []
axis_y1 = []
axis_x2 = []
axis_y2 = []
for i in range(m):
axis_x1.append(dataArr1[i,0])
axis_y1.append(dataArr1[i,1])
axis_x2.append(dataArr2[i,0])
axis_y2.append(dataArr2[i,1])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(axis_x1, axis_y1, s=50, c='red', marker='s')
ax.scatter(axis_x2, axis_y2, s=50, c='blue')
plt.xlabel('x1'); plt.ylabel('x2');
plt.savefig("outfile.png")
plt.show()
5.测试
def main():
datafile = "data.txt"
XMat = read_txt_to_array(datafile)
k = 2
return pca(XMat, k)
if __name__ == "__main__":
finalData, reconMat = main()
plotBestFit(finalData, reconMat)
