Glophy Zone Backup

Python 的sparse matrix 入門範例 PDF 列印 E-mail
作者 曾正男   
2013/05/22, 週三

Sparse matrix 的意思是稀疏矩陣,多數的元素是0只有少數非零的元素。在作數值偏微分方程時很容易造出這樣的矩陣,現今很流行的Web service 推薦系統運算中也很容易遇上Sparse matrix。由於矩陣多數的元素是0,所以把整個矩陣先開在記憶體中來儲存這些少量的非零的元素是很不划算也沒有效率的。特別是在Web service 的應用上,矩陣的大小隨時會隨著服務的項目和人數而增加,所對應的矩陣大小會隨時改變。為了讓儲存Sparse matrix 的資料更有效率,多半的作法是儲存(i,j,v) 元組的方式,其中(i,j) 是矩陣元素位置的index ,v 是元素所對應的值。這樣的儲存方式會大幅增加儲存的效率,然而在矩陣運算上便必需相對應改變。  

Python 當然也會有關於Sparse matrix 的套件,在http://docs.scipy.org/doc/scipy/reference/sparse.html  有2D sparse matrix 基本的介紹。這一個套件是依附在 scipy 這個運算套件底下,有作Python 數值計算的人通常都會安裝到scipy 和numpy 這兩個計算套件。以下是最基本的中文使用範例說明:

要宣告一個矩陣是sparse matrix 我們會利用到 scipy.sparse 底下的 lil_matrix 函數。 用法如下

from scipy.sparse import lil_matrix

from numpy.random import rand

A = lil_matrix((1000, 1000))

A[0, :100] = rand(100)

A[1, 100:200] = A[0, :100]

A.setdiag(rand(1000)) 

在lil_matrix 函數裡的tuple (1000,1000) 是控制矩陣的大小,當A=lil_matrix((1000,1000)) 指令下好後,我們可以從A.shape 看到A 矩陣已經是一個1000 x 1000 的矩陣。要小心lil_matrix 內不要用[1000,1000] 的方式來作為控制矩陣大小的方法。因為lil_matrix 兼具矩陣轉換型別的功能,B=lil_matrix([1000,1000]) 這個指令,Python 會視為要將[1000,1000] 這一個1 x 2 的矩陣從標準矩陣型別轉換成稀疏矩陣型別,因此B.shape 的結果是 (1,2) 而不是 (1000,1000)。開出稀疏矩陣之後,我們可以用一般在Python 中的矩陣操作方法對矩陣進行控制,A[0,:100] 就是A 矩陣第一行的前100列元素,對個別元素的操作是A[i,j] 的方式來控制。A.setdiag 是將串列直接放入A 的對角線,這在解數值偏微分方程時會很經常使用。

在解sparse matrix 時會用到scipy.sparse.linalg  套件底下的spsolve 函數,並且在求解之前要先將sparse matrix 轉換成特定的sparse format 。例如,

from scipy.sparse.linalg import spsolve

A = A.tocsr()

b = rand(1000)

x = spsolve(A, b)

sparse matrix 當然也可以轉換成一般的矩陣型別,我們用A.todense() 這個函數來達成,程式如下,

B = A.todense()

在前言裡我們提到用sparse matrix 的儲存方式可以避免矩陣大小隨時擴增的困擾,COO 格式的儲存法正好符合我們的需要,我們只要將對應的(i,j,v) 元素放進一個專門紀錄個別i,j,v 的陣列I,J,V 中,就可以很容易的開出對應的sparse matrix,範例如下

from scipy.sparse import coo_matrix

from numpy import array

I = array([0,0,1,3,1,0,0])

J = array([0,2,1,3,1,0,0])

V = array([1,1,1,1,1,1,1])

B = sparse.coo_matrix((V,(I,J)),shape=(4,4)).tocsr()

這個範例中I,J 存放的是index,而V 存放的是對應元素的值,我們看到有許多位置是重複的,在最後一行的B 矩陣的定義,我們在最後做了一個tocsr() 格式的轉換,這是把B 矩陣用CSR format 來表示(我們可以選擇使用CSR 或CSC 格式),對於重複的位置他所對應的值會用累加的方式來顯示。這樣的特性很適合Web service 來紀錄網頁活動。 

 

最後更新 ( 2013/05/22, 週三 )
 
下一個 >