为啥标题会有TensorFlow？因为它的API设计叕把我惊喜(吓)到了。这部分放到后面来说，先上点前菜。
1. 主成分分析
首先来讲什么是主成分分析（Principal Component Analysis），关于这个概念的介绍网上有很多，推荐最大方差解释和最小平方误差解释这两篇博客。
用最概况的话来说：主成分是一种用简单矩阵变换来近似模拟复杂矩阵变换的方法。
稍啰嗦一点的说法是，向量表示的是某个高维空间的位置（或状态），矩阵的现实意义是表示线性坐标空间的变换（参见孟岩老师经典的“理解矩阵”三部曲，链接：【1】/【2】/【3】）。所以向量与矩阵相乘即为位置的移动（变换到另一个位置），矩阵与矩阵相乘相当于对原矩阵所有空间坐标的一次映射，是一种对空间坐标轴方向和各方向单位长度的变换。对于同一种变换结果，通常可以使用几种不同的变换方式获得，这些变换所代表的矩阵互为相似矩阵，它们具有相同的维度。或者，我们也可以找出一些变换，使得变换后的空间降维，却保持与原矩阵空间所能表示的位置集合尽可能接近（这类变换矩阵通常不是方阵，因为一个矩阵只能乘以长度相同，宽度更小的矩阵才会变得更小）。
我们可以将矩阵中的每一行的向量看做该矩阵的一个坐标轴方向，大部分矩阵代表的坐标轴之间并没有正交。假如每个坐标表示一个影响最终状态的因素，则一种因素的变化将引起该空间中的向量（表示一种状态）在两个或以上坐标方向上移动。因此若想用最少的坐标表达原有空间的所有状态，就是找到一个矩阵将这个空间变换为正交空间。
2. 特征值分解
特征值分解和奇异值分解是实现空间坐标正交变换的方法，关于这两种数学方法的详细解释，推荐一篇写得很好的文章：强大的矩阵奇异值分解及其应用。简要的说明一下结论，特征值分解是指任意一个可对角化的方阵M都可以被分解为以下形式：
M = Q x Σ x I(Q)
其中Σ是一个对角矩阵，对角上的每个值称为特征值。Q是一个方阵，上面每一行称为特征向量，这些向量都是正交的。I(Q)表示对矩阵Q取逆矩阵。
TensorFlow提供了求特征值分解的方法tf.self_adjoint_eig，效果上和Numpy的np.linalg.eig差不多（返回的特征值也都没排序）。Numpy的array写起来方便一点，同时为了向等下要说的另一种网上改进版算法致敬，这里用Numpy举例：
>>> import numpy as np
>>> M = np.array([[ 5, 8, 6, 7],
[ 8, 5, 7, 4],
[ 6, 7, 5, 3],
[ 7, 4, 3, 5]])
>>> s, Q = np.linalg.eig(M)
# 打印所有特征值
>>> print(np.diag(s))
[[ 22.7986798 0. 0. 0. ]
[ 0. 2.8476628 0. 0. ]
[ 0. 0. -3.86490633 0. ]
[ 0. 0. 0. -1.78143627]]
# 打印所有特征向量
>>> print(Q)
[[ 0.56359593 0.17366844 0.7480347 -0.3043731 ]
[ 0.53290936 -0.31019456 -0.55620238 -0.55715874]
[ 0.47049019 -0.53931934 0.05405539 0.69631289]
[ 0.42072107 0.76338278 -0.35799584 0.33478277]]
# 验证性质
>>> print(M)
[[5 8 6 7]
[8 5 7 4]
[6 7 5 3]
[7 4 3 5]]
>> print(Q.dot(np.diag(s).dot(Q.I)))
[[ 5. 8. 6. 7.]
[ 8. 5. 7. 4.]
[ 6. 7. 5. 3.]
[ 7. 4. 3. 5.]]
在有的地方经常会见到下面这个公式：
M x v = λ x v
这里的λ是任意的一个特征值，v是任意的一个特征向量。同样可以用Numpy验证这个特性（注意Numpy返回的特性向量是列向量）：
>>> print(np.matmul(m, v[:,0]))
[[ 12.84924318]
[ 12.14962988]
[ 10.72655529]
[ 9.59188488]]
>>> print(l[0] * v[:,0])
[[ 12.84924318]
[ 12.14962988]
[ 10.72655529]
[ 9.59188488]]
后面这个公式我们现在不关心，回到M = Q x Σ x I(Q)这里来，我们把两边同时乘以矩阵Q，得到：
M x Q = Q x Σ
这个等式的意义在于，右侧的内容变为一个矩阵Q乘以对角矩阵Σ，根据矩阵乘法，矩阵乘以对角矩阵等同于，将被乘矩阵的每一行分别乘以对角矩阵相应行上的数值。同时由于矩阵Q的每一行都是一个正交空间坐标轴的向量，所以为了最大限度的保留变换结果，我们应该将数值绝对值最小的特征值所对应的行去除。
假设原本矩阵M表示的是n个数据样本（每行是一个样本），每个样本有n个特征（每列是一个特征），降维后的每个特征维度为r（即将原本n x n的矩阵用n x r的矩阵表示），则上述过程可表示为：
import numpy as np
def pca(M, r):
s, Q = np.linalg.eig(M) # 求特征值和特征向量
idx = np.argsort(np.fabs(s)) # 对特征值绝对值进行排序
top_idx = idx[:-(r + 1):-1] # 取出特征值绝对值最大的下标
rs = s[top_idx] # 取出排序后下标对应的特征值
rQ = Q[:, top_idx] # 取出排序后下标对应的特征向量
rM = rQ * np.diag(rs) # 得到降维后的矩阵
bM = rM * rQ.I # 使用降维后的矩阵还原原始矩阵
return rM, bM # 返回降维后的矩阵和还原结果
测试一下效果：
# 这是一个能够包含四维数据特征的矩阵
>>> M = np.array([[ 5, 8, 6, 7],
[ 8, 5, 7, 4],
[ 6, 7, 5, 3],
[ 7, 4, 3, 5]])
# 降到三维，查看还原结果
>>> rM, bM = pca(m, 3)
>>> print(bM)
[[ 5.16503758 8.30210333 5.62244433 6.81847366]
[ 8.30210333 5.5530039 6.30887965 3.66771378]
[ 5.62244433 6.30887965 5.86373231 3.41527695]
[ 6.81847366 3.66771378 3.41527695 5.19966249]]
# 降到两维，查看还原结果
>>> rM, bM = pca(m, 2)
>>> print(bM)
[[ 5.07914999 8.45550979 5.88916425 6.44094335]
[ 8.45550979 5.27899989 5.83248297 4.3420323 ]
[ 5.88916425 5.83248297 5.03544588 4.58767992]
[ 6.44094335 4.3420323 4.58767992 3.5401777 ]]
# 降到一维，查看还原结果
>>> rM, bM = pca(m, 1)
>>> print(bM)
[[ 7.24178118 6.84748197 6.04544292 5.40594729]
[ 6.84748197 6.4746515 5.71628173 5.11160524]
[ 6.04544292 5.71628173 5.04673908 4.51288778]
[ 5.40594729 5.11160524 4.51288778 4.03550804]]
我们忽略降维后的矩阵rM值的输出，单看还原效果。可以看出，在下降维度不多时，从低维矩阵可以近似的还原回原本矩阵的数值。
除了照搬原始公式实现的PCA降维方法，在网上还有一个不明觉厉的版本，例如主成分分析(PCA)—基于python+numpy和Python实现PCA等博客都介绍了这种算法，实测效果更佳。
3. 奇异值分解
看完特征值分解，再来看它的扩展形式：奇异值分解。
两者的区别在网上一搜一大把。一句话概括，特征值分解很好用，但只能用于长宽相等的方阵，对于非方阵的矩阵呢，就要使用奇异值分解了。实际上，特征值分解只是奇异值分解的一种特例。
直接上公式：
M = U x Σ x T(V)
同样的，M是原始矩阵。U和V是两个特征向量矩阵，分别称为“左奇异向量”和“右奇异向量”（其实分别是M x T(M)和T(M) x M方阵的特征向量，所以形状一大一小），T(V)即矩阵V的转置。Σ又是一个奇异值组成的对角矩阵（长宽不一致，所以会有一些空的行/列），而且还是已经由顶向下依据数值绝对值大小降序排列的。
关于这几个矩阵的大小，一般教科书是这样写的（网上查到的版本也都是这个）：
- 若M是m行n列矩阵
- 则U是m行m列矩阵
- 且Σ是m行n列矩阵
- 且V是n行n列矩阵
这样右侧的等式符合矩阵乘法要求的大小，最后会得到一个m行n列矩阵，与左边一致。如果在Numpy里验证这个数值也是如此（注意Numpy里返回的V是已经转置过的，无需再次转置）：
>>> import numpy as np
>>> M = np.mat([[1,2,3,4],[5,6,7,8],[2,3,4,5]])
>>> u, s, v = np.linalg.svd(M)
>>> print(M.shape)
(3, 4)
>>> print(u.shape)
(3, 3)
>>> print(s.shape)
(3,) # 根据矩阵乘法，这里实际表示的对角阵shape是(3, 4)
>>> print(v.shape)
(4, 4)
# 验证一下公式
>>> print(u.dot(np.column_stack((diag(s), np.zeros(3))).dot(v)))
[[ 1. 2. 3. 4.]
[ 5. 6. 7. 8.]
[ 2. 3. 4. 5.]]
此时TensorFlow再一次奇葩了！首先请注意它返回值中的u、s、v矩阵顺序，这是新手很容易被坑的地方。其次，它返回的矩阵大小十分的另类：
>>> import tensorflow as tf
>>> tf.InteractiveSession()
>>> M = tf.constant([[1,2,3,4],[5,6,7,8],[2,3,4,5]], dtype=tf.float32)
>>> s, u, v = tf.svd(M)
>>> print(M.shape)
(3, 4)
>>> print(u.shape)
(3, 3)
>>> print(s.shape)
(3,) # 根据矩阵乘法，这里实际表示的对角阵shape是(3, 3)
>>> print(v.shape)
(4, 3)
带进公式里计算一下，证实你没有看错。
>>> print(tf.matmul(tf.matmul(u, tf.diag(s)), tf.transpose(v)).eval())
[[ 1.00000107 2.00000095 3.00000143 4.00000191]
[ 5.00000191 6.00000191 7.00000286 8.00000286]
[ 2.00000048 3.00000095 4.00000095 5.00000095]]
不论怎样，我们现在得到了一个与公式结构一致的奇异值分解式。怎样提取特征矩阵的主成分呢？继续照猫画虎，不难发现，等式右侧的U是一个由正交向量组成的矩阵，而且Σ是对角矩阵，所以U x Σ是和特征值分解中的Q x Σ如出一辙。然后还是多出一个矩阵V，将两边同时乘以T(V)的转置，相当于把V挪到等号左边，等式变成：
M x I(T(V)) = U x Σ
接下来等号左边的部分就不用看了，等号右边的就是我们要的主成分矩阵。假设将原始矩阵的n维特征减少到r维，由于返回的s奇异值是已经排序的，直接取出最前面的r个值即可，同时相应裁剪u和v矩阵，得到：
u2 = u[:, :r]
s2 = s[:r]
v2 = v[:, :r]
M ≈ u2 x tf.diag(s2) x tf.transpose(v2) # 注意是约等于
代码写出来就是下面这个样子（忽略还原原始矩阵的计算）：
def pca_via_svd(M, r):
'''
M: 原始矩阵
r: 新的矩阵长度
'''
# 将输入矩阵做奇异值分解
s, u, v = tf.svd(M)
# 使用奇异值将矩阵的第二维压缩到指定长度
return tf.matmul(u[:, :r], tf.diag(s[:r]))
结论非常简单，两行代码的事情，只是为了摸索清楚相关的概念和最后这个API的用法，的确没少绕路。且行且分享，将踏过的路记录下来，让更多的人避免踩同样的坑。
网友评论