読者です 読者をやめる 読者になる 読者になる

MATHGRAM

主に数学とプログラミング、時々趣味について。

RustでForward自動微分を実装してみた

rust 数学 最適化

Rustの理解がまだまだすぎて. モジュールとして全く使えない実装になっているので悪しからず. Forward自動微分そのものを理解することを目標とします.

ちなみに僕が参考にしたのは以下のサイトです.

ありがとうございました.

自動微分とは

自動微分そのものの説明は上にあげたサイトがとても詳しく説明しています. 数値微分を理解している人ならば, 簡単に違いを理解できると思います.

自分なりに自動微分を簡単にまとめると,
導関数を先に定義して計算効率&精度を上げようって感じです.
例えば,  \sin{x}微分 \cos{x}である, と先に定義しておくわけですね.

参考にしたサイトと同じようなことを説明しても意味はないと思うので, ここでは例を複数あげることで理解の助けになれればと思います.

とりあえず, Forward自動微分の実装で重要なのは, 以下の3点だと自分は理解しました.

  • 二重数と呼ばれる, Dual型の定義
  • Dual型に対する基本演算の定義
  • Dual型の初期値の設定

1つずつ詳しく見ていくことで, どの言語でも実装できるような説明を心がけてみます.

Dual型の定義

Dual型はある変数xとその微小量を表すdxをもつ型です.
Rustでは以下のように実装しています.

#[derive(Debug, Copy, Clone)]
pub struct Dual {
    var: f32,
    eps: f32,
}

みてわかる通り, 普通の変数に微小量がくっついただけなので, int型やfloat型の拡張と考えていいと思います.

基本演算の定義

では, 次にDual型に対する基本演算を定義します.
というより, 微小量 dxに関する基本演算を定義します.

例として, まずはを定義してみましょう.

こんな式があったとします.

 y = x + \sin{x}

突然ですが, この y xに関して微分して見てください.

・・・はい. 恐らくみなさんの頭の中では第1項, 第2項をそれぞれ微分して最後に和をとる, というような暗算をしただろうと思います.

より詳しくいうのであれば, まず第1項を計算する.

\displaystyle \frac{d}{dx} x = 1

次に第2項も微分する.

\displaystyle \frac{d}{dx} \sin{x} = \cos{x}

そして最後に和をとる.

\displaystyle \frac{d}{dx}y = 1 +  \cos{x}

以上のようなステップを踏むことで, y xに関する微分を導いたと思います.

ここで重要なのは最後に和をとったことです.

つまりyの微小量は, それぞれの項の微小量を単に足すことで求めることができます.

より簡潔に言うならば, 和の微分微分の和ということです.

以上のルールに則りRustでDual型の和を定義するとこのようになります.

// + 演算子のオーバーロード
impl Add for Dual {
    type Output = Dual;
    fn add(self, r: Dual) -> Dual {
        Dual {
            var: self.var + r.var,
            //ある変数の和に対する微小量の和を定義する.
            eps: self.eps + r.eps
        }
    }
}

それでは次にを定義をしてみましょう.
先ほどと同じように具体例を出してみます.

 y = x  \sin{x}

和の時と同様にx微分してみてください.

・・・はい. 今度は高校の時に呪文のように覚えた微分そのまま, そのまま微分のルールに基づいて暗算したのではないでしょうか. (この呪文は僕だけかもしれませんが…)

具体的に書いてみると,

 \displaystyle
\frac{d}{dx}y = \sin{x}\frac{d}{dx}x + x \frac{d}{dx}\sin{x}

こうですね. それではこのルールに則りDual型の積, もといDual型の微小量に対する積を定義しましょう.

// 積のオーバーロード
impl Mul for Dual {
    type Output = Dual;
    fn mul(self, r: Dual) -> Dual {
        Dual {
            var: self.var * r.var,
            // ある変数の積に対する微小量の積を定義する.
            eps: self.eps*r.var + self.var*r.eps
        }
    }
}

以上のように差や商に関してもDual型の演算を定義してあげることで, 勝手に微小量が計算されちゃうよっていう寸法です.

ここで簡単な式に実践してみましょう.

 y = x^2 + x

