ベイズ推論による機械学習入門を読んだので軽くメモ

はじめに

タイトルの通りです。以下の本を一通り読んだので、再度読む際の手助けになるようメモを残します。

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

まだ理解が怪しい部分も多々ありますが、取りあえず忘れる前に。。。

全体を通してのポイント

ベイズ推論の枠組み

ベイズ推論においては、観測データも分布のパラメータ、回帰等の重み係数も潜在変数も確率変数として扱う。観測されたデータをとし、未知の確率変数（分布のパラメータ、重み係数、潜在変数等々）をとした時、以下の2ステップを行うのがベイズ推論の枠組み。

1. モデルの構築（モデリング）
2. 事後分布の計算

この事後分布を求めることがベイズ推論における「学習」。ベイズ推論は全てこの枠組みの中で扱われるので、後は如何にモデリングするかと、如何に求めたい確率変数の事後分布を計算するか、という話になる。モデルが複雑な場合、事後分布は解析的に求められない場合もあるので、その時は変分推論やギブスサンプリングなどを使って近似的に事後分布を求める。

予測分布

更に、求めた事後分布を使って、未観測のデータに対し、何らかの知見を得たい場合もある（知見の内容はタスクによって異なる。例えば、回帰の場合は、入力に対応する出力を求める等）。その時は以下の予測分布を計算する。

1章

基本的にこの章は後の章で語られる内容のアウトラインが載っている。この章で語られる内容で特に重要なのは、上記のベイズ推論のポイントと以下の関係式。

• 周辺化
• 条件付き分布

また、独立の諸概念も重要。

• 変数間の独立性（これが成り立つ時、とは独立という）
• 条件付き独立（を観測した下で、は独立だが、を観測しなければとは独立でないことに注意）

2章

各確率密度関数については書籍を参照するとして、取りあえず期待値が大事。エントロピーとかKLダイバージェンスとかも大事だけど、期待値がわからないと後々お話にならなくなる。特に、

と、期待値を取る関数に出てこない確率変数は、期待値計算から外せることは重要。
また、同時分布の期待値とか、条件付き期待値とかは込み入ってくると混乱するので注意。迷った時は期待値の定義、

に立ち返るのが吉。

各確率密度関数の対数を取った関数、基本的な期待値（とかとか）も、後々の章（特に4章、5章の変分推論辺り）で何度も使う。本当に何度も使うので、何度も見返す。

3章

2章で挙がっている確率分布の各パラメータおよび線形回帰の重み係数をベイズ推論の枠組みで推定する章。事前分布として共役事前分布を使うので、（この章で示されている問題については）事後分布が解析的に計算でき、しかも事前分布として同じ形の分布で求まるので、データを得るたびに逐次的に更新することも可能。考え方としては、

1. 学習したいパラメータの事前分布として共役事前分布を設定
2. こうすると事後分布は事前分布と同じ形になり、　としては定数として扱うことが可能（後から帳尻を合わせる）。
3. 尤度はと観測データが独立だと仮定（後の章で観測データ間にも独立ではないモデルもあり）

また、計算のテクニックとして、以下のようなものが使われている。

• 事後分布が事前分布と同じ形になることから、をうまく展開し、2章で計算済みの事後分布の対数をとった式のパラメータと対応を取ることで、事後分布のパラメータを計算
• 予測分布を求める際、積分の計算が難しい時は、ベイズの定理に立ち返りとなることから、上記の事後分布計算と同じく対数をとって分布の式とパラメータの対応を取ればよい
• 予測分布計算時、いきなり事後分布のパラメータで予測分布求めようとすると、ややこしくなるケースもあり、その場合はまず事前分布を用いて予測分布を計算し、後から対応するパラメータを事後分布のものに置き換えればよい

4章

生成モデル

観測データの生成過程を仮定して、確率分布をモデリングしたものが生成モデル。例えば、観測されたデータに多峰性が見れられる場合、これを単一の分布のみで表現しようとしても、うまく表現しきれない。このような多峰性のあるデータは、

1. 個のクラスタがあり、観測データはそのいずれかのクラスタ（例えば、を持つ番目のクラスタ）から生成されたと仮定
2. が割り当たっているクラスタを（割り当たっているクラスタの要素のみ1、それ以外0）で表現すると、の確率分布は
3. このクラスタの割り当ても割り当て比率パラメータに支配されるカテゴリ分布の確率変数と仮定
4. 各クラスタのパラメータの事前分布を仮定（ここでも共役事前分布使うと幸せになれる）
5. クラスタ割り当てのカテゴリ分布の事前分布を仮定（こちらもカテゴリ分布の共役事前分布であるディリクレ分布を使うと幸せ）
6. 以上のことから、データ生成のモデルをと仮定

等のようにモデリングする。
このモデルから、観測データに対するクラスタ割り当ての事後分布　を計算したいが、モデルが複雑すぎてが解析的に求まらない。そのため、近似的に事後分布の推定を行う。

ギブスサンプリング

上のモデルを例にすると、

1. が既に観測されているものとして（2でサンプリングしたを使う）、とサンプリング
2. が既に観測されているものとして（1でサンプリングしたを使う）、とサンプリング
3. この操作を繰り返す

この繰り返し回数が十分に多い場合、真の分布に近づくことが知られているため、繰り返し続けると所望のデータをサンプリングできるようになる。もちろん、サンプリングする対象となる分布、は十分にサンプリングしやすい分布である必要があるので、式を展開して2章で挙がっているような基本的な分布を導かないといけない（モデルによってはそれができないこともあるので、共役事前分布を使うと幸せになるというのはここに繋がる）。

変分推論

等のように、各確率変数間の独立性を仮定した分布で近似するのが変分推論。真の分布とのKLダイバージェンスが最小となるを求めればよく、それは

として計算できる（あくまで、上記モデルに対する近似式。もう少し一般的な表現は書籍参照）。こちらも一方の結果を用いで、もう一方の計算を行う形になっているため、更新を十分な回数繰り返す必要がある。また、実際には、下記のような式展開して期待値を解析的に求める。

1. 変分推論を基に、各近似分布の形を求める。2章に述べられている期待値の性質（上記の無関係な確率変数が期待値計算から外せる性質含む）をうまく使い、が何分布の関数かがわかるまで式を展開する。
2. 全ての近似分布の形がわかったら、それぞれの期待値を計算する。これまた2章で計算済みの各分布の基本的な期待値を流用する。

ここでもギブスサンプリングと同様、期待値を解析的に計算するために、2章で挙がっているような基本的な分布を導く必要があるので、やはり共役事前分布を使った方がハッピー。

崩壊型ギブスサンプリング

上記のギブスサンプリングでは、パラメータについてもサンプリングで求めたが、

として周辺化すれば、からサンプリングするだけでよくなる。ただし、周辺化した結果が十分にサンプリングしやすい確率分布になることが条件。このモデルの例では、周辺化した結果が互いに依存し合うようになり（元々はを与えられた下での条件付き独立）、このままでは個の潜在変数を同時にサンプリングしなければならず計算量が莫大になる。この問題をと、サンプリングしたい以外の潜在変数は全て観測済みとして、個別にサンプリングするアイデアで解決。

5章

各種応用に対して、モデルの構築と変分推論（一部ギブスサンプリング）を導出している。この中のいくつかは実際に実装してみようと思うので、ここでは考え方のみメモ。

線形次元削減

観測されたデータが低次元の潜在変数を元に、という線形変換で生成されると考える。また、実際にはノイズが乗っているはずなので、ガウス分布で生成されると仮定。更に他のパラメータについても事前分布がガウス分布に従うと仮定（は既知のものとして話を進めているが、ガンマ事前分布など使って推定してもよい）。後は変分推論の期待値計算をすれば低次元の潜在変数や変換行列を推測することができる。

非負値行列因子分解

観測された行列が全ての要素非負であるとし、という非負値の行列に近似分解することを考える。ここで、となる潜在変数を考えることで、効率的な変分推論が導き出せる。また、非負値なのではポアソン分布、の事前分布はガンマ分布を考える。

隠れマルコフモデル

これまで観測データ間はたがいに独立して得られるものとして扱っていたが、時系列データ等ではその仮定は成り立たない。隠れマルコフモデルでは、

1. 観測されるデータは現在の状態に基づいてデータが観測される（状態をクラスタ割り当ての潜在変数として持つ混合分布）。
2. データを観測したら、遷移確率に基づいて次の状態に遷移する（）。

というモデルを構築する。変分推論は状態の潜在変数の分布をと分離して近似する完全分解変分推論とと分離せずに近似する構造化変分推論がある（その中間であるミニバッチを使った近似もあり）。前者は特別なことをせずに、変分推論の期待値を導き出せるが、後者はうまいこと式変形してフォワードバックワードアルゴリズムと呼ばれる動的計画法が使える形に持っていかないと、計算量が増大する。

トピックモデル

複数の文書を観測データとして、単語の出現頻度を元に文書を解析。各文書中には複数の潜在トピックが存在すると仮定。文書中のn番目に出てくる単語には、どのトピックから発生したかというトピック割り当て潜在変数が紐づき、各単語は対応するトピックのカテゴリ分布に従うと仮定。また、トピック割り当て自体もカテゴリ分布に従うと仮定し、これらのカテゴリ分布のパラメータはディリクレ事前分布に従うモデルを構築。後は頑張って変分推論の期待値計算。

テンソル分解

協調フィルタリングに時間軸を導入し、その時のトレンドを踏まえた推薦を行うようにしたモデル。時間方向の軸には前回の状態を平均とするガウス分布を仮定（線形動的システム）。ここでは、完全分解変分推論を考えているので、推論計算自体はこれまでと同様。

ロジスティック回帰、ニューラルネットワーク

ロジスティック回帰、ニューラルネットワークのベイズ推論版。これまでの近似（平均場近似）とは違い、ガウス分布で事後分布を近似。変分パラメータを用いた勾配法で真の分布とのKLダイバージェンス最小化を行う。この勾配法を使うため、再パラメータ化トリックを使うが、これが良く分からない。

おわりに

兎にも角にもわかりやすいです。式展開がわからなくなって何度も行ったり来たりしましたが、PRML10章で卒倒した僕でも割りと短期間で最後まで行けたので、ベイズ推論の入り口として、適切な書籍の1つだと思います。ただ機械学習自体の入門としては少ししんどいので、機械学習全体を俯瞰したのち、その1つの枠組みとしてベイズ推論があると捉えられるようになった上で、読むべきかなと思います。

競プロとかに使うアルゴリズム実装メモ（二分探索、2次元累積和、しゃくとり法）

はじめに

アルゴリズムメモ第３段です。今回は二分探索法、２次元累積和、しゃくとり法と様々な問題に使える汎用的なアルゴリズムを書いていきます。
今回は勉強のため、アルゴリズムの本質的な部分を記述した抽象クラスと実際の問題を解く具象クラス（関数オブジェクト）を分け、テンプレートメソッドパターンもしくはストラテジパターンを用いてコーディングしました。

二分探索

用途

• 広義単調増加関数または広義単調減少関数が条件を満たし、かつ、条件成立と不成立の境界となる点を探索（条件成立と不成立の境界が１つのみ）
• ソート済み配列から条件が成立する値が存在するかの判定

アルゴリズムのポイント

1. 変数を探索範囲の最小値、変数を最大値で初期化
2. をの中点（例えば、で求める）として、が条件を満たすか否かを判定
   • このとき、見つけたい点（条件成立と不成立の境界）がより小さいか、大きいかはが広義単調増加（減少）関数なのでわかる
3. 見つけたい点が以上なら、そうでなければとする
4. この2～3の操作をを満たすまで繰り返す（は収束判定基準の定数で、問題に応じて適切な値を定める）
5. 収束したら、その時点のから最終的な出力値を決める

