奇异值分解实在太好玩了

最近有个师姐想找我帮她把一张Arcgis地图中的颜色处理一下，突然勾起了我探索一下R中图像处理的兴趣，在网上看了几个图像处理的包，有个帖子写了奇异值分解在数据降维中的应用，突然发现自己好像学过奇异值分解这个东西，但又忘了一干二净，于是索性找来资料看看，却发现自己线性代数的东西完全看不懂了，幸好R代码还是看的懂的，于是用R中的奇异值分解做了个小东西玩玩，还挺有意思。

奇异值分解其实是对一个m*n矩阵进行的分解，对于任何一个这样的矩阵，它满足

svdformula

其中U是一个正交阵，Σ是一个对角阵，V也是一个正交阵。

在图像中运用奇异值分解十分直观，图像本身就是一个矩阵，在R中raster数据类型就是一种图像的栅格数据，稍微做一下转化就能变成一个大型矩阵。因此直接对图像矩阵做奇异值分解是十分简单的。

奇异值分解后形成的三个矩阵中，Σ对角线上的元素就是A的奇异值，对于Σ中的r个奇异值，满足如下关系：

svd2

很牛逼，只用r个奇异值就能大致还原一个矩阵，这么一来图像存储占的空间就少了许多，假如原图要500x500=250000的空间，奇异值分解后，用50个奇异值来存储，就只需要500x50+50x50+50x500=52500的空间了。实际上，我通过R模拟发现，大概50个左右的奇异值就能得到比较好的还原效果，真是神奇！

多说无用，在R里用代码实际操作一下，就知道奇异值分解在数据降维中的威力了，R中自带了svd()函数专门用来做矩阵的奇异值分解，所以整个模拟十分简单，几行代码就可以搞定，代码如下：

    pic_svd <- function(sourcename, howmany){
    pic.ras <- raster(sourcename) #获得raster数据
    n <- howmany
    pic.flip <- flip(pic.ras, direction = "y")
    pic.m <- t(as.matrix(pic.flip))
    pic.svd <- svd(pic.m)
    d <- diag(pic.svd$d)
    v <- as.matrix(pic.svd$v)
    u <- pic.svd$u
    u1 <- as.matrix(u[, 1:n])
    d1 <- as.matrix(d[1:n, 1:n])
    v1 <- as.matrix(v[, 1:n])
    pic_after <- u1 %*% d1 %*% t(v1)
    windows()
    image(pic_after, col = grey(seq(0, 1, length = 256)), main = n) 
    }