この式の微分

 \displaystyle
\frac{d}{dx}y = 2x + 1

ですね. 簡単です. この式によると x = 2の点での傾きは5です. さてDual型を使ってこの5は計算できるのか.

// 式の定義
fn example1(x: Dual) -> Dual {
    x*x + x
}

fn main(){
    // x = 2 なので varは 2 とします.
    let x = Dual{var: 2f32, eps: 1f32};

    println!("{:?}", example1(x));
}

出力

Dual { var: 6, eps: 5 }

おーちゃんと計算できてますね!

先ほど説明した四則演算のように,  \sin \expも定義してあげればこれらの演算が入っている式でも問題なく傾きを求めることができます.

こんな感じでざっと定義してあげて…

impl Dual {
    fn sin(self) -> Dual {
        Dual {
            var: self.var.sin(),
            eps: self.eps*self.var.cos()
        }
    }

    fn exp(self) -> Dual {
        Dual {
            var: self.var.exp(),
            eps: self.eps*self.var.exp()
        }
    }
}

こいつらを含んだ式を適当に作ってあげて…

 y = \sin{x} + xe^{x}
fn example2(x: Dual) -> Dual {
    x.sin() + x*x.exp()
}

 x = 0の傾きを求める!

 \displaystyle \frac{d}{dx} y(0) = \cos{0} + e^{0} + 0e^{0} = 2
fn main(){
    // x = 0
    let x = Dual{var: 0f32, eps: 1f32};
    
    println!("{:?}", example2(x));
}

出力じゃ!!

Dual { var: 0, eps: 2 }

おもしれェェェエエエ!

仕組みは単純なのに簡単に傾きが求められる!やばい!

初期値の設定

ここでより一般的な拡張を考えてみると,

おいおい, 多変数関数の時はどうすんだい?

ってなりますよね. 実はForward modeは多変数関数に弱いんです. でもDual型の初期値の設定によって一応求めることができます.

例えばこんな式を用意しましょう.

 \displaystyle z(x, y) = \sin{y} + xy + ye^{x}

そして z xによる偏微分を求めてみる. 手計算するとこうです.

 \displaystyle
\frac{\partial}{\partial x} z(x, y) = y + ye^{x}

 xに関する傾きを求めるには次のように初期値を設定してあげましょう.

let x = Dual{var: 0f32, eps: 1f32};
let y = Dual{var: 2f32, eps: 0f32}; // epsを0にする!

つまり y xに対する微小量は0と, 定義してあげます. 上のように定義すると出力は

Dual { var: 2.9092975, eps: 4 }

こうなります. 実際,

 \displaystyle
\frac{\partial}{\partial x} z(0, 2) = 2 + 2e^{0} = 4

なのであってますね!逆に y偏微分を求めるには,

let x = Dual{var: 0f32, eps: 0f32}; // epsを0にする!
let y = Dual{var: 2f32, eps: 1f32}; 

こうすればOKですね.

1度の演算で1つの変数に対する傾きしか求められないのが残念ですね.

全ての変数に対して1度の演算で勾配を求められるのがBackward modeなのですが, ちょっとまだ理解できてません. と言うか仕組みはわかるけど実装ができない・・・. もうちょっと勉強してみます.

一応全コード載せとく

use std::ops::{Add, Sub, Mul, Div};
use std::f32;

#[derive(Debug, Copy, Clone)]
struct Dual {
    var: f32,
    eps: f32,
}

// + 演算子のオーバーロード
impl Add for Dual {
    type Output = Dual;
    fn add(self, r: Dual) -> Dual {
        Dual {
            var: self.var + r.var,
            //ある変数の和に対する微小量の和を定義する.
            eps: self.eps + r.eps
        }
    }
}

impl Sub for Dual {
    type Output = Dual;
    fn sub(self, r: Dual) -> Dual {
        Dual {
            var: self.var - r.var,
            eps: self.eps - r.eps
        }
    }
}

// 積のオーバーロード
impl Mul for Dual {
    type Output = Dual;
    fn mul(self, r: Dual) -> Dual {
        Dual {
            var: self.var * r.var,
            // ある変数の積に対する微小量の積を定義する.
            eps: self.eps*r.var + self.var*r.eps
        }
    }
}

