2018-04-09

Chainer ImageDatasetに与える画像ファイルの一覧を作成する

ChainerのImageDatasetを使用する際，引数pathsにはデータセットに含める画像のファイル名の一覧を渡す必要がある．
画像ファイルの一覧(.txt形式)を作成するために，以下のlinuxコマンドを使用した：

ls image_dir | grep jpg > list.txt

image_dirは画像が含まれているディレクリ，jpgは画像の拡張子，list.txtは出力される画像ファイル一覧を意味する．
ここで，実行環境はmacOS Sierra バージョン10.12.6とする

後学のために，失敗例も記載する：

ls image_dir | grep jpg > list.txt

この場合，list.txtも一覧に含まれて記述されてしまうため，期待される出力を得ることができない．

2018-02-15

ヒストグラムで頻度順にソートする

メモ

matplotlibで頻度分布を作図するときに，頻度順にソートしたかった
結論：mattplotlib.pyplot.hist()のxには均一な頻度のデータを与え，ソートした頻度をweightsとして与えれば良い

参考

matplotlib.pyplot.hist — Matplotlib 2.1.2 documentation

2018-01-26

CuPyで配列を部分配列に切り分けて結合してから平均

メモ WIP

メモ．後でちゃんと書く
結論：次が最速

cupy.mean(cupy.stack(cupy.split(data, split_size)), axis=axis)

ダメなケース：cupy.split(data, split_size)がリストで返ってくるからといって，cupy.array()でcupy配列にしてから平均を取ろうとすると，GPUとの無駄な通信が発生して，実行速度が遅くなる．GPU・CPU間の通信は極力減らすのが鉄則．

cupy.mean(cupy.array(cupy.split(data, split_size), dtype='f'), axis=axis) # かなり遅い
cupy.array(numpy.mean(numpy.split(cupy.asnumpy(data), split_size), axis=axis), dtype='f') #少しマシになったけどまだ遅い

2017-12-09

共分散行列適応進化戦略(CMA-ES)と自然勾配法

進化計算最適化 kosen10s AdventCalender CMA-ES 情報幾何

この記事は #kosen10s Advent Calendar 2017 の 9日目の記事です．前回は marin72_com さんの「#kosen10s だより 2017」でした． kosen10s，自由気ままにだけど秩序持ってのんびりやってていいところだと思います．僕みたいなプログラミング雑魚勢でもちゃんと居場所があるので，少しでもいいなと思った2010年卒の高専生はslackチームへおいでなさいまし．

はじめに

今回は自分の研究テーマの1つである進化計算について，ちょっと紹介しようかと思います．まだまだ勉強中であるため，ところどころぼかしますがご容赦を．

遺伝的アルゴリズム(Genetic Algorithm; GA)や蟻コロニー最適化(Ant Colony Optimization; ACO)など進化計算にも様々な種類がありますが，対象とする問題は基本的に同じです．進化計算が取り扱う問題においては，最適値を求めたい関数，つまり目的関数について次の1や2のような制約が課されている場合が多いです．

目的関数を明に書き表すことができない（ただし，入力値に対応する出力値を得ることはできる）
目的関数の導関数が不明である

このような制約の中で行う最適化は「Black-box最適化」という問題に属します．つまり，進化計算はBlack-box最適化を行うための最適化アルゴリズムと言うことができます．ちなみに，Black-box最適化アルゴリズムはいくつか提案されており，例えばベイズ的最適化(Bayesian Optimization)もこれに該当します．

Covariance Matrix Adaptation Evolution Strategies

今回は，数あるBlack-box最適化アルゴリズムの中から，性能が比較的良好なCMA-ES(Covariance Matrix Adaptation Evolution Strategies)*1について紹介します．CMA-ESは現在でも改良が続けられており様々な手法*2 *3が提案されています．

CMA-ESとは，目的関数への入力値（以後，サンプルと呼びます）を平均ベクトルと共分散行列をパラメータに持つ正規分布からサンプリングし，得られた出力値を基に正規分布のパラメータを更新，更新された正規分布からまた新たにサンプルをサンプリングし,...という操作を繰り返すことによって，最適なサンプルを得るという手法です．CMA-ESにおける最適化の概要図を次のFig.1に示します．オレンジの点線の形状を決めているのが分散共分散行列で，平均ベクトルが中心を決めています．

