这篇博文简要总结关于数据降维的另外两种方法,PCA与SVD。之前讲过用JL变换和随机投影的方式做降维,那种方式的出发点是从计算效率出发的,所以计算方式会更快同时也是数据透明(data-oblivious)的,也就是不关心数据的特点,而今天要讲的两种方式都关心数据的特点。同时这里给出了一些PCA的相关直观性解释

PCA

首先我们说降维的目标始终都是用一个低维空间($\mathcal{R}^m$) 来近似的表示高维空间 $\mathcal{R}^n,m\ll n$ 中的向量。

也就是说 PCA的目标就是

将 $m$ 个 $n$ 维空间中的向量 $\mathbf{x}_1,\dots,\mathbf{x}_m$ 用 $k$ 个$n$ 维空间中的向量来近似表示:

$$ \mathbf{x}_i\approx \sum_{j=1}^k a_{ij}\mathbf{v}_j,\forall l,m\le k,l\neq m,\mathbf{v}_l^T\mathbf{v}_m=0,||\mathbf{v}_j||=1 $$ 也就是说将 每个n维向量往 $k$ 个正交向量张成的空间投影。

因此PCA的目的就是找到这 $k$ 正交向量。怎么找我们留在后面(注意在随即投影中,这几个向量都是random选择的)

与linear regression 的关系

也就是说最大的不同在于 线性回归找到的线是使得说有点到线的竖直距离最小,而pca是让找到线到说有点的垂直距离(perpendicular) 距离最小。

问题的严格定义

1维情况

我们的目标就是要:

$$ \mathbf{v} = \arg \min_{||v||=1} \sum_i dist(\mathbf{x}_i\leftrightarrow v)^2 $$

很显然这个目标等价于

$$ \mathbf{v} = \arg \max_{||v||=1} \sum_i <\mathbf{x}_i,v>^2\tag{1} $$

detail

  1. 很显然这个投影距离是shift,和scale敏感的,所以做pca前通常需要0均值化,以及去scale

多维

定义 $n$ 维空间中 $k$ 个正交向量,所张成的子空间 $S=span\{\mathbf{v}_1,\dots,\mathbf{v}_k\}$,很显然,多为情况就是

$$ S = \arg \max \sum_i \sum_j (<\mathbf{x}_i,\mathbf{v}_j>)^2 $$

compute PCA

k=1

首先考虑 $k=1$ 的情形,

$$ \begin{aligned} &\sum <\mathbf{x}_i,\mathbf{v}> ^2\\
=&(Xv)^T (Xv)\\
=&v^TX^TXv
\end{aligned} $$

其中 $X$ 是行向量矩阵,记 $A = X^TX$,(NOTE 这个矩阵一定是semi-positive的矩阵)

因此,原来的问题可以简化为

$$ v=\arg \max v^TAv=v^TQ^TDQv\tag{2} $$

由于 $A$ 是半正定的,我们一定可以用特征值分解成 $A=Q^TDQ$, 其中 $D=diag\{\lambda_1,\dots,\lambda_m\}$

矩阵变换

注意对于正交矩阵 $Q$, 任意一个向量 $\mathbf{v}$ ,$Qv$ 其实是做了一个旋转变换,是保长度和角度的 $||Qv||^2 = vQ^TQv=v^Tv=||v||$

同时对角矩阵 $D$, $Dv$ 是拉伸变换。

也就是说对于$v=\{v_1,\dots,v_d\}$

$Dv=\{\lambda_1 v_1,\dots,\lambda_d v_d\}$

同时

$$ v^TDv=\sum \lambda_i v_i^2\tag{3} $$

由于 $||v||=\sum v_i^2 =1$, 所以最大化(3) 式,其实是在仅需令 $v_1=1$ 就行了(特征值由大到小排序)

而对于 (2) 式来说,仅需让 $v=Q_1$ 就行了,也就是第一个最大特征值

更大的 $k$

对于更大的 $k$ 来说由于特征向量之间是正交的,所以仅需让 $v_i=Q_i$ 就行了,也就是一般教科书中的定义

第 $k$ 个主成分为$X^TX$ 的第 $k$ 个特征向量

power iteration

特征值的计算,这里有一个迭代式的最大特征值计算方法

import numpy as np


