Sized Linear Algebra Package のチュートリアル

Sized Linear Algebra Package のチュートリアル 2016/7/9 ML勉強会 2016/7/7 A.J.A.勉強会

自己紹介 ● ● ● 名前：阿部晃典 Twitter：@ackey_65535 興味：プログラミング言語理論、機械学習 ○ ○ ○ 幽霊型を使った線形代数ライブラリ Sized Linear Algebra Package (SLAP) ■ http://akabe.github.io/slap/ 幽霊型トリック集 @ML Advent Calendar 2015 ■ http://akabe.github.io/2015/12/PhantomTypeTricks EvilML (a compiler from ML to C++ template) ■ http://akabe.github.io/evilml/

• http://akabe.github.io/slap/
• http://akabe.github.io/2015/12/PhantomTypeTricks
• http://akabe.github.io/evilml/

SLAP (Sized Linear Algebra Package) ● ● ● ● 行列演算の次元の整合性を幽霊型で保証 ○ 例：2次元と3次元のベクトルを足すと型エラー ○ ほとんどの行列演算を検査可能（行列積、固有値分解、最小二乗法など） ○ ファイルから読み込んだ行列などもサイズ検査可能 割と実用的 (cf. Abe & Sumii, EPTCS’15 post-proccedings of ML/OCaml) ドキュメントも（研究で作ったライブラリにしては）多い ○ チュートリアルやデモもあるよ！ 今日は、SLAP のチュートリアルと宣伝が目的 ○ 幽霊型？ブログ読んでね

おしながき ● ● ● ● ● 環境構築 簡単なベクトル・行列演算の紹介 ファイルから読み込んだベクトルのサイズ検査 LAPACK 関数の紹介 実用的なプログラムの紹介

ことはじめ ● お好みの BLAS, LAPACK をインストール ○ ● オーソドックスな本家も良し！ATLAS 等で最適なコードを錬成するも良し！ OPAM のインストール（OCaml のパッケージマネージャ） $ wget https://raw.github.com/ocaml/opam/master/shell/opam_installer.sh -O - | sh -s /usr/local/bin $ eval `opam config env` ● SLAP をインストール（ついでに utop もあると良い） ○ utop：OCaml の高級対話環境 $ opam install slap utop

対話環境で SLAP を読み込む プロンプト $ utop utop # utop # utop # utop # utop # utop # #use “topfind”;; #use “slap”;; #use “slap.top”;; #use “slap.ppx”;; open Slap.D;; open Slap.Common;;

ベクトルとか行列を作ってみる

ベクトル utop # open Slap.Size;; utop # open Slap.D;; utop # let x = [%vec [1.0; 2.0; 3.0; 4.0]];; val x : (four, 'a) vec = R1 R2 R3 R4 1 2 3 4 ベクトルの中身 ベクトルのサイズを表す型 utop # let y = Vec.cons 5.0 x;; val x : (four s, 'a) vec = R1 R1 R2 R3 R4 5 1 2 3 4 4 + 1 = 5 次元ベクトル utop # dot x x;; - : float = 30.

違う次元のベクトルの内積を取ると型エラー val x : (four, ‘a) vec val y : (four s, ‘a) vec val dot : (‘n, _) vec → (‘n, _) vec → float utop # dot x y;; Error: This expression has type (four s, 'a) vec = (four s, float, rprec, 'a) Slap_vec.t but an expression was expected of type (four, 'b) vec = (four, float, rprec, 'b) Slap_vec.t Type four s is not compatible with type four = z s s s s Type four = z s s s s is not compatible with type z s s s Type z s is not compatible with type z 「4 次元 ≠ 5 次元 だからエラー」 と言っている

