SVMの学習用アルゴリズムSMOを実装してみる

SVMは２次最適化問題になるので、それを勉強してみてはということだったのですが、SVMに特化したSMO(Sequential Minimal Optimisation)アルゴリズムがあるということなので、そちらをやってみました。

SVMの制約条件に

というのがあって、yiは正例なら1、負例なら-1となる値なのですが、そうすると、ようするにこの条件は、正例のαの合計と負例のαの合計が等しくなるということを示してるわけです。
この条件をつかうと、ひとつαを操作したときには、ほかのαを操作して、正例と負例のバランスを取る必要があることがわかります。
で、このことを利用して、同時に２つのαを操作することにすると、解析的に一つ目のαが求められて、２つ目のαはそこから足し算引き算で求められてお徳かも、というのがSMOの考え方です。
問題は、いかに効率よく更新する２つのαを決めるかということになります。

で、サポートベクターマシン入門に書いてある通りに実装してみたのですが、あんまりうまくいってない。線形分離可能なほうでは、うまくいったりいかなかったりします。線形分離不可能なほうは、実行するたびに同じ分離面が出るのですが、これ、もっと速く計算できていいはず。じゃないとSMOの意味ない。
ひとつ、etaが0以下のときの処理を行っていないので、それが敗因でしょうか。

とりあえずは、ソースを。
実行にはid:nowokay:20080327のGraph.java、id:nowokay:20080326のLerningMachine.javaが必要です。

サポートベクターマシンの本

うちにある本で、サポートベクターマシン(SVM)について書いてある本をあげてみます。

まずは、これ。機械学習ってなんなの？という人におすすめ。パーセプトロンからSVM、ニューラルネットワークときて、そうやってできた学習機械の評価方法についても書いてあります。

フリーソフトでつくる音声認識システム パターン認識・機械学習の初歩から対話システムまで

SVM・カーネルに関しては記述が軽く、これも機械学習の入門書。けど、学習機械を組み合わせるブースティングなどの話が書いてある。

パターン認識と学習の統計学―新しい概念と手法 (統計科学のフロンティア 6)

SVMについて勉強したいなら、これ。読みやすく、SVMについての全体が書いてます。

サポートベクターマシン (知の科学)

SVMのバイブルと呼べるかも。章の並びからして、いまからSVMってなに？ってときにはおすすめできないけど、SVMについての詳しい話が書いてあります。とくに、SVMの実装の話が書いてあるのは日本語の本ではこの本だけなので、実装したいならこの本を読むしかないと思います。

サポートベクターマシン入門

まだよく見てないけど、カーネル法に関する、数学的な位置づけが書いてあるような本だと思う。SVMに関してはあまり書いてなくて、そのあとのRVMとかそういう話につながっている。この本でSVMの概要を知るというよりは、他の本を読んだ後で、数学的な根拠が知りたければ読むといいのかも。

パターン認識と機械学習 下 - ベイズ理論による統計的予測

RでSVMをやる方法が書いてある。ただ、RでSVMをやるためにこの本を読むなら、Ｒとカーネル法・サポートベクターマシンがほぼ同じ内容です。

Rによるデータサイエンス データ解析の基礎から最新手法まで

バイオインフォマティクスで使えるカーネルについて書かれています。カーネルの実例について書いてあるのは、この本だけかも。

バイオインフォマティクスの数理とアルゴリズム (アルゴリズム・サイエンスシリーズ―適用事例編)

ソフトマージンサポートベクターマシン

いままでやってたSVMは、動きがそうなってるかどうかは別として、学習データが誤判別されないように識別面を選ぶのですが、学習データの中にありうる外れデータを考慮できるようにするのがソフトマージンSVMです。
いままでやってたのはハードマージンSVMということになります。
で、ソフトマージンSVMやろうと思ったのですが、最適化部分や収束判定ちゃんとできてないときにやってもあんまり意味なさそうだし、SVMの最適化部分のアルゴリズムであるSMOを実装すると自動的にソフトマージンSVMになるようなので、一回休み。

非線形サポートベクターマシン

とりあえず最適化の問題は置いておいて、ここを参考に非線形分離できるようにしてみました。
http://www.neuro.sfc.keio.ac.jp/~masato/study/SVM/SVM_3_1.htm

