Breezeで行列式

大規模データに対する本格的な分析を行うためには線形代数は避けて通れません。

そこでScalaで広く使われている線形代数ライブラリBreezeを調べてみることにしました。

• https://github.com/scalanlp/breeze

最終的には大規模データに対して分散演算を行いたいのでSparkのMLlibが有力候補ですが、MLlibもBreezeを使っているようなのでMLlibを使うときにもそのまま役に立ちそうです。

一応電気工学科を出ているので線形代数は習ったはずなのですが、すっかり忘れてしまったので「ゼロから学ぶ線形代数」(以下、ゼロ線)を参考にしています。以下に出てくる用語はこのゼロ線をベースにします。

目的

Breezeの基本的な使い方として以下の項目を調べます。

• ベクトルの表現
• 行列の表現
• 行列式の計算
• 連立一次方程式の解

準備

build.sbtは以下になります。

view plainprint?

1. name := "breeze-sample"  
2.   
3. version := "1.0"  
4.   
5. scalaVersion := "2.11.7"  
6.   
7. scalacOptions ++= Seq("-Xlint", "-deprecation", "-unchecked")  
8.   
9. val scalazVersion = "7.2.0"  
10.   
11. libraryDependencies ++= Seq(  
12.   "org.scalanlp" %% "breeze" % "0.12",  
13.   "org.scalanlp" %% "breeze-natives" % "0.12",  
14.   "org.scalanlp" %% "breeze-viz" % "0.12"  
15. )  
16.   
17. initialCommands in console := "import breeze.linalg._"  

name := "breeze-sample"

version := "1.0"

scalaVersion := "2.11.7"

scalacOptions ++= Seq("-Xlint", "-deprecation", "-unchecked")

val scalazVersion = "7.2.0"

libraryDependencies ++= Seq(
  "org.scalanlp" %% "breeze" % "0.12",
  "org.scalanlp" %% "breeze-natives" % "0.12",
  "org.scalanlp" %% "breeze-viz" % "0.12"
)

initialCommands in console := "import breeze.linalg._"

sbtをconsoleモードで起動する場合は以下になります。

$ sbt console

ベクトルの表現

BreezeではベクトルはDenseVectorオブジェクトで記述します。

ベクトル$\vec a = \begin{pmatrix}1 \cr2 \cr3\end{pmatrix}$をDenseVectorで作成すると以下になります。

scala> val a = DenseVector(1, 2, 3)
a: breeze.linalg.DenseVector[Int] = DenseVector(1, 2, 3)

行列の表現

行列式を求めるために、2つのベクトル$\vec{a} = \begin{pmatrix}3 \cr2\end{pmatrix}$, $\vec{b} = \begin{pmatrix}7 \cr{-}5\end{pmatrix}$を作成します。

scala> val a = DenseVector(3, 2)
a: breeze.linalg.DenseVector[Int] = DenseVector(3, 2)
scala> val b = DenseVector(7, -5)
a: breeze.linalg.DenseVector[Int] = DenseVector(7, -5)

次に$\vec{a}, \vec{b}$から行列を作成します。

scala> val m = DenseMatrix(a, b)
m: breeze.linalg.DenseMatrix[Int] =
3  2
7  -5

ここで気になることがあります。

ゼロ線の記法では$\vec{a} = \begin{pmatrix}3 \cr2\end{pmatrix}$, $\vec{b} = \begin{pmatrix}7 \cr{-}5\end{pmatrix}$から$\begin{pmatrix}3 & 7 \cr2 & {-}5\end{pmatrix}$という行列を作っていますが、DenseMatrixの表示形式では$\begin{pmatrix}3 & 2 \cr7 & {-}5\end{pmatrix}$になっています。

ゼロ線の表記方法と比べると転置した形になっています。これはDenseVectorがカラムベクトルであるものの要素を横に並べた表記になっているので、これに合わせてカラムベクトルを横に並べる表記にしているのではないかと推測します。

ここではこの解釈で確認を先にすすめることにします。

行列式の計算

行列の1次従属/1次独立の判定式は以下になります。(ゼロ線 P.26)

