fastmathを用いたClojureでのガウス過程の実装

機械学習技術人工知能技術デジタルトランスフォーメーション技術確率的生成モデル本ブログのナビ深層学習技術ノンパラメトリックベイズとガウス過程サポートベクトルマシン Clojureと関数プログラミング

ガウス過程とは

ガウス過程は関数形をランダムに出力する箱(確率過程)のようなものであり、例えばサイコロが1,2,3,4,5,6の自然数を生成する過程がサイコロのゆがみに依存するのと考えた場合、ガウス過程もパラメータ(この例ではサイコロの歪み具合)次第で関数の出現の様子(サイコロの目が出る確率を表した関数)が変化すると考える。

ガウス過程回帰の例として、入力xと出力yの間の関数関係(赤点線)と、その一部に関する観測データ(青丸印)が得られているとき、この関数関係をガウス過程を使うと以下のように表され、事前分布ではランダムな雲のような分布が実データが与えられると収束していく。

このようなガウス過程(関数の雲)を用いることで、(1)関数の雲が得られると、わからなさの程度がわかる、(2)関数の雲が得られると、自信のある領域とない領域の違いがわかる、(3)関数の雲が得られると、モデル選択や特徴選択ができる等の他の機械学習にない利点が得られる。

ガウス過程の実装

ガウス過程回帰は、データ間の相関係数を用いて解析することから、カーネル法を用いたアルゴリズムが用いられたり、ベイズ解析的手法と組み合わせたMCMCを用いたアルゴリズム等が適用される。

これらの解析に用いられるツールとして、MatlabやPython、R、Clojure等様々な言語でのオープンソースがある。今回はClojureでのアプローチについて述べる。

これは、Clojureの高速数学計算ライブラリーであるfastmathをベースに構築されている。Fastmathは、MKLおよび/またはOpenBlas経由でBLAS/LAPACKに依存するSMILE 2.6.0に依存している。

MKLはMath Kernel Libraryの略で、Intelが2003年に公開したもので、科学・高額・金融アプリケーションに提供される最適化(高速化)された数学ルーチンを含むライブラリで、提供される中心的な数学関数にはBLAS、LAPACK、ScaLAPACK、スパースソルバー(疎行列)、FFT、ディープラーニング、ベクトル演算が含まれており、CPUはインテルプロセッサーのみをサポートしているものとなる。

OpenBlasとは、行列演算ライブラリで、BLASのオープンソース実装となる。BLASとは、Basic Linear Algebra Subprogramsの略で、数値線形代数の基礎的演算に必要な関数を定義するAPIで、ベクトル・行列演算を含む38の関数からなるLevel1 BLASが1979年に発表された後、Level2及びLevel3まで拡張された、この分野におけるデファクトスタンダートとなっており、深層学習でお馴染みのtensorflowやpythonの総額的ライブラリであるnumpy,scikit-learn,pandas,scipyなどにも用いられているものとなる。

SMILEはStatistical Machine Intelligence and Learning Engineの略でJava系(Java、Scala、Koltlin、Clojure、Grrovy)の言語で利用することができる教師あり学習、教師なし学習、NLP(自然言語処理)、各種数学関数等を持った高速の機械学習ライブラリーとなる。

ガウス過程のClojureによる実装

以下にfastmathを用いたガウス過程の実装について述べる。

Clojureの開発環境立ち上げに関しては以下のページを参照のこと、まず用いるライブラリとしてはfastmathの他に、グラフ出力ライブラリとしてcljplotとclojure2dを用いる。よってproject.cljの依存関係は以下のように記述される。

(defproject GP-test01 "0.1.0-SNAPSHOT"
:description "FIXME: write description"
:url "http://example.com/FIXME"
:license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0"
:url "https://www.eclipse.org/legal/epl-2.0/"}
:dependencies [[org.clojure/clojure "1.10.1"]
[generateme/fastmath "1.4.0-SNAPSHOT"]
[clojure2d "1.2.0-SNAPSHOT"]
[cljplot "0.0.2-SNAPSHOT"]]
:repl-options {:init-ns fastmath-test01.core})

