Blog

2014.12.22

オープンソースライブラリbreezeを導入したときちょっと困ったこと

こんにちは、AMoAdの清水です。AMoAdでは広告配信ネットワークを提供しており、その配信システムを最適化するために広告毎のCTR(click率), CPM(1000回表示あたりの課金額)やimpression数(広告表示回数:imp)など様々な予測を行っています。
予測バッチ他、多くのバッチをAMoAdではScalaで実装しています。
なぜScalaを採用しているのか、それは

元々のシステムのコードはJavaで書かれていたため資産を引き継ぎやすい
Javaよりコードが簡潔に書くことができ、可読性/保守性が高くなる
Scala大好き

という理由です＾＾

今日は予測バッチにScalaで行列計算ができるオープンソースライブラリ「scalanlp/breeze」を導入してみてちょっと困った話を紹介したいと思います！

予測モデルについてざっくりと

まず、一定期間のデータを元に重回帰分析を用いて予測モデルを作りますが
例えばZの値を予測したい場合、Zに相関のあるA, B, C, D,… といった値とベクトルβを掛けて計算します。

Z = impとして予測する場合、A, B, C, D,…にはimpに相関関係のあると思われる指標を入れるというイメージです。
相関しそうなものなんてCPMとかCTRとか色々ありますがそこは工夫次第…。

nは当日をn = 0, 前日をn = 1のように表した値です。
ここまでは特別なことではありませんね。
秘伝のタレβをどれだけ上手に作るかが予測精度に大きく関係してきます。
予測モデルの話をしたいわけではないのでこの辺はさっと飛ばします。~~というかできません笑~~

ではβを求めます。

X:過去のimpやCPMやCTRやイベント変数等の行列
Y:impやCPMのベクトル
これをバッチ化して計算するのは非常に大変です。
まずXの行列のデータを何で持つか、二次元配列？List？
転置と逆行列に相当する関数を自分で実装するのは非常にめんどくさい上にバグの保証はどうするのかと…？

そんなときに便利なのが冒頭でも説明したScalaで行列が書けてしまうオープンソースライブラリ、ScalaNLPの「breeze」です！

基本的な使い方はgithubのドキュメントとこちらのブログで勉強しました。(ブログ主様にリンク許可は頂きました)

ScalaでRやMatlabみたいに行列の操作をおこなう

ではbreezeで実装していきましょう

breezeのバージョンはv0.9を使用

式はもうあるのであとは実装していくだけです。

XとYにDBから取ってきた値を詰めていざ計算してみます＾ω＾

(めっちゃ端折ってますが笑)

val X: DenseMatrix[Double] = ...
val Y: DenseMatrix[Double] = ...
val beta = inv(X.t * X) * X.t * Y

// X.tはXの転置行列
// invは逆行列を求める関数

val X: DenseMatrix[Double] = ...

val Y: DenseMatrix[Double] = ...

val beta = inv(X.t * X) * X.t * Y

// X.tはXの転置行列

// invは逆行列を求める関数

結果：MatrixSingularException

あれ？？

原因を考えると、日によって行列Xに縦に0が並ぶことがある -> 1行全部0は逆行列の定義的に計算できない！
はい。高校のときに習ったやつを思い出します。

|a b|
|c d|

1 2	\|a b\| \|c d\|

2 × 2 の正則行列で逆行列の存在条件は「 ad-bc ≠ 0 」でした。

4 × 4 でも10 × 10でも縦に0が並んだら存在条件を満たさなくなってしまいます。
よってエラるのは納得です＾ω＾

そうゆうわけで疑似逆行列を使いたいところですが幸いなことにbreezeにはムーア・ペンローズの疑似逆行列関数も実装されていました。
pinvってやつです。

では早速…

val beta = pinv(X.t * X) * X.t * Y

1	val beta = pinv(X.t * X) * X.t * Y

結果：MatrixSingularException