f:id:satuma77:20170623230903p:plain

Fig.1 CMA-ESによる最適化過程の概要図*4

また，最適化をアニメーションにしたものを次のFig.2に示します．最適解は左上の箇所となります．

f:id:satuma77:20171208235245g:plain

Fig.2 2次元Rastrigin関数に対するCMA-ESの最適化過程*5

Fig.2からわかるように，探索の過程でサンプルの広がりが大きく広がり，最適解で収束するという現象が見られます．この現象は，CMA-ESの工夫のひとつであるCommurative Step-Size Adaptation(CSA)によるものです．CMA-ESにおける共分散行列の更新には，CSA, Rank-1 Update, Rank-μ Updateの3種類の更新規則が適用されています．それぞれ説明すると非常に長くなりますので，今回はRank-μ Updateに絞って説明します．また，平均ベクトルの更新もRank-μ Updateで行われています．そこで，次節ではRank-μ Updateにおける平均ベクトルと共分散行列の更新方法について説明します．

Rank-μ Updateと確率的勾配降下法

ある時点 $t$ において， $n$ 番目のサンプルを $x_{n}$ ，平均ベクトル $m^{(t)}$ ，共分散行列を $C^{(t)}$ としたとき，CMA-ESのRank-μ Updateは次式で定義されます*6．

$m^{(t+1)} = m^{(t)} + \eta_{m} \sum_{i=1}^{\lambda} w_{\mathrm{rk}}(x_{i}) (x_{i} - m^{(t)})$

$C^{(t+1)} = C^{(t)} + \eta_{C} \sum_{i=1}^{\lambda} w_{\mathrm{rk}}(x_{i}) \left( x_{i}-m^{(t)})(x_{i}-m^{(t)})^{\mathrm{T}}-C^{(t)} \right)$

$\eta_{m}$ および $\eta_{C}$ は学習率（ユーザが指定するパラメータ）で，一回の更新に使用するサンプル数を $\lambda$ としています．また， $w_{\mathrm{rk}}(x_{i})$ はサンプル $x_{i}$ が $\lambda$ 個のサンプルの中で何番目に良い出力値かという「ランキング」に基づいた値を意味します．CMA-ESについて，これらの更新式はハンドメイドによって作成されたものとなっています．

これらの更新式は2010年に自然勾配と深い関連があることが明らかになっています*7．自然勾配*8は，確率分布モデルが形成する多様体が作る幾何的構造を分析する情報幾何学の中で発展した手法となります．

確率分布パラメータを「勾配」を用いて更新することを考えます．このとき，パラメータ間の擬距離にユークリッド距離ではない距離を採用した場合，ユークリッド空間における勾配は最急方向とは異なります．この場合の最急方向は，勾配にリーマン計量テンソルの逆行列をかけたものとなります．これが自然勾配と呼ばれるものです．情報幾何学において，確率分布モデルに対するパラメータ間の擬距離にはKL-Divergenceを採用する場合が多く，その場合のリーマン計量テンソルはフィッシャー情報行列(Fisher Information Matrix; FIM)となることが知られています．

話を戻してPure Rank-μ Updateの場合を考えます．多次元正規分布 $\pi(x; m, C)$ の対数尤度関数 $\ln \pi(x; m, C)$ に対して自然勾配を求めると

$\tilde{∇}_{m} \ln \pi(x; m, C) = (x_{n} - m)$

$\tilde{∇}_{C} \ln \pi(x; m, C) = (x_{i}-m)(x_{i}-m)^{\mathrm{T}}-C$

となります．この自然勾配を見るとPure Rank-μ Updateの更新式に類似した式がそのまま現れていることがわかります．

