遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

使ってみる〜疎行列計算、SuiteSparse for R

R Matrix SuiteSparse 疎行列
  • さて、うんちくは前の記事に任せて、使ってみることにしよう
  • まずは、密行列と疎行列のうちの、密行列の扱い
  • 密行列
    • 形に応じて、７種類の密行列オブジェクトを用意している
      • dgeMatrix:普通の実数行列を普通に保管
      • dsyMatrix:実対称行列を圧縮しないで保管
      • dspMatrix:実対称行列を圧縮保管(三角部分の保管)
      • dtrMatrix:実三角行列の非圧縮保管
      • dtpMatrix:実三角行列の圧縮保管(三角部分のみ)
      • dpoMatrix:Positive semi-definite symmetric realを非圧縮保管
      • dppMatrix:Positive semi-definite symmetric realを圧縮保管
    • この形式を、計算負荷が小さくなるように、適宜強制変換しながら演算(行列の積、cross積、行列norm、reciprocal condition number(数値計算精度を保つための指標)、LU factorization,これスキー分解、逆行列)を実施するように作ってある
    • 分解計算などでは、分解成分を使い捨てにせず、再利用することを前提に保管しながら計算を進める
  • 疎行列
    • 疎行列保管の基礎
      • 非ゼロのセルを、行・列・値の３つ組で保管
      • 要素の並べ方は、行でのつながりを優先するか(csr:compressed sparse row)、列でのそれを優先するか(csc:compressed sparse column)、で二通り定められる
        • RとMatlabは密行列の保管において、column優先システムを採用
    • 形に応じて、４種類の疎行列オブジェクトを用意している
      • dgTMatrix: ３つ組の原則での保管形式(順序は適当？圧縮していない？)
      • dgCMatrix: csc順に並び替えて(csc形式)の疎行列圧縮保管
      • dsCMatrix: 実対称疎行列をcsc形式で圧縮保管
      • dtCMatrix: 三角実疎行列をcsc形式で圧縮保管
      • ddiMatrix: 対角行列の型もある
  • 使ってみる
  • Matrix()関数に適当に0がちなデータを与えると、dgeMatrix形式(密行列の形式)としたり、dgCMatrix(疎行列の形式としたりする)
for(i in 1:6){
print(Matrix(sample(0:1,28,replace=TRUE,prob=c(0.5,0.5)),4,7))
}

4 x 7 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    1    1    0    1    0
[2,]    1    0    1    1    1    1    1
[3,]    1    1    0    1    0    0    0
[4,]    0    1    1    1    1    0    0
4 x 7 sparse Matrix of class "dgCMatrix"
                  
[1,] 1 1 . . 1 1 1
[2,] 1 1 1 . 1 . .
[3,] . . . 1 1 1 .
[4,] . . . . 1 . .
4 x 7 sparse Matrix of class "dgCMatrix"
                  
[1,] . 1 1 1 1 1 .
[2,] . . 1 1 1 1 .
[3,] . 1 . 1 1 . .
[4,] . . . 1 . . .
4 x 7 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    1    1    1    0    0
[2,]    1    0    0    1    1    1    1
[3,]    0    1    0    0    1    0    1
[4,]    1    1    0    0    0    1    0
4 x 7 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    0    1    1    0    1
[2,]    1    0    0    1    0    1    0
[3,]    0    1    0    0    0    1    1
[4,]    0    1    1    1    1    0    0
4 x 7 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    1    0    1    1    0
[2,]    1    0    0    0    0    1    0
[3,]    0    0    1    0    1    1    1
[4,]    0    1    0    1    1    1    1
  • 対称行列なら、dsyMatrix,dsCMatrix形式になる
for(i in 1:6){
	m <- matrix(sample(0:1,25,replace=TRUE,prob=c(0.8,0.2)),5,5)
	m <- m+t(m)

	print(Matrix(c(m),5,5))
}

5 x 5 sparse Matrix of class "dsCMatrix"
              
[1,] . . . 1 1
[2,] . . 1 . .
[3,] . 1 . . .
[4,] 1 . . . .
[5,] 1 . . . 2
5 x 5 sparse Matrix of class "dsCMatrix"
              