$\vec{a} = \begin{pmatrix}a \cr c\end{pmatrix}$, $\vec{b} = \begin{pmatrix}b \cr d\end{pmatrix}$に対して、

$$\begin{matrix}ad {-} bc = 0 & ⇔ & \vec{a}と\vec{b}は1次従属 \cr ad {-} bc \ne 0 & ⇔ & \vec{a}と\vec{b}は1次独立\end{matrix}$$

この$ad {-} bc$の式を行列式(determinan)と呼び以下のように定義します。(ゼロ線 P.28)

$$\rm{det}\begin{pmatrix}a & b \cr c & d\end{pmatrix} = ad {-} bc$$

Breezeでは行列式の計算を行う関数detを提供しています。2つのベクトルからdet関数で行列式を計算する関数は以下になります。

view plainprint?

1. def 行列式(a: DenseVector[Double], b: DenseVector[Double]): Double =  
2.     det(DenseMatrix(a, b))  

def 行列式(a: DenseVector[Double], b: DenseVector[Double]): Double =
    det(DenseMatrix(a, b))

この関数を使って行列式の各種法則が動くことを以下のプログラムで確認しました。法則名はゼロ線によります。(ゼロ線 P.30)

view plainprint?

1. val a = DenseVector(3.0, 2.0)  
2.     val b = DenseVector(7.0, -5.0)  
3.     val c = DenseVector(-5.0, 4.0)  
4.     val k = 10.0  
5.     println(s"行列式: det(a,b) = ${行列式(a, b)}")  
6.     println(s"一致退化法則: det(a, a) = ${行列式(a, a)}")  
7.     println(s"交代法則: det(a,b) = ${行列式(a, b)}, det(b,a) =${行列式(b, a)}")  
8.     println(s"分配法則: det(a, b + c) = ${行列式(a, b + c)}")  
9.     println(s"分配法則: det(a, b) + det(a, c) = ${行列式(a, b) + 行列式(a, c)}")  
10.     println(s"分配法則: det(a + b, c) = ${行列式(a + b, c)}")  
11.     println(s"分配法則: det(a, c) + det(b, c) = ${行列式(a, c) + 行列式(b, c)}")  
12.     println(s"スカラー倍法則: det(ka, b) = ${行列式(k * a, b)}")  
13.     println(s"スカラー倍法則: det(a, kb) = ${行列式(a, k * b)}")  
14.     println(s"スカラー倍法則: kdet(a, b) = ${k * 行列式(a, b)}")  

val a = DenseVector(3.0, 2.0)
    val b = DenseVector(7.0, -5.0)
    val c = DenseVector(-5.0, 4.0)
    val k = 10.0
    println(s"行列式: det(a,b) = ${行列式(a, b)}")
    println(s"一致退化法則: det(a, a) = ${行列式(a, a)}")
    println(s"交代法則: det(a,b) = ${行列式(a, b)}, det(b,a) = ${行列式(b, a)}")
    println(s"分配法則: det(a, b + c) = ${行列式(a, b + c)}")
    println(s"分配法則: det(a, b) + det(a, c) = ${行列式(a, b) + 行列式(a, c)}")
    println(s"分配法則: det(a + b, c) = ${行列式(a + b, c)}")
    println(s"分配法則: det(a, c) + det(b, c) = ${行列式(a, c) + 行列式(b, c)}")
    println(s"スカラー倍法則: det(ka, b) = ${行列式(k * a, b)}")
    println(s"スカラー倍法則: det(a, kb) = ${行列式(a, k * b)}")
    println(s"スカラー倍法則: kdet(a, b) = ${k * 行列式(a, b)}")

結果は以下になります。

5 07, 2016 10:44:10 午前 com.github.fommil.jni.JniLoader liberalLoad
情報: successfully loaded /var/folders/s_/qy3xzz_s22g99nh0x8dgrfjm0000gn/T/jniloader5473535094494407732netlib-native_system-osx-x86_64.jnilib
行列式: det(a,b) = -28.999999999999996
一致退化法則: det(a, a) = 0.0
交代法則: det(a,b) = -28.999999999999996, det(b,a) = 28.999999999999996
分配法則: det(a, b + c) = -6.999999999999999
分配法則: det(a, b) + det(a, c) = -6.9999999999999964
分配法則: det(a + b, c) = 25.0
分配法則: det(a, c) + det(b, c) = 25.0
スカラー倍法則: det(ka, b) = -290.00000000000006
スカラー倍法則: det(a, kb) = -289.99999999999994
スカラー倍法則: kdet(a, b) = -289.99999999999994

