奇異值分解和圖像壓縮

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

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

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

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

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

M=UDVTM=UDVT

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

我們?cè)倩仡櫼幌?#xff0c;主成分分析有什么作用？答曰，降維。換言之，就是用幾組低維的主成分來(lái)記錄原始數(shù)據(jù)的大部分信息，這也可以認(rèn)為是一種信息的（有損）壓縮。在 SVD 中，我們也可以做類似的事情，也就是用更少項(xiàng)的求和?MkMk?來(lái)近似完整的?nn?項(xiàng)求和。為什么要這么做呢？我們用一個(gè)圖像壓縮的例子來(lái)說(shuō)明我們的動(dòng)機(jī)。

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

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

為了直觀地來(lái)看看 SVD 壓縮圖像的效果，我們拿一幅 1000×622 的圖片來(lái)做實(shí)驗(yàn)（圖片來(lái)源：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)取一個(gè)成分時(shí)，景物完全不可分辨，但還是可以看出原始圖片的整體色調(diào)。取 5 個(gè)成分時(shí)，已經(jīng)依稀可以看出景物的輪廓。而繼續(xù)增加?kk?的取值，會(huì)讓圖片的細(xì)節(jié)更加清晰；當(dāng)增加到 100 時(shí)，已經(jīng)幾乎與原圖看不出區(qū)別。

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

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

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

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

在R中，可以使用?svd()?函數(shù)來(lái)對(duì)矩陣進(jìn)行 SVD 分解，但考慮到 SVD 是一項(xiàng)計(jì)算量較大的工作，我們使用了?rARPACK?包中的?svds()?函數(shù)，它可以只計(jì)算前?kk?項(xiàng)的分解結(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維基頁(yè)面

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

總結(jié)

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

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