[1,] . . . . .
[2,] . . 1 . 1
[3,] . 1 . 1 .
[4,] . . 1 . 1
[5,] . 1 . 1 .
5 x 5 Matrix of class "dsyMatrix"
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    1    1    1
[2,]    0    0    0    0    1
[3,]    1    0    2    1    1
[4,]    1    0    1    2    0
[5,]    1    1    1    0    0
5 x 5 sparse Matrix of class "dsCMatrix"
              
[1,] . 1 2 . .
[2,] 1 . . . 1
[3,] 2 . . . .
[4,] . . . 2 1
[5,] . 1 . 1 2
5 x 5 sparse Matrix of class "dsCMatrix"
              
[1,] 2 . 1 . .
[2,] . 2 1 . .
[3,] 1 1 . . .
[4,] . . . 2 1
[5,] . . . 1 .
5 x 5 Matrix of class "dsyMatrix"
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    1    0    0
[2,]    1    2    1    1    1
[3,]    1    1    0    0    0
[4,]    0    1    0    0    1
[5,]    0    1    0    1    0
  • 同様に三角行列ならdtrMatrix,dtCMatirx形式になる
for(i in 1:6){
	m <- matrix(sample(0:1,25,replace=TRUE,prob=c(0.2,0.8)),5,5)
	m[lower.tri(m)] <- 0

	print(Matrix(c(m),5,5))
}

5 x 5 sparse Matrix of class "dtCMatrix"
              
[1,] 1 . 1 1 1
[2,] . . . 1 1
[3,] . . 1 . 1
[4,] . . . 1 1
[5,] . . . . 1
5 x 5 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    1    1
[2,]    .    1    1    1    1
[3,]    .    .    1    1    1
[4,]    .    .    .    1    1
[5,]    .    .    .    .    1
5 x 5 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    .    1    1    0    1
[3,]    .    .    1    1    1
[4,]    .    .    .    1    0
[5,]    .    .    .    .    1
5 x 5 sparse Matrix of class "dtCMatrix"
              
[1,] 1 1 1 1 .
[2,] . 1 1 . 1
[3,] . . 1 1 1
[4,] . . . 1 .
[5,] . . . . 1
5 x 5 Matrix of class "dtrMatrix"
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    .    1    1    1    0
[3,]    .    .    0    1    1
[4,]    .    .    .    1    1
[5,]    .    .    .    .    1
5 x 5 sparse Matrix of class "dtCMatrix"
              
[1,] 1 1 . . 1
[2,] . . 1 . .
[3,] . . 1 1 1
[4,] . . . 1 1
[5,] . . . . 1
  • 実際、igraphパッケージの隣接行列はdgCMatrix形式も使って返ってくるようだ
> library(igraph)
> n.e <- 20
> n.v <- 10
> el <- cbind(sample(1:n.v,n.e,replace=TRUE),sample(1:n.v,n.e,replace=TRUE))
> g <- graph.edgelist(el)
> m <- get.adjacency(g)
> m
10 x 10 sparse Matrix of class "dgCMatrix"
                         
 [1,] 1 . 1 . 2 . . . . 1
 [2,] . . 1 . . 1 . . . .
 [3,] 1 . . . . 1 1 . . .
 [4,] . . . . . 2 . . . .
 [5,] . . 1 . . . . . . .
 [6,] . . . . . 1 . 1 . .
 [7,] . . . . . . . . . .
 [8,] . 1 . . 1 . . . . .
 [9,] . . 1 . . . . . . .
[10,] 2 . . . . . . . . .

> n.e <- 1000
> n.v <- 10
> el <- cbind(sample(1:n.v,n.e,replace=TRUE),sample(1:n.v,n.e,replace=TRUE))
> g <- graph.edgelist(el)
> m <- get.adjacency(g)
> m
10 x 10 sparse Matrix of class "dgCMatrix"
                                   
 [1,]  6 11 11  4 13 12 15  7  9  8
 [2,] 11  8 13 11  6  6  6 14 11 16
 [3,] 11  8 14  6  9  7 16 13 17 11
 [4,]  8 14  9 14 11 12 14 14 15 10
 [5,] 14  7 10  5  7  7 11 10 14  9
 [6,] 10  9  5 10  6  9 13 11  5  5
 [7,]  7  7 13  8  9  7  9 12  4 11
 [8,] 10  4 11  8  8  9 10  8 12  7
 [9,] 11  7 14  8 11 13 11 12 14  7
