---
title: Sized Linear Algebra Package のチュートリアル
tags: 
author: [Akinori Abe](https://www.docswell.com/user/akabe)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/2EVL53VMJQ.jpg?width=480
description: Sized Linear Algebra Package のチュートリアル by Akinori Abe
published: July 09, 16
canonical: https://www.docswell.com/s/akabe/KYWQJ9-2016-07-09
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/2EVL53VMJQ.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/57GX63LXEL.jpg)

自己紹介
●
●
●
名前：阿部晃典
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/


# Page. 3

![Page Image](https://bcdn.docswell.com/page/4EQ845Y5JP.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/KJ42Q9WVE1.jpg)

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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/LE1MPKY47G.jpg)

ことはじめ
●
お好みの 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


# Page. 6

![Page Image](https://bcdn.docswell.com/page/GEW32DGZE2.jpg)

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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/47ZVDGLLE3.jpg)

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


# Page. 8

![Page Image](https://bcdn.docswell.com/page/YJ651YWM7V.jpg)

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


# Page. 9

![Page Image](https://bcdn.docswell.com/page/GJ5Y4GMQJ4.jpg)

違う次元のベクトルの内積を取ると型エラー
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, &#039;a) vec = (four s, float, rprec, &#039;a) Slap_vec.t
but an expression was expected of type
(four, &#039;b) vec = (four, float, rprec, &#039;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 次元 だからエラー」
と言っている


# Page. 10

![Page Image](https://bcdn.docswell.com/page/LE35GRWZJ5.jpg)

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


# Page. 11

![Page Image](https://bcdn.docswell.com/page/8EDD1YK4EG.jpg)

行列
utop # let a = [%mat [1.0, 2.0, 3.0; 4.0, 5.0, 6.0]];;
val a : (two, three, &#039;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, &#039;a) mat =
C1 C2
R1 14 32
= A AT
行列 ✕ 行列演算
R2 32 77
左項の行列は非転置
右項の行列は転置


# Page. 12

![Page Image](https://bcdn.docswell.com/page/V7PXQ6KVE8.jpg)

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


# Page. 13

![Page Image](https://bcdn.docswell.com/page/2JVL53VR7Q.jpg)

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


# Page. 14

![Page Image](https://bcdn.docswell.com/page/5EGX63L6JL.jpg)

ファイルからベクトルを読み込む
(* ファイルからリスト or 配列を読む関数は自作してください *)
val load_list : string -&gt; float list
utop # module X = Vec.Of_list(struct
let value = load_list “hoge.txt”
end);;
module X : sig type n val value : (n, &#039;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


# Page. 15

![Page Image](https://bcdn.docswell.com/page/4JQ845Y27P.jpg)

違うファイルから読み込んだベクトルは違う次元
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, &#039;a) vec = (Y.n, float, rprec, &#039;a) Slap_vec.t
but an expression was expected of type
(X.n, &#039;b) vec = (X.n, float, rprec, &#039;b) Slap_vec.t
Type Y.n is not compatible with type X.n
※ 仮にファイル名が同じでも、読込中に更新されるかもしれないので、型エラーになる


# Page. 16

![Page Image](https://bcdn.docswell.com/page/K742Q9WPJ1.jpg)

LAPACK 関数の紹介


# Page. 17

![Page Image](https://bcdn.docswell.com/page/LJ1MPKYXEG.jpg)

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, &#039;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, &#039;a) mat =
C1 C2
C3
これが逆行列
R1
2.66667 -3
1.33333
R2 -0.666667 1 -0.333333
R3
5.66667 -7
3.33333


# Page. 18

![Page Image](https://bcdn.docswell.com/page/GJW32DGKJ2.jpg)

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, &#039;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, &#039;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


# Page. 19

![Page Image](https://bcdn.docswell.com/page/4EZVDGLNJ3.jpg)

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


# Page. 20

![Page Image](https://bcdn.docswell.com/page/Y7651YW9EV.jpg)

もう少し実用的な例


# Page. 21

![Page Image](https://bcdn.docswell.com/page/G75Y4GMD74.jpg)

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 : (&#039;a, &#039;b, &#039;c) mat → (&#039;d, &#039;a, &#039;e) mat → (&#039;a, &#039;f) vec →
(&#039;d, &#039;g) vec → (&#039;b, &#039;h) vec → (&#039;d, &#039;i) vec = &lt;fun&gt;
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;;


# Page. 22

![Page Image](https://bcdn.docswell.com/page/9J2GVYRMJR.jpg)

準ニュートン法 (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 → (&#039;a, &#039;a, &#039;b) mat → (&#039;a, &#039;c) vec → (&#039;a, &#039;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 &lt;= loops then begin
let s = symv ~alpha:(-1.) h df_dx in (* search direction *)
scal (wolfe_search df f s x) s;
let x&#039; = Vec.add x s in
let df_dx&#039; = df x&#039; in
let y = Vec.sub df_dx&#039; df_dx in
if dot y s &gt; 1e-9 then begin update_h h y s; aux (i + 1) x&#039; df_dx&#039; end
end
in
aux 1 x0 (df x0);;


# Page. 23

![Page Image](https://bcdn.docswell.com/page/DEYV3GDPEM.jpg)

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


