求資料前n個主成分並進行高維資料對映為低維資料的操作

求資料前n個主成分

先前的將多個樣本對映到一個軸上以求使其降維的操作，其中的樣本點本身是二維的樣本點，將其對映到新的軸上以後，還不是一維的資料，對於n維資料來說，他應該有n個軸，第一個軸是方差最大的，第二個軸次之，以此類推，可以將主成分分析法看做是將資料從一個座標系轉換到另一個座標系中

那麼在求出第一主成分以後，如何求出下一個主成分呢？我們可以對資料進行改變來達到這個效果，即將資料在第一主成分上的分量給去掉

先前的Xi點乘上w以後是等於Xprojecti的模，那麼我們讓Xprojecti乘上w，最後得到的結果就是Xproject這個向量，如果我們想要將X這個樣本將Xproject上相應的分量給去除掉，我們只要用X減去Xproject就好了，其幾何意義就是將X對映到和Xproject相垂直的一個軸上，這個軸就是我們將資料的第一主成分上的分量去除以後的結果

這樣我們要求下一個主成分需要的操作也很明瞭了，即在新的資料上重新求一下第一主成分，那麼此時得到的第一主成分就是原來資料的第二主成分，以此類推，就可以得到後續的主成分

使用程式碼進行實現

（在notebook中）

以下為求第一個主成分的程式碼

簡單來說就是載入好包，建立一個虛假用例，進行均值歸零以後繪製一下圖形看看整體的分佈情況，詳細可以翻上一個關於第一主成分的求解，這裡不再贅述，詳情可見【筆記】求資料的對應主成分PCA（第一主成分）

  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] + 3. + np.random.normal(0. ,10. ,size=100)

  def demean(X):
      return X - np.mean(X,axis=0)

  X_demean = demean(X)

  plt.scatter(X[:,0],X[:,1])
 

影象如下：

與求第一主成分的程式碼不同的是下面的程式，有一些改動，公式方面去除掉數學公司，只使用df，相應的，去除掉first_component原來要選擇方法的輸入

  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)
      cur_iter = 0

      while cur_iter < n_iters:
          gradient = df(w, X)
          last_w = w
          w = w + eta * gradient
          w = direction(w)
          if (abs(f(w,X) - f(last_w,X)) < epsilon):
              break

          cur_iter += 1

      return w
 

在求解過程中，先初始化一個w的值，注意不能為0，然後設定eta為0.01，使用first_component來求出第一主成分

  initial_w = np.random.random(X.shape[1])
  eta=0.01
  w = first_component(X,initial_w,eta)
  w

結果如下

求解第二主成分

我們將原始資料X中的相應的對於第一主成分的分量給進行去除，設去掉以後的結果為X2，初始化的時候先設定成與X行列相同的空矩陣，然後我們對X中的每一個樣本都要進行是Xi減去Xi點乘w再乘上w這個方向（Xproject），這樣我們就得到了X2i

  X2 = np.empty(X.shape)
  for i in range(len(X)):
       X2[i] = X[i] - X[i].dot(w) * w

當然這個過程是可以向量化的，X2實際上可以直接使用整體的X進行計算，即對X減去X這個矩陣點乘上w這個向量進行變換，使其變成一個m行1列的矩陣以後再乘以w

  X2 = X - X.dot(w).reshape(-1,1) * w

然後我們看一下X2這個資料點是什麼樣子的

  plt.scatter(X2[:,0],X2[:,1])

影象如下

從影象我們可以發現由於我們只有兩個維度，去除掉第一個維度的主成分以後，剩下的第二維度就是所有的內容了。相應的資料就全部分佈在這一個直線上

那麼我們看一下這個對應的軸，同樣適用first_component這個函式

  w2 = first_component(X2,initial_w,eta)
  W2

結果如下

驗證一下w和w2之間是不是垂直關係

  w.dot(w2)

結果如下

結果趨近於0，說明是垂直關係

將上方的內容合併為一個函式first_n_component，給定一個n，給定一個x，我們就求出來X的前n個主成分分別是誰，因為X一直在變換，所以我們先設定一個副本X_pca進行相應的操作，對其進行一個demean操作，我們迴圈n次進行n次操作，將這個n個主成分送到res這個列表中

每次求主成分，首先要初始化一個initial_w，隨機對其生成，然後呼叫first_component這個方法來求出此時的相應的第一主成分對應得軸的方向向量，得出以後將其append進res列表中，然後我們就可以讓X_pca減去剛剛得到的主成分方向上的所有的分量，然後基於新的X_pca再去求新的主成分，最後返回res，其中是X中的前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_pca,initial_w,eta)
          res.append(w)
    
          X_pca = X_pca - X_pca.dot(w).reshape(-1,1) * w
    
      return res

這裡我們設定n為2，注意的是，這裡的X本來就兩個維度，最多隻有兩個主成分

  first_n_component(2,X)

結果如下

高維資料對映為低維資料的操作

高維資料對映為低維資料

先前的操作中的資料集仍然是n維的，並沒有降維，那麼我們如何將X的m維轉換為W的k維（前k個更加重要的方向）呢

我們將這一個樣本和這k個w分別做點乘，得到的就是這一個樣本對於這k個方向上相應的每個方向上的大小，那麼這k個元素合在一起就能表示這一個樣本對映到k個軸所代表的座標系的相應的樣本的大小，對每個wk進行計算以後得到的陣列成的向量就是樣本對映到新的座標系上得到的k維的向量