impl Div for Dual {
    type Output = Dual;
    fn div(self, r: Dual) -> Dual {
        Dual {
            var: self.var / r.var,
            eps: self.eps/r.var - r.eps*self.var/r.var/r.var
        }
    }
}

impl Dual {
    fn sin(self) -> Dual {
        Dual {
            var: self.var.sin(),
            eps: self.eps*self.var.cos()
        }
    }

    fn cos(self) -> Dual {
        Dual {
            var: self.var.cos(),
            eps: -self.eps*self.var.sin()
        }
    }

    fn tan(self) -> Dual {
        Dual {
            var: self.var.tan(),
            eps: self.eps/(self.var.cos()*self.var.cos())
        }
    }

    fn exp(self) -> Dual {
        Dual {
            var: self.var.exp(),
            eps: self.eps*self.var.exp()
        }
    }
}

fn newton_sqrt(var: Dual) -> Dual {
    let mut y = Dual{var:2f32, eps:0f32};
    let two = Dual{var:2f32, eps:0f32};
    for i in 0..10 {
        y = (y + var/y) / two;
        println!("{:?}", y);
    }
    y
}

// 式の定義
fn example1(x: Dual) -> Dual {
    x*x + x
}

fn example2(x: Dual) -> Dual {
    x.sin() + x*x.exp()
}

fn example3(x: Dual, y: Dual) -> Dual {
    y.sin() + x*y + y*x.exp()
}

fn main(){
    // x = 0 の傾きを求めてみる
    let x = Dual{var: 0f32, eps: 1f32};
    let y = Dual{var: 2f32, eps: 0f32};

    //println!("{:?}", example2(x));
    println!("{:?}", example3(x, y));
    //println!("{:?}", newton_sqrt(var));
}

まとめ

1年前に僕の尊敬する師匠から自動微分の話を聞いたのですが, 当時は全く意味がわかりませんでした. そもそもコンピュータがどのように微分しているかすら考えたことなかったですからね.

Backward modeも頑張って実装したいところ・・・.

以上です.

chainerでYahoo!の占有率推定を実装してみた

Deep Learning chainer python 機械学習 画像処理

こんにちは. 今回は僕も参加したYahoo! JAPAN データ&サイエンス カンファレンス 2016で取り上げられていた占有率推定を実装してみました.

前書き

Yahoo! Japanは働き方改革のロールモデルとして, 頻繁にメディアに取り上げられていますね. 特に週休三日制度のインパクトは大きく, 各所で話題になっていることだろうと思います. 業務効率を落とさずにこの制度を現実なものとするためには単純作業をAIに行わせるということが不可欠となって来ます.

そんな中, つい先日こんな記事を見かけました. itpro.nikkeibp.co.jp

この技術は僕がカンファレンスに参加した時に一番興味を持った(画像系を主としている僕にはちゃんとわかるものがこれくらいしかなかった)ものでした. しかしその時はトリミングで使われるなどの説明は恐らくしていなかったように思えます. 僕の記憶が正しければ, 主に検索の最適化で使われると説明していました. 何はともわれ, この記事のおかげでカンファレンスのことやこの論文を実装しようとしていたことなどを思い出したので, 今回実装してみた次第です.

論文の概要

実装した論文はこちらにあります.

論 文 - Yahoo! JAPAN研究所 - ヤフー株式会社

モチベーション

  • 画像中に特定の物体が占める割合(占有率) を高速に求めたい.

既存手法との差異, 本手法の利点

  • 既存手法ではsegmentationを行なった後に占有率を算出していた.
  • 本手法ではsegmentationをせず, 占有率を直接算出できる.
  • 既存のImageNetの学習済みパラメータを利用できる.

アルゴリズム

  • 出力層をクラス数のunitからなる全結合にする.
  • 出力層の活性化関数をsoftmaxにし出力を [0, \,1]とする.
  • 占有率を以下のように定め教師データにすることで回帰問題に落とし込める.
 \displaystyle
t_n( \boldsymbol{x}) =  \frac{\sum_p \delta_{n, \,j_p}}{ \sum_i \sum_p \delta_{i, \,j_p}}