さぁ，これで君も非線形SVMのコーディングができちゃうのだ．素晴らしき哉．

ほんとにできた。
うん、ぼくにも非線形SVMのコーディングができちゃいましたよ！
すばらしきかな

SVMも基本は線形分離なので、非線形分離に対応するにはパーセプトロンでやったようにデータの次元を増やしてそこで線形分離します。
で、SVMがすごいのはそこでの計算をごにょごにょして、データの次元を実際には増やさずに高次元で計算したことにしてしまうのです。
SVMでは全データ同士の内積の計算をしていたのですが、その代わりにカーネル関数と呼ばれる関数を使います。カーネル関数には、高次元で内積を計算したことになるようなものを選びます。そうすると、データを高次元に飛ばすことなく、高次元に持っていったのと同じ結果がえられるのです。これをカーネルトリックといいます。
カーネルトリック！

まず、多項式カーネルというカーネル関数を使ってやってみました。こんなカーネルです。

次にガウシアンカーネル。

なんか、もう、まさにサポートベクターマシンという分離面ができました。
2σ^2=1.2だとこんな感じに。

2σ^2=2だとこう。

ソースはこれで。
実行にはid:nowokay:20080327のGraph.java、id:nowokay:20080326のLerningMachine.javaが必要です。

線形サポートベクターマシン失敗

とりあえずの目標はサポートベクターマシンてことで、次のサイト参考に組んでみました。*1
http://www.neuro.sfc.keio.ac.jp/~masato/study/SVM/index.htm

なんか、なぜかちゃんと角度的にはいい感じになってるけど、ぜんぜんマージン最大化していません。見てみたら、λは0以下にならないはずなのに、-4000とかになってるし。

で、よく考えたら、どうにか無理やり

の条件にあてはめてやってたけどλ≧0っていう条件を考えていないし、といろいろ試行錯誤して、どうも単純な最急降下法でやるのがまずいんじゃないかなと思ったのです。
上記のサイトでは「それは考えてみてよ．」と軽く書いてあったけど、これって、研究分野になるじゃないの？
結局、そういうのをどうにかするために、SMO法とかSVMLightとかいう分割法があって、それを使わないといけないということまでわかりました。
SMO法はサポートベクターマシン入門に載ってて & いまのところ実装方法が載ってる唯一の日本語の本なのだけど、福岡においてきてしまった*2。なので福岡に戻るまでおあづけかなぁ。もう一冊買っちゃおうかな。

ただ、上記のサイトは2003年と古くて、そのころはサポートベクターマシンに関する本はまったくでてないし、恐らくWebでも情報がなかったと思うので、しかたないと思われます。
bの求め方も上記サイトの式ではうまくいかなかったので、次のようにしてます。

とりあえず、ちゃんと動いてないソース。
実行にはid:nowokay:20080327のGraph.java、id:nowokay:20080326のLerningMachine.javaが必要です。

*1:サポートベクターマシンについては、はてなキーワード作って書いておきました

*2:今、東京にいる

バックプロパゲーションでニューラルネットの学習

ニューラルネットというのは、入力があって、複数の階層を経て出力を得るようなグラフ構造のことです。通常は、入力層・中間層・出力層のように層構造になっているようなものを差します。中でも、中間層が1層の、３層構造になっているものが多くとりあげられます。バックプロパゲーションは、誤差逆伝播法とも言って、ニューラルネットワークのパラメータを学習するための手法です。

ニューラルネットについてのサイトや本では、中間層を多層に対応した一般的な表現で説明されることが多いのですが、なかなか式を読み解くのが難しかったりするので、今回は3層で入力が2パラメータ、出力は1つ、中間層のニューロンは2つという、単純なものを取り上げます。

では、3層ニューラルネットワークでの判定時のデータの流れを見てみます。
3層ということになっていますが、実際の処理は2層になっています。実装するときには2層だと考えたほうがわかりやすいです。

それでは入力層から中間層g1までの流れを見てみます。まず、それぞれの入力x1、x2に係数w11、w21を掛けて足します。また、バイアスとしてw01という値も足します。
こうして集計した結果に対して、正なら1、負なら0という関数fを適用したものがg1の値になります。
式として書くと、次のようになります。