抽象クラスの実装

上記アルゴリズムを実現する抽象クラスを作成し、見つけたい点が現時点以上かの判定、条件判定、中点計算、戻り値計算を純粋仮想関数とし、上記アルゴリズムをテンプレートメソッドパターンで実装。

#include <vector>

template <typename T>
class abstract_binary_search
{
private:
    virtual bool is_in_right(T x) = 0; //求めているのがxより大きいか
    virtual bool is_converge(T L, T R) = 0; //収束判定
    virtual T ret_val(T L, T R) = 0; //見つからなかった場合に返す値
    virtual T get_middle(T L, T R) = 0; //中央値を返す

protected:
    T search_rec(T L, T R);

public:
    virtual T search() = 0;
    virtual bool is_match(T x) = 0;  //それが求めている値か

    virtual ~abstract_binary_search();
};

template <typename T>
T abstract_binary_search<T>::search_rec(T L, T R) {
    if (is_converge(L, R) == true) {
        return ret_val(L, R);
    }

    T mid = get_middle(L, R);

    if (is_in_right(mid) == true) {
        return search_rec(mid, R);
    }

    return search_rec(L, mid);

}


template <typename T>
abstract_binary_search<T>::~abstract_binary_search()
{
}

具体的な問題を解く実装

Lower bownd

要素数の数列が与えられて、指定された値以上の値が入っているindexを求める問題。もし数列内に以上の値がなければを返す。

#include <iostream>

class lower_bound :
    public abstract_binary_search<int>
{
private:
    std::vector<int> data_list;//検索する数列
    int N; //数列のサイズ
    int search_data;//検索対象データ
    virtual bool is_in_right(int x); //求めているのがxより大きいか
    virtual int ret_val(int L, int R); //見つからなかった場合に返す値
    virtual bool is_converge(int L, int R); //収束判定
    virtual int get_middle(int L, int R);
public:
    lower_bound(const std::vector<int>& data_list);
    void set_search_data(int a_search_data);
    virtual int search();
    virtual bool is_match(int x);  //それが求めている値か
    ~lower_bound();
};

bool lower_bound::is_in_right(int x) {
    return data_list[x] <= search_data;
}


int lower_bound::ret_val(int L, int R) {
    if (search_data <= data_list[L])
        return L;
    else
        return L + 1;
}

int lower_bound::get_middle(int L, int R) {
    return (L + R) / 2;
}

lower_bound::lower_bound(const std::vector<int>& data_list):
    data_list(data_list)
{
    N = data_list.size();
}

int lower_bound::search() {
    return search_rec(0, N);
}

bool lower_bound::is_match(int x) {
    if (x >= N) return false;
    return data_list[x] == search_data;
}

bool lower_bound::is_converge(int L, int R) {
    return R - L <= 1;
}

void lower_bound::set_search_data(int a_search_data) {
    search_data = a_search_data;
}

lower_bound::~lower_bound()
{
}

int main()
{

    int n, q;
    std::cin >> n;
    std::vector<int> S(n);
    
    for (int i = 0; i < n; i++) {
        int  s;
        std::cin >> s;
        S[i] = s;
    }
    std::cin >> q;

    int count = 0;
    
    lower_bound lbs(S);
    for (int i = 0; i < q; i++) {
        int  t;
        std::cin >> t;
        lbs.set_search_data(t);
        
        int res = lbs.search();
        std::cout << res << std::endl;
        if (lbs.is_match(res)) 
            count++;
    }
    std::cout << count << std :: endl;
    return 0;
}

POJ 1064：Cable master

蟻本に二分探索の例として載っている問題。問題はリンクを参照。が切り取る紐の長さ、が作れる同じ長さ（長さ）の数だとすると、は広義単調減少関数なので二分探索法が使える。
つまり、が条件を満たしつつ、が最大となる点を二分探索する（条件を満たさなければを小さくし、満たすならまだ長くできると考えを大きくする）。

#include <vector>
#include <algorithm>
#include <functional>
#include <iostream>
#include <iomanip>
#include <math.h>

class cable_master :
    public abstract_binary_search<double>
{
private:
    double tol;
    std::vector<double> len_list;
    int K;
    int N;
    bool is_condition(double x);
    virtual bool is_in_right(double x); //求めているのがxより大きいか
    //virtual bool is_in_left(double x); //求めているのがxより小さいか
    virtual double ret_val(double L, double R); //見つからなかった場合に返す値
    virtual double get_middle(double L, double R);
    virtual bool is_converge(double L, double R); //収束判定
public:
    cable_master(const std::vector<double>& len_list, int K, double tol);
    virtual double search();
    virtual bool is_match(double x);
    ~cable_master();
};

bool cable_master::is_condition(double x) {
    int count = 0;
    for (int i = 0; i < N; i++) 
        count += int(len_list[i] / x);

    if (count >= K)
        return true;

    return false;
}

bool cable_master::is_in_right(double x) {
    return is_condition(x);
}

double cable_master::ret_val(double L, double R) {
    return R;
}

double cable_master::get_middle(double L, double R) {
    return (L + R) / 2.0; 
}

bool cable_master::is_converge(double L, double R) {
    return R - L < tol;
}

cable_master::cable_master(const std::vector<double>& len_list, int K, double tol) :
    len_list(len_list),
    K(K),
    tol(tol)
{
    N = len_list.size();
}

double cable_master::search(){
    return search_rec(0.0, *std::max_element(len_list.begin(), len_list.end()));

}

bool cable_master::is_match(double x) {
    return false;
}

cable_master::~cable_master()
{
}



int main()
{

    std::vector<double> cable_list;
    int N, K;
    std::cin >> N >> K;

    for (int i = 0; i < N; i++) {
        double len;
        std::cin >> len;
        cable_list.push_back(len);
    }
    cable_master cm(cable_list, K, 0.0001);
    std::cout << std::fixed << std::setprecision(2) << floor(cm.search() * 100) / 100.0 << std::endl;
}

2次元累積和

用途

• の計算

アルゴリズムのポイント

1. あらかじめ、全てのについてを求めておく
   • これは、で効率的に計算できる
2. 目的の値はで求まる

抽象クラスの実装

このアルゴリズムはストラテジーパターンで実装する。具体的には、を計算する関数オブジェクト受け取り、それに応じた部分和を求める。

#include <vector>
#include <functional>

template <typename T>
class func_sum
{
private:
    std::vector<std::vector<T>> sum_table;
protected:
    virtual T acc_sum(int i, int j);
public:
    func_sum(std::function<T(int, int)> func, int nx, int ny) ;
    ~func_sum();
    //目的の区間挿話を求める関数
    T calc_section_sum(int px, int qx, int py, int qy);
};

template <typename T>
T func_sum<T>::acc_sum(int i, int j) {
    return sum_table[i][j];
};

template <typename T>
func_sum<T>::~func_sum() {}

template <typename T>
func_sum<T>::func_sum(std::function<T(int, int)> func, int nx, int ny) {
    sum_table = std::vector<std::vector<T>>(nx + 1, std::vector<T>(ny + 1, T()));
    //j=0の場合の総和の計算
    for (int i = 1; i <= nx; i++) {
        sum_table[i][0] += func(i - 1, 0);
    }
    //i=0の場合の総和の計算
    for (int j = 1; j <= ny; j++) {
        sum_table[0][j] += func(0, j - 1);
    }
    //総和の計算
    for (int i = 1; i <= nx; i++) {
        for (int j = 1; j <= ny; j++) {
            sum_table[i][j] += func(i - 1, j - 1) +
                sum_table[i - 1][j] +
                sum_table[i][j - 1] -
                sum_table[i - 1][j - 1];

        }
    }
}

template <typename T>
T func_sum<T>::calc_section_sum(int px, int qx, int py, int qy) {
    return acc_sum(qx, qy) - 
           acc_sum(qx, py - 1) - 
           acc_sum(px - 1, qy) + 
           acc_sum(px - 1, py - 1);
}

具体的な問題を解く実装

の計算

として、上記アルゴリズムに関数オブジェクトを渡す。

int main()
{
    int N, M;
    int px, py, qx, qy;
    std::cin >> N >> M ;
    std::cin >> px >> qx >> py >> qy;
    std::vector<double> x_list(N, 0);
    std::vector<double> y_list(M, 0);
    for (int i = 0; i < N; i++) {
        double x;
        std::cin >> x;
        x_list[i] = x;
    }
    for (int i = 0; i < M; i++) {
        double y;
        std::cin >> y;
        y_list[i] = y;
    }
    auto func = [](int i, int j, 
                   const std::vector<double>& x_list,
                   const std::vector<double>& y_list) { return x_list[i] / y_list[j]; };

    std::function<double(int, int)>  func_curried = std::bind(func, 
                                                              std::placeholders::_1, 
                                                              std::placeholders::_2, 
                                                              x_list, 
                                                              y_list);

    func_sum<double> ds(func_curried, x_list.size(), y_list.size() );
    std::cout<<ds.calc_section_sum(px, qx, py, qy)<<std::endl;
    return 0;
}

ABC106 D：AtCoder Express 2

問題はリンクを参照。この問題はを, のときそのペアの存在する個数を返す関数とすれば2次元累積和で解ける。

#include <iostream>
int main()
{
    int N, M, Q;
    std::cin >> N >> M >> Q;
    std::vector<std::vector<int>> table(N, std::vector<int>(N, 0));
    for (int i = 0; i < M; i++) {
        int L, R;
        std::cin >> L >> R;
        table[L - 1][R - 1]++;
    }

    auto func = [](int i, int j,
        const std::vector<std::vector<int>>& table) { return table[i][j]; };

    std::function<int(int, int)>  func_curried = std::bind(func,
        std::placeholders::_1,
        std::placeholders::_2,
        table);

    func_sum<int> ds(func_curried, N, N);

    for (int i = 0; i < Q; i++) {
        int p, q;
        std::cin >> p >> q;
        std::cout << ds.calc_section_sum(p, q, p, q) << std::endl;
    }
    return 0;
}

しゃくとり法

用途

• 解きたい問題が以下のいずれかの性質を満たす時、条件を満たすの数え上げ、もしくは評価関数の最大化・最小化
  • として、が条件を満たすならばも条件を満たす（以降、これを上位区間条件と呼ぶ）
    • 例えば、整数の数列の中からその和が以上となる連続する部分数列を数え上げる問題
  • として、が条件を満たすならばも条件を満たす（以降、これを部分区間条件と呼ぶ）
    • 例えば、整数の数列の中からその和が以下となる連続する部分数列を数え上げる問題

アルゴリズム

上位区間条件を満たす問題と部分区間条件を満たす問題とで若干異なる（うまいことやれば、共通にできるかもしれない）

• 上位区間条件の場合

1. 条件を満たさない間をインクリメント（右端まで来てインクリメントできない場合はそのまま）
2. 条件を満たしたところで、インクリメントを止め、その時点のを記録（この時点で, を満たすは全て条件を満たす）
3. 条件が満たさなくなるまでをインクリメント
4. この1～3の操作をがインクリメントできなくなり、かつ、条件を満たさなくなるまで繰り返す

• 部分区間条件の場合

1. 条件を満たす間をインクリメント（右端まで来てインクリメントできない場合はそのまま）
2. 条件を満たさなくなったところで、インクリメントを止め、その時点のを記録（この時点で, を満たすは全て条件を満たす）
3. 条件を満たすまでをインクリメント
4. この1～3の操作をがインクリメントできなくなり、かつ、がと同じになるまで繰り返す

抽象クラスの実装

こちらもテンプレートメソッドパターンで実装する。上述の通り上位区間条件と部分区間条件を分けて考えるので、それぞれ抽象クラスを作る。

#include <functional>
 
