第十五讲：PCA的奇异值分解、独立成分分析¶

回顾一下上一讲的内容——PCA算法，主要有三个步骤：

1. 将数据正规化为零期望以及单位化方差；
2. 计算协方差矩阵$\displaystyle\varSigma=\frac{1}{m}x^{(i)}\left(x^{(i)}\right)^T$；
3. 找到$\varSigma$的前$k$个特征向量。

在上一讲的最后，我们还介绍了PCA在面部识别中的应用。试想一下，在面部识别中$x^{(i)}\in\mathbb R^{10000}$，那么$\varSigma\in\mathbb R^{100000000}$，而要求这种亿级别的矩阵的特征向量并不容易，这个问题我们先放一放。

来看另一个例子，在前面（第五讲）我们讨论过关于垃圾邮件分类的问题，我们会构造一个词汇表向量，其每个分量所对应的单词如果在邮件中出现则置为$1$，反正则置为$0$。当我们对这种数据应用PCA时，产生的算法换了个名字，叫做潜在语义索引（LSI: latent semantic analysis）。不过在LSI中，我们通常跳过预处理阶段，因为正规化数据使它们具有相同的方差可能会大幅增加不参加单词的权重。比如说我们拿到了一个文档的内容$x^{(i)}$，现在希望知道这份文档与我们已有文档中的哪个相似度最高，也就是衡量两个高维输入向量所表示的文档的相似性$\mathrm{sim}\left(x^{(i)},x^{(j)}\right), x\in\mathbb R^{50000}$。我们通常的做法是衡量这两个向量直接的夹角，如果夹角越小则认为他们越相似，于是$\displaystyle\mathrm{sim}\left(x^{(i)},x^{(j)}\right)=\cos\theta=\frac{\left(x^{(i)}\right)^Tx^{(j)}}{\left\lVert x^{(i)}\right\rVert\left\lVert x^{(j)}\right\rVert}$。再来看上面这个分子，$x^{(i)}\left(x^{(j)}\right)^T=\sum_kx_k^{(i)}x_k^{(j)}=\sum_k1\{\textrm{文档i和j都包含词语k}\}$，如果文档中没有词语重复，则该式结果为$0$。但是，如果文档$j$包含单词“study”而文档$j$包含单词“learn”，那么如果一篇介绍“study strategy”的文章和一篇介绍“method of learning”的文章在这种算法下就是无关的，我们现在想要使它们相关。于是一开始，“learn”向量与“study”向量是相互正交的，它们的内积为零，我们在这两个向量之间再找一个向量$u$，然后将“learn”与“study”投影在$u$上，此时两个向量的投影点在$u$上将相距很近，那么它们做内积时将会得到一个正数，表示它们相关。于是，如果算法再遇到一篇关于“政治”的文章和另一篇包含很多政治家名字的文章时，它会判断这两篇文章是相关的。

奇异值分解（SVD: Singular Value Decomposition）¶

我们引入奇异值分解（可以参考线性代数笔记中的奇异值分解）来解决一开始遇到的大矩阵求特征向量的问题。比如在潜在语义索引中的输入，是一个$50000$维的向量，那么其对应的$\varSigma\in\mathbb R^{50000\times50000}$，这种规模的矩阵太大了，我们需要使用另一种方法实现PCA。

对矩阵$A\in\mathbb R^{m\times n}$，总有$A=UDV^T,U\in\mathbb R^{m\times m},D\in\mathbb R^{m\times n},V^T\in\mathbb R^{n\times n}$，其中$D=\begin{bmatrix}\sigma_1&&&\\&\sigma_2&&\\&&\ddots&\\&&&\sigma_n\end{bmatrix}$是一个对角矩阵，而$\sigma_i$称为矩阵的奇异值。分解之后的矩阵为$\begin{bmatrix}A\\m\times n\end{bmatrix}=\begin{bmatrix}U\\m\times m\end{bmatrix}\begin{bmatrix}D\\m\times n\end{bmatrix}\begin{bmatrix}V^T\\n\times n\end{bmatrix}$。

