大きめのSparse行列をnormalizeしたい時の書き方

R

掲題の件、あると思います。

私は今まで行列っぽいデータを normalize するとき*1は以下のように書いていました。

最後の結果を見たらわかるように行ごとに和をとると１になるようになっています。

> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667
> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667

適当に作ったSparse行列でもこれで動くと言えば動いている。

> library("Matrix")
> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x
3 x 3 sparse Matrix of class "dgCMatrix"
          
[1,] 1 4 7
[2,] 2 . 8
[3,] . . 9
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

でも、行列のサイズがでかくなるとアウト、なんか一回Denseな行列にしてるんだろうと思われる。

ここの例は恣意的 データに見えるがもともと俺が使ってハマったデータがこんなんだったのでこうなっている。

> size <- 10^5
> i <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> j <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> x <- sparseMatrix(i,j,x=rnorm(3*10^7 + 1), dims=rep(size + 1, 2))
> x <- as(x + t(x), "symmetricMatrix")
> str(x)
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:29910604] 2 23 14 6 13 28 18 24 0 19 ...
  ..@ p       : int [1:100002] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ Dim     : int [1:2] 100001 100001
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:29910604] -1.511 -0.878 0.777 1.258 1.185 ...
  ..@ uplo    : chr "U"
  ..@ factors : list()
> x/Matrix::rowSums(x)
 エラー:  ベクトルのメモリを使い切りました (上限に達した？)

こういうときは等価な計算のこれを使おう、これならいける

> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
100001 x 100001 sparse Matrix of class "dgCMatrix"
                                                                                                                              
 [1,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [2,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [3,] . . . . . . . . . . . . . . . . . . . . . . -0.08169018 . . . . . . . . . . . . . . . . .  .        . . . . . . . ......

ちゃんと計算も合う（確認）

> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000
> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

元ネタ

*1:各行ごとに足したら１になる、でnormalizeを定義。他L2ノルムなどによるnormalizeも同様