ここで,  n n番目の画像データ,  pピクセル,  iはクラス,  \deltaクロネッカーのデルタ,  j_pピクセル pのラベルである.

実験

環境

  • chainer 1.20.1
  • AWS: Bitfusion Boost Ubuntu 14 Torch 7

方法

  • 以前実装したFCNによるSegmentationと速度比較する. FCNについてはこちら.

www.mathgram.xyz

  • モデルはVGG16を使用する.
  • データセットPASCALVOC 2012のSegmentationClassを用いる.
  • 画像の大きさは224*224に固定. (全結合を用いるので固定せざるを得ない)
  • FCNは学習済みの重みを使用していないので今回も使っていない.
  • 論文にはLossに関する記述がないが回帰問題なのでMSEを用いた.
  • 最適化手法はAdamで学習率は 1.0*10^{-5}.

コードの説明

部分的にコードの説明をしておきます.

  • 占有率について
def calc_occupancy(x, n_class=21):
    v = xp.empty((n_class))
    h, w = x.shape

    for i in range(n_class):
        v[i] = np.sum(x == i)
    v /= h*w
    if v.sum != 1.0:
        diff = 1.0 - v.sum()
        v[0] += diff
    return v

教師データにはピクセルごとにクラスが記録されているので, それぞれのクラスごとに和をとり, 画像の大きさで割るだけです. また境界値は背景と同じクラス0でマスクしました.

しかし教師データによっては合計が1にならなかったりしたので, その差分は全て背景クラスに加算しています. 恐らく境界データの部分がうまくマスクできていないと思われます. 以下の部分がちょうどその演算にあたる部分です.

    if v.sum != 1.0:
        diff = 1.0 - v.sum()
        v[0] += diff
  • VGGモデル

PASCAL VOCのSegmentationClassは背景を含めて21のクラスで構成されていまので出力層のfull connectを21unitにしました.

class VGG(chainer.Chain):
    def __init__(self, n_class=21):
        super(VGG, self).__init__(
            conv1_1=L.Convolution2D(3, 64, 3, stride=1, pad=1),
            conv1_2=L.Convolution2D(64, 64, 3, stride=1, pad=1),

            conv2_1=L.Convolution2D(64, 128, 3, stride=1, pad=1),
            conv2_2=L.Convolution2D(128, 128, 3, stride=1, pad=1),

            conv3_1=L.Convolution2D(128, 256, 3, stride=1, pad=1),
            conv3_2=L.Convolution2D(256, 256, 3, stride=1, pad=1),
            conv3_3=L.Convolution2D(256, 256, 3, stride=1, pad=1),

            conv4_1=L.Convolution2D(256, 512, 3, stride=1, pad=1),
            conv4_2=L.Convolution2D(512, 512, 3, stride=1, pad=1),
            conv4_3=L.Convolution2D(512, 512, 3, stride=1, pad=1),

            conv5_1=L.Convolution2D(512, 512, 3, stride=1, pad=1),
            conv5_2=L.Convolution2D(512, 512, 3, stride=1, pad=1),
            conv5_3=L.Convolution2D(512, 512, 3, stride=1, pad=1),

            fc6=L.Linear(25088, 4096),
            fc7=L.Linear(4096, 4096),
            fc8=L.Linear(4096, n_class)
        )
        self.train = False

以上のようにほとんどVGGと変化はありません.

結果

申し訳ないのですが, サーバー側で色々実験をするのが面倒だった(画像を移動させる等)ので予測の部分はCPUの演算のみになっています. 速度比較という点では問題はないと思うので許してください.

環境は以下です.

また一番注意しなければならない点なのですが, FCNも含めて学習が不十分なため"精度“の点にはあまり注目しないで欲しいです. (論文では精度も本手法の方が良くなったという記述がある.) とにかく速度に注目して比較してみよう.

1枚目

まずはこの画像から.
占有率は以下.

label: occupancy
background: 38.18%
person: 61.82%

邪魔なのでほぼ0%のクラスは出力してません. 足しても100%にならないかもです.

FCNによる推定

label: occupancy
background: 40.68%
person: 59.22%