ここで処理上x0=1という入力値を仮においておくと、実装しやすくなります。

gについても一般化すると、次のようになります。

中間層から出力層までも同様に、g0=1と置いておくと次のようになります。

ここで、関数fに、次のようなしきい値関数を適用すると、fが微分できなくなるので、あとのバックプロパゲーション処理で不都合です。

そこで、tanhや、次のようなシグモイド関数を使います。

このグラフは次のようになります。

シグモイド関数の微分は次のようになるので、都合がよいです。

ニューラルネットワークの学習では、出力uと実際の値bとの誤差から、wやhの値を変更して行きます。
このとき、出力uと学習データbとの誤差から中間層→出力層の係数hを修正します。そして、hを修正した量によって入力層→中間層の係数wを修正します。このように、誤差を逆に伝播させることからバックプロパゲーション、誤差逆伝播法といいます。

具体的には、中間層→出力層の係数hの場合、次のように学習係数k、修正量eとすると

この修正量e1は、出力u、学習データbとして

のようになります。(u-b)が誤差、u(1-u)というのはシグモイド関数の微分になっています。
入力層→出力層の係数wの修正量cとすると

のようになります。前の層の修正量に伝播係数を掛けた分を修正しています。
学習係数kは、修正量をどの程度反映させるかという係数です。

ということで、実装してみたバックプロパゲーションによる学習の状況をみてみると次のようになりました。線形分離可能な場合にもうまく識別できています。*1

また、id:nowokay:20080330でやったようなパーセプトロンでの分離に比べて、両データの中間点に識別面ができています。これは単純なしきい値関数ではなく、シグモイド関数を使ったためともいえます。
実行にはid:nowokay:20080327のGraph.java、id:nowokay:20080326のLerningMachine.javaが必要です。

//BackPropergation.java
import java.util.*;

public class BackPropergation implements LearningMachine {

    List<Map.Entry<Integer, double[]>> patterns =
            new ArrayList<Map.Entry<Integer, double[]>>();
    double[][] w;//入力→中間層の係数
    double[] hidden;//中間層→出力の係数
    int dim;//入力パラメータ数
    int hiddendim;//中間層の数+1

    public BackPropergation(int dim, int hiddendim) {
        this.dim = dim;
        this.hiddendim = hiddendim + 1;
    }

    public static void main(String[] args) {
        new Graph("バックプロパゲーション評価") {

            @Override
            public LearningMachine createLearningMachine() {
                return new BackPropergation(2, 2);
            }
        };
    }

    public void learn(int cls, double[] data) {
        int yi = cls == 1 ? 1 : 0;

        patterns.add(new AbstractMap.SimpleEntry(yi, data));

        final double k = .3;//学習係数
        w = new double[hiddendim - 1][dim + 1];
        for(int i = 0; i < w.length; ++i){
            for(int j = 0; j < w[i].length; ++j){
                w[i][j] = Math.random() * 2 - 1;
            }
        }
        hidden = new double[hiddendim];
        for(int i = 0; i < hiddendim; ++i){
            hidden[i] = Math.random() * 2 - 1;
        }

        for (int t = 0; t < 10000; ++t) {
            //学習を繰り返す
            boolean fin = false;
            for (Map.Entry<Integer, double[]> entry : patterns) {
                double[] pattern = new double[entry.getValue().length + 1];
                for (int i = 0; i < entry.getValue().length; ++i) {
                    pattern[i + 1] = entry.getValue()[i];
                }
                pattern[0] = 1;

                int pcls = entry.getKey();//正解
                double[] hiddenvalue = new double[hiddendim];//中間層の出力値
                //入力層→中間層
                for (int j = 0; j < w.length; ++j) {
                    double in = 0;
                    for (int i = 0; i < pattern.length; ++i) {
                        in += pattern[i] * w[j][i];
                    }
                    hiddenvalue[j + 1] = sigmoid(in);
                }
                hiddenvalue[0] = 1;
                //中間層→出力層
                double out = 0;//出力
                for (int i = 0; i < hiddenvalue.length; ++i) {
                    out += hidden[i] * hiddenvalue[i];
                }
                out = sigmoid(out);
                //出力層→中間層
                double p = (pcls - out) * out * (1 - out);
                double[] e = new double[hiddendim];//中間層の補正値
                double[] oldhidden = hidden.clone();//補正前の係数
                for(int i = 0; i < hiddendim; ++i){
                    e[i] = p * hiddenvalue[i];
                    hidden[i] += e[i] * k;
                }
                //中間層→入力層
                for(int i = 1; i< hiddendim; ++i){
                    double ek = e[i] * oldhidden[i] * hiddenvalue[i] * (1 - hiddenvalue[i]);
                    for(int j = 0; j < dim + 1; ++j){
                        w[i - 1][j] += pattern[j] * ek * k;
                    }
                }
            }
            if (fin) {
                break;
            }
        }

    }

