谈谈SVD和LSA(zz)

原文地址：http://www.isnowfy.com/introduction-to-svd-and-lsa/

首先SVD和LSA是什么呢，SVD全称是singular value decomposition，就是俗称的奇异值分解，SVD的用处有很多，比如可以做PCA(主成分分析)，做图形压缩，做LSA，那LSA是什么呢，LSA全称Latent semantic analysis，中文的意思是隐含语义分析，LSA算是topic model的一种，对于LSA的直观认识就是文章里有词语，而词语是由不同的主题生成的，比如一篇文章包含词语计算机，另一篇文章包含词语电脑，在一般的向量空间来看，这两篇文章不相关，但是在LSA看来，这两个词属于同一个主题，所以两篇文章也是相关的。

特征值特征向量

要谈到SVD，特征值和特征向量是需要首先交代的。具体内容可以在wiki上看，这里我做个简单的介绍。对于方阵M如果有

M∗v=λ∗v

v是个向量， λ 是个数，那么我们称v是M的特征向量， λ 是M的特征向量，并且我们可以对M进行特征分解得到

M=Q∗Λ∗Q−1

其中Q是特征向量组成的矩阵， Λ 是对角阵，对角线上的元素就是特征值。对于特征的几何理解就是矩阵M其实是一种线性变换，而线性变换对于向量的影响有两种，旋转和拉伸，而特征向量就是在这种线性变换下方向保持不变的向量，但是长度还是会作相应的拉伸，特征值就是拉伸的程度。

从另一个角度说如果我们取特征值比较大的几项，那么就是对原矩阵做了一种近似。

M≈Q1..k∗Λ1..k∗Q−11..k

这样我们就可以用更少的元素去近似的表示原矩阵，但是特征分解的限制比较多，比如要求矩阵必须是方阵

奇异值分解

wiki是个好东西，你要想深入了解的话，建议还是去看wiki。奇异值分解是将矩阵变成了这样的形式

M=U∗Σ∗VT

其中 Σ 依旧是对角阵，而U和V是正交矩阵正交矩阵是说 U∗UT=I 。

我们还是先回到矩阵是线性变换这个思路上。

如果我们用M去作用空间里的一组基，那么我们就会得到另一组基，如上图那样。那么我们旋转一下最初的一组基。

这样我们经过M的变换由一组正交基变换到了另一组正交基上面。也是也就是下面这样。

也就是我们有

M∗v1=σ1∗u1
M∗v2=σ2∗u2

并且对于任意一个向量x，我们有

x=v1∗(vT1∗x)+v2∗(vT2∗x)

于是我们可以得到

M∗x=M∗v1∗(vT1∗x)+M∗v2∗(vT2∗x)
M∗x=σ1∗u1∗(vT1∗x)+σ2∗u2∗(vT2∗x)
M=σ1∗u1∗vT1+σ2∗u2∗vT2
M=U∗Σ∗VT

恩，我们得到了和特征值和特征向量相似的东西，SVD分解出来的就是在M的线性变换下，正交基变换仍是正交基，而奇异值就是拉伸的程度。其实SVD和特征值和特征向量的关系还是很大的。

M∗MT=U∗Σ∗VT∗V∗ΣT∗UT
M∗MT=U∗Σ2∗UT

也就是说SVD求出的是 M∗MT 和 MT∗M 的特征向量。同样的得到这SVD分解这种形式后我们就可以利用他来对原数据进行降维操作。

这里我们分别将RBG矩阵进行SVD，左上角的是原图，其他的依次是取最大的100个，50个，20个，10个，5个奇异值做的近似图像。

1. # -*- coding: utf-8 -*-
2.  
3. from scipy import linalg, dot
4. from PIL import Image
5.  
6. def main (num= 5 ):
7.     im = Image. open ( 'ai.jpg' )
8.     pix = im. load ( )
9.     ma = [ [ ], [ ], [ ] ]
10.     for x in xrange (im. size [ 0 ] ):
11.         for i in xrange ( 3 ):
12.             ma [i ]. append ( [ ] )
13.         for y in xrange (im. size [ 1 ] ):
14.             for i in xrange ( 3 ):
15.                 ma [i ] [ -1 ]. append (pix [x, y ] [i ] )
16.     for i in xrange ( 3 ):
17.         u, s, v = linalg. svd (ma [i ] )
18.         u = u [:, :num ]
19.         v = v [:num, : ]
20.         s = s [:num ]
21.         ma [i ] = dot (dot (u, linalg. diagsvd (s, num, num ) ), v )
22.     for x in xrange (im. size [ 0 ] ):
23.         for y in xrange (im. size [ 1 ] ):
24.             ret = [ ]
25.             for i in xrange ( 3 ):
26.                 tmp = int (ma [i ] [x ] [y ] )
27.                 if tmp < 0:
28.                     tmp = 0
29.                 if tmp > 255:
30.                     tmp = 255
31.                 ret. append (tmp )
32.             pix [x, y ] = tuple (ret )
33.     im. show ( )
34.     im. save ( 'test.jpg' )
35.  
36. if __name__ == '__main__':
37.     main ( )

