这篇博文简要总结关于数据降维的另外两种方法，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

很显然这个投影距离是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$

我们知道任意一个矩阵右乘一个向量$Av$，其实是对 $v$ 做一个线性变换 $\mathcal{R}^n\rightarrow \mathcal{R}^m$, 同时上面的那个问题相当于是找到一组正交向量 $v$,使得其在线性变换A下仍旧是正交的.而什么是线性变换呢？不就是拉伸和旋转嘛。
旋转不会改变正交性质，拉伸可能会wikipedia
如果拉伸的向量的方向正好是线性空间的基向量的方向，那拉伸就不会了。因此问题就变为能否在线性空间A中找到$n$ 个基向量了，这是显然的因为 $rank(A)\le min(m,n)$ 将前 $r$ 个向量设置为 $A$ 的基向量，其他在进行扩充就行了。
无法在线性空间A中表示的向量一定垂直与 A的基向量，因此投影一定为 $0$, 即变换后$Av$ 一定是 0，对应的 $\sigma_i$ 取0就好了，其他进行扩充。(这个挺直观的，想象在二维平面上的y轴上的点投影到x轴，一定是0)

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

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

SVD 的另外一种表示

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

SVD 还有很多有趣的性质，比如：

矩阵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

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

版权声明

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

作者: taotao

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

Dimension Reduction: PCA and SVD

PCA