[10,] 10 10 16 10 13  7  9 12 10 13
  • 密行列形式と疎行列形式との強制変換
library(igraph)
n.e <- 20
n.v <- 10
el <- cbind(sample(1:n.v,n.e,replace=TRUE),sample(1:n.v,n.e,replace=TRUE))
g <- graph.edgelist(el)
m <- get.adjacency(g)
m
m.dens <- as(m,"denseMatrix")
m.dens
m2 <- as(m.dens,"sparseMatrix")
m2

> m
10 x 10 sparse Matrix of class "dgCMatrix"
                         
 [1,] 2 . 2 . . . . . . .
 [2,] . . . . . . . . . .
 [3,] . . 1 1 . . . . . .
 [4,] . 1 2 1 . 1 . . . .
 [5,] . 1 . . 1 . . . 1 .
 [6,] . . . . 1 . . . . .
 [7,] . . . . . . 1 1 . .
 [8,] . . . . . . . 1 . .
 [9,] . 1 . . . . . . . .
[10,] . 1 . . . . . . . .
> m.dens <- as(m,"denseMatrix")
> m.dens
10 x 10 Matrix of class "dgeMatrix"
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    2    0    2    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    1    1    0    0    0    0    0     0
 [4,]    0    1    2    1    0    1    0    0    0     0
 [5,]    0    1    0    0    1    0    0    0    1     0
 [6,]    0    0    0    0    1    0    0    0    0     0
 [7,]    0    0    0    0    0    0    1    1    0     0
 [8,]    0    0    0    0    0    0    0    1    0     0
 [9,]    0    1    0    0    0    0    0    0    0     0
[10,]    0    1    0    0    0    0    0    0    0     0
> m2 <- as(m.dens,"sparseMatrix")
> m2
10 x 10 sparse Matrix of class "dgCMatrix"
                         
 [1,] 2 . 2 . . . . . . .
 [2,] . . . . . . . . . .
 [3,] . . 1 1 . . . . . .
 [4,] . 1 2 1 . 1 . . . .
 [5,] . 1 . . 1 . . . 1 .
 [6,] . . . . 1 . . . . .
 [7,] . . . . . . 1 1 . .
 [8,] . . . . . . . 1 . .
 [9,] . 1 . . . . . . . .
[10,] . 1 . . . . . . . .
  • 疎行列形式の値の格納方法はどうなっているのか？
    • 以下に示すように、m2@iには、各列について、第何行に0以外の値が入っているのかが格納されている(0始まり。第0列は第0行、第1列は第3,4,8,9列…。またm2@pには、0以外の値が通しでベクトル化してあるのだが、そのベクトルの番地のうち、各列の最初の値に相当するのはいくつか、が格納されている。最初の0は、ベクトルの第0番地が、第0列の最初の値、と宣言してあり、次の1はベクトルの第1番地が、第1列の最初の値と宣言してある、最初の値が0で次の値が1だから、第1列には１個しか非0の値が入っていないこともそこからわかる。第2列には4個の非0があるので、@pの第3番の値は、1+4=5
> str(m2)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:17] 0 3 4 8 9 0 2 3 2 3 ...
  ..@ p       : int [1:11] 0 1 5 8 10 12 13 14 16 17 ...
  ..@ Dim     : int [1:2] 10 10
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:17] 2 1 1 1 1 2 1 2 1 1 ...
  ..@ factors : list()
> m2
10 x 10 sparse Matrix of class "dgCMatrix"
                         
 [1,] 2 . 2 . . . . . . .
 [2,] . . . . . . . . . .
 [3,] . . 1 1 . . . . . .
 [4,] . 1 2 1 . 1 . . . .
 [5,] . 1 . . 1 . . . 1 .
 [6,] . . . . 1 . . . . .
 [7,] . . . . . . 1 1 . .
 [8,] . . . . . . . 1 . .
 [9,] . 1 . . . . . . . .
[10,] . 1 . . . . . . . .
> m2@i
 [1] 0 3 4 8 9 0 2 3 2 3 4 5 3 6 6 7 4
> m2@p
 [1]  0  1  5  8 10 12 13 14 16 17 17
  • ３つ組表示。疎行列は行番号、列番号、値の３つ組が基本だが、その形式で表示するには