现在来观察协方差矩阵的定义式$\displaystyle\varSigma=\sum_{i=1}^mx^{(i)}\left(x^{(i)}\right)^T$，在前面的章节中（第二讲）我们介绍过一种叫做“设计矩阵”的构造方式，也就是将每一个样本向量作为矩阵$X$的一行拼凑一个矩阵：$X=\begin{bmatrix}—\left(x^{(1)}\right)^T—\\—\left(x^{(2)}\right)^T—\\\vdots\\—\left(x^{(m)}\right)^T—\end{bmatrix}$，则我们可以将协方差矩阵用设计矩阵来表示：$\varSigma=\begin{bmatrix}\mid&\mid&&\mid\\x^{(1)}&x^{(2)}&\cdots&x^{(m)}\\\mid&\mid&&\mid\end{bmatrix}\begin{bmatrix}—\left(x^{(1)}\right)^T—\\—\left(x^{(2)}\right)^T—\\\vdots\\—\left(x^{(m)}\right)^T—\end{bmatrix}$。

最后就是计算$\varSigma$的前$k$个特征向量了，我们选择对$X$做奇异值分解$X=UDV^T$，而$X^TX=VDU^TUD^TV^T=VD^2V^T=\varSigma$，于是在$D$中从大到小排列奇异值并在$V$中取前$k$个奇异值对应的特征向量即可。

容易看出$X\in\mathbb R^{m\times50000}$，做这种规模矩阵的奇异值分解会比直接计算$\varSigma\in\mathbb R^{50000\times50000}$的特征向量快很多。这就是使用SVD实现PCA算法的计算过程。

（不过，值得注意的是，在不同的计算软件，甚至是在同一种软件的不同版本中，对SVD的计算可能遵循某种默认的维数约定，因为SVD经常会得到带有很多零元素的$U$和$D$，而软件可能会按照某种约定舍弃这些零元素。所以，在使用SVD时，需要注意这样的维数约定。）

无监督学习各算法的对比¶

我们在前面介绍因子分析模型时指出，它是对每个因子$z^{(i)}$进行高斯建模，是一种对概率密度进行估计的算法，它试图对训练样本$X$的概率密度进行建模；而在其后介绍的PCA则有所不同，它并不是一个概率算法，因为它并没有使用任何概率分布拟合训练集，而是直接去寻找子空间。从这里我们可以大致的看到，如何在解决问题时在因子分析与PCA间做取舍：如果目标就是降低数据维数、寻找数据所在的子空间，我们就更倾向于使用PCA；而因子分析会假设数据本来就在某子空间内，如果有维度很高的数据$X$，而我又想对$X$建模，那么就应该使用因子分析算法（比如做异常检测，我们可以建立关于$P(X)$的模型，如果有一个低概率事件，就可以将这个事件分解在各因子分布中，进而估计其异常情况）。这两种算法的共同点就是，它们都会假设数据位于或靠近某个低维的子空间。

再来看回顾一开始介绍的两种无监督学习方法：混合高斯模型以及$k$-means算法。这两种算法的共同点是——它们都会假设数据聚集在位于某几个簇内。不同点是混合高斯模型是一种对概率密度进行估计的算法，而$k$-means则不是。所以，如果我们需要将数据分成簇并对每一个簇建模，那么我们就倾向于使用高斯混合模型；而如我我们只想将数据分簇，并不要求确定每个簇的概率系统，则就更倾向于使用$k$-means算法。

综合上面的观点可以得到表格，便于记忆：

$$ \begin{array} {c|c|c} &\textbf{Model }P(x)&\textbf{Not probabilistic}\\\hline \textbf{Subspace}&\textrm{Factor Analysis}&\textrm{PCA}\\\hline \textbf{Cluster}&\textrm{Mixtures of Gaussians}&k\textrm{-means} \end{array} $$

第十二部分：独立成分分析（Independent components analysis）¶

接下来，我们将介绍独立成分分析（ICA: Independent components analysis）。类似于PCA，独立成分分析也会寻找一组新的基，用来重新表示训练样本，然而这两个算法的目标截然不同。

举一个实际问题作为例子，在第一讲，我们介绍了一个在鸡尾酒会上从嘈杂的背景音中分离发言者音频的应用。假设有$n$个发言者在一个舞会上同时说话，而在屋子里的各话筒仅捕捉到这种$n$个发言者的声音叠加在一起的音频。但是，假设我们有$n$个不同的话筒，由于每个话筒到个发言者的距离都不相同，则这些话筒记录下来的是不同形式的发言者声音叠加。那么，通过使用这些话筒的音频记录，我们能否分离出这$n$个发言者各自的音频信号？

为了正式的描述这个问题，我们假设数据$x\in\mathbb R^n$，这些数据由$n$个相互独立的来源生成，而我们能够观测到的数据为：

$$ x=As $$

这里的$A$是一个未知的方阵，通常被称为混合矩阵（mixing matrix）。通过重复观测，我们得到一个数据集$\left\{x^{(i)};i=1,\cdots,m\right\}$，而我们的目标是还原那一组生成“观测到的数据集（$x^{(i)}=As^{(i)}$）”的声音源$s^{(i)}$。