如果对矩阵先进行归一化，再SVD就是PCA的形式了，这种形式可以用方差最大化或者误差最小化来求得，具体可以去看PCA相关的东西。所以和scturtle讨论了下直接SVD的意义，但是最后也没得出什么结论。。。

隐含语义分析

终于讲到最后的隐含语义分析了，首先我们构造文本和词语的矩阵，也就是对于矩阵来说每一个向量表示一篇文章，每个向量里就是单词的出现次数(更好的是每个是单词的tf/idf值，tf/idf不在赘述，具体可以看wiki)。那么SVD分解之后，我们就得到了降维的矩阵，就是下面这个样子

就是说原来我们有1000000篇文章，总共有500000个单词，我们保留最大的100个来做降维，于是现在我们可以这样理解，我们保留了100个主题，其中U是文章对应的主题分布，而V则是主题对应的词语的分布，这样，我们可以减少噪音，并且这样计算文章间的相关性也更加合理，并且可以把相关的单词聚合到一起。代码如下

1. # -*- coding: utf-8 -*-
2.  
3. import os
4. import re
5. import heapq
6. import codecs
7. from math import log
8. from scipy import linalg
9.  
10. import unigram_good_turing as seg
11.  
12. seg. init ( )
13.  
14. def tfidf (docs ):
15.     doclen = len (docs ) +1.0
16.     for doc in docs:
17.         wordtotal = sum (doc. values ( ) ) +0.0
18.         for word in doc:
19.             tf = doc [word ]/wordtotal
20.             idf = log (doclen/ ( sum ( [word in tmp for tmp in docs ] ) +1 ) )
21.             doc [word ] = tf *idf
22.     return docs
23.  
24. def solve (data ):
25.     re_zh, re_other = re. compile (ur "([\u4E00-\u9FA5]+)" ), re. compile (ur "[^a-zA-Z0-9+#\n]" )
26.     blocks = re_zh. split (data )
27.     for item in blocks:
28.         if re_zh. match (item ):
29.             for i in seg. solve (item ):
30.                 yield i
31.         else:
32.             tmp = re_other. split (item )
33.             for x in tmp:
34.                 if x != '':
35.                     pass
36.  
37. def show (dic, p ):
38.     p = heapq. nlargest ( 10, enumerate (p ), key= lambda x:x [ 1 ] )
39.     print ' '. join ( map ( lambda x:dic [x [ 0 ] ], p ) )
40.  
41. def main ( ):
42.     names = os. listdir ( 'text' )
43.     dic = { }
44.     cnt = 0
45.     ma = [ ]
46.     for name in names:
47.         data = codecs. open ( 'text/'+name, 'r', 'utf-8' ). read ( )
48.         doc = { }
49.         for word in solve (data ):
50.             if not word in dic:
51.                 dic [word ] = cnt
52.                 cnt += 1
53.             tmp = dic [word ]
54.             if tmp not in doc:
55.                 doc [tmp ] = 0
56.             doc [tmp ] += 1
57.         ma. append (doc )
58.     ma = tfidf (ma )
59.     ret = [ ]
60.     for item in ma:
61.         tmp = [ ]
62.         for i in xrange (cnt ):
63.             if i in item:
64.                 tmp. append (item [i ] )
65.             else:
66.                 tmp. append ( 0 )
67.         ret. append (tmp )
68.     u, s, v = linalg. svd (ret )
69.     for i in xrange ( 10 ):
70.         show ( dict ( zip (dic. values ( ), dic. keys ( ) ) ), list (v [i ] ) )
71.  
72. if __name__ == '__main__':
73.     main ( )

原文地址：http://www.isnowfy.com/introduction-to-svd-and-lsa/