> summary(m2)
10 x 10 sparse Matrix of class "dgCMatrix", with 17 entries 
    i j x
1   1 1 2
2   4 2 1
3   5 2 1
4   9 2 1
5  10 2 1
6   1 3 2
7   3 3 1
8   4 3 2
9   3 4 1
10  4 4 1
11  5 5 1
12  6 5 1
13  4 6 1
14  7 7 1
15  7 8 1
16  8 8 1
17  5 9 1
  • 疎行列ならばimage()関数でプロットするのも軽い
    • さすがに、以下の例(１辺の長さが１の16次元超立方格子グラフの隣接行列を疎行列で表したもののimage()プロット)

library(Matrix)
library(igraph)
g.lattice <- graph.lattice(length=2, dim=16)
m.lattice <- get.adjacency(g.lattice)
image(m.lattice)
  • この規則的な多次元単位格子グラフの隣接行列は、1次元の場合をスタートして、それを疎行列として、単位行列と列バインド・行バインドすることを再帰的に行うことによって作ることができる
    • k次元を一つ上げるときには、k次元格子２つを上下に重ねて、対応する点ペア間にエッジを張ればよい。k次元格子の行列を対角に置くのが、２つのk次元格子を持ってくることに相当し、残りのスペースに単位行列を置くのが、２つのk次元格子の対応点間エッジを張ることに相当する
    • cBind(),rBind()という疎行列用の連結関数を使っていることに注意
library(Matrix)
library(igraph)
# 地道に作る
dfs <- 1:8
Lattice.seq.jimichi <- list()
for(i in 1:length(dfs)){
	g.lattice <- graph.lattice(length=2, dim=dfs[i])
	Lattice.seq.jimichi[[i]] <- get.adjacency(g.lattice)
}
g.lattice <- graph.lattice(length=2, dim=8)
m.lattice <- get.adjacency(g.lattice)
image(m.lattice)

# 再帰的に作る
my.serial.lattice2 <- function(X){
	tmp <- cBind(X,Matrix(diag(rep(1,nrow(X)))))
	rBind(tmp, cBind(Matrix(diag(rep(1,nrow(X)))),X))
}
g.lattice <- graph.lattice(length=2, dim=1)
m.lattice1 <- get.adjacency(g.lattice)

Lattice.seq <- list()
Lattice.seq[[1]] <- m.lattice1
for(i in 2:length(dfs)){
	Lattice.seq[[i]] <- my.serial.lattice2(Lattice.seq[[i-1]])
}
  • さらに言えば、入れ子になっているので、各次元の行列をいちいち作る必要はなくて、最大次元の行列を作り、その部分を取り出せばよいことになる
library(Matrix)
library(igraph)
# 地道に作る
dfs <- 1:8
Lattice.seq.jimichi <- list()
for(i in 1:length(dfs)){
	g.lattice <- graph.lattice(length=2, dim=dfs[i])
	Lattice.seq.jimichi[[i]] <- get.adjacency(g.lattice)
}
g.lattice <- graph.lattice(length=2, dim=8)
m.lattice <- get.adjacency(g.lattice)
image(m.lattice)

# 再帰的に作る
my.serial.lattice2 <- function(X){
	tmp <- cBind(X,Matrix(diag(rep(1,nrow(X)))))
	rBind(tmp, cBind(Matrix(diag(rep(1,nrow(X)))),X))
}
g.lattice <- graph.lattice(length=2, dim=1)
m.lattice1 <- get.adjacency(g.lattice)

Lattice.seq <- list()
Lattice.seq[[1]] <- m.lattice1
dfs <- 1:8
for(i in 2:length(dfs)){
	Lattice.seq[[i]] <- my.serial.lattice2(Lattice.seq[[i-1]])
}

# 再帰的に作るけれど、保管は最大次元のものだけにして
Lattice.seq. <- m.lattice1

for(i in 1:10){
	Lattice.seq. <- my.serial.lattice2(Lattice.seq.)
}
# 必要に応じて取り出す
Lattice.seq.3 <- list()
for(i in 1:8){
	Lattice.seq.3[[i]] <- Lattice.seq.[1:2^i,1:2^i]
}

for(i in 1:8){
	print(identical(Lattice.seq[[i]],Lattice.seq.3[[i]]))
}