//ある区間が条件を満たすならば、
//それを包含する上位区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_superset
{
private:
    virtual void move_left_pointer() = 0;//lを右に動かす処理
    virtual void move_right_pointer() = 0;//rを右に動かす処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_superset<T>::next(){
    while (is_cond() == false && is_right_end() == false) {
            move_right_pointer();
    }
    if (is_cond() == false)
        return no_next();
    std::pair<T, T> ret = get_LR();
    move_left_pointer();
    return ret;
}
 
//ある区間が条件を満たすならば、
//それが包含する部分区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_subset
{
private:
    virtual void move_left_pointer() = 0;//lを右にインクリメントする処理
    virtual void move_right_pointer() = 0;//rを右にインクリメントする処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
    virtual bool is_left_eq_right() = 0;//lとrが一致しているかを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_subset<T>::next() {
    std::pair<T, T> ret;
    while (is_right_end() == false && is_cond() == true) {
        move_right_pointer();
    }
    ret = get_LR();
    if (is_left_eq_right() == true) {
        return no_next();
    }
    
    move_left_pointer();
    return ret;
}

具体的な問題を解く実装

POJ 3061：Subsequence

蟻本に例として載っている和が以上になる部分数列を数え上げる問題。上位区間条件型のしゃくとり法で解ける。

#include <iostream>
#include <algorithm>
#include <vector>

class subsequence :
    public abstract_inchworm_superset<int>
{
private: 
    int Left, Right;
    long long S;
    int N;
    long long LRsum;
    std::vector<long long> sequence;
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
public:

    void set_S(long long S);
    subsequence(std::vector<long long>& sequence, int N, int start);
};

subsequence::subsequence(std::vector<long long>& sequence, int N, int start):
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    S(0),
    LRsum(0)
{
}

void subsequence::set_S(long long S) {
    this->S = S;
}

void subsequence::move_left_pointer() {
    LRsum -= sequence[Left];
    Left++;
}

void subsequence::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        Right++;
    }
}

bool subsequence::is_cond() {
    if (LRsum >= S)
        return true;
    return false;
}

bool subsequence::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}

std::pair<int, int> subsequence::no_next() {
    return std::make_pair(-1, -1);
}

std::pair<int, int> subsequence::get_LR() {
    return std::make_pair(Left, Right);
}

int main()
{

    int N;
    std::cin >> N;
    for (int k = 0; k < N; k++) {
        std::vector<long long> sequence;
        int n;
        long long S;
        std::cin >> n >> S;
        for (int i = 0; i < n; i++) {
            long long a;
            std::cin >> a;
            sequence.push_back(a);
        }

        subsequence seq(sequence, sequence.size(), 0);
        seq.set_S(S);
        int min_len = INT_MAX;
        while (1) {
            std::pair<int, int> LR = seq.next();
            if (LR.first == -1) break;
            min_len = std::min(min_len, LR.second - LR.first);
        }
        if (min_len == INT_MAX)
            std::cout << 0 << std::endl;
        else
            std::cout << min_len << std::endl;
    }
    return 0;
}

ARC098 D：Xor Sum 2

部分区間条件型のしゃくとり法で解けると気付くのが難しい問題。気づいてしまえば、普通にしゃくとり法を適用するだけ（部分区間を数え上げる問題なので、何も考えずしゃくとり法を適用してACしちゃったって人も結構いるかも）。なぜ、しゃくとり法で解けるかは解説参照。
ここで注意するのは、関数is_cond()のtrueを返す条件が現在指しているの更に一つ先のにおいて、条件が成立している場合としている。つまり、right pointer（）を動かしても条件成立するなら、実際にright pointerを動かすという処理になっている。なお、この実装では区間をという半区間で扱っているため、一つ先のはが指示しているインデックスの値を足すこと（XORを取ること）に当たる。

#include <vector>
#include <stdio.h>
class xor_sum :
    public abstract_inchworm_subset<int>
{
private:
    int Left, Right;
    
    int N;
    long long LR_XOR;
    long long LRsum;
    std::vector<long long> sequence;
 
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
    virtual bool is_left_eq_right();
public:
 
    xor_sum(std::vector<long long>& sequence, int N, int start);
};
 
xor_sum::xor_sum(std::vector<long long>& sequence, int N, int start) :
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    LR_XOR(0),
    LRsum(0)
{
}
 
 
void xor_sum::move_left_pointer() {
    LRsum -= sequence[Left];
    LR_XOR ^= sequence[Left];
    Left++;
 
}
 
void xor_sum::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        LR_XOR ^= sequence[Right];
        Right++;
    }
}
 
bool xor_sum::is_cond() {
    if (Left >= N) 
        return false;
    //Rightが1つ先に行っても条件を満たすか
    if ((LRsum + sequence[Right]) == (LR_XOR ^ sequence[Right]))
        return true;
    return false;
}
 
bool xor_sum::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}
 
std::pair<int, int> xor_sum::no_next() {
    return std::make_pair(-1, -1);
}
 
std::pair<int, int> xor_sum::get_LR() {
    return std::make_pair(Left, Right);
}
 
bool xor_sum::is_left_eq_right() {
    return (Left == Right);
}
 
 
int main()
{
    int N;
    std::vector<long long> A;
    scanf("%d", &N);
 
    for (int i = 0; i < N; i++) {
        long long a;
        scanf("%lld", &a);
        A.push_back(a);
    }
    xor_sum xs(A, N, 0);
    long long count = 0;
    while (1) {
        std::pair<int, int> LR = xs.next();
        if (LR.first == -1) break;
        count += LR.second - LR.first ;
    }
    printf("%lld\n", count);
    return 0;
}

おわりに

今回は勉強のため、本質的なアルゴリズムを記述した抽象クラスと実際の問題を解く具象クラスを分けましたが、実際のコンテスト時はそんなことせずに、普通に書いた方が速いと思います。

競プロとかに使うアルゴリズム実装メモ（幅優先・深さ優先探索、union-find、最小全域木）

はじめに

以前の記事にて最短経路問題を解くアルゴリズムの実装を書きましたが、今回はその続きとしてグラフアルゴリズムの中でも幅優先探索、深さ優先探索、union-find、最小全域木問題を解くアルゴリズム2種について実装を書いていきます。
例によって、自分の脳内整理メモの側面が強い上、検索すればそれぞれのアルゴリズムについて素晴らしい解説記事が出てくるので、まじめな解説を見たい場合はそちらを当たった方が良いと思います。

幅優先探索

用途

• コストが全て同じグラフ（重みなしグラフ）の最短経路探索。
• 連結成分の列挙。
• 重みなしグラフの最小全域木。

アルゴリズムのポイント

• 隣接する頂点から順々に探索していく。
• 距離を求める場合は、遷移元の距離に+1をする。

実装のポイント

• キューを使う。
• まずスタートの頂点をキューに入れ、キューの先頭から頂点を取り出し、それに隣接する頂点を探索。
• 未探索の全て頂点はキューに追加して、またキューの先頭から頂点を取り出して・・・ということを目的の頂点が見つかるか、キューが空になるまで繰り返していく。
• 既に探索済みの頂点かというフラグを持たせておき、既に見た頂点だった場合はキューに追加しない（下記実装では距離が初期値-1でなければ探索済みと判定）。

実装

以前と同様、メイン関数は省略。

#include <unordered_map>
#include <vector>
#include <queue>

class breadth_first_search
{
private:
    int node_num;
    int abs_time;
    std::unordered_map<int, std::vector<int>> adj_list;
    std::vector<int> dists;//始点からの距離

public:
    breadth_first_search(const int node_num,
        const std::unordered_map
        <int, std::vector<int>>& adj_list);
    void exec_search(int start);
    int get_dist(int v);
};

breadth_first_search::breadth_first_search(const int node_num,
    const std::unordered_map
    <int, std::vector<int>>& adj_list) :
    node_num(node_num),
    adj_list(adj_list),
    dists(node_num + 1, -1)
{
}

void breadth_first_search::exec_search(int start) {
    std::queue<int> v_queue;
    v_queue.push(start);
    dists[start] = 0;
    int d = 0;
    while (1) {
        if (v_queue.empty() == true)
            break;
        int v = v_queue.front();
        v_queue.pop();
       
        for (int &adj : adj_list[v]) {
            if (dists[adj] == -1) {
                dists[adj] = dists[v] + 1;
                v_queue.push(adj);
            }
        }
    }
}

int breadth_first_search::get_dist(int v) {
    return dists[v];
}

AOJのALDS1_11_Cが通ることを確認。

深さ優先探索（再帰呼び出し版）

用途

• コストが全て同じグラフ（重みなしグラフ）において、与えられたコスト以内でたどり着ける頂点の列挙。
• 現実における迷路探索（現実だと幅優先探索は分身かワープでもできないと無理）。
• 深さを時系列的な情報と関連づけて探索するケース。

（基本、幅優先探索と同じことができるが、それぞれ解きたい問題に応じて使い分ける）

アルゴリズムのポイント

• 行けるところまでどんどん深くまで探索していく。
• これ以上探索できないというところまで探索したら根まで戻り、更に別なルートの探索を進める。

実装のポイント

• 幅優先探索と同じノリでスタックを使うパターンと再帰呼び出しを使うパターンがある（今回の実装は再帰版）。
• 再帰呼び出しで行う場合、分岐して最深部まで探索した後、分岐点まで戻って来たときに既に探索済みか否かという情報をクリアしたいケースの実装が楽になる（上記の深さと時系列情報を関連付ける場合、分岐点に戻る=時系列的にも戻る=探索済み情報をクリアする、とするべき問題もある。これとか）。
• 今回の実装では、AOJの問題に対応するべく、探索で移動したらカウントアップ（後戻り時もカウントアップ）するabs_time変数を用意し、exec_search_sub関数に入った時と出る時にカウントアップさせている。また、find_abs_timesには各頂点を見つけたときのabs_timeを保持しておく。
• find_abs_timeが初期値（-1）でなければ、確定頂点として探索しない。
• このタイムスタンプのつけ方は、問題に応じて適宜カウントアップするタイミング、回数、場合によっては初期値にクリアする等の変更が必要。

実装

#include <unordered_map>
#include <vector>
class depth_first_search
{
private:
    int node_num;
    int abs_time;
    std::unordered_map<int, std::vector<int>> adj_list;
    std::vector<int> depths;//ルートからの深さ（距離ではない）
    std::vector<int> find_abs_times;//頂点に到達した時刻
    std::vector<int> finish_abs_times;//サーチが完了した時刻
    
    void exec_search_sub(int start, int depth);

public:
    depth_first_search(const int node_num, 
                       const std::unordered_map
                        <int, std::vector<int>>& adj_list);
    
    void exec_search(int start);
    int get_depth(int v);
    int get_find_abs_time(int v);
    int get_finish_abs_time(int v);
};

depth_first_search::depth_first_search(const int node_num, 
                                       const std::unordered_map
                                        <int, std::vector<int>>& adj_list):
    node_num(node_num),
    abs_time(0),
    adj_list(adj_list),
    depths(node_num + 1, -1),
    find_abs_times(node_num + 1, -1),
    finish_abs_times(node_num + 1, -1)
{
}

void depth_first_search::exec_search_sub(int start, int depth) {
    
    if (find_abs_times[start] != -1) return;
    depths[start] = depth;
    abs_time++;//探索進みのカウントアップ
    find_abs_times[start] = abs_time;
    for (auto &v : adj_list[start]) {
        exec_search_sub(v, depth + 1);         
    }
    abs_time++;//探索後戻りのカウントアップ
    finish_abs_times[start] = abs_time;
}

void depth_first_search::exec_search(int start) {
    for (int i = start; i < node_num + 1; i++) {
        exec_search_sub(i, 0);
    }
}

int depth_first_search::get_depth(int v) {
    return depths[v];
}
int depth_first_search::get_find_abs_time(int v) {
    return find_abs_times[v];
}