time: 1.88s

本手法による推定

background: 27.52%
cat: 1.85%
chair: 1.24%
person: 65.3%

time: 1.19s

2枚目

次はこの画像.
占有率は以下.

label: occupancy
background: 73.42%
bus: 26.58%

FCNによる推定

label: occupancy
background: 80.55%
bus: 19.45%

time: 1.81s

本手法による推定

label: occupancy
background: 37.96%
aeroplane: 4.65%
boat: 7.58%
bottle: 1.12%
bus: 24.83%
car: 6.87%
diningtable: 1.32%
motorbike: 4.79%
person: 1.23%
sofa: 2.14%
train: 1.72%
tv/monitor: 1.08%

time: 1.29s

ちょっとノイズが多いですね. でもバスはちゃんと24%です.

3枚目


占有率は以下.

label: occupancy
background: 87.16%
dog: 12.84%

FCNによる推定

label: occupancy
background: 90.6%
bird: 2.05%
dog: 5.58%
person: 1.32%

time: 1.89s

本手法による推定

background: 87.83%
aeroplane: 3.45%
dog: 8.41%

time: 1.30s

ラスト


占有率は以下.

label: occupancy
background: 20.0%
car: 80.0%

ぴったり?まじかw

FCNによる推定

label: occupancy
background: 25.92%
aeroplane: 2.89%
car: 70.61%

time: 1.99s

本手法による推定

label: occupancy
background: 61.32%
aeroplane: 3.54%
bird: 1.95%
boat: 12.49%
bus: 4.01%
car: 1.4%
motorbike: 2.11%
person: 4.12%
train: 2.68%

time: 1.26s

うーん, これは全然ダメだ.

考察

まず間違いなく学習が足りてませんね. もっと厳密に実験したいのですが, 今月のAWSの請求2万くらい行きそうなんでもう無理っすw

速度だけで見ると本手法の圧勝です. 精度に関しても学習が進めばしっかりと占有率を出してくれそうですね.

一番個人的に気になっているのが, トリミングの実用的にはどうしているのかって部分です. このモデルはFCNと違い, 全結合を用いているので任意の大きさに対応できません. つまり固定の大きさの画像しか占有率を計算できないということになります.

今のところ, 僕が運用として思いつくのは, 様々な比に対応したモデルを複数用意するくらいです. ここはYahoo!ニュースなどがどのようなルールでトリミングを行なっているかによりますが, トリミングの比がある程度定まっているのであれば, そこまで大変なことではないと思います. AWSを借りてるような僕にはできませんがw

まぁとりあえず僕的には満足の行く実験&結果でした.

既存モデルのたった1層に工夫をするだけで, 会社に大きな貢献をすることができる.

とっても夢のある話じゃないでしょうか.

まとめ

  • CNNによる占有率推定を実装しました
  • 速度向上は測れたが画像によっては精度が悪い
  • 学習が不十分なのがとても心残り

一応レポジトリあげときますが, 出力の部分とか手元でざっとスクリプト書いたのでここに含まれてないです. まぁほんの参考程度に.

GitHub - k3nt0w/occupancy_estimate

以上です.

chainerに復帰したくてFCN実装した

Deep Learning OpenCV chainer python 機械学習 画像処理

もともと, chainerユーザーだった僕ですが, 5月くらいにKerasを使い出してからchainerにはあまり触れてきませんでした. しかし先日の分散処理の発表なども含め, バージョンアップが早く, 常に速度の向上を図っているchainerにはもう1度触れておきたいと思い, 今回ハンドリングとしてFCNを実装した次第です. (ちょうどFCNを使ってやりたいこともあったしね.)

レポジトリはここに置いときます.

GitHub - k3nt0w/chainer_fcn: Implementation FCN via chainer

あと学習済みパラメータもここに置いときます. chainerってついてる方です. お間違いのないように.

ぱらめーたー - Google ドライブ

FCNは以前Kerasでも実装していて, FCNそのものの説明もこちらに詳しくしているので, もし"FCNってなんじゃい“って方がおりましたら, こちらを参照してください.

www.mathgram.xyz