実際は，解くべき最適化問題を目的関数の期待値最小化問題 $\mathrm{argmin}_{x} \mathbb{E}_{\pi(x; m, C)}f(x)$ とした上で，自然勾配を用いるとPure Rank-μ Updateの式を導くことができます．詳しくは文献*9を参照していただければと思います．

追記(2017.12.09)

時間の許す限り，導出過程を書いていきます．目的関数の期待値最小化 $\mathrm{argmin}_{x} \mathbb{E}_{\pi(x; m, C)}f(x)$ を行うために， $\mathbb{E}_{\pi(x; m, C)}f(x)$ の勾配を計算します．文献*10に従って $\theta = \lbrack m^{\mathrm{T}}, \mathrm{vec}(C)^{\mathrm{T}} \rbrack ^{\mathrm{T}}$ とすることで，2つのパラメータに関する勾配をまとめて求めます． $\frac{\partial \log f(x)}{\partial x} = \frac{\partial f(x)}{\partial x} \frac{1}{f(x)}$ というテクニックを使えば，

$∇_{\theta} \mathbb{E}_{\pi(x; m, C)}f(x)$ $= ∇_{\theta} \int f(x) \pi(x; m, C) \mathrm{d}x$ $= \int f(x) ∇_{\theta} \pi(x; m, C) \mathrm{d}x$ $= \int f(x) \frac{∇_{\theta} \pi(x; m, C)}{\pi(x; m, C)} \pi(x; m, C) \mathrm{d}x$ $= \int \left( f(x) ∇_{\theta} \log \pi(x; m, C) \right) \pi(x; m, C) \mathrm{d}x$ $= \mathbb{E}_{\pi(x; m, C)} \left\lbrack f(x) ∇_{\theta} \log \pi(x; m, C) \right\rbrack$

とこのように勾配は期待値の形で書き表すことができます．今回 $f(x)$ が未知であるため，この期待値を代数的に求めることは難しいです．そこで，モンテカルロ法によって勾配を推定します． $\lambda$ 個のサンプルを使って勾配を近似すれば，

$∇_{\theta} \mathbb{E}_{\pi(x; m, C)}f(x) \approx \frac{1}{\lambda} \sum_{i=0}^{\lambda} f(x) ∇_{\theta} \log \pi(x; m, C)$

と求めることができます．このとき，FIMは次の式となります（文献*10）． $F = \lbrack \lbrack C^{-1}, 0 \rbrack, \lbrack 0, \frac{1}{2}(C^{-1} \bigotimes C^{-1}) \rbrack \rbrack$

FIMの逆行列を先程の勾配の式と掛け合わせると，自然勾配は

$\tilde{∇}_{m} \mathbb{E}_{\pi(x; m, C)}f(x) = \sum_{i=1}^{\lambda} f(x_{i}) (x_{i} - m^{(t)})$ $\tilde{∇}_{C} \mathbb{E}_{\pi(x; m, C)}f(x) = \sum_{i=1}^{\lambda} f(x_{i}) \left( x_{i}-m^{(t)})(x_{i}-m^{(t)})^{\mathrm{T}}-C^{(t)} \right)$

となります．ここで悪スケールな関数に対してうまく最適化できるように，自然勾配を目的関数 $f(x)$ のスケールに依存しないように変更します．その変換方法として目的関数値 $f(x)$ のかわりに，解のランキングに基づいた値 $w_{\mathrm{rk}}(x_{i})$ を使用します．これで，Pure Rank-μ Updateが自然勾配によって更新するアルゴリズムということが確認されました．

Pure Rank-μ UpdateのPython実装

先程のPure Rank-μ Updateについては，私が現在作成しているPython用の進化計算ライブラリEvoltierにて実装済みですので，軽く紹介します．次にEvoltierにおけるPure Rank-μ Updateの実行方法を示します．（リポジトリにあるmain.pyを一部改編）

肝心のPure Rank-μ Updateの部分はevoltier.optimizers.GaussianNaturalGradientOptimizerにあります．実際に動かす際には，最適化したい関数などをevoltier.optimizers.updater.Updaterに入れてから最適化を行うように作成しています．Updaterの初期化におけるパラメータについて重要な部分のみ触れていきます．pop_sizeは $\lambda$ に相当し，thresholdは停止条件に関するパラメータです．具体的には目的関数の値がこの値を下回るサンプルを得られれた時点で停止するというものになります．max_iterも停止条件に関するパラメータであり，1試行における最大更新回数を設定するものとなります．