在鸡尾酒舞会问题中，$s^{(i)}$是一个$n$维向量，$s_j^{(i)}$表示发言者$j$在第$i$次音频采集时发出的声音。而$x^{(i)}$也是$n$维向量，$x_j^{(i)}$表示话筒$j$在第$i$次音频采样时记录下来的音频。

令$W=A^{-1}$作为分离矩阵（unmixing matrix）。我们的目标就是求出$W$，进而使用$s^{(i)}=Wx^{(i)}$从话筒收集的音频中还原出各独立的音源。按照一贯的标记法，我们通常使用$w_i^T$表示$W$矩阵的第$i$行，则有$W=\begin{bmatrix}—w_1^T—\\\vdots\\—w_n^T—\end{bmatrix}$。那么，对$w_i\in\mathbb R^n$，有第$j$个发言者音源可以使用$s_j^{(i)}=w_j^Tx^{(i)}$表示。

1. ICA二义性¶

$W=A^{-1}$能够做出什么程度的还原？如果没有关于音源和混合矩阵的先验经验，不难看出，矩阵$A$存在固有二义性使得它不可能只通过$x^{(i)}$就还原并对应出每一个$s^{(i)}$。

• 令$P$为$n$阶置换矩阵，也就是$P$的每行每列都只有一个$1$，其余元素均为零。举个例子：$P=\begin{bmatrix}0&1&0\\1&0&0\\0&0&1\end{bmatrix},\ P=\begin{bmatrix}0&1\\1&0\end{bmatrix},\ P=\begin{bmatrix}1&0\\0&1\end{bmatrix}$。其作用是，对向量$z$，$Pz$会产生一个将$z$的各分量重新排列得到的的置换后的向量$z'$。对于给定的$x^{(i)}$，我们无法分辨出$W$和$PW$（也就是无法确定分离矩阵中的每一行向量对应哪一位发言者）。不难预料，音源也存在这种置换二义性，不过这种二义性对于大多数应用而言并不是重要问题。

• 此外，我们也无法确定$w_i$的大小。比如当$A$变为$2A$，而每个$s^{(i)}$变为$0.5s^{(i)}$时，对于我们的观测值$x^{(i)}=2A\cdot(0.5)s^{(i)}$没有任何影响。更广泛的，如果$A$的某列被加上了缩放因子$\alpha$，而相应的音源又被缩放因子$\frac{1}{\alpha}$调整了大小，那么我们将无法从$x^{(i)}$这一单一条件中还原这次缩放调整。因此，我们无法还原音源的原有大小。不过，对于我们关心的问题（包括鸡尾酒舞会问题），这个二义性也并不重要。特别是在本例中，使用正缩放因子$\alpha$调整发言者音源$s_j^{(i)}$的大小只会影响到发言者音量的大小。而且即使改变音源的符号也没关系，$s_j^{(i)}$和$-s_j^{(i)}$在扩音器上听起来是完全一样的。综上，如果算法得出的$w_i$被某非零缩放因子影响，那么使用$s_i=w_i^Tx$得到的相应音源也会受到这个缩放因子的影响，然而这种二义性也并不重要。（这个ICA二义性同样适用于后面介绍的脑MEG中。）

我们不禁要问，上面提到的这两种情况是ICA涉及的所有可能的二义性吗？只要源$s_i$不服从高斯分布，那么答案就是肯定的。

