奇異值分解和圖像壓縮

from: http://cos.name/2014/02/svd-and-image-compression/

【2.18更新】：楠神寫了一個非常gelivable的Shiny應(yīng)用，用來動態(tài)展示圖片壓縮的效果隨k的變化情況。謝大大把這個應(yīng)用放到了RStudio的服務(wù)器上，大家可以點進(jìn)去玩玩看了。

=====================代表正義的分割線=====================

今天我們來講講奇異值分解和它的一些有意思的應(yīng)用。奇異值分解是一個非常，非常，非常大的話題，它的英文是 Singular Value Decomposition，一般簡稱為 SVD。下面先給出它大概的意思：

對于任意一個?m×nm×n?的矩陣?MM，不妨假設(shè)?m>nm>n，它可以被分解為

M=UDVTM=UDVT

其中

• UU?是一個?m×nm×n?的矩陣，滿足?UTU=InUTU=In，InIn?是?n×nn×n?的單位陣
• VV?是一個?n×nn×n?的矩陣，滿足?VTV=InVTV=In
• DD?是一個?n×nn×n?的對角矩陣，所有的元素都非負(fù)

先別急，我看到這個定義的時候和你一樣暈，感覺信息量有點大。事實上，上面這短短的三條可以引發(fā)出 SVD 許多重要的性質(zhì)，而我們今天要介紹的也只是其中的一部分而已。

前面的表達(dá)式?M=UDVTM=UDVT?可以用一種更容易理解的方式表達(dá)出來。如果我們把矩陣?UU?用它的列向量表示出來，可以寫成

U=(u1,u2,…,un)U=(u1,u2,…,un)

其中每一個?uiui?被稱為?MM?的左奇異向量。類似地，對于?VV，有

V=(v1,v2,…,vn)V=(v1,v2,…,vn)

它們被稱為右奇異向量。再然后，假設(shè)矩陣?DD?的對角線元素為?didi?（它們被稱為?MM的奇異值）并按降序排列，那么?MM?就可以表達(dá)為

M=d1u1vT1+d2u2vT2+?+dnunvTn=∑i=1ndiuivTi=∑i=1nAiM=d1u1v1T+d2u2v2T+?+dnunvnT=∑i=1ndiuiviT=∑i=1nAi

其中?Ai=diuivTiAi=diuiviT?是一個?m×nm×n?的矩陣。換句話說，我們把原來的矩陣?MM?表達(dá)成了?nn?個矩陣的和。

這個式子有什么用呢？注意到，我們假定?didi?是按降序排列的，它在某種程度上反映了對應(yīng)項?AiAi?在?MM?中的“貢獻(xiàn)”。didi?越大，說明對應(yīng)的?AiAi?在?MM?的分解中占據(jù)的比重也越大。所以一個很自然的想法是，我們是不是可以提取出?AiAi?中那些對?MM?貢獻(xiàn)最大的項，把它們的和作為對?MM?的近似？也就是說，如果令

Mk=∑i=1kAiMk=∑i=1kAi

那么我們是否可以用?MkMk?來對?Mn≡MMn≡M?進(jìn)行近似？

答案是肯定的，不過等一下，這個想法好像似曾相識？對了，多元統(tǒng)計分析中經(jīng)典的主成分分析就是這樣做的。在主成分分析中，我們把數(shù)據(jù)整體的變異分解成若干個主成分之和，然后保留方差最大的若干個主成分，而舍棄那些方差較小的。事實上，主成分分析就是對數(shù)據(jù)的協(xié)方差矩陣進(jìn)行了類似的分解（特征值分解），但這種分解只適用于對稱的矩陣，而 SVD 則是對任意大小和形狀的矩陣都成立。（SVD 和特征值分解有著非常緊密的聯(lián)系，此為后話）

我們再回顧一下，主成分分析有什么作用？答曰，降維。換言之，就是用幾組低維的主成分來記錄原始數(shù)據(jù)的大部分信息，這也可以認(rèn)為是一種信息的（有損）壓縮。在 SVD 中，我們也可以做類似的事情，也就是用更少項的求和?MkMk?來近似完整的?nn?項求和。為什么要這么做呢？我們用一個圖像壓縮的例子來說明我們的動機(jī)。