いやーしかし・・・, chainer使いやすすぎだろ!! 誰だよ,
あれ...?chainerよりkerasの方が書きやすくね? - MATHGRAM
なんて言ったの!死んで詫びろ!
昔使っていた貯金はありますが, 今回の実装には3時間もかかってないと思います. しかもほとんどは画像の出力方法の部分. ただ僕が使っていた時には無かったtrainerの部分は全く触れずに実装しちゃいました. 必要性を感じたら勉強しよう.

keras2も気になりますが, chainerもいっぱい使っていこうと思います.

目次

環境

実験方法

実験にはPASCAL VOC2012のSegmentationClassを使いました. クラス数は背景を含めて21. trainデータは, VOCのtrain&valを全部使っちゃってて約2500枚, epochは100, batchsizeは5です.

境界線には-1を当てて計算から除外しました. 任意のサイズの画像をinputできますが, 速度向上のため256に統一しています. 予測は任意のサイズいけます.

また今回もこの論文からそのまま実装しました. 厳密にはFCN-8sの実装となっています.

実験結果

lossはまだまだ下がりきってないけど, AWSの請求が半端なくなりそうなので, とりあえずここで区切る.

あとtrainデータにのみ予測を行なっています. お叱りを受けそうですが, 今回の目的はFCNの精度を上げることではないので許してください. むしろいいGPUマシン持っている方, fine-tuningしてくれたらくそ喜びます.

人間

f:id:ket-30:20170218120353p:plain

f:id:ket-30:20170218115535p:plain

バス

f:id:ket-30:20170218115050p:plain

クルマ

f:id:ket-30:20170218120620p:plain

Kerasの時もそうだったけど, 犬って難しいのかな? f:id:ket-30:20170218120815p:plain

おまけ

VGG16のpretrain modelを使う

今回はなぜか1から学習させてしまいましたが, FCN−8sはほとんどVGGと同じ構成なので学習済みのパラメータを使うことができます. モデルの定義をVGGとその出力を受け取る部分に分ければいいだけですね.

インデックスカラーについて

前回の実装では, PASCALVOCの画像をRGBに変換し, 色ごとに自分でクラスを分けるというクッソだるいことしていましたが, なんとインデックスカラーなるものがこの世には存在するようです.

インデックスカラー - Wikipedia

VOCのSegmentationClassに存在する画像はこのインデックスカラーによって着色されたもので, pngで画像を読み込めばそのままクラスを表すarrayになってくれるんです.

画像で説明するとこんな感じ. f:id:ket-30:20170218124136p:plain

当たり前な人にとっては当たり前なのかもしれませんが, 多分1人で勉強してるとなかなか分からなかったりすると思います. この情報だけで実装にかかる時間は何倍も短縮できますからね. この情報が誰かの役に立つことを望みます.

アルファブレンドについて

今回の出力は元画像にsegmentationした画像を重ねているのですが, こういうのをアルファブレンドって言ったりするそうです.

アルファブレンド - Wikipedia

こういう画像そのものの学習を自分は怠っていると感じました, pngとかjpgとかちゃんとわかってないしなぁ. もっと基本的なところから勉強して, いずれ記事にまとめたいところ.

PILでアルファブレンドする方法がわからなかったので, そこだけopencv使っています.

o = cv2.imread("out/original.jpg", 1)
p = cv2.imread("out/pred.png", 1)

# 0.6, 0.4 は比です.
pred = cv2.addWeighted(o, 0.6, p, 0.4, 0.0)

参考

今回もたくさんのブログ, レポジトリを参考にさせていただきました.
ありがとうございました.

chainerそのものとFCNの参考
https://github.com/yusuketomoto/chainer-fast-neuralstyle https://github.com/pfnet/PaintsChainer
memo: Fully Convolutional Networks 〜 Chainerによる実装 〜

PILでは透化の合成ができなかったけど, このへん試しました.
Python + Pillow(PIL)で、透過したpng画像を作成する - Symfoware 画像処理についてあれこれ: Python Imaging Libraryを使用して画像を重ね合わせる

結果的にはこちらの方法です. ありがとうございました.
Python3 OpenCV3で平均値画像(Alpha Blending)を作成 | from umentu import stupid