【笔记】求数据前n个主成分以及对高维数据映射为低维数据

求数据前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)

图像如下

以上就是使用虚拟的数据集来测试时用的全部了