奇異值分解和圖像壓縮

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

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

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

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

對于任意一個?m×nm×n?的矩陣?MM，不妨假設?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?的對角矩陣，所有的元素都非負

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

為了直觀地來看看 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

可以看出，當取一個成分時，景物完全不可分辨，但還是可以看出原始圖片的整體色調。取 5 個成分時，已經依稀可以看出景物的輪廓。而繼續增加?kk?的取值，會讓圖片的細節更加清晰；當增加到 100 時，已經幾乎與原圖看不出區別。

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

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

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

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

在R中，可以使用?svd()?函數來對矩陣進行 SVD 分解，但考慮到 SVD 是一項計算量較大的工作，我們使用了?rARPACK?包中的?svds()?函數，它可以只計算前?kk?項的分解結果。完整的 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); }

參考文獻

Image Compression with the SVD in R

Foundations of Data Science

SVD維基頁面

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。