int depth_first_search::get_finish_abs_time(int v) {
    return finish_abs_times[v];
}

AOJのALDS1_11_Bが通ることを確認。

union-find

用途

• 集合を扱うデータ構造。
• 値xとyが同じ集合に属するか否かの判定や2つの集合の結合が高速で実行可能。
• 後に説明する最小全域木問題を解くアルゴリズムであるクラスカル法でも使われる。

アルゴリズムのポイント

• 各集合、木構造でデータを保持する。
• 値xとyが同じ集合か否かを判断するには（以下、find操作）、木の根を見て同じだったら同じ集合、異なれば別の集合に属すると判断すればよい。
• そのため、なるべく木の高さが低くなるようにデータを保持した方が、高速に判断できるようになる。
• このことから、2つの集合を結合する操作（以下、unite操作）では木の高さが低い方の根の親を高い方の根になるように連結させる。こうすれば、連結後の木の高さは元の集合のうち高い方の高さと同じになるので、高さは増えない（両方の木の高さが同じ場合は、木の高さが1だけ増える）。

実装のポイント

• 木は親ノードの値を保持する1次元配列と所属する木の高さを保持する1次元配列で表現する。
• こうすれば、unite操作時、低い方の木の根の親を書きかえるように親を保持する配列を操作だけで良くなる。
• find操作時に全てのノードを直接親ノードとつなぐことで木の縮約を行う。
• 木の縮約時に木の高さの更新を行わないので実際の木の高さと配列に保持している高さが合わなくなるが、あまり気にしなくてよい（重要なのはノードの親子関係であり、それさえ合っていればunion-findは問題なく動作する。木の高さの配列に入っている値がおかしいとunite操作時に、接続する関係が逆になり木が高くなる恐れがあるが、それよりもfind操作で縮約する方がメリットが大きい）。

実装

#include <vector>

class union_find
{
private:
    std::vector<int> parent;
    std::vector<int> rank;
    int node_num;
public:
    union_find(int node_num);
    void unite(int x, int y);
    int find(int x);
};

union_find::union_find(int node_num):
    parent(node_num + 1, 0),
    rank(node_num + 1, 0)
{
    for (int i = 0; i <= node_num; i++) {
        parent[i] = i;
    }
}

void union_find::unite(int x, int y) {
    int x_root = find(x);
    int y_root = find(y);
    //既に同じ根なら、同じ所属なので何もしない
    if (x_root == y_root) return;
    //ランクが大きい方に小さい方をくっつける
    //結合後のランクは大きい方になる。
    if (rank[x_root] > rank[y_root]) {
        parent[y_root] = x_root;
    }
    else {
        parent[x_root] = y_root;
        if(rank[x_root] == rank[y_root]){
            //ランクが等しい場合は、必ずランクが1上がる。
            rank[y_root]++;
        }
    }
}

int union_find::find(int x) {
    if (x != parent[x]) {
        // findの結果を各parentに入れることで経路を圧縮。
        //ただし、変更コストが大きいのでランクは変えてない。
        parent[x] = find(parent[x]);
    }
    return parent[x];
}

後に説明するクラスカル法の中で一緒に確認。

プリム法

用途

• 最小全域木問題を解く。
• 最小全域木問題自体は与えられた無向グラフから、ループがなくなるように枝を取り除き、辺のコストの合計が最小になる全域木を求める問題。
• 例えば、全ての街間で行き来できるように最小コストで道を作りたい的な問題に使える（完全グラフを構成して最小全域木を求める等*1）。

アルゴリズムのポイント

• ダイクストラ法と似ており、まず探索の初期位置となる頂点を決め、そこから出ている辺を探索中リストに追加する。
• 探索中リストの中で最小のコストの辺を取り出す。この辺は最小全域木の構成要素として確定する。
• 上記で取り出した辺につながっている頂点を次の探索対象の頂点とし、その頂点を確定頂点とする。また同様にそこから出ている辺を探索中リストに追加する（ただし、確定頂点へつながる辺は登録しない）。
• これらの操作を全ての頂点が確定するまで繰り返す。

実装のポイント

• ダイクストラ法と同様に優先度付きキューを使う。
• 各頂点に対して現在どのコストで優先度付きキューに登録しているかを保持する外部変数keyを用意する。
• 優先度付きキューから最小コストの辺（と接続先の頂点）を取り出し、この辺を最小全域木の構成要素として確定させる。
• 同時に接続先の頂点も確定頂点として記録し、次はその頂点につながっている辺を見ていく。
• 辺のコストを見て、接続頂点のkeyの値より小さいコストの辺だったらkeyの値を更新し、優先度付きキューに追加する。もし、大きいコストの辺だったら、その辺は最小全域木の構成要素にはなり得ないので、優先度付きキューには追加しない。これ結構大事。これをやらないと、どんどん優先度つきキューに登録されていくので、サイズがばかでかくなり、確定ノードかの判定回数が大幅に増えるので遅くなる可能性がある。

実装

#include <vector>
#include <unordered_map>
#include <algorithm>
#include <vector>
#include <unordered_map>
#include <queue>
#include <functional>
#include <climits>


class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }

class prim
{
private:
    //unordered_mapを使って隣接リストを構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    std::vector<int> pred;
    std::vector<long long> key;
    int node_num;
    std::vector<bool> is_done;//確定ノードか否かを保持
    long long MST_cost;
public:
    prim(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    long long get_MST_cost() const;
    //std::unordered_map<int, std::vector<edge>> calc_MST();
    std::unordered_map<int, std::vector<edge>> get_MST_list();
    void calc_MST();
};

prim::prim(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list):
    adj_list(adj_list),
    node_num(node_num),
    is_done(node_num + 1, false),
    pred(node_num + 1, INT_MAX),
    key(node_num + 1, LLONG_MAX),
    MST_cost(0)
{

}

void prim::calc_MST() {
    MST_cost = 0;
    key = std::vector<long long>(node_num + 1, LLONG_MAX);
    pred = std::vector<int>(node_num + 1, INT_MAX);
    //優先度付きキューpairで頂点番号と最短距離を保持
    //firstの要素で比較されるので、firstが距離、secondを遷移先の頂点とする
    std::priority_queue<std::pair<long long, int>,
        std::vector<std::pair<long long, int>>,
        std::greater<std::pair<long long, int>>> p_queue;
    p_queue.push(std::make_pair(0, adj_list.begin()->first));
    key[adj_list.begin()->first] = 0;
    pred[adj_list.begin()->first] = -1;

    while (1) {
        if (p_queue.empty() == true) break; //キューが空なら終了
        std::pair<long long, int> cost_node = p_queue.top();
        p_queue.pop();
        if (is_done[cost_node.second] == false) {
            MST_cost += cost_node.first;
        }
        
        int from_node = cost_node.second;
        is_done[from_node] = true; //キューから取り出した頂点を確定頂点とする
        for (auto &it : adj_list[from_node]) {
            int adj_node = it.get_to();
            if (is_done[adj_node] == false) {
                long long adj_cost = it.get_cost();
                if (adj_cost < key[adj_node]) {
                    p_queue.push(std::make_pair(adj_cost,adj_node));
                    key[adj_node] = adj_cost;
                    pred[adj_node] = from_node;
                }
            }
        }
    }
}


long long prim::get_MST_cost() const {
    return  MST_cost;
}

std::unordered_map<int, std::vector<edge>> prim::get_MST_list() {
    std::unordered_map<int, std::vector<edge>> MST_list;
    for (int i = 0; i < node_num + 1; i++) {
        if (pred[i] >= 0 && pred[i] < node_num + 1) {
            MST_list[i].push_back(edge(pred[i], key[i]));
            MST_list[pred[i]].push_back(edge(i, key[i]));
        }
    }
    return MST_list;
}

AOJのGRL_2_Aが通ることを確認。

クラスカル法

用途

• プリム法と同じく最小全域木問題を解く。

アルゴリズムのポイント

• 最初は全ての頂点間が未連結の状態（辺がない状態）からスタート。
• 辺をコストでソート後、コストが小さい辺から順に全域木の構成要素に追加していく。ただし、ループができる場合はその辺は追加しない。
• 辺を追加してループができるか否かは、その辺で結んでいる両頂点が既に連結関係にあるか否かで判断可能（既に連結関係にあれば、その辺を追加すると必ずループができる）。

実装のポイント

• 連結されている頂点同士が同じ集合になるようunion-findで管理し、頂点間の連結関係を効率よく判定および操作できるようにする。
• 実装では隣接リストをコンストラクタで渡しているが、辺を管理する変数を用意する方が楽なので、隣接リストから辺のリストに変換させている（最初から辺のリストを構築してもよい）。
• 辺はstd::pair<int, std::pair<int, int>>型で(コスト, (頂点, 頂点))という変数にし、リストにしたとき第1要素のコストでソートされるようにする。

実装

上述のunion-findを使用する。

#include <algorithm>
#include <vector>
#include <unordered_map>
#include <queue>
#include <functional>
#include <climits>

class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }

class kruskal
{
private:
    //unordered_mapを使って隣接リストを構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    //最小全域木の隣接リスト
    std::unordered_map<int, std::vector<edge>> MST_list;
    int node_num;
    long long MST_cost;
public:
    kruskal(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    long long get_MST_cost() const;
    std::unordered_map<int, std::vector<edge>> get_MST_list();
    void calc_MST();
};


kruskal::kruskal(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list):
    adj_list(adj_list),
    node_num(node_num),
    MST_cost(0)
{
}

void kruskal::calc_MST() {
    //辺のリスト。firstがコスト、secondがつながっている両頂点
    std::vector<std::pair<long long, std::pair<int, int>>> edge_list;

    //隣接リストから辺のリストへ変換
    for (auto &v : adj_list) {
        for (auto &e : v.second) {
            edge_list.push_back(std::make_pair(e.get_cost(), 
                        std::make_pair(v.first, e.get_to())));
        }
    }

    std::sort(edge_list.begin(), edge_list.end());

    union_find uf(node_num + 1);
    MST_cost = 0;
    MST_list.clear();
    for (auto &ed : edge_list) {
        int v1 = ed.second.first;
        int v2 = ed.second.second;
        long long cost = ed.first;
        //同じ根を持つなら、既に両頂点は同じ木の中
        //=辺を追加するとループができるのでなにもしない。
        if (uf.find(v1) != uf.find(v2)) {
            //根が異なる場合、両頂点間に辺を追加
            uf.unite(v1, v2);
            MST_cost += cost;
            MST_list[v1].push_back(edge(v2, cost));
            MST_list[v2].push_back(edge(v1, cost));
        }
    }
}

long long kruskal::get_MST_cost() const {
    return MST_cost;
}

std::unordered_map<int, std::vector<edge>> kruskal::get_MST_list() {
    return MST_list;
}

プリム法同様、AOJのGRL_2_Aで確認。

おわりに

前回、今回とグラフアルゴリズムについて見てきましたが、これら以外にもネットワークフロー問題やマッチング問題など様々な問題とそれらを解くアルゴリズムがあります。また、前回紹介した最短経路問題や今回扱った最小全域木問題、他にもマッチング問題など多くの問題は離散凸解析の枠組みで抽象化され、その枠組みの中で最適化アルゴリズムが提案されています。この離散凸解析についてはまだ全然理解していませんが、いつか勉強して記事を書きたいところです。
あと、このシリーズはもう少し続きます。次回は尺取り法や累積和法などを書いていこうと思います。

*1:競プロだと、完全グラフを考えると大抵の場合TLEになるので問題の性質から元のグラフをどのように構築するかというところが肝になる。

競プロ関係の雑メモ 2018/9/9

はじめに

以前と同様、atcoderの400～500点問題穴埋め時のメモに加え、ABC109参加のメモです。

ARC088-D Wide Flip

区間[l - 1, r]を反転させた後、区間[l, r]を反転させると、l - 1番目の1文字のみ反転させることができる。同様に区間[l, r + 1]を反転させた後、区間[l, r]を反転させると、r + 1番目の1文字のみ反転させることができる。ということは、選べる区間の最小の長さKが与えられたとして、左端からN-K個、右端からもN-K個、合計2(N-K)個の要素は自由に0,1を選べる、ということになる。
じゃあ、左端もしくは右端から何番目の要素まで自由に変えられるようにする必要があるかと考えると、真ん中の要素がいくつ連続で同じ文字かを見ればよいという結論に至る。つまり、中央N/2番目の要素から左右にC個が全て同じ要素だとすると、そのC個は変化させる必要がなく、左端から(N-C)/2番目までの要素もしくは右端から(N-C)/2番目までの要素のどれかが自由に0, 1を選べるようKを定めればよい。よって、N - K = (N - C)/2という方程式が出てくるので、これをKについて解いてK = N - (N - C)/2となる。
Cについては、普通に真ん中から見て同じ要素をカウントすればいいのだけど、文字列の長さが偶数か奇数かで話が変わってくるので、少し注意が必要（実際、1度for文回す回数間違えてWAになった）。

ARC064-D An Ordinary Game

まず個数に着目。1個ずつ文字を削除していくゲームなので、負けるときの文字数の法則性を見つけ出せれば、後は逆算して元の文字数から勝者がわかると考えた。じゃあ、負けるときの文字列とは？ってところについては、2パターン考えられる。ひとつはabのような文字数が偶数での負け、もうひとつはabaのような文字数が奇数での負け（この辺り、大体の予想は立ったが、じゃあ本当？っていうところが解説見るまで分からなかった）。では、どういうときに偶数負けパターンになるか？奇数負けパターンになるか？
ポイントはゲームの性質上、先頭の文字と後尾の文字は必ず最後まで残るということ。このことから、先頭の文字と後尾の文字が同じであれば奇数負けパターン、異なるのであれば偶数負けパターンになる。なぜなら、例えばa**b***a*b（*はaとb以外の文字）のように先頭と後尾が異なるなら、必ず*は取り除くことができ、最終的な負けパターンとしてはabかababになる（そうなるように、両プレイヤーがコントロール可能）。どちらになるかは、一意に定まらないが偶数の文字数になることはかわらない。これと同じような考えで、a*b**a*b**a等のように、先頭と後尾が同じ場合は、abaかa*aかababaかのいずれかが最終形になるので必ず奇数の文字数になる。ここまで来れば、上述の通り逆算して元の文字列数からどちらが勝つかはわかる。
どうもゲーム系は苦手意識があったが（未だに「互いが最適に行動を取る」ってことのイメージがつかめてない）、最終的に詰み（負け）となる形を考えて、そこから逆算していくっていう方法が有効そう。
この問題では一発でACを取れたが、根拠の薄い予想がたまたま当たったっていうパターンなので、あまりよろしくない*1。

ARC067-D Walk and Teleport

まず、いくつか移動パターンを考える。

1. 西から東へ順々に移動するパターン（テレポートしても町は飛び越えない）。
2. テレポートを使って、まだ移動していない町を飛び越え、その飛び越えた町へ行くため東から西へ逆走をするパターン。
3. 全てテレポートで移動するパターン。

この中で2の移動パターンは使う必要がない。なぜなら、最終的には全ての町を訪れる必要があるため、逆走しても順方向に移動しても距離は変わらず、町を飛び越すメリットがないから。また、1のパターンはB<A(X[i + 1] - X[i])であればテレポートを使うし、そうでなければ徒歩で移動する。3のパターンは1の条件で、全ての町間の距離がテレポート使用条件に合うなら自動的に適用される。ということで、min(A(X[i + 1] - X[i]), B)を足しこんでいけばよい。
500点問題の中では明らかに簡単なような・・・

ABC 109参加メモ

2018/9/8開催のABC109に参加。4完達成もD問題の実装に時間がかかりすぎたのがよろしくない。

A問題

A*Bが奇数か偶数か。

B問題

setを使って既出かの判定と、前回の文字列の最後の文字と今見てる文字列の頭文字が同じかという判定を行う。

C問題

Xも数列に組み込んでソートしてx[i + 1] - x[i]の最大公約数を求めたが、解説読むとソートなんかしなくても|X - x[i]|の最大公約数でよかった。

D問題

着目しているマスのコイン数が奇数なら隣のマスにコインを1枚ずらし、着目するマスもそのずらした先に移動。逆に、着目しているマスのコイン数が偶数なら何もせず、着目するマスのみ隣にずらす。この操作を左端からジグザグでもグルグルでも何でもいいから順々に全てのマスに対して行えば、目的が達成される。
その理由としては、まず奇数のマスから奇数のマスへコインを一つずらすと両方偶数になる。その上で、ずらし元が奇数、ずらし先が偶数の場合、操作後はずらし先が奇数になるが、直後に奇数になったずらし先に対して操作を行えば、そのマスも偶数になる。この操作を続けていくと、どこかかでずらし元、ずらし先ともに奇数という個所が出てきて、両方偶数になる。もしそのような個所が出てこない場合、最後のマスまでコインが伝搬され最後のマスのみが奇数になる。従って、最適な最終形としては、コインの総枚数が偶数の場合は全てのマスが偶数になり、コインの総枚数が奇数の場合、1つのマスのみが奇数でそれ以外全手のマスが偶数になる。
ただ、いざ実装ってなると、配列のジグザグ操作またはグルグル操作の実装が面倒だなーと億劫になり、結局、既に操作済みのマスか？という情報を記録する配列を別途用意。ずらせるマスを探し、ずらした先の座標を次のずらし元とする作戦を取る。ACだったが実装に時間がかかりすぎた。

おわりに

最近、競プロ関連チラ裏ブログと化しているので、次回はもう少しまともな内容を書きます。

*1:もう少し正確に言うと、予想した負けパターンは、abのような2文字か、abaのような3文字、もしくは最初からabababのようなFirstプレイヤーが取れず負けるというパターンの3種類を考え、最後のパターンだけ（無駄な）例外処理をする実装でACになった。

競プロ関係の雑メモ 2018/9/1

はじめに

競プロの問題解いてて、解けなかった問題は何故解けなかったのか、何故その発想が出なかったのかを記していく自分用メモです。atcoderの400、500点くらいの問題なら、解けないなりにもボチボチ近いところまでは行けてるっていうパターンが多いので（それ以上の問題は今のところ大体歯が立たない）、あと一歩何が足りなかったのかを探って書いていきます。

ARC074-D 3N Numbers

とりあえず、境界を決めて左側の数列から大きい順にN個、右側の小さい順にN個取ればよく、後はスコアが最大になる境界位置を全探索すればよい。ただし、これだとO(N^2logN)だからTLEになる。じゃあ、どうする？優先度付きキューか？っていうところまでは発想が至った。
ということで、左側数列用優先度付きキューと右側数列用優先度付きキューを用意して、境界位置を左から順にずらし、両優先度付きキューを更新していくってことを行おうとしたところ、右側数列用優先度付きキューから指定要素を取り除く操作が出てきて、その操作には1回O(N)かかるぞ、どうすればいいんだーとなる。
境界位置探索について1回で済まそうと、左側と右側の更新を同時（同じループ内）で行おうと考えたのが間違い。
まず境界位置を左からずらして、左側数列優先度付きキューを更新していく。それが終わった後、今度は境界位置を右からずらして、右側数列の優先度付きキューを更新していく。こうすれば優先度付きキューの操作は、先頭要素取り出しと要素追加だけになるので1回の操作はO(logN)で可能。最後にどの境界位置が最適かを見ればよい。

ARC083-D Restoring Road Network

与えられている最短距離表を各頂点間を結ぶ隣接行列だとみなしてワーシャル・フロイドで解いた後、矛盾がないか探る（矛盾があれば、ワーシャルフロイドで求めた最短距離と最短距離表の値が異なる部分が出てくる）。矛盾がなければ、頂点i,j間を直接結ぶ辺が必要かを考えればよい。つまり、i,j間を直接結ばずともどこかを経由して最短距離となる経路があるかを探す。
っていうところまでは出てきたけど、ここで躓く。i, j間の最短経路で直接結ぶ以外のパスがあるかを求める方法がわからなかった（i, jを直接結ぶパスのみを削除したグラフG'を用意して、これに対してダイクストラ法でi, j間の最短距離を求め、元のグラフで求めた最短距離と変わらないかってことをやろうとしたけどTLE）。
これを求めるには、頂点kを経由して、最短距離が等しくなる経路が存在するかを調べればよい。つまり、i-k間の最短距離+k-j間の最短距離が最短距離表にあるAijと等しくなるkが存在するかみればよい。

ARC076-D Built?

x座標、y座標それぞれでソートして、隣り合う街同士にのみ道を作ることを考えればOK。x座標、y座標ごった煮にして、コストの小さい道から順に追加していけばいいってところまでは考えた。ところが、やってみると各頂点をが全て連結されているのか効率よく判断する方法がわからず、言葉だけは知っていたUnion-Findを使うのかな？って思ったところで、なんとなく敷居が高い気がして断念。
というか、ここまで来れば、まんまクラスカル法。もし、クラスカル法知らなくても、ソートして隣り合う街同士のみ考えればいいってことを鑑みると、考慮する辺数が2Vだけなので、O(VlogV)で解ける最小全域木問題だって事気付けばプリム法で解いても良かった（これ解いた時、プリム法は知ってた）。いずれにしても、問題そのものが最小全域木の典型問題にもかかわらず、ソートとか隣り合う街同士のみ考慮とか色々考えていくうちに、最小全域木問題を解くアルゴリズム使うという発想がすっぽり頭から抜けていたのが反省点。

おわりに

もう少しレベルを上げたいので、時々こういう反省ポエムも書いていきます。

競プロとかに使うアルゴリズム実装メモ（最短経路探索系）

はじめに

ここ１年くらい、ちまちまとatcoder中心に競技プログラミングに参加してたりします。ABCのD問題を解くのがやっとなのにAGCに突撃して爆死するってことを繰り返し続け、未だに緑コーダーです（パフォーマンスも1200前後がやっと）。こりゃまずいということで、もう少しまじめにアルゴリズムの勉強をしようと思い立ったので、今回から数回にわたって基本的なアルゴリズムのメモと実装を書いて、自分の理解を深めていこう思います。
今回はその中でも最短経路探索系（ダイクストラ法、ベルマン・フォード法、ワーシャル・フロイド法）のアルゴリズムについて書いていきます。これらのアルゴリズムについては、ネットや探せばいくらでも解説記事が出てくるので、ここでは主に自分が実装していく中で理解にしにくかった、ハマった点をポイントとして記していきます。そのため、ほぼ自分用のメモ&脳内整理用の記事なので、まじめな解説は他を当たった方が良いかもしれません。
ちなみに、今まで競プロも機械学習もpythonで実装してきましたが、今後本業でC++を使う機会が増えそうなので、今回はリハビリの意味も込めてC++での実装を行っていきます。

ダイクストラ法

用途

• 負のコストがない場合の単一始点最短経路探索。

アルゴリズムのポイント

• 始点に隣接している頂点からの順次最短経路を確定させていくアルゴリズム。
• 最短経路が確定した頂点に隣接する頂点と始点間のコストを計算・更新する処理を全部の頂点の最短経路が確定するまで繰り返す。
• 負のコストがない前提なので、現時点の始点からコスト計算されている頂点の中で最もコストが小さい頂点について、そのコストを最短経路として確定させてよい。
• 具体的には、始点の頂点のコストを0、それ以外の頂点までのコストを∞として初期化後、下記の処理を終了条件を満たすまで繰り返す。
  1. 最短経路が確定確定していない頂点の中からコストが最も小さい頂点を選択。
  2. その選択した頂点を最短経路確定頂点とする。
  3. その選択した頂点に隣接する頂点の現時点のコストが、選択した頂点へのコスト+辺のコストより大きければ、選択した頂点へのコスト+辺のコストに置き換える。
• 最短経路の復元のため、どの頂点から来たかも記憶しておき、最後にバックトラックで辿れるようにする。

実装のポイント

• 優先度付きキューを使って、コスト最小の頂点から取り出して最短経路確定頂点としていく。
• 優先度付きキューの要素はSTLのpairを使う。pairの比較演算子は第1要素で比較するので、pairの第1要素にコスト、第２要素に頂点番号を入れる。
• 同じ頂点でコストが異なる要素が優先度付きキューに残るけど気にしない（プライオリティキューから取り出した時に確定頂点か否かを確認するから、下手に操作せずに残しておく）。
• 隣接頂点をすぐに取り出せるようにするため、グラフの表現にはunordered_mapを使った隣接リストを採用。

実装

メイン関数（呼び出し元）は省略。

#include <algorithm>
#include <vector>
#include <unordered_map>
#include <queue>
#include <functional>
#include <climits>

class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }

class dykstra {
private:
    //unordered_mapを使って隣接を構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    int node_num;
    int start_point;
    std::vector<int> pred;//経路（遷移元）を保持
    std::vector<long long> costs;//距離を保持
    std::vector<bool> is_done;//確定ノードか否かを保持

public:
    dykstra(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    void calc_min_cost(int start);
    long long get_cost(int end) const;
    std::vector<int> get_min_path(int end) const;
};

dykstra::dykstra(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list) :
    adj_list(adj_list),
    start_point(0),
    node_num(node_num),
    pred(node_num + 1, INT_MAX),
    costs(node_num + 1, LLONG_MAX),
    is_done(node_num + 1, false)
{}

void dykstra::calc_min_cost(int start) {
    this->start_point = start;
    this->pred = std::vector<int>(node_num + 1, INT_MAX);
    this->costs = std::vector<long long>(node_num + 1, LLONG_MAX);
    this->is_done = std::vector<bool>(node_num + 1, false);

    //優先度付きキューpairで頂点番号と最短距離を保持
    //firstの要素で比較されるので、firstが距離、secondを頂点番号とする
    std::priority_queue<std::pair<long long, int>,
        std::vector<std::pair<long long, int>>, std::greater<std::pair<long long, int>>> p_queue;
    p_queue.push(std::make_pair(0, start));
    pred[start] = -1;
    costs[start] = 0;

    while (1) {
        if (p_queue.empty() == true) break; //キューが空なら終了
        std::pair<long long, int> cost_node = p_queue.top();
        p_queue.pop();
        is_done[cost_node.second] = true; //キューから取り出した頂点を確定頂点とする
        for (auto &it : adj_list[cost_node.second]) {
            int adj_node = it.get_to();
            if (is_done[adj_node] == false) {
                long long adj_cost = it.get_cost() + cost_node.first;
                if (adj_cost < costs[adj_node]) {//計算された隣接頂点のコストが現在のコストより小さいなら
                    costs[adj_node] = adj_cost; //隣接頂点のコストを更新
                    pred[adj_node] = cost_node.second;//隣接頂点の遷移元を更新
                    p_queue.push(std::make_pair(adj_cost, adj_node));//キューに隣接頂点の情報を突っ込む
                }
            }
        }
    }
}

long long dykstra::get_cost(int end) const {
    return costs[end];
}

std::vector<int> dykstra::get_min_path(int end) const {
    int node = end;
    std::vector<int> vec;
    vec.push_back(node);
    //終点から始点までの経路をたどる
    while (1) {
        //始点から辿れない頂点or
        //始点であれば終了
        if (pred[node] >= INT_MAX || pred[node] == -1) break;
        node = pred[node];
        vec.push_back(node);
    }
    std::reverse(vec.begin(), vec.end());
    return vec;
}

AOJのALDS1_12_Cが通ることを確認。

ベルマン・フォード法

用途

• 負のコストがある場合の単一始点最短経路探索および負のサイクル検出。

アルゴリズムのポイント

• 全ての辺のコストを見て、各頂点までの最短経路（コスト）を更新するという処理を所定の回数繰り返すアルゴリズム。
• 具体的には、始点の頂点のコストを0、それ以外の頂点までのコストを∞として初期化後、下記の処理を所定の回数繰り返す（回数については後述）。
  1. 各辺に対して、接続元の頂点のコスト+辺のコストが接続先の頂点のコストより小さければ、接続先頂点のコストを接続元の頂点のコスト+辺のコストに置き換える。
  2. この処理を全ての辺に対して行う。

• 負のサイクルがあれば、始点からの最短経路は全ての頂点に対して求まらない（いくらでもコストを小さくできるので）。
• 頂点数をVとして、負のサイクルがグラフ上にない場合、始点から各頂点までの最短経路の経由頂点数は高々V-1個（始点は除く）。
• これらのことから負のサイクルがなければ、V-1回の繰り返しで全頂点までの最短経路が求まる。逆にV回目に更新が発生するのならばグラフに負のサイクルがあると判定可能。
• 更に、V～V*2回目の繰り返しにおいて、最短経路が更新される頂点は負のサイクルに含まれると判断される。
• 加えて、負のサイクルに含まれていると確定している頂点から辿れる頂点についても、負のサイクルに含まれる頂点。
• 上記のことより、特定終点までの最短経路と負のサイクルの有無を求める場合は、全頂点の最短経路更新および負のサイクル情報の伝搬処理をV*2回繰り返す（グラフ全体に負のサイクルがあるか否かを求めるだけの場合はV回の繰り返しでよい）。

実装のポイント

• for文使って、V*2回のループと辺数E回のループで素直にアルゴリズムを実装すればよい。
• ダイクストラ法みたいに優先度付きキューは不要。
• 今回は上記ダイクストラ法に合わせて、グラフの表現に隣接リストを使ったが、辺（fromとtoとcostを持つクラス）のリストの方が楽かも。
• 負のサイクルが見つかった後もコスト配列の更新を怠ってはいけない（自戒）。

実装

#include <algorithm>
#include <vector>
#include <unordered_map>
#include <functional>
#include <climits>

class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }


class bellman_ford {
private:
    //unordered_mapを使って隣接グラフを構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    int node_num;
    int start_point;
    std::vector<int> pred;//経路（遷移元）を保持
    std::vector<long long> costs;//距離を保持
    bool is_negative_graph;//グラフ内に負のサイクルをもつか否か
    std::vector<bool> is_negative_pass;//経路上に負のサイクルを持つか否か

public:
    bellman_ford(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    void calc_min_cost(int start);
    long long get_cost(int end) const;
    bool get_is_negative_graph() const;
    std::vector<int> get_min_path(int end) const;
};

bellman_ford::bellman_ford(const int node_num,
    const std::unordered_map<int, std::vector<edge>>& adj_list) :
    adj_list(adj_list),
    start_point(0),
    node_num(node_num),
    pred(node_num + 1, INT_MAX),
    costs(node_num + 1, LLONG_MAX),
    is_negative_graph(false),
    is_negative_pass(node_num + 1, false)
{}

void bellman_ford::calc_min_cost(int start) {
    this->start_point = start;
    this->pred = std::vector<int>(node_num + 1, INT_MAX);
    this->costs = std::vector<long long>(node_num + 1, LLONG_MAX);
    this->is_negative_pass = std::vector<bool>(node_num + 1, false);

    pred[start] = -1;
    costs[start] = 0;
    for (int i = 0; i < 2 * node_num; i++) {
        for (auto &node : adj_list) {
            if (costs[node.first] == LLONG_MAX) continue;
            for (auto &adj : adj_list[node.first]) {
                int adj_node = adj.get_to();
                long long adj_cost = adj.get_cost() + costs[node.first];
                if (adj_cost < costs[adj_node]) {//計算された隣接頂点のコストが現在のコストより小さいなら
                    costs[adj_node] = adj_cost;//隣接頂点のコストを更新
                    pred[adj_node] = node.first;//隣接頂点の遷移元を更新
                    if (i >= node_num - 1) {
                        //頂点数回以上繰り返しても、
                        //まだ更新が発生するなら、その頂点への経路には負のサイクルあり
                        is_negative_pass[adj_node] = true;
                        is_negative_graph = true;
                    }
                }
                if (is_negative_pass[node.first] == true) {//負のサイクル情報伝搬
                                                           //経路に負のサイクルを含む頂点に連結されているなら、
                                                           //その頂点への経路も負のサイクルを含む
                    is_negative_pass[adj_node] = true;
                }
            }
        }
    }
}

bool bellman_ford::get_is_negative_graph() const {
    return is_negative_graph;
}

long long bellman_ford::get_cost(int end) const {
    if (is_negative_pass[end] == false)
        return costs[end];
    else
        return LLONG_MIN;
}

std::vector<int> bellman_ford::get_min_path(int end) const {
    int node = end;
    std::vector<int> vec;
    vec.push_back(node);
    //終点から始点までの経路をたどる
    while (1) {
        //負のサイクルを含む頂点or
        //始点から辿れない頂点or
        //始点であれば終了
        if (is_negative_pass[node] == true ||
            pred[node] >= INT_MAX || pred[node] == -1) break;
        node = pred[node];
        vec.push_back(node);
    }
    std::reverse(vec.begin(), vec.end());
    return vec;
}

AOJのGRL_1_Bが通ることを確認。加えて、グラフそのものに負のサイクルを含むか否かではなく、目的の経路に負サイクルがあるかという判定が正しく動作するか検証するために、atcoderのABC061Dでも確認。

ワーシャル・フロイド法

用途

• 全頂点ペアの最短経路探索。
• 負のコストがあっても良く、負のサイクル検出も可能。

アルゴリズムのポイント

• 頂点iから頂点jへ頂点0～k-1のみを経由する最短経路が既に得られているとして、それに頂点kを追加して最短経路を更新していくという動的計画法の一種。
• もし頂点kを追加してi→j間の最短経路更新が発生するということは、更新後のi→jの最短経路はi→kの最短経路+k→jの最短経路。
• 従って、i→kの最短経路のコスト+k→jの最短経路のコストが、既に得られているi→jへの最短経路のコストよりも小さければ更新。
• 具体的には頂点iから頂点jへの最短経路のコストをcost[i][j]とし、隣接行列で初期化した後、以下の処理を行う。
  1. k =0～Vに対し、以下の処理を繰り返す。
  2. 全頂点のペアi, jについて、既に得られている0～k-1のみを経由する最短経路cost[i][j]よりもcost[i][k] + cost[k][j]が小さければ、頂点kを追加することで経路が更新されるということなので、cost[i][j]をcost[i][k] + cost[k][j]に置き換える。
• 上記処理完了後、costの対角成分（自分から自分への距離）が負ならば、そのグラフには負のサイクルがあり、その頂点は負のサイクルに含まれる。
• 対角成分が負でなくても、負のサイクルに含まれる頂点から辿れる頂点も負のサイクルに含まれる。
• 上記の方法で、各頂点が負のサイクルに含まれるか判断可能なので、頂点iから頂点jへの経路に負のサイクルを含むか否かは、終点jが負のサイクルに含まれるか否かで判断可能（この辺り、あまり文献が見つからず自分で考えたので少し怪しい）。やっぱり間違えてる（というより処理が足りない）。以下のようなグラフで、4→5は負のサイクルを含まないのに、負のサイクルありと判定されてしまう。ちゃんと伝搬元の負のサイクルに含まれる頂点が始点iから辿れるか？も判断しないといけない（時間がないのでコードの修正は別途）。

• ダイクストラ法などと同様、経路復元のためどの頂点から来たかについても、cost[i][j]が更新されたときに合わせて記憶する（どう記憶するかは、少しややこしいのでソースコード内のコメント参照）。

実装のポイント

• k=0～V、i=0～V、j=0～Vの3重ループを回す。
• アルゴリズムの性質上、グラフの表現には隣接行列を採用。
• 隣接行列の未接続の表現にLLONG_MAXを使っているので、呼び出し元は注意が必要。

実装

#include <algorithm>
#include <vector>
#include <functional>
#include <climits>
 
class warshall_floyd
{
private:
    int node_num;
    bool is_negative_graph;
    std::vector<std::vector<long long>> cost_matrix;//隣接行列を構築
    std::vector<std::vector<int>> pred;//経路（遷移元）を保持
   
public:
    warshall_floyd(const int node_num, const std::vector<std::vector<long long>> &adj_matrix);
    bool get_is_negative_graph() const;
    bool get_is_negative_pass(int start, int end) const;
    long long get_cost(int start, int end) const;
    std::vector<int> get_min_path(int start, int end) const;
};
 
warshall_floyd::warshall_floyd(const int node_num, 
        const std::vector<std::vector<long long>> &adj_matrix):
    cost_matrix(adj_matrix),
    node_num(node_num),
    is_negative_graph(false),
    pred(node_num + 1, std::vector<int>(node_num + 1, INT_MAX))
{
    //predの初期化
    for (int i = 0; i <= node_num; i++) {
        for (int j = 0; j <= node_num; j++) {
            if (i == j) pred[i][j] = -1;
            else if (adj_matrix[i][j] < LLONG_MAX) {
                pred[i][j] = i;
            }
        }
    }
    //ワーシャルフロイドの更新式
    for (int k = 0; k <= node_num; k++) {
        for (int i = 0; i <= node_num; i++) {
            for (int j = 0; j <= node_num; j++) {
                if (cost_matrix[i][k] < LLONG_MAX && cost_matrix[k][j] < LLONG_MAX) {
                    if (cost_matrix[i][k] + cost_matrix[k][j] < cost_matrix[i][j]) {
                        cost_matrix[i][j] = cost_matrix[i][k] + cost_matrix[k][j];
                        //経路復元用
                        //pred[i][j]にはi→jの最短経路におけるjの1つ前の頂点が入る。
                        //pred[k][j]にはk→jの最短経路におけるjの1つ前の頂点が入っている。
                        //この処理を通るということは最短経路がkを使った経路に変わる
                        //ということでありi→jへの経路の中にk→jの最短経路を含むことになる。
                        //従って新しいi→jにおけるjの１つ前の頂点は
                        //k→jにおけるjの１つ前の頂点、すなわちpred[k][j]となる。
                        pred[i][j] = pred[k][j];
                    }
                }
            }
        }
    }
 
    //負のサイクル検出
    //[注意!]間違っているので修正必要!（アルゴリズムのポイントを参照）
    for (int i = 0; i <= node_num; i++) {
        //自分への距離が負になるならその頂点は負のサイクルに含まれる
        if (cost_matrix[i][i] < 0) {
            for (int j = 0; j <= node_num; j++) {
                //負のサイクルに含まれる頂点から辿れる頂点も負のサイクルになる
                if (i == j) continue;
                if (cost_matrix[i][j] < LLONG_MAX && cost_matrix[j][j] >= 0)
                    cost_matrix[j][j] = -1;
            }
            is_negative_graph = true;
            //break;
        }
    }
}
 
bool warshall_floyd::get_is_negative_graph() const {
    return is_negative_graph;
}
 
bool warshall_floyd::get_is_negative_pass(int start, int end) const {
    //startからendにたどり着けない場合はfalseを返す
    if (cost_matrix[start][end] < LLONG_MAX) {
        //startからend辿りつける場合、
        //endが負のサイクルに含まれるならtrueを返す
        if (cost_matrix[end][end] < 0) return true;
    }
    return false;
}
 
 
long long warshall_floyd::get_cost(int start, int end) const {
    return cost_matrix[start][end];
}
 
std::vector<int> warshall_floyd::get_min_path(int start, int end) const {
    int node = end;
    std::vector<int> vec;
    vec.push_back(node);
    //終点から始点までの経路をたどる
    while (1) {
        //負のサイクルを含む頂点or
        //始点から辿れない頂点or
        //始点であれば終了
        if (cost_matrix[node][node] < 0 || 
            pred[start][node] >= INT_MAX || 
            pred[start][node] == -1) break;
        node = pred[start][node];
        vec.push_back(node);
    }
    std::reverse(vec.begin(), vec.end());
    return vec;
}

AOJのGRL_1_Cが通ることを確認。後、負サイクル検出ロジックの検証のため、上記同様atcoderのABC061Dも確認（ワーシャル・フロイドで解くには少しきつい問題設定だけどC++なら通る）。

おわりに

今回は、競プロレベルアップに向け最短経路探索系アルゴリズムに関するメモと実装を書きました。
ただ、ここに書いたコードがそのまま競プロに使えるケースはあまりないと思います。例えば、ダイクストラ法を使うにしても、辺のコストがあらかじめ与えられておらず（もしくは与えようとすると膨大なメモリ使用量になる等）キューから取り出して、コスト計算時に初めて辺のコストがわかるなんて問題もあります。そんな問題だとここに書いた実装はそのまま使えないので、カスタマイズする必要があります。
あとそのまま使えるにしても、問題文からはそのことがパッとわからない、または、アルゴリズムを使うのはわかるけどどうやって適用するかを考える必要がある（問題で与えられたコストを負値にする、スタートとゴール両方から計算する、等々）なんてこともあります。その辺りの考察力が勝負になってくる世界なので、強化するべきはそこなんだと思いますが、まずは基本を理解しないと考察力も鍛えられないかなと
ということで、次回以降はその他のグラフアルゴリズム（幅優先、深さ優先、最小全域木とか）、探索系（二分探索、尺取り法とか）、数学系（素数、modとか）、動的計画法辺りを書いていこうと思います。

交互方向乗数法による最適化と画像ノイズ除去への応用

はじめに

これまでの記事で近接勾配法と、それによるスパース解や低ランク解に導く正則化項を付随した最適化問題の解法、そしてその応用を見てきました。正則化項に変数間の絡みがなく各変数が独立に扱える場合は、正則化項のproximal operatorが解析的に求まるため近接勾配法は有効な手段です。ところが、各変数間に絡みがあり分離できない場合、一般的にproximal operatorの計算は容易ではありません（解析的に求まらない）。そこで今回は、変数間に絡みがある正則化項が付いていても最適解を導出することができる交互方向乗数法（Alternating Direction Method of Multipliers : ADMM）のアルゴリズムを見ていきます。また、それを用いた画像のノイズ除去をPythonで実装したので紹介します。

交互方向乗数法（ADMM）とは

いきなりですがADMMのアルゴリズムを記載します。まずADMMの対象となる最適化問題は下記のようなものになります。

この問題を解くために、下記に示す「拡張ラグランジュ関数」というものを定義します。

この拡張ラグランジュ関数に関する詳細は補足の項や参考文献をご参照ください。とりあえず、ここでは通常のラグランジュ関数に等式制約が満たされないことに対する罰則項（第3項）がついたものと考えて問題ありません。この拡張ラグランジュ関数を用いると式（1）を解くADMMアルゴリズムは下記のようになります。

1. 、、を適当な値で初期化する
2. 以下の操作を収束するまで繰り返す

これだけです。アルゴリズムとしては非常に簡単です。特に式（3）と式（4）ではとそれぞれに関する最適化問題を独立で解いているので、元々の最適化問題が関数同士の和が含まれており解くのが難しい場合でも、式（1）の形式に変換できればADMMにより簡単に解くことができます。このことを次の章で見ていきます。が、その前に補足として、上記アルゴリズムがどのように出てきたかを理解するために、双対上昇法、拡張ラグランジュ関数について簡単に説明します（細かいことはかなり省略するので、詳しい説明は参考文献をご参照ください）。

【補足】双対上昇法、拡張ラグランジュ関数について

ADMMは最適化問題の強双対性という性質を利用しています。ここでは、

という最適化問題を考えます。この問題のラグランジュ関数は、

と定義されます。が真凸関数である場合は、となるとが存在すれば、が式（6）の最適解です。
更にが真凸関数であるとき式（7）について下記も成り立ちます。

これは強双対性と呼ばれる性質で、左辺は式（6）と同値の問題です。この強双対性を利用して最適化問題を解くのが双対上昇法や拡張ラグランジュ関数法、そして今回の主題であるADMMです。

双対上昇法

双対上昇法は式（6）の問題を強双対性から双対問題という問題に変換し、そちらを解こうという方針の手法です。まず、式（8）の右辺に問題を当てはめ変形すると、

という問題になります。ここで、は以前説明した共役関数でです。この式（9）の問題は双対問題と呼ばれ、これを解く方法を双対上昇法と言います。更にこれを勾配法に基づいて解く場合は宇沢の方法といい、

という更新式になります。で、この勾配を求めるためには、

を計算する必要があります。結局のところこれは、式（7）にあるラグランジュ関数のに関する最小化です。つまり、主問題であるに関する最適化と双対問題であるに関する最適化を交互に行うアルゴリズムとなります。
この双対上昇法では共役関数の勾配を求めるため、共役関数が微分可能でなければ使えません。しかし、一般的に共役関数は微分可能ではなく、しかも関数が強凸という性質*1を持っていなければアルゴリズムの収束も保証されません*2。そこで、次に説明する拡張ラグランジュの登場です。

拡張ラグランジュ関数法

まず式（6）の問題について、

という式に書きなおします。目的関数の第２項としてを付け加えましたが制約条件よりになるため、式（12）は式（6）と等価です。次に式（12）のラグランジュ関数は

となります。これは拡張ラグランジュ関数と呼ばれる関数で、ラグランジュ関数に対し制約を満たさないことに対する罰則項を加えたものとなっています。実はこの罰則項があることにより主問題は強凸となり、たとえの共役関数が微分不可能だったとしても双対問題が滑らかになるため双対上昇法が適用できます（その理由については、参考文献をご参照ください）。このように、ラグランジュ関数に制約条件を満たさないことに対する罰則項を加えることで問題を解きやすくした上で、 との更新を交互していくのが拡張ラグランジュ関数法です。

ADMM再考

さて、式（1）は式（6）において変数がとに分解可能という特別なケースです。式（1）を拡張ラグランジュ関数法で解こうとした時にとに関する同時最適化が困難なことが多いです。そこで1回の更新でとの同時最適化はあきらめ、個別に最適化していこうというのがADMMです。

変数間に絡みのある正則化項付き最適化問題のADMMによる解法

L1ノルム正則化は変数間が独立しており、proximal operatorが解析的に求まります。L2ノルム正則化（２乗しない）はルートを取る段階で変数間の絡みが出てきますがMoreau decompositionをうまく使ってproximal operator計算することができます（近接勾配法の記事参照）。p=2の重複なしのグループ正則化については、重複がないため各グループ独立したL2ノルムとしてproximal operatorとして計算できます。
しかし、上述のように変数間に絡みがあるとproximal operatorは容易には計算できません。例えば、重複ありのグループ正則化（p=2）の場合、同じ変数が複数のグループにまたがって現れるため、各グループ独立してproximal operatorを計算できなくなってしまいます。そこで、そのような正則化がつく最適化問題をうまく式（1）の形式かつがproximal operatorが計算しやすい形に変形し、ADMMで解くことを考えます。重複ありグループ正則化を例にとって考えてみると、

という問題を解くことになるのですが、これを下記のような問題に変形します。

ここでをうまく定義しがグループ間の重複がないようにします。具体的にはを

などのように、それぞれと左からかけると、グループに所属する要素のみを抽出する行列とします。そして、とすると、式（15）より とは各グループに対応する部分ベクトルを並べたベクトルとなります。こうすることで、正則化項は重複なしのグループ正則化になります。
ということで、この式（15）に対してADMMを適用することを考えます。この中でに関する更新式は、

となります。ここで、 は、グループ正則化項のproximal operatorです（こちらの式（21）参照）。このように、元々の問題が正則化項のproximal operaorを求めるのが困難であり近接勾配法を利用できない場合でも、うまく式変形を行いADMMを適用すれば解けるようになることもあるため、ADMMはかなり強力なアルゴリズムだと言えます。

TV正則化によるノイズ除去への応用

冒頭でも述べたように、今回はこのADMMを応用して画像のノイズ除去を行ってみます。今、入力画像をとして、このからのノイズ除去は、全変動正則化（TV正則化）と呼ばれる正則化項を用いて以下のような定式化が提案されています。

式（17）の第１項は出力画像がなるべく入力画像に近する働きがあります。第２項がTV正則化項で下記のような式となります。

ここではピクセル位置を表します*3。基本的に画像は滑らかな性質を有しているはずなので、隣り合うピクセルの差は小さいはずです。そのため、このTV正則化は隣り合うピクセルの差が大きい場合の罰則項となり、画像を滑らかにしノイズ成分を除去する働きがあります。
さて、この式（17）をADMMで解くことを考えます。TV正則化項も重複ありグループ正則化と同じく変数間に絡みがあり、直接proximal operatorを計算するのは困難です。よって、重複ありグループ正則化を考えたときと同じく、として、重複なしグループ正則化項に変換することを考えます。これは以下のように、

と、となるように定義すれば、式（18）はのL2ノルムの和なので、 とすればTV正則化を重複なしグループ正則化に変換することができます。よって、このを用いれば式（16）の更新式でを更新することができます。一方、についての更新式は、式（3）についてで微分してとなる点を求めればよいので解析的に求まり、以下のようになります。

実装と実験

ノイズ除去アルゴリズムの実装

ということで、TV正則化によるノイズ除去を実装し、実験してみました。ここで、全てのピクセルを１つのベクトルとして扱ってしまうと画像サイズによってはメモリが足りなくなってしまうので、超解像を行ったときのように画像パッチを切り出し、それぞれのパッチに対してノイズ除去を行った後、そのパッチをつなぎ合わせて画像を再構築するという方法を取りました。なお、この画像切り出しやパッチのつなぎ合わせなどは、超解像の記事で記した「おれおれ画像・ベクトル・行列変換関数群」も利用してるので、そちらもご参照ください。

import os
import sys
import numpy as np
import scipy as sci
from scipy import sparse
from numpy.random import *
from PIL import Image

#ブロックソフト閾値計算
def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b, 2)) * b

