Rで疎な多次元配列を扱う（Matrix,slamパッケージ）

概要

Rで大規模テンソル解析がしたい。
そのために、slamパッケージを使う。
slamパッケージとは、疎な行列、配列を扱うためのパッケージである。
Matrixパッケージ（疎行列）やffパッケージ（ストレージの利用）とは異なり、疎な多次元配列を扱うことができる。
保持する値が2^31以下なら、配列自体の大きさはいくらでも大きくできる。
たぶんRで疎なarrayが使えるのはslamだけなので、使ってみる。

結論

配列の密度次第ではあるが、
これくらい軽く

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(300,400,300)) 
> object.size(as.array(data))
288000208 bytes
> object.size(data)
1096 bytes

これくらい速い

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(3000,400,300))
> tmp <- Matrix(data[,,1],nrow=dim(data[,,1])[1])
> dim(tmp)
[1] 3000  400
> system.time(hoge<-tmp %*% t(tmp))
   user  system elapsed 
  0.002   0.000   0.004 
> tmp <- matrix(data[,,1],nrow=dim(data[,,1])[1])
> system.time(hoge<-tmp %*% t(tmp))
   user  system elapsed 
  6.482   0.073   6.724

ただ、slamは行列積すらまともに利用できないので、疎なarrayとして保持するためだけに利用して、最終的にMatrixパッケージで積を処理することにした。
テンソルの演算のために行列化が必要な場所は、都度その関数を書いている。
果たしてこれで良かったのかは謎。

基本関数

simple_sparse_array

疎配列の作成関数
i:配列の要素がある位置を行列で指定する
v:配列の要素の値をベクトルで指定する
dim:配列の次元数をベクトルで指定する

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3),nrow=3,ncol=3,byrow=T),v=c(1,2,3),dim=c(3,4,5))
> data
A simple sparse array of dimension 3x4x5.
> as.array(data)
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    2    0    0
[3,]    0    0    0    0

, , 3

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    3    0

（以下省略

> str(data)
List of 4
 $ i       : int [1:3, 1:3] 1 2 3 1 2 3 1 2 3
 $ v       : num [1:3] 1 2 3
 $ dim     : int [1:3] 3 4 5
 $ dimnames: NULL
 - attr(*, "class")= chr "simple_sparse_array"

abind_simple_sparse_array

配列を結合させる関す
...:結合させる配列（いくつでも可）
MARGIN:結合させる方向(row=1,col=2,...)

> data <- abind_simple_sparse_array(array(1:8,c(2,2,2)), array(9:12, c(2,2)),MARGIN=3) 
> data
A simple sparse array of dimension 2x2x3.
> as.array(data)
, , 1

     [,1] [,2]
[1,]    1    3
[2,]    2    4

, , 2

     [,1] [,2]
[1,]    5    7
[2,]    6    8

, , 3

     [,1] [,2]
[1,]    9   11
[2,]   10   12

extend_simple_sparse_array

配列のインデックスを増加させる関数
x:配列
MARGIN:インデックスを追加する位置

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(10,11,12))
> data
A simple sparse array of dimension 10x11x12.
> extend_simple_sparse_array(data,MARGIN=0)
A simple sparse array of dimension 1x10x11x12.

rollup

疎配列用のapply関数
x:配列
MARGIN:applyの方向
FUN:関数

まともな計算がしたいなら、このFUNを自作せざるをえない？

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(3,4,3)) 
> as.array(data)
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    2    5    0
[3,]    0    0    0    0

, , 3

     [,1] [,2] [,3] [,4]
[1,]    0   10    0    0
[2,]    0    0    0    0
[3,]    0    0    3    0

> as.array(rollup(x=data,MARGIN=3,FUN=sum))
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1   10    0    0
[2,]    0    2    5    0
[3,]    0    0    3    0

#行名がつく
> as.array(rollup(x=data,MARGIN=1,INDEX=c('a','b','c')))                                                                                                                             
, , 1

  [,1] [,2] [,3] [,4]
a    1    0    0    0
b    0    0    0    0
c    0    0    0    0

, , 2

  [,1] [,2] [,3] [,4]
a    0    0    0    0
b    0    2    5    0
c    0    0    0    0

, , 3

  [,1] [,2] [,3] [,4]
a    0   10    0    0
b    0    0    0    0
c    0    0    3    0

as.sparse.array

普通の配列を疎配列化する関数
x:データの配列

> A <- matrix(round(runif(3*4)),3,4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    0    1    0
[2,]    1    0    0    1
[3,]    0    0    0    1
> as.simple_sparse_array(A)
A simple sparse array of dimension 3x4.

tcrossprod_simple_triplet_matrix

二つの疎行列の内積を取る関数
x:行列
y:行列

> A <- as.simple_triplet_matrix(A)
> A
A 3x4 simple triplet matrix.
> as.matrix(A)
     [,1] [,2] [,3] [,4]
[1,]    1    0    1    0
[2,]    1    0    0    1
[3,]    0    0    0    1
> tcrossprod_simple_triplet_matrix(A)
     [,1] [,2] [,3]
[1,]    2    1    0
[2,]    1    2    1
[3,]    0    1    1
> tcrossprod_simple_triplet_matrix(t(A))
     [,1] [,2] [,3] [,4]
[1,]    2    0    1    1
[2,]    0    0    0    0
[3,]    1    0    1    0
[4,]    1    0    0    2

Matrix

疎行列に変換する関数（Matrixパッケージ）
x:配列
nrow:行数

cBind

疎行列を列方向に結合する関数（Matrixパッケージ）
...=:行列

> data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(10,11,12))
> dim(data)
[1] 10 11 12
> tmp <- apply(data,MARGIN=3,FUN=Matrix)                                                                                                                                           
> dim(tmp[[1]])
[1] 10 11
> length(tmp)
[1] 12
> dim(do.call(args=tmp,what=cBind))
[1]  10 132

apply系はRの配列数は2^31未満ルールが適応されるので、arrayが大きすぎる場合はfor文でやらざるをえない。
したがって並列計算はdoMCパッケージが良いと思われる。

registerDoMC(10)
data <- simple_sparse_array(i=matrix(c(1,1,1,2,2,2,3,3,3,1,2,3,2,3,2),ncol=3,byrow=T),v=c(1,2,3,10,5),dim=c(3000,4000,3000)) 
tmp <- foreach(i = 1:3000) %dopar% {Matrix(data[,,i],nrow=3000)}

なんか色々汚い。Cで書きたい。