PCA与ICA
关于机器学习理论方面的研究，最好阅读英文原版的学术论文。PCA主要作用是数据降维，而ICA主要作用是盲信号分离。在讲述理论依据之前，先思考以下几个问题：真实的数据训练总是存在以下几个问题：
①特征冗余情况，比如建立文档-词频矩阵过程中，"learn"和"study"两个特征，从VSM（计算文档向量间的相似度，Lucene评分机制由此推导而来）角度来看，两者独立，但是从语义角度看，是冗余的……
②特征强相关性，两个特征间具有很强的相关性，需要去除其中一个……
③训练样本数远小特征数，容易造成过拟合……
还有一类情况：在鸡尾酒宴会上，有n个嘉宾，在会场的不同角落里，布置n个麦克风。源信号是每个人的声音，n个麦克风接收到的信号混合声音，如何把声音还原分离出来？
以上两类大的情况，第一大类，涉及到数据降维，第二类涉及到盲信号分离。先谈第一类，关于理论上的深入论述，不探讨(能够搞理论创新的人，都是很少的)，只说明个人理解。PCA的计算步骤：①计算各个特征维度的均值；②更新x = x - u，就是去均值化，让均值为0；③计算各个特征的variance④归一化处理：x = x / variance；⑤计算协方差矩阵(方阵)⑥计算协方差矩阵的特征值和特征向量，把特征值按照降序排列，依次对应特征vector,取前k个列特征向量；⑦把归一化处理后的矩阵数据投影到列特征向量，就是最终结果。整个过程非常简单，求协方差矩阵的特征向量就是最后理想的结果。理论依据：①最大方差理论②最小平方误差理论③坐标轴相关度理论。说白了，PCA就是找出所有两两相互正交的坐标轴(投影轴的方向向量，对应协方差矩阵的特征向量)，这些坐标轴具有这样的特性，样本点在坐标轴上的投影(由于均值已经去零化，所以投影在轴上的点均值也为零)距离最远(均值为0的情况下，就是具有最高的方差，对应协方差矩阵的特征值)。这是从几何角度来理解，从数学角度理解，就是找出协方差矩阵的特征值(对角阵)，按降序排列，取前k个对应的特征向量，剔除不重要的特征。关于PCA最经典的解读，可以观看吴恩达的公开课，或者下载讲义和学习笔记。按照第一种理论来理解：
在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好 。因此我们认为，最好的k为特征是将n为样本点转换成k为后，每一维上的样本方差都很大。比如下图有5个样本点：(已经做过预处理，均值为0,特征方差归一)
下面将样本点投影到某一维上，这里用一条过原点的直线表示(前处理的过程实质是将原点移到样本点的中心点)。
假设我们选择两条不同的直线作投影，那么左右两条中哪个好呢？根据我们之前的方差最大化理论，左边的好，因为投影后样本点之间的方差最大。这里先解释一下投影的概念：
红色表示样例
,蓝色表示在u上的投影，u是直线的斜率也是直线的方向向量，而且是单位向里。蓝色点是在u上的投影点，离原点的距离是<,u>.由于这些样本点的每一维度均值都为0，因此投影到u上的样本点(只有一个到原点的距离)的均值仍然是0。由于投影后均值为0，因此方差为：
中间的部分正是协方差矩阵(一般协方差矩阵除以m-1，这里用m)！。
今天试着用Python写一个PCA算法，然后加载到numpy模块中，以后就可以直接使用了：
# encoding: utf-8
#PCA.py
from numpy import *;
import numpy as np;
import numpy.linalg as nl;
'''用PCA对数据降维，实现步骤：
1.零均值化；
2.归一化；
3.计算协方差矩阵；
4.对上述特征值分解；
5.按照特征值降序排列选取特征向量组成矩阵；
6.把第一步得到的矩阵投影到5中返回最后的结果.
'''
def zeromean(dataMat):
meanVal = np.mean(dataMat,axis=0);
newData = dataMat - meanVal;
return meanVal,newData;
def pca(dataMat,k):
meanVal,newData = zeromean(dataMat);
varData = np.std(newData,axis=0);
normData = newData / varData;
covMat = np.cov(normData,rowvar=0);
eigVals,eigVector = nl.eig(np.mat(covMat));
eigValIndice = np.argsort(eigVals);
eigValIndice = eigValIndice[-1:-(k+1):-1];
n_eigVect = eigVector[:,eigValIndice];
lowDDataMat = normData * n_eigVect;
return lowDDataMat;
def main():
data = np.mat(np.arange(6).reshape(3,2));
returnMat = pca(data,1);
print(returnMat);
if __name__ == "__main__":
main();
把这个文件加载到Python35\Lib\site-packages\numpy目录下，然后在修改numpy的__init__.py文件，加上如下语句：
就是第172行的语句，意思是把pca函数变为numpy包下可以直接使用的函数。使用时这样用：import numpy as np;np.pca(data,k);如果用from numpy import *的话，需要在__init__.py中修改__all__变量，把PCA加入进去就行了。
重点研究ICA，他的难度大于PCA。本人之前研究了FastICA算法，花了好长时间算是弄明白了。为了不误人子弟，直接上传原始的学术论文。写这篇博客的目的，更多的是为了方便更多的人学习，关于这方面理论的研究，本人没有什么建树，因此不能像之前写的kmeans算法那篇博客那样详细。
以下为matlab实现的源代码：
function Z = ICA(X)
%----------去均值--------------
[M,T] = size(X);
average = mean(X')'
for i = 1:M
X(i,:) = X(i,:) - average(i,:)*ones(1,T);
end
%--------白化/球化---------
Cx = cov(X',1)%----协方差矩阵------
[eigvector,eigvalue] = eig(Cx);%----计算Cx的特征值和特征向量------
W = eigvalue^(-1/2)*eigvector';
Z = W*X %正交矩阵
%-------迭代--------
MaxCount = 10000;
W = rand(m);
Critical = 0.00001
m = M;%-----分量-------
for n = 1:m
WP = W(:,n)
LastWP = zeros(m,1)
count = 0;
W(:,n) = W(:,n) / norm(W(:,n))
while abs(WP-LastWP)&abs(WP+LastWP) > Critical
LastWP = WP;
count = count + 1;
for i = 1:m
WP(i) = mean(Z(i,:)*(tanh(LastWP*Z))) - mean(1-(tanh(LastWP*Z))^2)*LastWP(i)
end
WPP = zeros(m,1)
for j = 1:n-1
WPP = WPP + (WP'*W(:,j))W(:,j)
end
WP = WP - WPP;
WP = WP / (norm(WP))
if count = MaxCount
fprintf('未找到相应的信号)
return;
end
end
W(n,:) = WP;
end
Z = W*Z
下面为python代码实现：
#!/usr/bin/env python
#FastICA from ICA book, table 8.4
import math
import random
import matplotlib.pyplot as plt
from numpy import *
n_components = 2
def f1(x, period = 4):
return 0.5*(x-math.floor(x/period)*period)
def create_data():
#data number
n = 500
#data time
T = [0.1*xi for xi in range(0, n)]
#source
S = array([[sin(xi) for xi in T], [f1(xi) for xi in T]], float32)
#mix matrix
A = array([[0.8, 0.2], [-0.3, -0.7]], float32)
return T, S, dot(A, S)
def whiten(X):
#zero mean
X_mean = X.mean(axis=-1)
X -= X_mean[:, newaxis]
#whiten
A = dot(X, X.transpose())
D , E = linalg.eig(A)
D2 = linalg.inv(array([[D[0], 0.0], [0.0, D[1]]], float32))
D2[0,0] = sqrt(D2[0,0]); D2[1,1] = sqrt(D2[1,1])
V = dot(D2, E.transpose())
return dot(V, X), V
def _logcosh(x, fun_args=None, alpha = 1):
gx = tanh(alpha * x, x); g_x = gx ** 2; g_x -= 1.; g_x *= -alpha
return gx, g_x.mean(axis=-1)
def do_decorrelation(W):
#black magic
s, u = linalg.eigh(dot(W, W.T))
return dot(dot(u * (1. / sqrt(s)), u.T), W)
def do_fastica(X):
n, m = X.shape; p = float(m); g = _logcosh
#black magic
X *= sqrt(X.shape[1])
#create w
W = ones((n,n), float32)
for i in range(n):
for j in range(i):
W[i,j] = random.random()
#compute W
maxIter = 200
for ii in range(maxIter):
gwtx, g_wtx = g(dot(W, X))
W1 = do_decorrelation(dot(gwtx, X.T) / p - g_wtx[:, newaxis] * W)
lim = max( abs(abs(diag(dot(W1, W.T))) - 1) )
W = W1
if lim < 0.0001:
break
return W
def show_data(T, S):
plt.plot(T, [S[0,i] for i in range(S.shape[1])], marker="*")
plt.plot(T, [S[1,i] for i in range(S.shape[1])], marker="o")
plt.show()
def main():
T, S, D = create_data()
Dwhiten, K = whiten(D)
W = do_fastica(Dwhiten)
#Sr: reconstructed source
Sr = dot(dot(W, K), D)
show_data(T, D)
show_data(T, S)
show_data(T, Sr)
if __name__ == "__main__":
main()