    private double sigmoid(double d) {
        return 1 / (1 + Math.exp(-d));
    }

    public int trial(double[] data) {
        double[] pattern = new double[data.length + 1];
        for(int i = 0; i < data.length; ++i){
            pattern[i + 1] = data[i];
        }
        pattern[0] = 1;
            
        double[] hiddendata = new double[hiddendim];
        //入力層→中間層
        for (int j = 0; j < w.length; ++j) {
            double in = 0;
            for (int i = 0; i < pattern.length; ++i) {
                in += pattern[i] * w[j][i];
            }
            hiddendata[j + 1] = sigmoid(in);
        }
        hiddendata[0] = 1;
        //中間層→出力層
        double out = 0;
        for (int i = 0; i < hiddendata.length; ++i) {
            out += hiddendata[i] * hidden[i];
        }
        return (sigmoid(out) > .5) ? 1 : -1;
    }
}

*1:収束判定してないので、学習に失敗する場合もあります

パーセプトロンで非線形分離するには

パーセプトロンでは線形分離できない場合に対応できないということでしたが、これを可能にする方法をid:nowokay:20080318でしましまさんに教えてもらっていました。
実装してみるとこんな感じで、非線形分離できるようになりました。

x^2、xy、y^2を新たにデータとして追加してるわけですが、要するにこれは２次元データを５次元データに拡張してることになります。2次元じゃ線形分離は無理だけど5次元だと線形分離できますよ、という話。

//NonlinearPerceptron.java
import java.util.*;

public class NonlinearPerceptron implements LearningMachine{
    List<Map.Entry<Integer, double[]>> patterns = 
            new ArrayList<Map.Entry<Integer, double[]>>();
    double b;
    double[] p;
    int dim;
    public NonlinearPerceptron (int dim) {
        this.dim = dim;
    }

    public static void main(String[] args) {
        new Graph("非線形パーセプトロン評価"){
            @Override
            public LearningMachine createLearningMachine() {
                return new Perceptron(2);
            }
        };
    }

    public void learn(int cls, double[] data) {
        int yi = cls == 1 ? 1 : -1;

        patterns.add(new AbstractMap.SimpleEntry(yi, kernel(data)));
        
        final double k = .01;
        b = 0;
        p = new double[dim + 3];
        
        for(int j = 0; j < 1000; ++j){
            //学習を繰り返す
            boolean fin = true;
            for(Map.Entry<Integer, double[]> entry : patterns){
                double[] pattern = entry.getValue();
                int pcls = entry.getKey();
                double in = 0;
                for(int i = 0; i < pattern.length; ++i){
                    in += pattern[i] * p[i];
                }
                in += b;
                if(in * pcls <= 0){
                    //誤判定で再学習
                    fin = false;
                    for(int i = 0; i < p.length; ++i){
                        p[i] += pattern[i] * k * pcls;
                    }
                    b += k * pcls;
                }
            }
            if(fin) break;
        }
        
    }
    public int trial(double[] data) {
        double in = 0;
        double[] d = kernel(data);
        for(int i = 0; i < d.length; ++i){
            in += d[i] * p[i];
        }
        in += b;
        return (in > 0) ? 1 : -1;
    }

    private double[] kernel(double[] data){
        double[] d = new double[data.length + 3];
        for(int i = 0; i < data.length; ++i){
            d[i] = data[i];
        }       
        d[d.length - 3] = d[1] * d[1];
        d[d.length - 2] = d[0] * d[0];
        d[d.length - 1] = d[1] * d[0];
        return d;
    }
}