以此類推，對每個樣本都進行操作以後，我們就完成了樣本由高維（n維）到低維（k維）的對映，其實就是做了一個矩陣的乘法，即將X（mn）和wk的轉置（nk）做了個乘法（行乘列），最後得到Xk（m*k，即m個樣本，每個樣本有k個維度）這個矩陣

這樣就完成了高維資料向低維資料的對映

其實我們還可以進行恢復，給恢復成原來n維的資料，當然，由於在降維中丟失了一些資料，因此恢復以後的資料和原先的資料是不一樣的，但是這在數學上是成立的，同樣的思想，使Xk和Wk相乘即可，得到Xm，但是這個和原先的X是不一樣的

程式實現

在python chame中，建立PCA.py這個檔案，將PCA整體封裝成一個新的類，在建構函式中設定一個引數n_components，這是為了得到想要多少個主成分，即k，我們使其必須大於等於1，然後傳入資料，起一個新的變數components_來儲存新的主成分是什麼，在最後設定一個repr這個函式，可以打印出相應的結果

然後就開始具體的實現方法，設定兩個函式，一個fit，一個transform，其中fit函式，初始的時候使用者需要傳來一個X，同時因為是梯度上升法，所以還要設定一個eta，預設為0.01，以及一個最大迴圈次數，同時需要讓特徵數大於主成分數，之後的過程就是先前的操作了，設定好均值歸零的函式，求目標函式，求梯度，求w的方向以及梯度上升法的求解函式，過程就是，首先，要對X進行均值歸零操作，之後進行一個初始化self.components_的操作，然後迴圈操作，最後將self返回

transform函式的作用就是將X對映到各個主成分分量中，注意，X中的特徵數必須等於每一個單位方向向量對應的元素數量，過程就是使X點乘上Wk的轉置就可以了

我們還可以設定一個函式inverse_transform，用來將低維資料返回成高維資料，注意，此時的X的列數應該等於self.components_的行數，然後將X點乘上self.components_就可以了

具體程式碼如下：

  import numpy as np


  class PCA:

def __init__(self, n_components):
    """初始化pca"""
    assert n_components >= 1, "n_components must be valid"
    self.n_components = n_components
    self.components_ = None

def fit(self, X, eta=0.01, n_iters=1e4):
    """獲得資料集X的前n個主成分"""
    assert self.n_components <= X.shape[1], \
        "n_components must not be greater than the feature number of X"

    def demean(X):
        return X - np.mean(X, axis=0)

    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=0.01, n_iters=1e4, epsilon=1e-8):

        w = direction(initial_w)
        cur_iter = 0

        while cur_iter < n_iters:
            gradient = df(w, X)
            last_w = w
            w = w + eta * gradient
            w = direction(w)
            if (abs(f(w, X) - f(last_w, X)) < epsilon):
                break

            cur_iter += 1

        return w

    X_pca = demean(X)
    self.components_ = np.empty(shape=(self.n_components, X.shape[1]))
    for i in range(self.n_components):
        initial_w = np.random.random(X_pca.shape[1])
        w = first_component(X_pca, initial_w, eta, 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，對映到各個主成分分量中"""
    assert X.shape[1] == self.components_.shape[1]

    return X.dot(self.components_.T)

def inverse_transform(self,X):
    """將給定的X，反向映射回原來的特徵空間"""
    assert X.shape[1] == self.components_.shape[0]

    return X.dot(self.components_)

def __repr__(self):
    return "PCA(n_components=%d)" % self.n_components

（在notebook中）

載入好相應庫以及虛構的資料集以後，呼叫封裝好的PCA，首先例項化一個pca，同時傳入2這個數，然後進行fit操作將X傳入計算

  from PCA import PCA

  pca = PCA(n_components=2)
  pca.fit(X)

結果如下

然後我們可以用pca.components_來看一下計算出的兩個座標軸是什麼
結果如下

感覺有點差距，重新例項化一下，傳入1這個數，然後再fit

  pca = PCA(n_components=1)
  pca.fit(X)

結果如下

進行降維，我們叫降維以後的矩陣為X_reduction，將其等於呼叫transform函式的結果即可

  X_reduction = pca.transform(X)
  X_reduction.shape

結果如下（其為一百個樣本，每個樣本只有一個元素的矩陣）

我們還可以使用inverse_transform這個方法將資料恢復成二維的情況

  X_restore = pca.inverse_transform(X_reduction)
  X_restore.shape

結果如下（可知其又變成了二維）

那麼我們來看一下這個恢復以後的樣子

  plt.scatter(X[:,0],X[:,1],color='b',alpha=0.5)
  plt.scatter(X_restore[:,0],X_restore[:,1],color='r',alpha=0.5)

影象如下（其實就是在高維的空間裡面表達出低維的樣本而已）

在sklearn中的PCA（並不是用的梯度上升法）

首先呼叫sklearn中的PCA這個類

  from sklearn.decomposition import PCA

其餘的操作一樣，例項化傳入想要的主成分個數以後進行fit

  pca = PCA(n_components=1)
  pca.fit(X)

結果如下

簡單的驗證一下其中的內容

  pca.components_

結果如下

這樣就可以進行降維

  X_reduction = pca.transform(X)
  X_reduction.shape

結果如下

同理，也可以對其進行恢復

  X_restore = pca.inverse_transform(X_reduction)
  X_restore.shape

結果如下

繪製出原始的X和現在的X的影象

  plt.scatter(X[:,0],X[:,1],color='b',alpha=0.5)
  plt.scatter(X_restore[:,0],X_restore[:,1],color='r',alpha=0.5)

影象如下

以上就是使用虛擬的資料集來測試時用的全部了