以下の部分はBreezeがネイティブの数値計算ライブラリを使用していることを情報として出力しているようです。

view plainprint?

1. 5 07, 2016 10:44:10 午前 com.github.fommil.jni.JniLoader liberalLoad  
2. 情報: successfully loaded /var/folders/s_/qy3xzz_s22g99nh0x8dgrfjm0000gn/T/jniloader5473535094494407732netlib-native_system-osx-x86_64.jnilib  

5 07, 2016 10:44:10 午前 com.github.fommil.jni.JniLoader liberalLoad
情報: successfully loaded /var/folders/s_/qy3xzz_s22g99nh0x8dgrfjm0000gn/T/jniloader5473535094494407732netlib-native_system-osx-x86_64.jnilib

行列式

行列式の計算結果は-28.999999999999996となりました。期待していた結果は-29ですから、少し丸め誤差が出るようです。

$ad {-} bc$を文字通り計算すると丸め誤差が出るはずはないので、行列演算の汎用ロジックとして内部で何か難しい計算を行っているのだと推測できます。

Breezeを使う場合(あるいは線形代数ライブラリの一般的な仕様として)、丸め誤差が出ることを前提で考えて、必要な有効桁数で四捨五入していく必要があるということのようです。

この前提を受け入れると、(適当なところで四捨五入すると)行列式の計算結果として-29を得ることができました。

一致退化法則

一致退化法則は0.0を得ることができました。

交代法則

交代法則はベクトルを入れ替えることで行列値の符号が変わることが確認できました。

分配法則

以下の2つの計算が丸め誤差を許容すると同じ値-7になることが確認できました。

$$\begin{matrix}\rm{det}(\vec{a}, \vec{b} + \vec{c}) & ⇒ & \rm{det}\begin{pmatrix}3 & 7 + ({-}5) \cr 2 & {-}5 + 4 \end{pmatrix} \cr\rm{det}(\vec{a}, \vec{b}) + \rm{det}(\vec{a}, \vec{c}) & ⇒ & \rm{det}\begin{pmatrix}3 & 7 \cr 2 & {-}5\end{pmatrix} + \rm{det}\begin{pmatrix}3 & {-}5 \cr 2 & 4\end{pmatrix} \cr\end{matrix}$$

また以下の2つの計算が同じ値25になることが確認できました。

$$\begin{matrix}\rm{det}(\vec{a} + \vec{b}, \vec{c}) & ⇒ & \rm{det}\begin{pmatrix}3 + 7 & {-}5 \cr 2 + ({-}5) & 4 \end{pmatrix} \cr\rm{det}(\vec{a}, \vec{c}) + \rm{det}(\vec{b}, \vec{c}) & ⇒ & \rm{det}\begin{pmatrix}3 & {-}5 \cr 2 & 4\end{pmatrix} + \rm{det}\begin{pmatrix}7 & {-}5 \cr {-}5 & 4\end{pmatrix} \cr\end{matrix}$$

スカラー倍法則

以下の3つ計算とも丸め誤差を許容すると同じ値-290になることが確認できました。

$$\begin{matrix}\rm{det}(k\vec{a}, \vec{b}) & ⇒ & \rm{det}\begin{pmatrix}10 × 3 & 7 \cr 10 × 2 & {-}5 \end{pmatrix} \cr\rm{det}(\vec{a}, k\vec{b}) & ⇒ & \rm{det}\begin{pmatrix}3 & 10 × 7 \cr 2 & 10 × ({-}5)\end{pmatrix} \cr {}k\rm{det}(\vec{a}, \vec{b}) & ⇒ & 10 × \rm{det}\begin{pmatrix}3 & 7 \cr 2 & {-}5\end{pmatrix}\end{matrix}$$