あ？(＾ω＾#)

色々試してダメなので簡単にpinvの検証をすることにします。

val m = DenseMatrix((0d,3d,6d), (0d,4d,7d), (0d,5d,9d))
val mi = pinv(m)
println(mi)

val m = DenseMatrix((0d,3d,6d), (0d,4d,7d), (0d,5d,9d))

val mi = pinv(m)

println(mi)

結果は同じく MatrixSingularException

これが既にだめですね。

一旦、本来期待する値をRの疑似逆行列関数ginvで出してみましょう

> m <- rbind(c(0,3,6), c(0,4,7), c(0,5,9))
> ginv(m)
[,1] [,2] [,3]
[1,] 0.000000 0.0000000 0.0000000
[2,] -2.526316 1.4210526 0.5789474
[3,] 1.421053 -0.7368421 -0.2631579

> m <- rbind(c(0,3,6), c(0,4,7), c(0,5,9))

> ginv(m)

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

[1,] 0.000000 0.0000000 0.0000000

[2,] -2.526316 1.4210526 0.5789474

[3,] 1.421053 -0.7368421 -0.2631579

これですよこれ。

なぜRなら計算できるのにbreezeだと計算できないのか？
それぞれの実装を見てみました。

R(3.1.0)の実装 (関数実装のコンソール出力より)

> ginv
function (X, tol = sqrt(.Machine$double.eps)) {
 if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X))) 
   stop("'X' must be a numeric or complex matrix")
 if (!is.matrix(X)) 
   X <- as.matrix(X)
 Xsvd <- svd(X)
 if (is.complex(X)) 
   Xsvd$u <- Conj(Xsvd$u)
 Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
 if (all(Positive)) 
   Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
 else if (!any(Positive)) 
   array(0, dim(X)[2L:1L])
 else 
   Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))
}

> ginv

function (X, tol = sqrt(.Machine$double.eps)) {

if (length(dim(X)) > 2L || !(is.numeric(X) || is.complex(X)))

stop("'X' must be a numeric or complex matrix")

if (!is.matrix(X))

X <- as.matrix(X)

Xsvd <- svd(X)

if (is.complex(X))

Xsvd$u <- Conj(Xsvd$u)

Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)

if (all(Positive))

Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))

else if (!any(Positive))

array(0, dim(X)[2L:1L])

else

Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))

}

breezeの実装 (doc略)
math/src/main/scala/breeze/linalg/functions/inv.scala: 54

object pinv extends UFunc {
  implicit def implFromTransposeAndSolve[T, TransT, MulRes, Result]
  (implicit numericT: T=>NumericOps[T],
   trans: CanTranspose[T, TransT],
   numericTrans: TransT => NumericOps[TransT],
   mul: OpMulMatrix.Impl2[TransT, T, MulRes],
   numericMulRes: MulRes => NumericOps[MulRes],
   solve : OpSolveMatrixBy.Impl2[MulRes, TransT, Result]):Impl[T, Result] = {
    new Impl[T, Result] {
      def apply(X: T): Result = {
        (X.t * X) \ X.t
      }
    }
  }
}