def power_iteration(A, num_simulations):
    """
    compute largest eigenvector and eigenvalue,
    code copy from wikipedia:https://en.wikipedia.org/wiki/Power_iteration
    Args:
        A:2-d array
            must be a semi-positive matrix
        num_simulationa: int
            the number of iterations
    return: (eigenvetor,eigenvalue)
    """
    # Ideally choose a random vector
    # To decrease the chance that our vector
    # Is orthogonal to the eigenvector
    b_k = np.random.rand(A.shape[1])

    for _ in range(num_simulations):
        # calculate the matrix-by-vector product Ab
        b_k1 = np.dot(A, b_k)

        # calculate the norm
        b_k1_norm = np.linalg.norm(b_k1)

        # re normalize the vector
        b_k = b_k1 / b_k1_norm

    return b_k,b_k.dot(A).dot(b_k)

这个的直观是有的,我们知道 $Av$ 相当于不停的往矩阵 $A$ 的特征子空间做投影,因此当向量变化到最大特征值对应的特征值方向的时候,这个投影就不会变化了,近似保护的了投影最大。然而这只是我们理解,这个算法背后有非常严格的数学证明,这里就掠过了

SVD

SVD 属于另外一种数据降维的方法,他提供了矩阵的一种低秩(low rank)近似。

SVD 简介

这部分属于线性代数中的类容,这里仅对SVD做个简单的直观上的证明。非严格的数学证明

SVD是说 任意的一个矩阵 $A$,均可以写成

$$ A=U_{m\mathbb{x}n}\Sigma_{m\mathbb{x}n} V_{n\mathbf{x}n}^T,U,V,s.t.||U||=1,U^TU=I $$

即 $U,V$ 都是正交阵,$\Sigma=diag\{\sigma_1,\dots,\sigma_r,\},r=rank(A)$ 是对角阵

注意这个形式和特征值分解特别像,不过特征值分解是对方阵而言的,而SVD是对任意矩阵而言的,它的证明有很多种方法,我这里提供一个简单的非常直观的不那么严格的证明,主要是提供一种直观上的理解

这个问题等价于,我们要在 $m$ 维空间中找到一组正交基 $\{v_1,\dots,v_m\}.s.t. Av_i=\sigma_i u_i$

  1. 我们知道任意一个矩阵右乘一个向量$Av$,其实是对 $v$ 做一个线性 变换 $\mathcal{R}^n\rightarrow \mathcal{R}^m$, 同时上面的那个问题相当于是找到一组正交向量 $v$,使得其在线性变换A下仍旧是正交的.而什么是线性变换呢?不就是拉伸和旋转嘛。
  2. 旋转不会改变正交性质,拉伸可能

    wikipedia

  3. 如果拉伸的向量的方向正好是线性空间的基向量的方向,那拉伸就不会了。因此问题就变为能否在线性空间A中找到$n$ 个基向量了,这是显然的因为 $rank(A)\le min(m,n)$ 将前 $r$ 个向量设置为 $A$ 的基向量,其他在进行扩充就行了。

  4. 无法在线性空间A中表示的向量一定垂直与 A的基向量,因此投影一定为 $0$, 即变换后$Av$ 一定是 0,对应的 $\sigma_i$ 取0就好了,其他进行扩充。(这个挺直观的,想象在二维平面上的y轴上的点投影到x轴,一定是0)

NOTE 以上证明并不严格,仅仅是有个直观理解,

同时SVD也可以理解为将线性变换分解为旋转拉伸再旋转的过程

SVD 的另外一种表示

$$ A=\sum \sigma_i u_i v_i^T $$

SVD 还有很多有趣的性质,比如:

  1. 矩阵A的 F-norm: $||A||_F^2=tr(A^TA)=||V\Sigma U^TU\Sigma V^T||_F^2=\sum \sigma_i^2$

low rank 近似

$$ A_k=\sum_{i=1}^k \sigma_i u_i v_i^T,k\le rank(A) $$

可以证明,对于任意的 $B,rank(B) =k$

$$ ||A-A_k||\le ||A-B|| $$

对于 F-norm 成立(这个证明很麻烦,我暂时没找到一个简单的证明,得复习一下代数了)

而这可以用来压缩矩阵。对于原来的矩阵你需要 $space=O(mn)$,现在就仅仅需要 $O(k(m+n))$ 了。

reference

  1. Gregory Valiant.2019 cs168 lecture note 7
  2. Gregory Valiant.2019 cs168 lecture note 8
  3. Gregory Valiant.2019 cs168 lecture note 9
  4. Anand Rajaraman, Jure Leskovec, and Jeffrey D. Ullman. mining massive data sets

版权声明

本作品为作者原创文章,采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议

作者: taotao

仅允许非商业转载,转载请保留此版权声明,并注明出处