実行速度については100次元のSphre関数で $\lambda=1000$ としたとき，評価値が $10^{-30}$ となる値を見つけるまでに約40秒ほどでした．まだ，並列化には対応していないこともあり個体数が多いときはmultiprocessingモジュールの利用によって高速化も期待されます．

おわりに

今回はBlack-box最適化アルゴリズムの1つであるCMA-ESについて，説明しました．また，Pure Rank-μ Updateによる分布更新が自然勾配との関連があることも紹介しました．最後に現在作成中のライブラリEvoltierにおいて，Pure Rank-μ Updateによる最適化を実行する方法を示しました．CMA-ESについて，その派生アルゴリズムやRank-1 Updateなどの他の更新ルール，制約付き問題に対する適用方法など，もし時間があればご紹介できればと思います．また，Evoltierについてもゆっくりですが開発を進めていき，時期を見てPypiに上げて実際に使えるようにできたらと思います．とりあえずCuPyによるGPU対応も盛り込むつもりで作っています．もしよければStarください．

次回はNKMR6194さんの「なぜpthread_mutexは機能するのか？」です!

*1:N. Hansen and A. Ostermeier, “Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation,” Proceedings of IEEE International Conference on Evolutionary Computation, pp. 312–317, 1996.

*2:N. Hansen, “Benchmarking a BI-population CMA-ES on the BBOB-2009 noisy testbed,” Proceedings of the Genetic and Evolutionary Computation Conference - GECCO ’09, pp. 2389-2396, 2009.

*3:Y. Akimoto and N. Hansen, “Projection-Based Restricted Covariance Matrix Adaptation for High Dimension,” Proceedings of the Genetic and Evolutionary Computation Conference - GECCO ’16, pp. 197-204, 2016.

*4:https://en.wikipedia.org/wiki/CMA-ES

*5:http://blog.otoro.net/2017/10/29/visual-evolution-strategies/

*6:S. Shirakawa, Y. Akimoto, K. Ouchi, and K. Ohara, “Sample Reuse in the Covariance Matrix Adaptation Evolution Strategy Based on Importance Sampling,” in Proceedings of the 2015 on Genetic and Evolutionary Computation Conference - GECCO ’15, 2015, pp. 305–312.

*7:Y. Akimoto, Y. Nagata, I. Ono, and S. Kobayashi, “Bidirectional relation between CMA evolution strategies and natural evolution strategies,” Parallel Problem Solving from Nature (PPSN XI), pp. 154–163, 2010.

*8:S. Amari, “Natural Gradient Works Efficiently in Learning,” Neural Computation, vol. 10, no. 2, pp. 251–276, 1998.

*9:Y. Ollivier, L. Arnold, A. Auger, and N. Hansen, “Information-Geometric Optimization Algorithms: A Unifying Picture via Invariance Principles,” Journal of Machine Learning Research, vol. 18, no. 18, pp. 1–65, 2017.

*10:Y. Akimoto, “Analysis of a Natural Gradient Algorithm on Monotonic Convex-Quadratic-Composite Functions,” Proceedings of the fourteenth international conference on Genetic and evolutionary computation conference - GECCO '12, pp. 1293–1300, 2012.

2017-09-20

TED Talks「Meg Jay: Why 30 is not the new 20」

雑記

TEDを眺めていたら，なんか感じるものがあったので，忘れないうちにブログへembedded．

大学院に進学して就職が遅くなっていて，確かにいろんな問題を先送りにしている感じはあるから，結構感じることは色々．

20代はあと7年弱しかないってマジですか…

2017-08-12

scikit-learnによる特徴選択

機械学習

お盆ですな昨日書いたの続き

satuma-portfolio.hateblo.jp

Motivation

特徴選択がやりたい
Wekaはお手軽
でも僕はPythonでやりたい