#グループ正則化のproximal operator
def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, np.max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

#グループ正則化計算
def gloup_l1_norm(x, groups, p):
    s = 0
    for i in range(1, np.max(groups) + 1):
        s += np.linalg.norm(x[groups == i], p)
    return s

#関数をメモ化する関数
def memoize(f):
    table = {}
    def func(*args):
        if not args in table:
            table[args] = f(*args)
        return table[args]
    return func

#ADMMを行う関数
#argmin_x, argmin_y, update_lam, objectiveはそれぞれ変数名に準ずる計算を行う関数
def ADMM(argmin_x, argmin_y, update_lam, p, objective, x_init, y_init, lam_init, tol= 1e-9):
    x = x_init
    y = y_init
    lam = lam_init
    result = objective(x)
    while 1:
        x_new = argmin_x(y, lam, p)
        y_new = argmin_y(x_new, lam, p)
        lam_new = update_lam(x_new, y_new, lam, p)
        result_new = objective(x_new)
        if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break
        x = x_new
        y = y_new
        lam = lam_new
        result = result_new
    return x_new, result_new


#TV正則化付きのスムージング
#N,M:画像のピクセル行数、列数
#v:入力画像のベクトル
#B:正則化の変換行列（スパースなので、scipyのsparse行列を渡す）
#groups:変換後ベクトルの各要素の所属グループ
#C:正則化の係数
#p：拡張ラグランジュの係数
def TV_reg_smoothing(N, M, v, B, groups, C, p) :
    #inv(2I + pB^T B)を計算する関数。アルゴリズムによってはpが可変なので、pを引数として受け取る。
    #関数をメモ化して同一のpが入力された場合、再計算不要としている。
    inv_H = memoize(lambda p:np.array(np.linalg.inv( 2.0 * np.eye(B.shape[1], B.shape[1]) + p * (B.T * B))))

    argmin_x = lambda y, lam, p:np.dot(inv_H(p), 2.0 * v - np.array(B.T * (-p * y + lam)))
    argmin_y = lambda x, lam, p:prox_group_norm((B * x) + lam/p, groups, 1.0/p, C)
    update_lam = lambda x, y, lam, p:lam + p*((B * x) - y)
    objective = lambda x:np.linalg.norm(v - x, 2)**2 + C * gloup_l1_norm((B * x), groups, 2)

    x_init = np.random.randn(B.shape[1])
    y_init = (B * x_init)
    lam_init = np.zeros(B.shape[0])

    (x, result) = ADMM(argmin_x, argmin_y, update_lam, p, objective, x_init, y_init, lam_init, 1e-9)

    return x, result

#グループ変換行列を計算する関数
#N,M:画像のピクセル行数、列数
def calc_group_trans_mat(N, M):
    B = sci.sparse.lil_matrix((2 * (M - 1) * (N - 1) + (M - 1) + (N - 1), M * N))
    groups = np.zeros(B.shape[0], 'int')
    k = 0
    for i in range(N):
        for j in range(M):
            base = i * M + j
            if i < N -1 and j < M -1:
                B[k, base] = 1
                B[k, base + 1] = -1
                B[k + 1, base] = 1
                B[k + 1, base + M] = -1
                groups[k] = int(k/2)  + int(k % 2) + 1
                groups[k + 1] = int(k/2)  + int(k % 2) + 1
                k += 2
            #一番下の行のピクセルは右隣のピクセルとの差分のみ計算
            elif i >= N - 1 and j < M - 1:
                B[k, base] = 1
                B[k, base + 1] = -1
                groups[k] = int(k/2) + int(k % 2) + 1
                k += 1
            #一番右の劣のピクセルは下のピクセルとの差分のみ計算
            elif i < N - 1 and j >= M - 1:
                B[k, base] = 1
                B[k, base + M] = -1
                groups[k] = int(k/2) + int(k % 2) + 1
                k += 1
    return B, groups

#グレースケール画像に対するノイズ除去
#img:画像データ（PILのimageクラス）
#patch_hight:切り出す画像パッチの高さ
#patch_width:切り出す画像パッチの幅
#shift_num:画像を切り出す際のずらし量（重複して切り出してもOK）
#C:正則化の係数
#p:拡張ラグランジュの係数
def denoise_gray_img(img, patch_hight, patch_width, shift_num, C, p):
    #グループ変換行列の計算
    [B, groups] = calc_group_trans_mat(patch_hight, patch_width)

    #画像パッチの切り出し
    patchs = gray_im2patch_vec(img, patch_hight, patch_width, shift_num)

    new_patchs = np.zeros((patch_hight * patch_width, patchs.shape[1]))
    for i in range(patchs.shape[1]):
        #各パッチに対してノイズ除去を施す
        new_patchs[:, i] = TV_reg_smoothing(patch_hight, patch_width, patchs[:, i], B, groups, C, p)[0]

    #パッチをつなぎ合わせて画像の再構成
    [new_img, img_arr] = patch_vecs2gray_im(new_patchs, img.size[0], img.size[1], patch_hight, patch_width, shift_num)
    return new_img

細かい点はソースコードとそのコメントをご参照ください。とりあえず、denoise_gray_imgにPILのimageクラスのオブジェクトと切り出すパッチのサイズ、パッチ切り出し時のずらし量を渡せばノイズ除去されたimageクラスのオブジェクトを返してくれるように実装しています。超解像の時と同様、パッチを切り出す際に領域の重複を許し、再構成時に重ね合わせて平均を取るようにしています。また、画像の最後の行および列は右隣のピクセル（および下のピクセル）とのみ差分を取るようにBを工夫しています（calc_group_trans_mat関数内）。
注意する点としては、TV_reg_smoothingが受け取るBはスパースな行列なので、scipyのsparse行列のクラス渡します。また、同関数内にあるinv_Hはpを引数にとり、式（19）の逆行列部分を計算する関数ですが、同一のpに対して再度同じ計算を行うと時間がかかるので、memoizeによりメモ化し同じpが入力されたときに再計算しないようにしています。
なお、更新式はADMM関数の中に直接記述するのではなく、ADMM関数はそれぞれの更新式の計算を行う関数を引数として受け取るようにし、目的関数の出力が収束するまでそれらの更新式関数を呼び出すという設計にしているため、更新式の定義もTV_reg_smoothing関数の中で行っています。

実験と結果

いつものように256×256のレナさん画像で実験を行います。

これに対し、標準偏差30のガウスノイズを乗せた画像に対してノイズ除去を行ってみました。

def main():
    test_img = Image.open("ADMM_input.png")
    img_arr = np.array(test_img.convert('L'),'float')
    img_arr = (img_arr  + 30 * (np.random.randn(test_img.size[0], test_img.size[1]) ))
    img_arr[img_arr >255] = 255
    img_arr[img_arr <0] = 0
    test_img = Image.fromarray(np.uint8(img_arr))
    patch_hight = 20
    patch_width = 20
    shift_num = 10
    C = 50
    p = 20
    img2 = denoise_gray_img(test_img, patch_hight, patch_width, shift_num, C, p)
    img2.save("ADMM_result.png")

画像パッチのサイズを20×20ずらし量10として切り出し、C=50、p=20で実験を行った結果は以下のようになります。

　　
左が入力画像（標準偏差30のガウスノイズを重畳）、右がノイズ除去後の画像

元画像レベルとまではいきませんが、それなりにノイズ除去ができているかと思います。ただ、パッチを切り貼りしている影響か、少しゆがんでしまっているようにも見えます。更に追実験として入力画像のノイズレベルを標準偏差50と上げ、その時のノイズ除去もやってみました（C=80、p=1で実験）。

　
左が入力画像（標準偏差50のガウスノイズを重畳）、右がノイズ除去後の画像

ノイズレベルの割にはある程度ノイズ除去できているかと思いますが、TV正則化項の罰則に引っ張られ画像がノッペリとしてしまいました。Cをもっと調整すればもう少し改善できるかも知れません。

まとめ

今回は近接勾配法よりも幅広い問題を解くことができるADMMについて説明し、それを用いて画像のノイズ除去をやってみました。今回紹介した重複ありグループ正則化やTV正則化以外にも、特徴間のグラフ構造や階層構造なんかを表す正則化項がついた最適化問題もADMMを使えば解けるようになり、加えて別途制約条件がある場合でも式変形をうまくやればADMMで解けるようになるケースもあるなど応用範囲が広い手法です。更に収束も割りと速く、の値も適当に決めてもより大きければ大体収束するので、使い勝手もなかなか良いです。

参考文献

今回はこちらの書籍の主に12章と15章を参考にしました。

機械学習のための連続最適化 (機械学習プロフェッショナルシリーズ)

*1:関数があって、が凸関数ならば、関数は強凸であるといいます。が凸関数でなくとも、強凸である可能性はあります。

*2:逆に強凸であれば共役関数は微分可能でありアルゴリズムの収束が保証されます。

*3:便宜上、と２次元の添え字で表していてが行列データのように見えますが、実際は画像のピクセル値を直列に連結したベクトルとして扱っています。