• 我们通过一个例子来看高斯分布下的数据会产生什么麻烦，考虑当$n=2$时的$s\sim\mathcal N(0,I)$，此处的$I$是一个二阶单位矩阵。注意到标准正态分布$\mathcal N(0,I)$概率密度的等高线图是一组圆心在原点的正圆，其概率密度具有旋转对称性。假设我们观测到某些$x=As$，$A$是混合矩阵，因为源源服从$\mathcal N(0,I)$，则混合后的$x$同样服从一个以$0$为期望、以$\mathrm E\left[xx^T\right]=\mathrm E\left[Ass^TA^T\right]=AA^T$为协方差的高斯分布。现在令$R$为任意正交矩阵（非正式的也叫作旋转矩阵或反射矩阵），则有$RR^T=R^TR=I$，并令$A'=AR$。那么如果源被$A'$混合（而不是被$A$混合），则我们会观测到$x'=A's$。而$x'$也服从高斯分布，同样以$0$围棋网、以$\mathrm E\left[x'\left(x'\right)^T\right]=\mathrm E\left[A'ss^T\left(A'\right)^T\right]=\mathrm E\left[ARss^T\left(AR\right)^T\right]=ARR^TA^T=AA^T$。可以看到，不论混合矩阵是$A$还是$A'$，我们观测到的值都会服从同一个高斯分布$\mathcal N\left(0,AA^T\right)$，于是，我们无法辨别观测到的随机变量是来自混合矩阵$A$还是$A'$。因此，混合矩阵中可以包含任意的旋转矩阵，而我们是无法从观测数据中看到这个旋转矩阵的痕迹的，所以我们也就无法完全还原出源了。

上面的这段论述是基于“多元标准正态分布是旋转对称的”这一性质的。尽管ICA对于服从高斯分布的数据存在这样的缺陷，但只要数据不服从高斯分布，在数据充足的情况下，我们还是可以分离出$n$个独立的源的。

2. 概率密度与线性变换¶

在展开ICA算法的推导之前，我们先来简要介绍一下线性变换对概率密度的影响。

假设有一个从$p_s(s)$中抽取的随机变量$s$，为了简约表达，令$s\in\mathbb R$为实数。现在按照$x=As$定义随机变量$x$（$x\in\mathbb R,A\in\mathbb R$），令$p_x$为$x$的概率密度。那么，什么是$p_x$？

令$W=A^{-1}$，要计算一个特定$x$的概率，也就是尝试计算$s=Wx$，之后再使用$p_s$估计再点$s$处的概率，从而得出$p_x(x)=p_s(Wx)$的结论——然而这是一个错误的结论。比如令$s\sim\mathrm{Uniform}[0,1]$服从均匀分布，则$s$的概率密度为$p_s(s)=1\{0\leq s\leq1\}$；再令$A=2$，则$x=2s$。显然$x$是一个在区间$[0,2]$上的均匀分布，则其概率密度函数为$p_x(x)=(0.5)1\{0\leq x\leq2\}$。而此时的$W=A^{-1}=0.5$，很明显，$p_s(Wx)$并不等于$p_x(x)$。此时，真正成立的是式子$p_x(x)=p_s(Wx)\lvert W\rvert$。

一般的，若$s$是一个向量值，来自以$p_s$为概率密度的某分布，且对于可逆矩阵$A$有$x=As$，$A$的逆矩阵为$W=A^{-1}$，那么$x$的概率密度为：

$$ p_x(x)=p_s(Wx)\cdot\lvert W\rvert $$

注意：如果曾经见过$A$将$[0,1]^n$映射到一组体积$\lvert A\rvert$的集合上（可以参考克拉默法则、逆矩阵、体积），则就有另一种记忆上面这个关于$p_x$的公式的方法，这个方法同样可以一般化前面$1$维的例子。令$A\in\mathbb R^{n\times n}$，令$W=A^{-1}$，再令$C_1=[0,1]^n$为$n$维超立方体，并定义$C_2=\{As:\ s\in C_1\}\subseteq\mathbb R^n$（即$C_2$是原像$C_1$在映射$A$作用下得到的像）。按照上面这些条件，在线性代数中有一个现成的标准公式可以使用（同时这也是行列式的一种定义方式）：$C_2$的体积为$\lvert A\rvert$。假设$s$服从在$[0,1]^n$上的均匀分布，则其概率密度函数为$p_s(s)=1\{s\in C_1\}$。很明显$x$为在$C_2$上的均匀分布，其概率密度为$\displaystyle p_x(x)=\frac{1\{x\in C_2\}}{\mathrm{Vol}(C_2)}$（因为它必须做从$C_2$到$1$的积分）。再利用“逆矩阵的行列式是原矩阵行列式的倒数”这一性质，则有$\displaystyle\frac{1}{\mathrm{Vol}(C_2)}=\frac{1}{\lvert A\rvert}=\left\lvert A^{-1}\right\rvert=\lvert W\rvert$。最终我们又得到$p_x(x)=1\{x\in C_2\}\lvert W\rvert=1\{Wx\in C_1\}\lvert W\rvert=p_s(Wx)\lvert W\rvert$。

3. ICA算法¶

现在我们可以开始推导ICA算法了。这个算法来自Bell和Sejnowski，我们这里给出的关于ICA的理解，会将其看做是一个用来最大化似然估计的算法（这与该算法原本的演绎不同，该演绎涉及到一个称为infomax原则的复杂的概念，不过在现代关于ICA的推导中已经没有必要提及了）。

我们假设每个源$s_i$来自一个由概率密度函数$p_s$定义的分布，而且$s$的联合分布为：

$$ p(s)=\prod_{i=1}^np_s(s_i) $$

注意到我们使用各源的边缘分布的乘积来对联合分布建模，也就是默认假设每一个源是独立的。使用上一节得到的公式能够得到关于$x=As=W^{-1}s$的概率密度：

$$ p(x)=\prod_{i=1}^np_s\left(w_i^Tx\right)\cdot\lvert W\rvert $$

接下来就剩为每一个独立的源$p_s$估计概率密度了。

以前在概率论中我们学过，对于给定的实随机变量$z$，其累积分布函数（CDF: Cumulative Distribution Function）$F$可以使用概率密度函数（PDF: Probability Density Function）来计算：$F(z_0)=P(z\leq z_0)=\int_{-\infty}^{z_0}p_z(z)\mathrm dz$。当然，也可以通过对CDF求导得到$z$的概率密度函数$p_z(z)=F'(z)$。

因此，要求出每个$s_i$的PDF，只需要求出它们对应的CDF的即可，而CDF是一个函数值从$0$到$1$的单调增函数。根据上一节的推导，我们知道ICA对服从高斯分布的数据无效，所以不能选择高斯分布的累积分布函数作为CDF。我们现在要选择一个合理的“默认”函数，其函数值也是从$0$缓慢的单调递增至$1$——这就是前面（第三讲）介绍的逻辑函数（即S型函数）$\displaystyle g(s)=\frac{1}{1+e^{-s}}$。于是有$p_s(s)=g'(s)$。

（如果我们有关于源数据集的先验知识，已经知道源的PDF的形式，那么就可以用该PDF对应的CDF来代替上面的逻辑函数。如果不知道PDF的形式，那么逻辑函数就是一个很合理的默认函数，因为在处理很多问题时，逻辑函数都有具有良好的表现。并且，在本例中我们使用的输入观测数据集$x^{(i)}$要么已经被预处理为期望为$0$的数据，要么$x^{(i)}$在自然状态下就是期望为$0$的数据集，如音频信号。而零期望是必须的，因为我们假设了$p_s(s)=g'(s)$，表明$\mathrm E[s]=0$——这里说明一下，对逻辑函数求导会得到一个对称函数，这也就是PDF，所以这个对称的PDF对应的随机变量必须保证期望为$0$——则有$\mathrm E[x]=\mathrm E[As]=0$。顺便再提一点，除了逻辑函数，PDF也经常选用拉普拉斯分布/Laplace distribution $\displaystyle p(s)=\frac{1}{2}e^{-\lvert s\rvert}$。）

方阵$W$是模型中的参数，对已给定训练集$\left\{x^{(i)};i=1,\cdots,m\right\}$，有对数似然函数为：

$$ \mathscr l(W)=\sum_{i=1}^m\left(\sum_{j=1}^n\log g'\left(w_j^Tx^{(i)}\right)+\log\lvert W\rvert\right) $$

我们的目标是用$W$最大化上面的函数，使用性质$\nabla_W\lvert W\rvert=\lvert W\rvert\left(W^{-1}\right)^T$（参见第二讲）对似然函数求导，便可以得到一个以$\alpha$为学习速率的随机梯度上升的更新规则。对于训练集$x^{(i)}$的更新规则是：

$$ W:=W+\alpha\left(\begin{bmatrix} 1-2g\left(w_1^Tx^{(i)}\right)\\ 1-2g\left(w_2^Tx^{(i)}\right)\\ \vdots\\ 1-2g\left(w_n^Tx^{(i)}\right) \end{bmatrix}\left(x^{(i)}\right)^T+\left(W^T\right)^{-1}\right) $$

在算法迭代收敛之后，通过$s^{(i)}=Wx^{(i)}$还原源即可。

注意：在构造似然函数时，我们默认假设了各样本$x^{(i)}$是相互独立的（要注意，这不是指样本$x^{(i)}$的各分量间是相互独立的）。于是能够得到训练集的似然函数为$\prod_ip\left(x^{(i)};W\right)$。这个假设在$x^{(i)}$表示演讲音频或其他基于时间序列的数据时是错误的，因为这种样本数据之间并不是相互独立的，不过这一点也可以证明在训练样本充足的前提下，存在相关关系的各样本并不会影响ICA的表现。不过，对于相关的训练样本，当我们使用随机梯度上升时，如果随机改变训练样本载入算法的次序，有时可能会加速算法的收敛过程。（建议准备多组经过乱序的训练集，然后用不同的顺序为模型加载样本，有时有利于快速收敛。）