ソースコード

さすがscikit-learn．めちゃシンプルにできた．対象とするデータはirisデータセット ( UCI Machine Learning Repository: Iris Data Set )

詳細

scikit-learnでの特徴選択に関する関数はsklearn.feature_selectionに含まれている．今回は特徴量にランキングを付けて，そのランキングの上からk個を採用するというsklearn.feature_selection.SelectKBestを採用*1．

次にランキングをつけるための評価関数を設定する必要がある．今回は相互情報量に基づいてランキングをつけるmutual_info_classifを採用*2．

あとはSelectKBest(選択した評価関数, k=ベストいくつか).fit(データ行列, ラベル)とすればOK．もし，どの特徴を採用したか知りたい場合は，get_support()関数を使うことでboolean配列として取得することができる．

最後に

wrapper法がない..?
- sklearn.feature_selection.SelectFromModelがwrapper法か？と思ったけど，なんか違う…
「wrapper法はそんな難しくないんだから自分で組め」というアレなんだろうか…?

*1:他のものは API Reference — scikit-learn 0.19.0 documentation を参考に．多分Selectとついてるものが選択関連

*2:他のものは API Reference — scikit-learn 0.19.0 documentation を参考に．多分小文字で始まっているものが評価関数

2017-08-11

機械学習ソフトウェアWekaの使い方のメモ

機械学習

なんだかんだ初投稿からかなり時間が経っていた（）

Motivation

特徴選択の話をかじり始め，とりあえず従来手法をサクッと試したかった
そこでJavaで作られたフリーのソフトウェアwekaを使うことを思いつく
インストール込みで使い方をメモる

Wekaとは

Wikipediaでの記述では

Weka (Waikato Environment for Knowledge Analysis) は、ニュージーランドのワイカト大学で開発した機械学習ソフトウェアで、Javaで書かれている。GNU General Public License でライセンスされているフリーソフトウェアである。 – Weka - Wikipedia

とのこと．公式のWebページを覗いてみると

Weka is a collection of machine learning algorithms for data mining tasks. The algorithms can either be applied directly to a dataset or called from your own Java code. Weka contains tools for data pre-processing, classification, regression, clustering, association rules, and visualization. It is also well-suited for developing new machine learning schemes. – Weka 3 - Data Mining with Open Source Machine Learning Software in Java

とのこと．前処理から結果の可視化まで一括してやりまっせって感じですね．

インストール

このへんを参考に．

Wekaのダウンロード – Wekaの日本語情報
Weka 3 - Data Mining with Open Source Machine Learning Software in Java
- 最新の安定版やdeveloper版が欲しい場合は公式Webページから取るのが安牌

最新ver.は3.8系．ただし，日本語化に対応しているのは3.6系の模様．とりあえず3.6系を入れました．

使い方

このへんを参考に

クイックスタート – Wekaの日本語情報 d.hatena.ne.jp d.hatena.ne.jp

特徴選択

このへんを参考に． http://www.weka-jp.info/archives/Weka_FeatureSelection-1.pdf

フィルタ法だと

Infomation Gain
ReliefF

という鉄板はもちろんありそう．

ラッパー法でも

Genetic Algorithm
Greedy Search(貪欲法)
全探索

がありそう．まず取っ掛かりとして試すのは十分．

最後に

鉄板のirisデータセットをはじめ，データセットもまあまああるのでサクッとはじめられそう
自作アルゴリズムをインポートできるかな？
(実際，提案手法の実装はPythonでやりたいので，wekaは従来手法のチェックで終わりなんだけど…)
次は本番で使うscikit-learnも漁る

そういえば

最近，CMA-ESのpure rank-μ updateだけ組んだ．main.pyでらくらくお試しできます． github.com

こっちの技術的背景というかInformation Geometric Optimizationの話も後でまとめたい． IGOフレームワークは今後，Bernoill分布にも対応させたいし，Gauss分布のほうでもCommulative Step-size Adaptationだったり，rank-1 updateだったり，学習率の自動調整だったり，やることはまだまだある．お星様付けてくれると悦びます．