SLAP の静的サイズ検査でできないこと ● 基本的に、サイズの等しさのみ検査 ○ ● 計算結果が等しいか検査するのは苦手 ○ ● 線形代数だと等しさだけでも結構しっかり検査できる 多少はできるけど、加算の交換則、分配則などは扱えない 大小比較は検査しない ○ ○ 検査する方法もあるが、複雑なので諦めている 具体例： ■ 添字アクセスなどの低級操作 ■ LAPACK 関数の作業用メモリサイズの検査 ■ orgqr (QR分解), syevr (上位固有値計算) utop # Vec.get_dyn x 2;; - : float = 2. 動的にサイズ検査する関数は _dyn がついている

行列 utop # let a = [%mat [1.0, 2.0, 3.0; 4.0, 5.0, 6.0]];; val a : (two, three, 'a) mat = C1 C2 C3 R1 1 2 3 行列の中身 2 ✕ 3 行列 R2 4 5 6 utop # open Slap.Common;; utop # gemm ~transa:normal ~transb:trans a a;; - : (two, two, 'a) mat = C1 C2 R1 14 32 = A AT 行列 ✕ 行列演算 R2 32 77 左項の行列は非転置 右項の行列は転置

行列演算でも次元の不整合を静的に検出できる やってみて！

ファイルから読み込んだって大丈夫

