處理 R 的 sparse matrix 問題

遇到 sparse matrix 是很棘手的事情。原因是一不小心它會使檔案變得肥大，佔盡了記憶體資源，同時拖垮了 script 的效率。

真的是很棘手的敵人。

幸虧已經有許多人研究、提出了一些解決辦法。比如說 R 的 Matrix 套件。



了解格式

---

比較推薦使用的是 MatrixMarket 格式。

5  5  8
1  1  1.000e+00
2  2  1.050e+01
3  3  1.500e-02

第一行記錄的分別是 rows, columns, entries (列、行、總數)。接著每一列表示每一個 entry。可加入註解，行首以 % 標示，放在檔案的開頭處。更詳細的說明可以見此網頁。


建立 sparse matrix

---

可參考底下的範例，試著在 R 重複一遍(鍵入 ?sparseMatrix 查更多說明)。

i <- c(1,3:8); j <- c(2,9,6:10); x <- 7 * (1:7)
(A <- sparseMatrix(i, j, x = x))

輸出結果

8 x 10 sparse Matrix of class "dgCMatrix"
                             
[1,] . 7 . . .  .  .  .  .  .
[2,] . . . . .  .  .  .  .  .
[3,] . . . . .  .  .  . 14  .
[4,] . . . . . 21  .  .  .  .
[5,] . . . . .  . 28  .  .  .
[6,] . . . . .  .  . 35  .  .
[7,] . . . . .  .  .  . 42  .
[8,] . . . . .  .  .  .  . 49

讀、寫外部 sparse matrix

---

我們可以藉著 read.table() 讀入 MatrixMarket 格式的檔案，再利用 sparseMatrix() 製作 sparse matrix。筆者相信應該還有更好的方法，讀者可以自行嘗試其他指令。底下的範例亦僅供參考。

讀

dd <- read.table("D:/test.mtx", header=FALSE)
mm <- sparseMatrix(i = dd$i, j = dd$j, x = dd$x, dims=c(row_dim,col_dim))

寫
如果你已經擁有了一個 sparse matrix，使用 writeMM() 可以幫您輸出成為外部檔案

data(KNex)
writeMM(KNex$mm, "mmMM.mtx")

這部份可鍵入 ?readMM 看更多說明。

閱讀更多 »

用 Jensen-Shannon divergence 計算相似度

今天認識到了 Jensen-Shannon divergence，它似乎是一個常用於計算兩個機率分佈間之相似程度的方法。

where  and

 is Kullback–Leibler divergence

根據維基百科的條目，它又被稱為 information radius (IRad) 或 total divergence to the average。儘管是筆者自己的猜測，但我想後者名字的由來應該與 M = 1/2(P+Q) 這個平均值脫不了干係。

2009 年 PNAS 刊了一篇充分運用 Jensen-Shannon divergence 方法以比較多種基因體的文章。對於想運用 Jensen-Shannon divergence 以計算 DNA 序列間相似程度的話，裡頭的公式與步驟應已經夠詳細與清楚了。

標題：Alignment-free genome comparison with featurefrequency profiles (FFP) and optimal resolutions



另外，筆者發現 R 有一個名為 textcat 的套件，裡頭含有可以計算 Jensen-Shannon divergence 的函式 ─ textcat_xdist()。這個函式的優點在於，給定 n-gram profiles，它還可以計算 Kullback-Leibler I-divergence, Kullback-Leibler J-divergence, the sum of the absolute differences in n-gram log frequencies, 以及 cosine dissimilarity 等等。

雖然有現成的 R 套件可以使用，但也不需要高興的太早。依照筆者這幾天的經驗，資料量太大(一千個維度以上) 的情況就開始可能會有「記憶體不足」的錯誤發生，差不多的意思就是這函式報廢了。

所以最好還是從頭了解，弄清楚每個算式，別擔心把手弄髒。如此一來最壞的情況下也有辦法實作演算法，把任務完成。



提醒

---

若要得到距離，記得將 Jensen-Shannon divergence 開平方根

The Jensen-Shannon divergence is not a distance (as it  does not obey the triangle inequality), but its square root is. ─ "What Makes a Query Difficult?" Carmel et al., ACM, 2006. 

閱讀更多 »