連立一次方程式の解

行列式の応用として連立一次方程式の解を求めてみます。

行列式を使って連立一次方程式を解くための公式としてクラメールの法則があります。(ゼロ線 P.35)

連立方程式$\left\{\begin{array}{ll}ax + by = e\\cx + dy = f\\\end{array}\right.$

を$\vec{a}=\begin{pmatrix}a \cr c\end{pmatrix}$と$\vec{b}=\begin{pmatrix}b \cr d\end{pmatrix}$を用いて、$x\vec{a} + y\vec{b} = \vec{p}$と表すとき、$\rm{det}(\vec{a}, \vec{b}) \ne 0$ならば、解$x$, $y$がただ1つ存在し、

$$x={\rm{det}(\vec{p}, \vec{b}) \over \rm{det}(\vec{a}, \vec{b})}, 　y={\rm{det}(\vec{a}, \vec{p}) \over \rm{det}(\vec{a}, \vec{b})}$$

となる。

クラメールの法則を使って連立一次方程式を解く関数を以下に定義します。

view plainprint?

1. def 連立一次方程式(a: DenseVector[Double], b: DenseVector[Double], p: DenseVector[Double]): (Double, Double) = {  
2.     val x = det(DenseMatrix(p, b)) / det(DenseMatrix(a, b))  
3.     val y = det(DenseMatrix(a, p)) / det(DenseMatrix(a, b))  
4.     (x, y)  
5.   }  

def 連立一次方程式(a: DenseVector[Double], b: DenseVector[Double], p: DenseVector[Double]): (Double, Double) = {
    val x = det(DenseMatrix(p, b)) / det(DenseMatrix(a, b))
    val y = det(DenseMatrix(a, p)) / det(DenseMatrix(a, b))
    (x, y)
  }

この関数を使って以下の2元連立1次方程式を解きます。

$$\left\{\begin{array}{ll}3x + 7y = 1\\2x {-} 5y = 20\\\end{array}\right.$$

この方程式を以下の3つのベクトルで記述し、これに先ほどの関数を適用します。

$$\begin{matrix}\vec{a} & = & \begin{pmatrix}3 \cr 2\end{pmatrix} \cr\vec{b} & = & \begin{pmatrix}7 \cr {-}5\end{pmatrix} \cr\vec{p} & = & \begin{pmatrix}1 \cr 20\end{pmatrix}\end{matrix}$$

view plainprint?

1. val a = DenseVector(3.0, 2.0)  
2.     val b = DenseVector(7.0, -5.0)  
3.     val p = DenseVector(1.0, 20.0)  
4.     val (x, y) = 連立一次方程式(a, b, p)  
5.     println(s"x = $x, y =$y")  

val a = DenseVector(3.0, 2.0)
    val b = DenseVector(7.0, -5.0)
    val p = DenseVector(1.0, 20.0)
    val (x, y) = 連立一次方程式(a, b, p)
    println(s"x = $x, y = $y")

実行結果は以下になりました。

x = 5.000000000000001, y = -2.0000000000000004

丸め誤差をならすと、結果は x = 5, y = 2 となります。

正しい解が得られたようです。

まとめ

線形代数ライブラリBreezeをお試しで使ってみました。

まず入り口ということで行列式と連立一次方程式を動かしてみましたが、法則通り動くことが確認できました。

実際に使ってみた結果の注意点として以下の2つが出てきました。

• 行列の表記方法(ベクトルから作成した時の縦横方向)
• 丸め誤差

行列の表記方法については、ドキュメントなどを調べてみた範囲ではよく分かりませんでした。連立一次方程式の解は正しくでてきたので、今回の解釈でよさそうなのではないかと考えています。継続して、使いながら解釈の補強をしていく予定です。

丸め誤差に関しては、確率計算のような用途ではほぼ問題にならないと推測しますが、今回応用例として用いた連立1次方程式では無視できない問題です。このあたりは使用時に注意する必要がありそうです。

諸元

• Java 1.7.0_75
• Scala 2.11.7
• Breeze 0.12