ファイルからベクトルを読み込む
(* ファイルからリスト or 配列を読む関数は自作してください *)
val load_list : string -> float list
utop # module X = Vec.Of_list(struct
let value = load_list “hoge.txt”
end);;
module X : sig type n val value : (n, 'cnt) vec end
utop # X.value;;
- : (X.n, ‘a) vec =
R1
R2 R3 …
R98
R99 R100
1.982 2.234 8 … 9.273 -1.293 3.202
次元が（コンパイル時に）
不明だった
※ X.n は生成的な幽霊型 cf. http://akabe.github.io/2015/10/GenerativePhantomTypes

• http://akabe.github.io/2015/10/GenerativePhantomTypes

違うファイルから読み込んだベクトルは違う次元 utop # module X = Vec.Of_list(struct let value = load_list “hoge.txt” end);; utop # module Y = Vec.Of_list(struct let value = load_list “fuga.txt” end);; utop # dot X.value Y.value;; Error: This expression has type (Y.n, 'a) vec = (Y.n, float, rprec, 'a) Slap_vec.t but an expression was expected of type (X.n, 'b) vec = (X.n, float, rprec, 'b) Slap_vec.t Type Y.n is not compatible with type X.n ※ 仮にファイル名が同じでも、読込中に更新されるかもしれないので、型エラーになる

LAPACK 関数の紹介

getri：逆行列 (LU 分解ベース) ● ● getrf (LU 分解) の結果を元に逆行列を計算 求めた逆行列は引数の行列に破壊的に代入 utop # let a = [%mat [3.0, 2.0, -1.0; 1.0, 4.0, 0.0; -3.0, 5.0, 2.0]];; val a : (three, three, 'a) mat = C1 C2 C3 R1 3 2 -1 R2 1 4 0 R3 -3 5 2 utop # getri a;; - : unit = () utop # a;; - : (three,three, 'a) mat = C1 C2 C3 これが逆行列 R1 2.66667 -3 1.33333 R2 -0.666667 1 -0.333333 R3 5.66667 -7 3.33333

syev：対称行列の固有値分解 ● ● 対称行列：A = AT、直交行列：V VT = E、対角行列：対角成分以外ゼロ 固有値分解：A = V D VT (V は直交行列、D は対角行列) utop # let a = [%mat [3.0, 2.0, -1.0; 2.0, 4.0, 0.0; -1.0, 0.0, -2.0]];; val a : (three, three, 'a) mat = C1 C2 C3 R1 3 2 -1 R2 2 4 0 R3 -1 0 -2 utop # syev ~vectors:true a;; - : (three, ‘a) vec = R1 R2 R3 固有値（D の対角成分） -2.21856 1.60627 5.6123 utop # a;; - : (three, three, 'a) mat = C1 C2 C3 直交行列 R1 0.213021 -0.750591 -0.625487 R2 -0.0685114 0.62713 -0.775896 R3 0.974643 0.208135 0.082168

その他、愉快な仲間たち ● ● ● ● ● ● ● LU 分解、コレスキー分解 QR 分解 連立方程式の求解、線形最小二乗法 ランク計算 行列・ベクトルのノルム 固有値分解（一般行列、正定値対称行列など） 特異値分解

もう少し実用的な例

2層ニューラルネット（順伝搬、逆伝搬）
●
●

機械学習のモデルの1つ（重要な部分だけ紹介）
http://qiita.com/akabe/items/b930e94543f2fa81570f

utop # let exec w1 w2 b1 b2 x =
順伝搬（予測）
let y = sigmv (gemv ~trans:normal w1 x ~beta:1.0 ~y:(Vec.copy b1)) in
let z = sigmv (gemv ~trans:normal w2 y ~beta:1.0 ~y:(Vec.copy b2)) in
z;;
val exec : ('a, 'b, 'c) mat → ('d, 'a, 'e) mat → ('a, 'f) vec →
('d, 'g) vec → ('b, 'h) vec → ('d, 'i) vec = <fun>
utop # let train eta w1 w2 b1 b2 x t =
逆伝搬（学習）
let y = sigmv (gemv ~trans:normal w1 x ~beta:1.0 ~y:(Vec.copy b1)) in
let z = sigmv (gemv ~trans:normal w2 y ~beta:1.0 ~y:(Vec.copy b2)) in
let delta2 = Vec.mul (Vec.sub z t) (sigmdv z) in
let delta1 = Vec.mul (gemv ~trans:trans w2 delta2) (sigmdv y) in
ignore (ger ~alpha:(~-. eta) delta2 y w2);
ignore (ger ~alpha:(~-. eta) delta1 x w1);
axpy ~alpha:(~-. eta) b2 ~x:delta2;
axpy ~alpha:(~-. eta) b1 ~x:delta1;;

• http://qiita.com/akabe/items/b930e94543f2fa81570f

準ニュートン法 (Quasi-Newton method)
●
●

最適化問題の数値解法の1つ（重要な部分だけ紹介）
http://akabe.github.io/slap/demo-gradopt.html

utop # let update_h ?up h y s =
let rho = 1. /. dot y s in
let hy = symv ?up h y in
ignore (ger ~alpha:(~-. rho) hy s h); (* h := h - rho * hy * s^T *)
ignore (ger ~alpha:(~-. rho) s hy h); (* h := h - rho * s * hy^T *)
ignore (syr ?up ~alpha:((1. +. dot y hy *. rho) *. rho) s h);;
val update_h : ?up:bool → ('a, 'a, 'b) mat → ('a, 'c) vec → ('a, 'd) vec → unit
utop # let quasi_newton ~loops df f x0 =
let h = Mat.identity (Vec.dim x0) in (* an approximated inverse Hessian *)
let rec aux i x df_dx = if i <= loops then begin
let s = symv ~alpha:(-1.) h df_dx in (* search direction *)
scal (wolfe_search df f s x) s;
let x' = Vec.add x s in
let df_dx' = df x' in
let y = Vec.sub df_dx' df_dx in
if dot y s > 1e-9 then begin update_h h y s; aux (i + 1) x' df_dx' end
end
in
aux 1 x0 (df x0);;

• http://akabe.github.io/slap/demo-gradopt.html

まとめ ● 結論「SLAP 使ってほしいなぁ〜」 今日の内容： ● 幾つかの関数の使い方を説明 ○ ● ● 他の関数は http://akabe.github.io/slap を参照してください ファイルから読み込んだベクトルでも静的サイズ検査できる ちょっと複雑な線形代数のプログラムを紹介した ○ ○ でも、型注釈は書いてない 行列のサイズは自動的に型推論される

• http://akabe.github.io/slap