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();

把这个文件加载到Python35Libsite-packages
umpy目录下，然后在修改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()