object pinv extends UFunc {

implicit def implFromTransposeAndSolve[T, TransT, MulRes, Result]

(implicit numericT: T=>NumericOps[T],

trans: CanTranspose[T, TransT],

numericTrans: TransT => NumericOps[TransT],

mul: OpMulMatrix.Impl2[TransT, T, MulRes],

numericMulRes: MulRes => NumericOps[MulRes],

solve : OpSolveMatrixBy.Impl2[MulRes, TransT, Result]):Impl[T, Result] = {

new Impl[T, Result] {

def apply(X: T): Result = {

(X.t * X) \ X.t

}

(X.t * X) \ X.t
あれ、なんだこれ？

公式ドキュメントに
「Breeze is a generic, clean and powerful Scala numerical processing library patterned after NumPy, Matlab and R and licensed under Apache Public License 2.0.」
て書いてあったのでてっきりbreezeはRやNumPyの機能を踏襲して作ってあるのかと思い込んでいました！
仮にも一部では有名なOSSでありプロダクションコードで実用されている例もいくつもあったので。。

これは僕のような一般人にはわからない思想の違いがあるのかと思い、breezeのメインコントリビューターのDavidさんに聞いてみたところ、
頂いた回答は本来は特異値分解(SVD)を使って実装しなければならない、修正する必要がある、との事でした。
つまりぶっちゃけpinvはほぼ未実装です。
またpinvだけでなく、breeze自体結構バグや未実装機能が多く存在するので使用するときは注意が必要です。

そうですよね、オープンソースってそうゆうことですよね…
そもそもまだバージョン1.0に達してませんし。

また機会があればbreezeで便利な機能やバグにあたるので注意すべきところを紹介したいと思います。(需要ないか…)

ということで詰みました

このままでは行列計算が使えません！
しかしbreezeをやめるのはもっと時間が掛かります。

よって一番手っ取り早い方法としてRのginvを自分で実装してしまえばいいわけです＾ω＾
ginvの実装はわかっているので簡単ですね！中学生でも書けちゃいますね！！

* 実装してみた結果

object ginv {
  def apply(m: DenseMatrix[Double]) = {
    //特異値分解
    val svd1 = svd(m)
    val u = svd1._1
    val d = svd1._2
    //Rでは特異値分解のvはbreezeのvの転置になっているため転置行列を取る
    val v = svd1._3.t
    //RのDoubleの最小の非負値浮動小数の平方根
    val tol = 1.490116E-08
    val positive = d.map(_ > max(tol * d(0), 0d))
    if (all(positive)) v * multipleVector(d.map(1 / _), u.t)
    else if (!any(positive)) m := 0d
    else dropMatrixColumns(v,positive) * multipleVector(dropVectorColumns(d, positive).map(1 / _), dropMatrixColumns(u, positive).t)
  }

  /**
   * ベクトル[positive]のfalseに対応する行列[dm]の列を削除する 
   * @param dm 対象行列
   * @param positive booleanのベクトル
   * @return DenseMatrix
   */
  def dropMatrixColumns(dm: DenseMatrix[Double], positive: DenseVector[Boolean]): DenseMatrix[Double] = {
    def go(dm2: DenseMatrix[Double], n: Int): DenseMatrix[Double] = {
      if (n <= 0 || positive(n)) dm2
      else go(dm2.delete(n, Axis._1), n - 1)
    }
    go(dm, positive.length - 1)
  }

  /**
   * ベクトル[positive]のfalseに対応するベクトル[dv]の要素を削除する
   * @param dv 対象ベクトル
   * @param positive booleanのベクトル
   * @return DenseVector
   */
  def dropVectorColumns(dv: DenseVector[Double], positive: DenseVector[Boolean]): DenseVector[Double] = {
    val list = dv.toArray
    val lb = ListBuffer(list: _*)
    for (i <- (list.length - 1) to 0 by -1) {
      if (!positive(i)) lb.remove(i)
    }
    DenseVector(lb: _*)
  }

  /**
   * ベクトル * 行列を行う
   * DenseVector(1,2,3) * DenseMatrix((1,2,3),(4,5,6),(7,8,9)) はできないため
   * DenseVector(1,2,3) -> DenseMatrix((1,1,1),(2,2,2),(3,3,3))に変換する必要がある
   * R での V * Mに当たる
   * @param dv ベクトル
   * @param dm 行列
   * @return DenseMatrix
   */
  def multipleVector(dv: DenseVector[Double], dm: DenseMatrix[Double]) = {
    val arr = for {
      d <- dv.toArray
    } yield Stream.fill(dm.cols)(d)
    val dm2 = DenseMatrix(arr.toSeq: _*)
    dm :* dm2
  }
}

object ginv {

def apply(m: DenseMatrix[Double]) = {

//特異値分解

val svd1 = svd(m)

val u = svd1._1

val d = svd1._2

//Rでは特異値分解のvはbreezeのvの転置になっているため転置行列を取る

val v = svd1._3.t

//RのDoubleの最小の非負値浮動小数の平方根

val tol = 1.490116E-08

val positive = d.map(_ > max(tol * d(0), 0d))

if (all(positive)) v * multipleVector(d.map(1 / _), u.t)

else if (!any(positive)) m := 0d

else dropMatrixColumns(v,positive) * multipleVector(dropVectorColumns(d, positive).map(1 / _), dropMatrixColumns(u, positive).t)

}

/**

* ベクトル[positive]のfalseに対応する行列[dm]の列を削除する

* @param dm 対象行列

* @param positive booleanのベクトル

* @return DenseMatrix

def dropMatrixColumns(dm: DenseMatrix[Double], positive: DenseVector[Boolean]): DenseMatrix[Double] = {

def go(dm2: DenseMatrix[Double], n: Int): DenseMatrix[Double] = {

if (n <= 0 || positive(n)) dm2

else go(dm2.delete(n, Axis._1), n - 1)

}

go(dm, positive.length - 1)

}

/**

* ベクトル[positive]のfalseに対応するベクトル[dv]の要素を削除する

* @param dv 対象ベクトル

* @param positive booleanのベクトル

* @return DenseVector

def dropVectorColumns(dv: DenseVector[Double], positive: DenseVector[Boolean]): DenseVector[Double] = {

val list = dv.toArray

val lb = ListBuffer(list: _*)

for (i <- (list.length - 1) to 0 by -1) {

if (!positive(i)) lb.remove(i)

}

DenseVector(lb: _*)

}

/**

* ベクトル * 行列を行う

* DenseVector(1,2,3) * DenseMatrix((1,2,3),(4,5,6),(7,8,9)) はできないため

* DenseVector(1,2,3) -> DenseMatrix((1,1,1),(2,2,2),(3,3,3))に変換する必要がある

* R での V * Mに当たる

* @param dv ベクトル

* @param dm 行列

* @return DenseMatrix

def multipleVector(dv: DenseVector[Double], dm: DenseMatrix[Double]) = {

val arr = for {

d <- dv.toArray

} yield Stream.fill(dm.cols)(d)

val dm2 = DenseMatrix(arr.toSeq: _*)

dm :* dm2

}

こんなことに／(^o^)＼
なんでこんなに内部関数作ってんだよと罵ってください。

Rにはbreezeにない機能があり

|1 2 3|
|4 5 6|
|7 8 9|

|1 2 3|

|4 5 6|

|7 8 9|

という行列と

|true true false|

1	\|true true false\|

というベクトルから

|1 2|
|4 5|
|7 8|

|1 2|

|4 5|

|7 8|

という行列を生成できる関数が標準で備わっていたり、

ベクトル×行列の計算はRではできるけどbreezeではできない等の理由から、補助関数を作る必要がありました。
残念ながら関数multipleVectorの中の実装が良くないため、この自家製ginvは正則行列以外には対応しておりません。
とりあえず今回は正則行列だけ対応できれば問題ありませんでした。
今後改良する機会が来るだろうか…
というかもうちょっと汎用的に実装してbreezeにプルリクエスト送れって話ですね笑

まとめ

ScalaNLPのbreezeの一機能を紹介しました。
breezeは行列計算ライブラリの中でもかなり直感的でわかりやすい仕様になっていて、数学が得意でない人も気軽に使えるライブラリになっています。
何かのデータをバッチで計算するときにbreezeの行列を使うことで計算過程のわかりやすいコードが書け、後から見た人にやさしい保守性の高い実装ができるかもしれません。
しかしまだまだ発展途上のライブラリであるため、今回紹介したpinvのように簡易実装になっている機能やバグがいくつかありますので、各機能を使う前に十分に検証してみたほうがいいかと思います＾ω＾

それでは良いお年を。

Author

shimizu