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