我們知道，電腦上的圖像（特指位圖）都是由像素點組成的，所以存儲一張 1000×622 大小的圖片，實際上就是存儲一個 1000×622 的矩陣，共 622000 個元素。這個矩陣用 SVD 可以分解為 622 個矩陣之和，如果我們選取其中的前 100 個之和作為對圖像數(shù)據(jù)的近似，那么只需要存儲 100 個奇異值?didi，100 個?uiui?向量和 100 個?vivi向量，共計 100×(1+1000+622)=162300個 元素，大約只有原始的 26% 大小。

【注：本文只是為了用圖像壓縮來介紹 SVD 的性質(zhì)，實際使用中常見的圖片格式（png，jpeg等）其壓縮原理更復(fù)雜，且效果往往更好】

為了直觀地來看看 SVD 壓縮圖像的效果，我們拿一幅 1000×622 的圖片來做實驗（圖片來源：http://www.bjcaca.com/bisai/show.php?pid=33844&bid=40）

SVD演示圖片，原圖SVD演示圖片，k=1SVD演示圖片，k=5SVD演示圖片，k=20SVD演示圖片，k=50SVD演示圖片，k=100

可以看出，當(dāng)取一個成分時，景物完全不可分辨，但還是可以看出原始圖片的整體色調(diào)。取 5 個成分時，已經(jīng)依稀可以看出景物的輪廓。而繼續(xù)增加?kk?的取值，會讓圖片的細(xì)節(jié)更加清晰；當(dāng)增加到 100 時，已經(jīng)幾乎與原圖看不出區(qū)別。

接下來我們要考慮的問題是，AkAk?是否是一個好的近似？對此，我們首先需要定義近似好壞的一個指標(biāo)。在此我們用?BB?與?MM?之差的 Frobenius 范數(shù)?||M–B||F||M–B||F?來衡量?BB?對MM?的近似效果（越小越好），其中矩陣的 Frobenius 范數(shù)是矩陣所有元素平方和的開方，當(dāng)其為 0 時，說明兩個矩陣嚴(yán)格相等。

此外，我們還需要限定?AkAk?的“維度”（否則?MM?就是它對自己最好的近似），在這里我們指的是矩陣的秩。對于通過 SVD 得到的矩陣?MkMk，我們有如下的結(jié)論：

在所有秩為?kk?的矩陣中，MkMk?能夠最小化與?MM?之間的 Frobenius 范數(shù)距離。

這意味著，如果我們以 Frobenius 范數(shù)作為衡量的準(zhǔn)則，那么在給定矩陣秩的情況下，SVD 能夠給出最佳的近似效果。萬萬沒想到啊。

在R中，可以使用?svd()?函數(shù)來對矩陣進(jìn)行 SVD 分解，但考慮到 SVD 是一項計算量較大的工作，我們使用了?rARPACK?包中的?svds()?函數(shù)，它可以只計算前?kk?項的分解結(jié)果。完整的 R 代碼如下：

library(rARPACK); library(jpeg);factorize = function(m, k) {r = svds(m[, , 1], k);g = svds(m[, , 2], k);b = svds(m[, , 3], k);return(list(r = r, g = g, b = b)); }recoverimg = function(lst, k) {recover0 = function(fac, k){dmat = diag(k);diag(dmat) = fac$d[1:k];m = fac$u[, 1:k] %*% dmat %*% t(fac$v[, 1:k]);m[m < 0] = 0;m[m > 1] = 1;return(m);}r = recover0(lst$r, k);g = recover0(lst$g, k);b = recover0(lst$b, k);m = array(0, c(nrow(r), ncol(r), 3));m[, , 1] = r;m[, , 2] = g;m[, , 3] = b;return(m); }rawimg = readJPEG("pic2.jpg"); lst = factorize(rawimg, 100); neig = c(1, 5, 20, 50, 100); for(i in neig) {m = recoverimg(lst, i);writeJPEG(m, sprintf("svd_%d.jpg", i), 0.95); }

參考文獻(xiàn)

Image Compression with the SVD in R

Foundations of Data Science

SVD維基頁面

《新程序員》：云原生和全面數(shù)字化實踐50位技術(shù)專家共同創(chuàng)作，文字、視頻、音頻交互閱讀

總結(jié)

以上是生活随笔為你收集整理的奇异值分解和图像压缩的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。