Rで最大部分列和問題 - もうカツ丼はいいよな

最大部分列和(Maximum Segment Sum、略してMSS)問題とは、与えられた整数列の部分列の和のうち最大のものを求めるという問題。
非常に簡単な例で言うと、

a = {-1, -1, 1, 1, 1, -1, -1}

という数列の最大部分列和mss(a)は{1, 1, 1}の和の3になる。
他にも、

## example
a1 <- c(8, -3, 65, -20, -45, 100, 8, -17, 4, -14)
## mss(a1) => 113 = sum(8, -3, 65, -20, -45, 100, 8)
a2 <- c(31, -41, 59, 26, -53, 58, 97, -93, -23, 84)
## mss(a2) => 187 = sum(59, 26, -53, 58, 97)

といった具合になる。

全ての部分列を求めて比較する方法

すぐ思いつくやり方は、全ての部分列を列挙し、それぞれの総和を計算して一番大きなものを探すというやり方だろう。

mss1 <- function(a){
  n <- length(a)
  s <- 0
  for(i in 1:n){
    for(j in i:n){
      s <- max(s, sum(a[i:j]))
    }
  }
  s
}

計算速度を測る関数を定義する。

test.mss <- function(fun, scale=1){
  n <- seq(100, 500, by=100) * 10^(scale-1)
  times <- numeric(length(n))
  names(times) <- as.character(n)
  for(i in 1:length(n)){
    s <- runif(n[i])
    times[i] <- system.time(fun(s))[[1]]
  }
  return(times)
}

測ってみる。

> test.mss(mss1, 1)
  100   200   300   400   500 
0.028 0.135 0.324 0.626 1.065

$O(n^2)$ の計算量で、長さ500の数列から最大部分列和を求めるだけで1秒程度かかる。

端っこから順次計算していくやり方

全ての部分列を求める必要はなくて、もうすこし上手いやり方がある。
a[1]、a[1]+a[2]、a[1]+a[2]+a[3]...というように頭から順に部分列を伸ばしながら和を計算していく。
和はsという変数に保持しておいて、sの最大値をmssという変数に保持しておく。もし途中でsが0を下回った場合、それまでの部分列を延長しても最大部分列和は得られないので、sを0にリセットして新たに部分列を計算していく。

mss2 <- function(a){
  s <- mss <- 0
  for(i in 1:length(a)){
    s <- s + a[i]
    if(s < 0) s <- 0
    if(mss < s) mss <- s
  }
  mss
}

> test.mss(mss2, 4)
1e+05 2e+05 3e+05 4e+05 5e+05 
0.293 0.582 0.872 1.157 1.448

このやり方だと $O(n)$ の計算量で済む。mss1よりかなり早い。
それでも長さ500000の配列の最大部分列和を求めるのに1秒強かかる。

そもそもRがいけない

Rcppを使えばRの遅さに煩わされることはない。
さらにinlineパッケージを使えば手軽にC++の恩恵にあずかれるので、使わない手はない。

library(Rcpp)
library(inline)
mss3 <- cxxfunction(signature(a = "integer"), "
IntegerVector aa(a);
int s = 0;
int mss = 0;
int n = aa.size();
for(int i=0; i<n; i++){
  s += aa[i];
  if(s < 0) s = 0;
  if(mss < s) mss = s;
}
return wrap(mss);
", plugin = "Rcpp")

> test.mss(mss3, 6)
1e+07 2e+07 3e+07 4e+07 5e+07 
0.079 0.163 0.252 0.342 0.472

5千万の配列から最大部分列和求めるのに0.5秒程度で、mss2と比較してざっと300倍くらい早い。

結論

Rcpp使おう。