それぞれのライブラリは最新のものでは異なるが、最新バージョンにするとコンフリクトを起こすため、上記のバージョンのままでの利用を推奨する。

次にsrcファイルファイルパス名はGP-test01/core.cljとした時、以下のようになる。

(ns GP-test01.core)

(require '[fastmath.core :as m]
         '[fastmath.kernel :as k]
         '[fastmath.random :as r]
         '[fastmath.regression :as reg]
         '[clojure2d.color :as c]
         '[cljplot.build :as b]
         '[cljplot.core :refer :all])

まずはSMILEバージョンでの回帰と点補間について見てみる。このバージョンでは「kernel」(カーネルオブジェクト:デフォルトではガウス分布)と「lambda」(ノイズ係数:デフォルトでは0.5)の2つのパラメータを受け取る。予測するには、predict/predict-all関数を呼び出すか、回帰オブジェクトそのものを呼び出す。

(def domain [-3.0 3.0]) ;; 1次元領域で解析
(def xs [0 1 -2 -2.001]) ;; 入力データは4点
(def ys [-2 3 0.5 -0.6]) ;; 出力データは4点
(def scatter-data (map vector xs ys)) ;; 散布図で表示

(def gp1 (reg/gaussian-process xs ys)) ;;デフォルト値でのガウス過程回帰
(def gp2 (reg/gaussian-process {:kernel (k/kernel :gaussian 0.8)
                                :lambda 0.000001} xs ys)) ;パラメータ変更後のガウス過程回帰

(reg/predict gp1 1.5);1.8472037067061506 ;;predict関数呼び出しでのデフォルトパラメータで入力点1.5で1.8472037067061506とガウス過程回帰予測
(gp1 1.5) ;1.8472037067061506 ;;回帰オブジェクトそのもの呼び出し

(gp2 1.5) ;13.289038788347785 ;;パラメータを変化させたときの回帰値
(reg/predict-all gp2 [0 1 -2 -2.001]) ; List(4)(-1.9999423944882437, 2.9999719139365197, 0.18904211011249572, -0.28904456237796694)  ;;全データでの回帰値

これらの値をプロットすると以下のようになる。

(defn draw-gp
  "Draw multi function chart with legend.
  
  * xs - domain points
  * ys - values
  * title - chart title
  * legend-name - legend title
  * labels - legend labels
  * gps - map of gaussian processes to visualize
  * h, w - optional height and width of the chart"
  [xs ys title legend-name labels gps & [h w]]
  (let [pal (c/palette-presets :category10) ;; use category10 palette
        leg (map #(vector :line %1 {:color %2}) labels pal) ;; construct legend
        scatter-data (map vector xs ys)] ;; prepare scatter data
    (-> (xy-chart {:width (or w 700) :height (or h 300)}
                  (-> (b/series [:grid]
                                [:scatter scatter-data {:size 8}])
                      (b/add-multi :function gps {:samples 400
                                                  :domain domain} {:color pal}))
                  (b/add-axes :bottom)
                  (b/add-axes :left)
                  (b/add-legend legend-name leg)
                  (b/add-label :top title)))))  ;;描画関数


(let [gp (reg/gaussian-process {:lambda 0.5} xs ys)] ;; create gaussian process
  (save (draw-gp xs ys "Default GP regression with Smile"
         "Lambda" [0.5] {:r1 gp})
      "results/default-smile.jpg"))

回帰曲線は全くデータ点を通過していない。これはラムダパラメータ(ノイズ係数)を0.5と大きめに取ったしたためとなる。

ラムダ(ノイズ)パラメータ

ラムダパラメータはデータのノイズレベルを測定し、回帰を調整するものとなる。ここで、ラムダパラメータをいくつか振ってみる。

(let [lambdas [0.00005 0.1 2.0]]
  (save (draw-gp xs ys "GP regression with different lambda"
                 "Lambda" lambdas 
                 (map vector lambdas 
                         (map #(reg/gaussian-process {:lambda %} xs ys) lambdas)))
      "results/smile-different-lambda.jpg"))

結果は以下のようになる。

このように、ラムダパラメータが小さいと関数補間はうまくいく(入力点に交わる)が、モデルのオーバーフィッティングにつながる。

ここで以下のノイズを加えたテスト用のデータを用いて

(defn toyfn [x] (+ (r/grand) (m/exp x) (* -5 (m/cos (* 3 x)))))

(save (xy-chart {:width 700 :height 300}
                (b/series [:grid] [:function toyfn {:domain domain}])
                (b/add-axes :bottom)
                (b/add-axes :left)) "results/toyfn.jpg")

ラムダパラメータの影響を評価する。

(let [lambdas [0.00005 0.1 2.0]
      xs (repeatedly 30 #(apply r/drand domain))
      ys (map toyfn xs)]
  (save (draw-gp xs ys "GP regression with different lambda"
                 "Lambda" lambdas 
                 (map vector lambdas 
                         (map #(reg/gaussian-process {:lambda %} xs ys) lambdas)))
      "results/toy-different-lambda.jpg"))

低いラムダパラメータの方が良いフィッティングが得られる。

カーネルパラメータ

次にもう一つの重要なパラメータであるカーネル関数について述べる。カーネルは相関を表す尺度となり、通常は、カーネルは近い点ほど高い相関を持つように作られる。

fastmath.kernel名前空間では多くのカーネルが定義されている。そのほとんどはGPで使うことができるが、そうでないものもある。GPで使うにはカーネルは対称で半正定値である必要がある。以下にfastmathが持つkernel関数を並べてみる。

(sort (keys (methods k/kernel))) ; :anova, :bessel, :cauchy, :chi-square-cpd, :chi-square-pd, :circular, :dirichlet, :exponential, :gaussian, :generalized-histogram, :generalized-t-student, :hellinger, :histogram, :hyperbolic-secant, :hyperbolic-tangent, :inverse-multiquadratic, :laplacian, :linear, :log, :mattern-12, 14 more…)

カーネルmultimethodは、ベクトルに対して操作する実際のカーネル関数を作成する。multimethodはキーワードとしてカーネル名を受け取り、さらに追加のパラメータ(主にスケーリングパラメータなど)を受け取る。

(def k1 (k/kernel :gaussian)) ;; デフォルトのスケーリングは1.0
(def k2 (k/kernel :mattern-32)) ;; mattern 3/2 kernel
(def k3 (k/kernel :mattern-32 0.5)) ;; mattern 3/2 kernelのスケーリングパラメータ設定

{:k1 (k1 [0] [1])
 :k2 (k2 [0] [1])
 :k3 (k3 [0] [1])}  ;{:k1: 0.6065306597126334, :k2: 0.4833577245965077, :k3: 0.1397313501923147}

ここでそれぞれのカーネル関数(gaussian, cauchy, anova, linear, inverse-multiquadratic, mattern-12, mattern-52, periodic,exponential)がどのように見えるかを示す。ブルーはスケーリング1.0で赤がスケーリング0.5となる。

(def kernels [:gaussian :cauchy :anova :linear :inverse-multiquadratic
              :mattern-12 :mattern-52 :periodic :exponential])

(let [fk (fn [f] #(f [0.1] [%]))
      make-data #(map (fn [k] [(str k)
                        (fk (k/kernel k %))]) kernels)
      cfg {:domain [-3 3] :samples 200 :stroke {:size 2}}]
  (-> (xy-chart {:width 700 :height 600}
                (-> (b/lattice :function (make-data 1) cfg {:grid true})
                    (b/add-series (b/lattice :function (make-data 0.5) 
                                             (assoc cfg :color (c/color 245 80 60)) {:label name})))
                (b/update-scales :x :ticks 5)
                (b/add-label :top "Various kernels around 0.1")
                (b/add-axes :bottom)
                (b/add-axes :left)
                (b/add-axes :right))
      (save "results/kernels.jpg")))

さらにそれぞれの共分散行列について示す。

(let [fk (fn [f] #(f [0] [%]))
      r (range -1 1 0.025)
      data (map (fn [k] [(str k)
                        (let [kernel (k/kernel k)]
                          (for [x r y r] [[x y] (kernel [x] [y])]))]) kernels)]
  (-> (xy-chart {:width 700 :height 600}
                (b/lattice :heatmap data {} {:label name :grid true})
                (b/add-label :top "Covariance matrices"))
      (save "results/covariances.jpg")))

それらのうちいくつかを選んで回帰を行った結果を以下に示す。

(save (draw-gp xs ys "GP regression with different kernel"
                 "Kernel" kernels 
                 (map vector kernels 
                         (map #(reg/gaussian-process {:lambda 0.01
                                                      :kernel (k/kernel %)} xs ys) kernels))) 
      "results/smile-different-kernels.jpg")

ここでいくつかのパラメータの例を示す。

(defn draw-gp-lattice
  [title gps]
  (xy-chart {:width 700 :height 500}
                (-> (b/lattice :function gps {:domain domain :samples 400} 
                               {:shape [-1 1] :grid true})
                    (b/add-series (b/lattice :scatter (zipmap (map first gps) (repeat scatter-data)) {:size 8} {:label name :shape [-1 1]})))
                (b/add-axes :bottom)
                (b/add-axes :left)
                (b/add-label :top title)))

(defn gen-gps
  "Generate gaussian processes for given parametrization."
  [lambdas kernels kernel-params]
  (map (fn [l k p]
         (let [n (str "Kernel: " k "(" p "), lambda=" l)
               kernel (k/kernel k p)]
           [n (reg/gaussian-process {:kernel kernel :lambda l} xs ys)])) (cycle lambdas) (cycle kernels) kernel-params))

まずガウスカーネルについて

(save (draw-gp-lattice "Various gaussian kernels"
                       (gen-gps [0.001] [:gaussian] [0.1 1 4])) "results/gk.jpg")

次にmatternカーネルについて

(save (draw-gp-lattice "Various mattern kernels"
                       (gen-gps [0.0001] [:mattern-12] [0.1 1 4])) "results/mk.jpg")

最後に周期的カーネルについて

(save (draw-gp-lattice "Various periodic kernels"
                       (gen-gps [0.0001] [:periodic] [0.1 1 4])) "results/pk.jpg")

gaussian-process + (GP+)

最後にSMILEバージョンのガウス過程を拡張させたgaussian-process+について述べる。これは以下のような追加情報を得ることができる。(1)事前分布、(2)事後分布、(3)標準偏差。前回と同様にガウス過程のパラメータであるカーネルとノイズ(ラムダ)パラメータに対してこれらがどのようになるか見てみる。

GP+では、「kernel」パラメータはSMILE版と同様で、その他のパラメータは「kscale」(分散巣スケーラ;カーネルの結果にこの値を掛ける(デフォルト1.0))、「noise」(ラムダと同じ、パラメータ名が異なるだけ(デフォルトは1.0e-6)、「normalize?」(ysの平均が0になるように正規化するかどうか(デフォルトはfalse)となる。

(def gp1+ (reg/gaussian-process+ xs ys))
(def gp2+ (reg/gaussian-process+ {:normalize? true} xs ys))

{:gp1+ (gp1+ 1.5)
 :gp2+ (gp2+ 1.5)}

事前分布

事前標本は何の知識(データがない状態)も無い状態から作られた関数の雲となる。

(let [gp (reg/gaussian-process+ xs ys)]
  (reg/prior-samples gp (range 0 1 0.1)))


(defn draw-prior
  ([gp] (draw-prior gp 20))
  ([gp cnt]
   (let [xs (range -3 3.1 0.1)
         priors (map #(vector % (map vector xs (reg/prior-samples gp xs))) (range cnt))]
     (xy-chart {:width 700 :height 300}
               (-> (b/series [:grid])
                   (b/add-multi :line priors {} {:color (cycle (c/palette-presets :category20c))})
                   (b/add-serie [:hline 0 {:color :black}]))
               (b/add-axes :bottom)
               (b/add-axes :left)
               (b/add-label :top "Priors")))))


(save (draw-prior (reg/gaussian-process+ xs ys)) "results/prior.jpg")

ノイズレベルが高くなるとすべての事前分布はよりランダムになる。

(save (draw-prior (reg/gaussian-process+ {:noise 0.5} xs ys)) "results/priorn.jpg")

異なるカーネルを使用した場合、カーネル形状を反映した関数構造になる

(save (draw-prior (reg/gaussian-process+ {:kernel (k/kernel :periodic 2)
                                          :noise 0.01} xs ys) 5)
      "results/priorp.jpg")

(save (draw-prior (reg/gaussian-process+ {:kernel (k/kernel :linear)} xs ys))
      "results/priorl.jpg")

事後分布

いくつかの点が分かっている場合、事後標本を使って事後関数を描くことができる。事後関数は、与えられたカーネルとノイズの値に対して、可能な関数の選択を表す。

(def xs [-4 -1 2])
(def ys [-5 -1 2])

(defn draw-posterior
  ([gp] (draw-posterior gp 20))
  ([gp cnt]
   (let [xxs (range -5 5.1 0.2)
         posteriors (map #(vector % (map vector xxs (reg/posterior-samples gp xxs))) (range cnt))]
     (xy-chart {:width 700 :height 300}
               (-> (b/series [:grid])
                   (b/add-multi :line posteriors {} {:color (cycle (c/palette-presets :category20c))})
                   (b/add-serie [:function gp {:domain [-5 5] :color :black :stroke {:size 2 :dash [5 2]}}])
                   (b/add-serie [:scatter (map vector xs ys) {:size 8}]))
               (b/add-axes :bottom)
               (b/add-axes :left)
               (b/add-label :top "Posteriors")))))


(save (draw-posterior (reg/gaussian-process+ {:kernel (k/kernel :gaussian 0.5)}
                                             xs ys))
      "results/posterior.jpg")

ガウス分布以外のカーネル関数を用いた結果は以下のようになる。

(save (draw-posterior (reg/gaussian-process+ {:kernel (k/kernel :periodic 0.2 6.5)}
                                             xs ys))
      "results/posteriorp.jpg")

標準偏差

事後標本を描く代わりに、標準偏差を使って信頼区間を描くことができる。これは、不確実性の領域に関する情報を与えることができる。この情報はベイズ最適化で重要になります。

下の図は、信頼度50％と95％、そして予測される関数（平均値）を示している。

(defn draw-stddev
  [gp]
  (let [xxs (range -5 5.1 0.2)
        pairs (reg/predict-all gp xxs true)
        mu (map first pairs)
        stddev (map second pairs)
        s95 (map (partial * 1.96) stddev)
        s50 (map (partial * 0.67) stddev)]
    (xy-chart {:width 700 :height 300}
              (b/series [:grid]
                        [:ci [(map vector xxs (map - mu s95)) (map vector xxs (map + mu s95))] {:color (c/color :lightblue 180)}]
                        [:ci [(map vector xxs (map - mu s50)) (map vector xxs (map + mu s50))] {:color (c/color :lightblue)}]
                        [:line (map vector xxs mu) {:color :black :stroke {:size 2 :dash [5 2]}}]
                        [:scatter (map vector xs ys) {:size 8}])
              (b/add-axes :bottom)
              (b/add-axes :left)
              (b/add-label :top "Confidence intervals"))))


(save (draw-stddev(reg/gaussian-process+ xs ys)) "results/predict.jpg")

既知点のない領域では、関数がゼロ値に向かって下がっていることがわかる。これは、正規化されていないことが原因となる。以下に正規化を行ったケースを示す。

(save (draw-stddev(reg/gaussian-process+ {:normalize? true} xs ys)) "results/predictn.jpg")

データの正規化により、線はより平坦になるが、これはより高い不確実性を生み出す。異なるカーネルとオプションを試してみる。

(save (draw-stddev(reg/gaussian-process+ {:kernel (k/kernel :mattern-32)
                                          :noise 0.2
                                          :normalize? true}
                   xs ys)) "results/predictnm.jpg")

(save (draw-stddev(reg/gaussian-process+ {:kernel (k/kernel :gaussian 0.2)
                                          :normalize? false}
                   xs ys)) "results/predictb.jpg")

(save (draw-stddev(reg/gaussian-process+ {:kscale 0.1}
                   xs ys)) "results/predict-low-scale.jpg")

次回はPythonを用いたガウス過程のフレームワークについて述べる。