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

はじめに

ここ１年くらい、ちまちまと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:便宜上、と２次元の添え字で表していてが行列データのように見えますが、実際は画像のピクセル値を直列に連結したベクトルとして扱っています。

近接勾配法応用編その3～トレースノルム正則化項付きロジスティック回帰による行列データの分類～

はじめに

前回の記事では、ベクトルデータの分類問題に対するスパースなロジスティック回帰を説明しましたが、今回はその拡張となる行列データに対するロジスティック回帰を見ていきます。この行列データ分類問題においても正則化項は重要な役割を持ちますが、中でもトレースノルムの正則化項を用いることで低ランクな解が得られ、多次元の時系列データの分類などに有効であることを見ていきます。

行列データの分類

まず、学習データ集合があるとします。ここで、はの行列、 はのラベルです。この学習データ集合から係数行列とバイアス項を学習し、以下に示すような線形モデルで分類することが今回のテーマです。

行列データ分類問題の特徴

では、行列データを直接扱うことにどのようなメリットがあるのでしょうか。確かに、行列データをベクトルに変換して（列ベクトルを縦に連結するなどして）、学習・分類するという方法も考えられます。しかしベクトルにしてしまうと行列構造がなくなり、行と列の意味合いが異なる場合に重要な情報が失われる可能性があります。
例えば、複数のセンサーから取得される多次元の時系列データの場合、行がセンサー、列が時間に対応する行列データとしてとらえることができます。このようなデータを連結してベクトルとして扱ってしまうと、空間方向と時間方向の特徴を一緒くたにしてしまい、重要な情報が失われてしまう可能性があります。特にこのような複数センサーのデータの場合、実際に観測できるデータはのように、潜在変数（信号源）の空間的な線形結合で表されることがしばしばあり、クラス分類において有用な情報はこの潜在変数に隠れていることも多いです（下図参照）。

一方で、行列データとして直接扱うことの意味は、と特異値分解して式（1）を以下のように変形すると見えてきます。

このようにすると、は空間方向の特徴を抽出する空間フィルタ、は時間方向の特徴を抽出する時間フィルタとして考えることができ、時空間両方の特徴を捉えられることがわかります。上述のように観測されるデータが潜在変数の空間的な線形結合である場合にも、観測データに空間フィルタをかけることでクラス分類に有用な潜在的な情報を抽出した上で時間方向に係数をかけるので、時空間両方の特徴をとらえることができます。

行列データロジスティック回帰の目的関数

行列データの分類問題においてもロジスティック回帰は使えます。具体的には下記の誤差関数を最小化するとを学習データ集合から求めればよいです。

トレースノルム正則化による低ランク化

さて、式（2）にあるように係数行列のランクが高いほど、かけ合わせる時空間フィルタの個数は増えます。もし、がフルランクなら個の時空間フィルタをかけて重み付き和を取ることになります。しかし、上述のような潜在変数の中でもクラス分類に有用な情報は、（空間的に）ごく一部の成分でありそれ以外はノイズというケースも多くあります。例えば、上の図の例だと、3つの信号源のうちクラス分類に有用な成分は1つです。この場合にがフルランクだと必要以上に複雑なモデルとなってしまい過学習となる恐れがあります。
ということで、を低ランクにしたいのですが、そのためにはこちらの記事で少し述べたトレースノルムの正則化項をつけると特異値がスパースになり、結果低ランクな解が得られます。具体的な式としては、

であり、これを近接勾配法で解いていきます。なお、今回はバイアス項は正則化の中に含めていないことに注意してください*1。は正則化項に含まれないので、は近接勾配法の更新式、は通常の勾配法の更新式で更新していきます。それぞれの更新式で必要な式（3）のによる微分および式（4）のによる微分は、


となるので、後はこれを元にとの更新を交互に行えばよいです。

実装と実験

アルゴリズムの実装

今回の行列データロジスティック回帰は下記のように実装しました。基本的に今までと同じですが、proximal operatorはトレースノルム正則化のproximal operator計算を行うようにしており、近接勾配法の関数はバイアス項の更新（勾配法計算）も行うように変更しています。

import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt
import sys

#ソフトしきい値作用素の計算
def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

#トレースノルム正則化のproximal operatorの計算
def prox_trace_norm(V, gamma, lam):
    [L, sig, R] = np.linalg.svd(V)
    sig_ = soft_thresh(sig, lam * gamma)
    if L.shape[1] > sig.shape[0]:
        L = L[:,:sig.shape[0]]
    if R.shape[1] > sig.shape[0]:
        R = R[:sig.shape[0], :]
    return np.dot(np.dot(L, np.diag(sig_)), R)

#近接勾配法とバイアス項の勾配法計算
def proximal_gradient_with_bias(grad_f, grad_b, prox, gamma, objective, init_x, init_b, tol = 1e-9):
    x = init_x
    b = init_b
    result = sys.maxint
    while 1:
        x_new = prox(x - gamma * grad_f(x, b), gamma)
        b_new = b - gamma * grad_b(x, b)
        result_new = objective(x_new, b_new)
        if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        b = b_new
        result = result_new
    return x_new, b_new, result_new

#トレースノルム正則化項付きロジスティック回帰
#X:学習用データベクトルを列ベクトルとして並べた行列
#y:ラベルを並べたベクトル
def trace_norm_LogisticRegression(X, y, lam):
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, W, b:sigma(np.trace(np.dot(Z.T, W).T) + b) #メモ化したい・・・
    objective = lambda W, b: - np.sum( y * np.log(p(X, W, b)) +
                        (1 - y) *np.log(( 1 - p(X, W, b)))) + lam * np.linalg.norm(W, 'nuc')
    grad_E = lambda W, b: np.dot(X, p(X, W, b) - y) #勾配
    grad_b = lambda W, b: np.sum(p(X, W, b) - y)

    (u, l, v) = np.linalg.svd(X)
    Gamma =1.0/max(np.average(l.real*l.real,0)) 
    prox = lambda V, gamma:prox_trace_norm(V, gamma, lam) #proximal operator
    W_init = np.zeros((X.shape[0], X.shape[1]))
    b_init = 0
    (W, b, result) = proximal_gradient_with_bias(grad_E, grad_b, prox,
                                                    Gamma, objective, W_init, b_init, 1e-4)
    return W, b, result

なお、停止判定の方法については収束性を考慮して元の方法に戻しています。また、p(Z, W, b)の計算がかなり重いですが１回の更新で4回計算しているので、pをメモ化するなどの工夫を行えばもう少し速くできます（pの引数ZとWはnumpyのarray型なので簡単に辞書のキーにできず断念）。
あと、学習のステップ幅は今回かなりテキトーに決めた実装になっているので、その辺りはもう少し検討が必要です。

実験用疑似データ

実験用疑似のデータは下記のコードで生成しました。

    #学習データ
    #30×100の行列データ負例500個、正例500個
    #信号源のデータZ
    Z_train1=np.random.randn(30,100,500)
    Z_train2=np.random.randn(30,100,500)
    #30個の信号源のうち2個のみクラス分類に有用とする
    Z_train2[0, 20:50, :] += 0.2
    Z_train2[29, 80:90, :] -= 0.5
    Z_train = concatenate((Z_train1,Z_train2),2)

    #信号源データZを観測データに変換する行列A
    A = np.exp(-np.abs(np.tile(np.linspace(0, 1,30), (30,1)) -
            np.tile(np.linspace(0, 1,30), (30,1)).T)/0.1) + np.random.randn(30, 30)

    #実際に観測されるデータX
    X_train = np.dot(Z_train.T, A).T

    y_train1 = ones((500,1))
    y_train2 = zeros((500,1))
    y_train = concatenate((y_train1,y_train2),0).squeeze()

    #テストデータ
    #学習データと同じく30×100の行列データ負例500個、正例500個
    Z_test1=np.random.randn(30,100,500)
    Z_test2=np.random.randn(30,100,500)
    Z_test2[0, 20:50, :] += 0.2
    Z_test2[29, 80:90, :] -= 0.5
    Z_test = concatenate((Z_test1,Z_test2),2)

    X_test = np.dot(Z_test.T, A).T

    y_test1=ones((500,1))
    y_test2=zeros((500,1))
    y_test=concatenate((y_test1,y_test2),0).squeeze()

学習、テストに使用するデータはX_train、X_testですが、上述の信号源からの線形結合を模擬するため、まず潜在変数（信号源）Z_train、Z_testの生成を行いました。この信号源のデータは30成分ありますが、0番目の成分と29番目の成分のみがクラス分類に有用となるようなデータとしています。これを行列Aで線形結合して観測データXを生成しています。なお、行列Aは近くの信号程大きな重み、離れるほど小さな重みで足し合わせるような線形結合としており、更にこのAにもランダム成分を加えています。以下、正例・負例におけるXとZの平均波形を表示します。

• y=1に対応する観測データ集合X_train1の平均波形

• y=0に対応する観測データ集合X_train2の平均波形

• y=1に対応する信号源データ集合Z_train1の平均波形

• y=0に対応する信号源データ集合Z_train2の平均波形

やはり、観測データはy=0クラスとy=1クラスの平均データを見比べても、ノイズにまみれて差異があまり分からなくなっていますが、信号源データの場合は、2信号のみy=0クラスとy=1クラスで差異があることがはっきりわかります。

実験結果

このように生成した疑似データに対し、以下のコードのようにと変化させて、その時の推定精度とのランクを求めてみました。

    num = 41
    scale = 50
    rank_of_W  = np.zeros(num)
    acc_arr = np.zeros(num)
    for i in range(0, num):
        (W, b, r) = trace_norm_LogisticRegression(X_train, y_train.squeeze(), i * scale)
        y_result = np.trace(np.dot(X_test.T, W).T) + b
        y_result[y_result>=0] = 1
        y_result[y_result<0] = 0
        rank_of_W[i] = np.linalg.matrix_rank(W)
        acc_arr[i] = np.sum(y_result == y_test)/float(y_test.shape[0])

その結果が以下の図です（赤が推定精度、青がランクです）。

L1ノルムロジスティック回帰と同様、を大きくしすぎても、制約がきつくなりすぎるのか精度が徐々に下がっていきます。実際今回の条件で推定精度が最大になったのはのときで、精度が68.0%、のランクが4でした。以下には、このでのを特異値分解して得られる空間フィルタを、y=0クラスのX_trainとy=1クラスのX_trainの平均行列に施した波形を示します。

• y=1に対応する観測データ集合X_train1の平均にをかけて得られる波形

• y=0に対応する観測データ集合X_train2の平均にをかけて得られる波形

この図の中で青の波形は最大特異値に対応する空間フィルタで抽出された波形、緑の波形は２番目に大きい特異値に対応する空間フィルタで抽出された波形です。この両者の注目すると、空間フィルタを施すことによって信号源に潜在しているクラス分類に有用な成分を抽出できていることがわかります。

まとめ

今回は行列データの分類アルゴリズムについて、ロジスティック回帰を例にとってトレースノルム正則化項によって低ランクな解が得られること、多次元の時系列データの分類問題などに有効であることを見てきました。トレースノルムの正則化は信号源の中でも一部の成分のみがクラス分類に有用であるとき特に威力を発揮します。
なお、時系列データの分類には今回紹介した方法以外にも隠れマルコフモデルを使った手法やリカレントニューラルネットワークの一種であるLong Short-Term Memory(LSTM)を使った手法もあります。特に最近は深層学習の流れもあってLSTMが主流になってきているとのことなので、そろそろニューラルネット辺りにも手を出していきたいところです。

参考文献

以前も紹介したこちらの本の第８章が非常に参考になります。

スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)

あと、同著者によるこちらのスライドもわかりやすいです。
行列およびテンソルデータに対する機械学習

*1:どうやら、前回のベクトルデータの分類においても、基本的にバイアス項は正則化には含めない方が良いようです

近接勾配法応用編その2 ～L1ノルム正則化項によるスパースなロジスティック回帰～

今回は前々回の記事で書いた近接勾配法の応用例第２段ということで、分類問題などで使われるロジスティック回帰の目的関数にL1ノルムの正則化項をつけて、スパースな解を導いてみます。スパースな解を導出するメリットとしては、クラス分類に寄与しない無駄な属性の係数が自動的に0となり、クラス間で有意差がある属性を見つけ出せ、かつ、分類モデルの複雑度も下げられるので過学習の防止も期待できるというメリットがあります。なお、ロジスティック回帰は多クラス分類にも利用できますが、今回は２クラス分類の例で説明します。

ロジスティック回帰

ロジスティック回帰について詳しいことはたくさん文献が出ているのでそちらを参照いただくとして、ここでは概略と最終的な目的関数を示します。
まず、学習データ集合があるとします。ここで、は次元のデータベクトル、はのラベルです。
ここで、が与えられてとなる事後確率は

と表せます。ここで、は

という式であり、ロジット関数とも呼ばれます。がクラスに所属する確率が大きければこのロジット関数は大きな値を取るので事後確率は1に近づきますし、逆にクラスに所属する確率が大きければロジット関数は小さな値を取るのでは0（つまり、クラスとなる事後確率は1）に近づきます。
そしてこのをの線形関数としてモデル化したのがロジスティック回帰です。具体的には、

ただし、, *1

となります。
さて、これまでの内容を踏まえて、学習データ集合に対する尤度関数を考えます。上記の議論より、とが与えられて、そのクラスがである尤度はです。一方で、である尤度は です。このように考えると尤度関数は、

となります。後は式（4）を最大化するを求めればよいです。ただ、この式（4）のままだと計算しにくいので、以下のような式（4）の負の対数を取った誤差関数を最小化することを考えます。

これがロジスティック回帰の目的関数です。

正則化項によるスパース化と近接勾配法の適用

ここまでで普通のロジスティック回帰の目的関数が得られました。このロジスティック回帰をスパースな解にするためには、L1ノルム正則化項をつけた式

を近接勾配法により最小化すればよいです。近接勾配法を解くためには、元の目的関数の勾配と正則化項 のproximal operatorが必要でした。
正則化項のproximal operatorは以前の記事で導出したとおりです。また、の勾配は、

となります。後は、これらに基づき近接勾配法を適用するだけです。

実装と実験

これまでの内容pyhonコードで実装すると下記のようになります。なお、近接勾配法のアルゴリズム部分についても、前々回記事で記したコードから停止条件部分を修正したので、改めて近接勾配法の関数も書きます（以前の記事から変更していない関数については下記に記していません）。

import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#近接勾配法
def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        #停止条件を説明変数差分のノルムに変更
        #if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
        if np.linalg.norm(x_new - x, 2) < tol:
            break;
        x = x_new
        result = result_new
    return x_new, result_new

#l1ノルム正則化ロジスティック回帰
#X:学習用データベクトルを列ベクトルとして並べた行列
#y:ラベルを並べたベクトル
def l1norm_LogisticRegression(X, y, lam):
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, w:sigma(np.dot(Z.T, w))  # ※1下記に補足あり
    X=concatenate((ones((1,size(X,1))), X), 0) #バイアス項の追加
    objective = lambda w: - np.sum( y * np.log(p(X, w)) +
                        (1 - y) *np.log(( 1 - p(X, w)))) + lam * np.sum(np.abs(w))
    grad_E = lambda w: np.dot(X, p(X, w) - y) #※2下記補足あり
    (u, l, v) = np.linalg.svd(X)
    gamma = 1.0/max(l.real*l.real)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam) #proximal operator
    w_init = np.zeros(X.shape[0])
    (w_hat, result) = proximal_gradient(grad_E, prox, gamma, objective, w_init, 1e-12)
    return w_hat, result

繰り返し構文を避けるため、行列・ベクトル演算で勾配計算などを行っています。これを少し補足しておくと、まずソースコード※1に当たる部分は、というベクトルを返す関数として定義しています。そして、※2の部分では、

というように、行列・ベクトルの演算で式（7）の勾配計算を行う関数を定義しています。
このL1ノルム正則化つきロジスティック回帰を疑似データで実験してみました。疑似データは以下のようなコードで生成しました。

def main():
    #学習データ
    #最初の２次元だけ分類に寄与する属性
    X_train1=np.random.randn(2,200)
    X_train2=np.random.randn(2,200) + ones((2,200))
    X_train = concatenate((X_train1,X_train2),1)
    #残り123次元は分類に寄与しないランダムな属性
    X_train = concatenate((X_train, np.random.randn(123, 400)), 0)
    y_train1 = ones((200,1))
    y_train2 = zeros((200,1))
    y_train = concatenate((y_train1,y_train2),0).squeeze()

    #テストデータ
    X_test1=np.random.randn(2,200)
    X_test2=np.random.randn(2,200)+ones((2,200))
    X_test = concatenate((X_test1,X_test2),1)
    X_test = concatenate((X_test, np.random.randn(123, 400)), 0)
    y_test1=ones((200,1))
    y_test2=zeros((200,1))
    y_test=concatenate((y_test1,y_test2),0).squeeze()

疑似データの次元は125次元ですが、その内最初の2次元だけクラス分類において意味があり、それ以外の123次元は全くのランダムでノイズとなる属性です。また、y=0のクラスにのみオフセットをはかせている形なので、バイアス成分もクラス分類に寄与します。そのため、今回のアルゴリズムを適用すると、最初の2次元の係数およびバイアス項のみが非0となり、それ以外の係数は0になることが期待されます。
このように生成した疑似データに対し、以下のようにと変化させて、その時の推定精度との非0要素の数を求めてみました。

    non_zero_cnt_arr = np.zeros(40)
    acc_arr = np.zeros(40)
    for lam in range(0, 40):
        (w, r) = l1norm_LogisticRegression(X_train, y_train.squeeze(), lam)
        y_result = np.dot(w[1:126], X_test) + w[0]
        y_result[y_result>=0] = 1
        y_result[y_result<0] = 0
        non_zero_cnt_arr[lam] = np.sum(np.abs(w)>1e-12)
        acc_arr[lam] = np.sum(y_result == y_test)/float(y_test.shape[0])

    fig, ax1 = plt.subplots()
    ax1.plot(acc_arr, "r.-")
    ax1.set_ylabel("accuracy")
    ax2=plt.twinx()
    ax2.plot(non_zero_cnt_arr, "bx-")
    ax2.set_ylabel("non zero element count")
    plt.legend()
    plt.show()

その結果は下記のようになりました（赤線が推定精度、青線が係数の非0要素数です）。

のときは、すべての要素が非0になりノイズである属性も判別に用いてしまっているので精度が60.5%とあまり良くありません。そこからを大きくしていくと非0要素数は減少していくことがわかります。一方で推定精度については、適切なであれば精度が上がっていますが、逆にが大きすぎると制約がきつくなりすぎるためか精度は下がっていきます。面白いのは、非0要素数が3のとき（クラス分類に寄与する２次元 + バイアス項のみが非0のとき）、推定精度がピークとなると予想していましたが、予想に反してで非0要素数8のときに推定精度が最大（75.5%）になる結果が得られました*2。これは、適切にクラス分類に寄与する属性を見つけられるまでを上げると、ロジスティック回帰の制約条件としては厳しくなりすぎてしまい結果、精度の低下を招くのだと考えられます。そのため、クラス分類問題とどの属性に有意差があるか見つける問題とは一緒くたにして解くのではなく、分けて解くべきなのだと思います。

まとめ

今回はロジスティック回帰にL1ノルム正則化項をつけることによりスパースな解が得られ、分類器としての精度も向上させられることを見てきました。また、分類に寄与する属性もこのL1ノルム正則化を使うことで見つけられますが、だからと言って制約条件が厳しすぎるを用いると分類精度低下を招くので注意が必要です。一旦、変数選択の目的でL1ノルム正則化項付きのロジスティック回帰を行って、そこで係数が非0となった要素だけを使って今度は制約なしのロジスティック回帰を再学習させる方法がよいのかも知れません。なお、今回はL1ノルムの正則化を試しましたが、もちろんL1ノルム以外にもグループL1ノルムやL1/L2の混合ノルムも適用できます。
今後の予定としては、もうひとつ近接勾配法の応用として行列データの分類について紹介したり、Alternating Direction Method of Multipliers (ADMM)を調べて書いたり、ICML2018の論文（そろそろ公開？）を読んだりなんかをしたいと思います。

*1:バイアス項をに統合し、次元のベクトルとして扱います。

*2:ランダムな疑似データなので試行毎に結果は変わりますが、傾向としては毎回同じで非0要素数が3よりも多い時点で推定精度のピークを迎えます。

近接勾配法応用編その1 ～スパースコーディング、辞書学習からの超解像～

はじめに

前回の記事で近接勾配法の基本と実装を見てきました。今回はこの近接勾配法を利用して、スパースコーディングの応用の１つである超解像技術を調べPythonで実装してみました*1。画像処理関係の実装はほぼ未経験だったので、その分野の人が見れば前処理等が足りなかったり、アンチパターン的なことをやってたりするかもしれませんが、その点は指摘頂けると幸いです。

スパースコーディングとは

スパースコーディングは観測した信号を、「辞書」と呼ばれる個の基底ベクトルの集合　のスパースな線形結合で表現する符号化手法です。
今、この線形結合の係数をとすると、スパースコーディングは下記のように定式化されます。
　
ただし、の要素の多くは（つまりはスパース）です。

このようなスパースコーディングを使うことによりデータの圧縮を行えるだけでなく、画像からのノイズ除去や今回取り上げる超解像など様々な技術へ応用が可能です。
スパースコーディングを成功させる上で重要となるのは辞書です。がテキトーだとうまくいきません*2。ということで、辞書は実データからの学習により構築します。ここで、学習データの個数を、学習データをまとめた行列、 各線形結合の係数をまとめた行列を定義すると、辞書学習の目的関数は下記のようになります。
　
ただし、はスパース。

この方法で辞書が得られたら、次はその辞書を用いて入力画像の符号化します。そして、その係数と辞書から画像を再構成し超解像画像を得るという流れになります。
ということで、スパースコーディングによる超解像生成は（1）辞書を学習するフェーズと（2）その辞書を元に画像を再構成するフェーズの2フェーズで実現できます。

画像再構成のための辞書学習

実装した手法の説明

では、今回実装した学習用の画像データを元に辞書を学習する方法を説明します。
まず、最終的に得たい解像度の画像データ（高解像画像データ）複数枚を学習データとして用意します。この時、超解像を行いたい画像のカテゴリが決まっているなら、そのカテゴリに合った画像を用意するのがよさそうです。
この学習データを用いて、辞書を学習する流れを下記に示します。

1. 学習用画像を小領域を取り出した「画像パッチ」を個用意し、それぞれ縦に並べたベクトルおよび、これらのベクトルを横に並べた行列をを作る。
2. の各行が平均、分散になるように標準化する。
3. をランダム行列で初期化する。
4. を固定して、次の最適化問題をについて解き、行列を得る。

5. を固定して、次の最適化問題を解きを得る。

6. 手順4,手順5の操作を収束するまで、もしくは一定回数繰り返す。

ここで、式（3）はlassoの式であるため解がスパースになることが期待され、更に近接勾配法で解くことができます。なお、今回、近接勾配法の応用編ということで手順4の方法でスパースな係数を求めていますが、この方法以外にもOrthogonal Matching PursuitやIterative Reweighted Least Squaresなどのアルゴリズムがあります。
また、手順5で式（4）の解を得る方法としてはK-SVDという手法がよく使われますが、今回は簡単のために下記のように解析的に解いた結果を使用します。

この手順4、手順5を繰り返しての同時最適化を行うわけですが、この両変数に対しては凸関数ではないため（式（3）、式（4）はそれぞれ凸関数ですが）、求まる結果は局所最適解です。
なお、この学習を収束まで待っているとかなり時間がかかるため、今回は適当な回数ループさせたら学習を完了させています。

ちなみに、最初は手順2の標準化を省いた実装をしていましたが、これだとうまく超解像ができませんでした。やはり、学習データの標準化は重要のようです。

実装したコード

ここまでの手順を実現するpythonコードです。

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

#dorectory以下のファイルを再帰的に検索する関数
def find_all_files(directory):
    for root, dirs, files in os.walk(directory):
        yield root
        for file in files:
            yield os.path.join(root, file)

#おれおれ画像・ベクトル・行列変換関数群
def color2gray_img(img):
    return img.convert('L')

def im2vec(img):
    img_arr = np.array(img.convert('L'),'float')
    return img_arr.reshape(-1)

def gray_im2float_vec(gray_img):
    return im_scale2float((im2vec(gray_img)))

def color_im2float_vec(color_img):
    return gray_im2float_vec(color2gray_img(color_img))

def vec2image(vec, H, W):
    tmp_vec = copy.deepcopy(vec)
    tmp_vec = tmp_vec.reshape((H, W))
    return Image.fromarray(np.uint8(tmp_vec))

def float_vec2image(vec, H, W):
    return vec2image(float2im_scale(vec), H, W)

def arr2image(arr):
    arr = copy.deepcopy(arr)
    return Image.fromarray(np.uint8(arr))

def float_arr2image(arr):
    return arr2image(float2im_scale(arr))

def im_scale2float(vec):
    return vec.astype(np.float)

def float2im_scale(vec):
    vec[vec > 255.0] = 255.0
    vec[vec < 0.0] = 0.0
    return vec

#画像データからパッチを切り出す関数
#shift_numで切り出す際のずらし量を指定
def gray_im2patch_vec(im, patch_hight, patch_width, shift_num):
    patchs = np.zeros((patch_hight*patch_width,
                        (1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)) *
                        (1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0))))
    im_arr = im_scale2float(np.array(im))
    k = 0
    for i in range(1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)):
        h_over_shift = min(0, im.size[0] - ((i* shift_num) + patch_hight)) #ずらしたことによって画像からはみ出る量(縦方向)
        for j in range(1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0)):
            w_over_shift = min(0, im.size[1] - ((j* shift_num) + patch_width)) #ずらしたことによって画像からはみ出る量(横方向)
            #パッチの切り出し　画像からはみ出たら、そのはみ出た分を戻して切り出し
            patchs[:, k] = im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
                                    (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift].reshape(-1)
            k += 1
    return patchs


#式（5）を計算する関数
def dictionary_learning_with_mod(W, Y):
    D = np.dot(np.dot(Y, W.T), np.linalg.inv(np.dot(W, W.T)))
    result = np.linalg.norm(Y - np.dot(D, W)) ** 2
    return D, result

#辞書クラス
class dictionary_with_learning():
    def __init__(self, weight_learner, dictionary_learner, base_num, max_loop_num = 8):
        self.Dic = np.array([])
        self.weight_learner = weight_learner
        self.dictionary_learner = dictionary_learner
        self.base_num = base_num

    #学習データ行列Yから辞書を学習するメソッド
    def learn(self, Y, eps = 1e-8):
        #手順3:ランダム行列で辞書の初期化
        D = np.random.randn(Y.shape[0], self.base_num)

        result = -1
        loop_num = 0
        while 1:
            loop_num += 1
            W = np.zeros((self.base_num, Y.shape[1]))
            #手順4:Dを固定して係数を求める
            for i in range(Y.shape[1]):
                [W[:, i], tmp] = self.weight_learner(D, Y[:, i])

            #手順5:Wを固定してDを求める
            [D, new_result] = self.dictionary_learner(W, Y)

            if result < 0 or np.abs(new_result - result) > eps or loop_num < max_loop_num:
                result = new_result
            else: #手順4、5を収束するまで繰り返す。
                break;
        self.Dic = D
        return D

    def get_dic_array(self):
        return self.Dic



#画像パッチから学習される辞書クラス
class image_patch_dictionary():
    def __init__(self, dictionary):
        self.dictonary = dictionary
        self.patch_hight = 0
        self.patch_width = 0

    #学習画像の集合から辞書を学習するメソッド
    def learn_images(self, image_list, patch_hight, patch_width, eps=2e-2, shift_num=1):
        self.patch_hight = patch_hight
        self.patch_width = patch_width
        Y = np.empty((patch_hight * patch_width, 0))
        #手順1:画像パッチの切り出し
        for img in image_list:
            gray_img = color2gray_img(img)
            Ytmp = gray_im2patch_vec(gray_img, patch_hight, patch_width,shift_num) 
            Y = np.concatenate((Y, Ytmp), axis=1)

        #手順2:学習データ行列の標準化
        Y_ = (Y - np.mean(Y, axis=1, keepdims = True))/np.std(Y, axis = 1, keepdims=True)
        self.dictonary.learn(Y_, eps)

    def get_dic(self):
        return self.dictonary.get_dic_array()

    def get_patch_hight(self):
        return self.patch_hight

    def get_patch_width(self):
        return self.patch_width


def main():
    lam = 5
    base_num = 144
    patch_hight = 12
    patch_width = 12
    train_img_list = []
    for file in find_all_files('C:\\data'):
        if os.path.isfile(file):
            try:
                train_img_list.append(Image.open(file))
                if len(train_img_list) >= 50:
                    break;
            except:
                print file

    weight_learning_method = lambda W, y:ISTA(W, y, lam) #近接勾配法を行う関数を使用
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)#解析的に辞書Dを求める方法を使用
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.learn_images(train_img_list, patch_hight, patch_width, shift_num = patch_hight)

    np.save('dic_array.npy', im_dic.dictonary.Dic)

本実装は、画像⇔行列・ベクトル変換を行う関数群、パッチ切り出しを行う関数、画像パッチ辞書クラスimage_patch_dictionary、辞書クラスdictionary_with_learningから構成されています。
image_patch_dictionaryは実際の辞書学習オブジェクトであるdictionary_with_learningをインスタンス変数として持ち、与えられた学習画像群からパッチの切り出しを行い、辞書を学習させる機能を持ちます。
そのため、辞書学習を行うのはdictionary_with_learningになりますが、係数の学習を行うweight_learnerと辞書の学習を行うdictionary_learnerをインスタンスとして持ち、生成時にどのアルゴリズムで係数推定、辞書推定を行うかを指定します。今回、係数の計算には近接勾配法を使用しますので、前回の記事記載のISTA関数をweight_learnerとして指定し、dictionary_learnerには新たに作成したdictionary_learning_with_mod関数を指定します（main関数参照）。もし係数推定として別のアルゴリズムを使いたければweight_learnerを変えればよいし、辞書学習もK-SVDなどを使いたければdictionary_learnerを変更すればよいです。

画像再構成による超解像

実装した手法の説明

ここからは、スパースコーディングによる超解像の手法の説明です。手法の大まかな考え方は「学習した辞書をダウンサンプリングした辞書で入力画像を符号化する係数を求めた後、その係数と元の辞書から画像を再構成する」というものです。以下、順を追って説明します。

1. 画像パッチの切り出し
今、入力画像（低解像画像）をとし、そこから辞書学習時と同様に画像パッチを切り出します。その切り出した画像パッチをとします。切り出す画像パッチのサイズは、辞書学習に使用した画像（高解像画像）のパッチサイズを、低解像画像から高解像画像への倍率をとすると、になります。ということで、は次元のベクトルになります。ここで画像パッチは、もれなくかつ領域を重複させて取得します。領域を重複させないと、画像再構成後パッチ間のつなぎ目が滑らかにならずブロック状のアーティファクト発生します。なるべく重複度が大きい方が精度は良くなりますが、取得する画像パッチも多くなるので、その分計算に時間がかかります。
2. 辞書のダウンサンプリング
次に辞書のダウンサンプリングです。ここの方法も色々考えられますが、今回は単純に学習したにダウンサンプリングの作用素をかけてを得る方法をとりました。 は学習画像パッチのサイズ 基底数の行列です。各基底はの画像を縦に並べたベクトルと見なせるので、これをの画像に戻し、近傍の画素を平均することでダウンサンプリング後の辞書を得ます（はの行列になります）。
3. 符号化係数の推定
画像パッチが切り出せたら次は、それぞれ符号化する係数を求めます。これは、ダウンサンプリング後の辞書を使い、下記のように式（3）と同様の最適化問題を解けばよいです。
4. 高解像画像パッチの生成
ここまでで、各パッチを符号化するスパースな係数が求まりました。次はこの係数を用いて高解像画像を構成する画像パッチを求めます。やることは簡単で、以下のように元々の辞書（高解像画像の辞書）と低解像辞書で推定した係数をかけるだけです。
5. パッチをつなぎ合わせて高解像画像を再構成
得られた高解像画像パッチから高解像画像を再構成します。手順1で画像パッチを領域を重複させて切り出したことに注意が必要です。今、手順1で画像を切り出す際ピクセルずらして切り出したとすると、再構成時はピクセルずらしてパッチを貼り合わせます。全てのパッチを貼り合わせた後、各ピクセルを加算回数分で割り平均をとります。この画像パッチのつなぎ合わせについても色々な工夫が提案されているようです。
6. バックプロジェクションによる画像の補正
最後に仕上げです。パッチを貼り合わせ高解像画像は、パッチごとに最適化されているため、画像全体の仮定が考慮されていません。その仮定というのは、「入力画像は、高解像画像に、劣化やダウンサンプリングを施す作用素をかけることによって得られたものである。」というものです。つまり、という仮定が前提にあり、それと推定したとのバランスをとることで、解の補正を行います。具体的には下記の最適化問題を解きます。 ここで、は画像全体の各要素を縦に並べたベクトルです。また、今回の実験では、はダウンサンプリングのみを施す行列として定義しています。 この最適化は解析的に求められそうですが、多くの論文では勾配法を使って求めています。逆行列計算が大変だからでしょうか？とりあえず、今回の実装でも論文に倣って勾配法で解きました。

実装したコード

以上の内容の実装が下記です。モジュールのimportやおれおれ画像変換関数群は辞書学習時のコードと同様のものを使っています。

#バックプロジェクションで使うダウンサンプラー行列を求める関数
#ダウンサンプラーはスパースかつサイズが大きい行列なのでscipy.sparseのスパース行列型として返す
def get_downsampler_arr(from_H, from_W, to_H, to_W):
    S = sci.sparse.lil_matrix((to_H * to_W, from_H * from_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            s = sci.sparse.lil_matrix((from_H, from_W))
            s[i : i + from_H/to_H, j:j + from_W/to_W] = (float(to_W)/float(from_W)) * (float(to_H)/float(from_H))
            S[k, :] = s.reshape((1, from_H*from_W))
            k += 1
    return  S.tocsr()

#ベクトルデータのダウンサンプリング関数
def vec_downsampling(vec, from_H, from_W, to_H, to_W):
    arr = vec.reshape((from_H, from_W))
    ret_arr = arr_downsampling(arr, from_H, from_W, to_H, to_W)
    return ret_arr.reshape(-1)

#行列データのダウンサンプリング関数
def arr_downsampling(arr, from_H, from_W, to_H, to_W):
    ret_arr = np.zeros((to_H, to_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            ret_arr[i/(from_H/to_H), j/(from_W/to_W)] = (np.average(arr[i : i + from_H/to_H, j:j + from_W/to_W]))
            k += 1
    return ret_arr

#画像パッチをつなぎ合わせて画像を再構成する関数
#shift_numがずらし量、パッチがオーバーラップしている個所は足し合わせて平均化
def patch_vecs2gray_im(patchs, im_hight, im_width, patch_hight, patch_width, shift_num):
    im_arr = np.zeros((im_hight, im_width))
    sum_count_arr = np.zeros((im_hight, im_width)) #各領域何回パッチを足したか表す行列
    i = 0
    j = 0
    w_over_shift = 0
    h_over_shift = 0
    for k in range(patchs.shape[1]):
        P = patchs[:, k].reshape((patch_hight, patch_width))
        im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
            (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift] += P #指定の領域にパッチの足しこみ
        sum_count_arr[(i* shift_num) + h_over_shift :(i* shift_num) + patch_hight + h_over_shift,
            (j* shift_num) + w_over_shift:(j* shift_num) + patch_width + w_over_shift] += 1 #当該領域のパッチを足しこんだ回数をカウントアップ
        if j < ((im_width - patch_width)/shift_num + (1 if (im_width - patch_width) % shift_num  != 0 else 0)):
            j += 1
            w_over_shift = min(0, im_width - ((j* shift_num) + patch_width)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数（横方向）
        else:
            j = 0
            i += 1
            h_over_shift = min(0, im_hight - ((i* shift_num) + patch_hight)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数（縦方向）
            w_over_shift = 0
    im_arr /= sum_count_arr
    return float_arr2image(im_arr), im_arr

#バックプロジェクションを行う関数
def back_projection(x0, y, L, c, mu, x_init = [], eps = 1e-7):
    if x_init == []:
        xt = sci.matrix(np.random.rand(x0.shape[0], x0.shape[1])) *255.0
    else:
        xt = x_init
    result = -1
    while 1:
        x_new = xt - mu * (L.T.tocsr()*(L * xt - y) + c* (xt - x0))
        result_new = sci.linalg.norm(L * x_new - y, 2) + c * sci.linalg.norm(x_new - x0, 2)
        if np.abs(result_new - result) < eps:
            break;
        xt = x_new
        result = result_new
    return x_new, result_new

#超解像計算クラス
class super_resolution:
    def __init__(self, dictionary, weight_learner, back_projection):
        self.Dh = dictionary.get_dic()
        self.patch_hight = dictionary.get_patch_hight()
        self.patch_width = dictionary.get_patch_width()
        self.weight_learner = weight_learner
        self.back_projection_method = back_projection

    def get_super_resolution(self, image, scaling, shift_num = 1):
        #手順1:画像パッチの切り出し
        image = color2gray_img(image)
        Y = gray_im2patch_vec(image, self.patch_hight/scaling, self.patch_width/scaling, shift_num)

        #手順2:辞書のダウンサンプリング
        Dl = np.zeros((self.Dh.shape[0]/(scaling**2), self.Dh.shape[1]))
        for i in range(self.Dh.shape[1]):
            Dl[:, i] = vec_downsampling(self.Dh[:, i], self.patch_hight, self.patch_width,
                                        self.patch_hight/scaling, self.patch_width/scaling)
       
        im_arr = np.zeros(image.size)
        i = 0
        j = 0
        Patchs = np.zeros(((self.patch_hight) * (self.patch_width), Y.shape[1]))
        W = np.zeros((self.Dh.shape[1], Y.shape[1]))
        for k in range(Y.shape[1]):
            #手順3:符号化係数の推定
            [w, result] = self.weight_learner(Dl, Y[:, k])

            #手順4:高解像画像パッチの生成
            Patchs[:, k] = np.dot(self.Dh, w)

        #手順5:パッチをつなぎ合わせて高解像画像を再構成
        [reconstructed_im, X0] = patch_vecs2gray_im(Patchs, image.size[0] * scaling , image.size[1] * scaling , self.patch_hight, self.patch_width, shift_num* scaling)
        L = get_downsampler_arr(image.size[0] * scaling, image.size[1] * scaling, image.size[0], image.size[1])

        #手順6:バックプロジェクションによる画像の補正
        x = np.array(self.back_projection_method(sci.matrix(X0.reshape(-1)).T, sci.matrix(np.array(image).reshape(-1)).T, L)[0])
        return float_arr2image(x.reshape((image.size[0] * scaling, image.size[1] * scaling)))


def main():
    test_img = Image.open('Lenna.png')
    scale = 3
    lam = 50
    base_num = 144
    patch_hight = 12
    patch_width = 12
    shift_num = 1
    c = 0.8
    mu = 0.9
    weight_learning_method = lambda W, y:ISTA(W, y, lam)
    back_projection_method = lambda x0, y, L:back_projection(x0, y, L, c, mu)
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.dictonary.Dic = np.load('dic_array.npy')
    im_dic.patch_hight = 12
    im_dic.patch_width = 12
    SR = super_resolution(im_dic, weight_learning_method, back_projection_method)
    result_img = SR.get_super_resolution(test_img, scale, shift_num)
    result_img.save("Lenna_SR.png")

超解像計算クラスsuper_resolutionは、画像パッチ辞書オブジェクト、係数推定関数オブジェクト、バックプロジェクション関数オブジェクトを生成時に指定し、get_super_resolutionメソッドで指定した辞書、関数に基づいた超解像画像の生成を行います。推定方法は上記で解説したとおりです。
画像パッチをつなぎ合わせて画像の再構成を行う関数patch_vecs2gray_imは、パッチのオーバーラップを許し、オーバーラップした部分は足し合わせた回数分割ることで平均をとっています。また、バックプロジェクションを行う関数back_projectionは、サイズが大きくスパースな行列Lを扱うので、scipy.matrixクラスで各引数を受け取る関数としています。ダウンサンプラー行列はget_downsampler_arr関数で求めています。

実験

まず辞書学習ですが、サイズ250×250の顔画像50枚から12×12のパッチ22000枚程度切り出して学習を行いました。また、辞書の基底数はパッチベクトルの次元と同じ144、学習時のはにしています。なお、辞書の学習繰り返し回数は8で打ち切っています（それ以上繰り返した学習辞書も作りましたが、生成される超解像画像の精度が逆に悪くなるという結果に・・・）。
超解像画像生成時は、、バックプロジェクション時のパラメータ、、パッチのずらし量で実験しています。実験に使用する画像データはおなじみのレナさんです。

83×83の画像を249×249の画像へ（拡大率3倍）

• 入力画像

• 左：通常の拡大　右：超解像アルゴリズムの出力

　

左の通常の拡大はPILのresizeで得られる画像です。少しわかりにくいですが、通常の拡大に比べて右の超解像画像は滑らかになっているのがわかると思います。

62×62の画像を248×248の画像へ（拡大率4倍）

• 入力画像

• 左：通常の拡大　右：超解像アルゴリズムの出力

　

もう少しわかりやすくするため、更に拡大率を上げた画像を生成してみました。通常の拡大に比べると滑らかになってはいますが、ブロック的なノイズが乗ってしまっています。アルゴリズム的な限界なのか、辞書やパラメータを洗練すればもっとうまくいくのか、その辺りの検証は今後の課題です。

まとめ

今回は近接勾配法の応用例として辞書学習および更にその応用である超解像画像生成をやってみました。
実際やってみると、パッチのつなぎ目が汚くなったりノイズまみれになったりと結構苦労しました。その過程で色々わかったこと・思ったことを下記にまとめます（その分野では当たり前のことばかりかも知れませんが・・・）。また、やり残したこと、今後やりたいことも併せて記します。

• やってみてわかったこと、思ったこと
  • 辞書学習時はパッチベクトルを平均、分散へ標準化した方がよい。
  • 画像再構成時はパッチ間の境界を滑らかにするため、元画像から切り出すパッチは領域を重複させて切り出し、再構成後のパッチを重ね合わせることで画像全体を再構成する方がよい。
  • 学習させすぎると精度が落ちるので、途中で学習を打ち切った方が超解像の結果がよい。
  • によって結構結果が変わる。係数推定結果の非零要素数がパッチベクトルの要素数と大体同じになる程度が妥当か。
  • バックプロジェクションはあまり効かない？（ダウンサンプリング以外にボケなどを考慮したすると有効か？）

• やり残したこと、今後やりたいこと
  • 学習用画像データのバリエーションを増やして実験。
  • 学習回数が増えると精度が落ちるのはなぜか数理的に考察。
  • 推定した画像の一部の要素が0~255の範囲を超えることがある。制約条件をつけられないか検討。
  • K-SVDなどの他の辞書学習手法の適用。
  • 独立成分分析やウェーブレットを用いた辞書行列の適用。
  • 他のダウンサンプリング辞書生成方法の適用。
  • 辞書学習以外の超解像手法との比較。

参考文献

• スパース表現の数理とその応用：スパースコーディング全般の解説。
• スパースコーディングを用いたマルチフレーム超解像：超解像手法の論文。
• 近傍パッチの情報に基づく重みを利用した スパース再構成による一枚超解像：画像パッチ再構成時近傍のパッチも利用する手法を提案。
• Image Super-Resolution via Sparse Representation：かなり有名な論文。今回使用した手法とは異なり、ダウンサンプリング辞書は辞書学習時に一緒に学習してる。
• 辞書学習のアルゴリズム：スライド資料。K-SVDなどの解説がわかりやすい。
• scikit-learnでSparse CodingとDictionary Learning -実践編：scikit-learnを用いた辞書学習の実装が書かれている。

*1:実務等で利用するときはscikit-learn辺りのライブラリを使うべきですが、今回は理解を深めるためnumpy、scipy、PILを使った実装を行っています。

*2:実はが正規分布に従ったランダムな行列でも、符号化としてはある程度うまくいくという研究もあります。ただ超解像やその他画像再構成への応用という観点ではどうなんでしょう？

近接勾配法とproximal operator

はじめに

前回、前々回とl1ノルム正則化項をつけた離散フーリエ変換により、スパースな解が得られること、そして基底ベクトルを並べた行列がユニタリ行列（直交行列）のため解析的に解けることを見てきました。それでは、が一般の行列の場合どうするか？今度こそ解析的に解けないので、反復法により最適解を求めます。しかし、目的関数にはl1ノルムという微分できない項があるため、最急降下法などのアルゴリズムは使えません。
このような微分できない項を含む最適化問題を効率的に解く手法の1つに近接勾配法（Proximal Gradient Method）があります。今回はこの近接勾配法と近接勾配法において重要な役割を果たすproximal operatorについて書いていきます。また、pythonで近接勾配法の実装を行いましたので、そのソースコードも載せておきます。

近接勾配法、proximal operatorとは

まずは、近接勾配法で解ける最適化問題の定義です。

ただし、は微分可能な凸関数、は微分可能とは限らない凸関数*1です。
l1ノルム正則化項付きの問題を式（1）に当てはめると、となります。

次にproximal operatorの説明です。とりあえず定義を見てみましょう。関数のproximal operatorは以下のように定義されます。

ここでは後述の更新式のステップ幅です。
このproximal operatorは何者か？というのは、として以下のような指示関数を考えると見えてきます。

ここでは閉凸集合です。この関数はが集合の要素なら0、の要素でないならに吹っ飛ばすという関数です。そして、この指示関数のproximal operatorは以下のようになります。

が集合の要素の場合、そのままとすればを最小化できます。一方で、が集合の要素ではない場合、集合の要素でとのユークリッド距離を最小とするものが選ばれます。すなわち、の集合への射影です（このような射影作用素をとして定義しています）。
よって、proximal operatorは射影作用素の関数への一般化と考えることができ、がによって定められる領域からはみ出していたら、その領域に入るように射影する作用素と見なせます。
proximal operatorの具体例を挙げると、がというl1ノルムであれば、前回述べたソフトしきい値作用素がproximal operatorになります。つまり、

です。そして、近接勾配法の更新式はこのproximal operatorを使って以下のようになります。

この更新式の収束性や正しさなどの説明は他の文献を参照いただくとして、直観的には次のような解釈ができると思います。
まず、は最急降下法の更新式です。また、上述のようにproximal operatorはによって定められる領域に射影する作用素です。従って、最急降下法の1ステップで得られる値がの領域からはみ出してたら、その領域へ射影するという操作を繰り返して最適解を探すアルゴリズムと見ることができます。
このように、関数の勾配と関数のproximal operatorがわかれば、近接勾配法により最適解が得られます。
なお、の決め方も収束のためには重要ですが、こちらについても別の文献を参照してください。

Moreau decompositionと双対ノルム

さて、ここまで見たように近接勾配法を使うためにはproximal operatorを求める必要があります。しかし、関数によっては、式（2）から直接的にproximal operatorを求めるのが難しいことがあります。例えば、とl2ノルムとした場合（2乗が付いていないことに注意）、式（2）を解こうとしてで微分しても、がかかる項が出てきてしまうので厄介です。
そこで便利なのが以下に示すMoreau decompositionが成り立つというproximal operatorの性質です。

ここで、はの共役関数のproximal operatorです。この共役関数とは下記の式で表される関数のことを言います。

よって、この共役関数のproximal operatorさえ分かれば、元の関数のproximal operatorも求められます。
では、というノルム関数の場合どうなるかを見ていきます。まずノルム関数の共役関数は以下の指示関数となります。

ここで、は双対ノルムと呼ばれるもので、以下の式で定義されます。

ということで、共役関数のproximal operatorは式（4）でみたように下記の射影作用素になります。

ただし、はで表される集合です。
つまり、をの範囲に射影する作用素がノルム関数の共役proximal operatorとなります。
おそらく式だけをザッと書いているので何故ノルムの共役関数が指示関数なのか？双対ノルムってそもそも何？どうやって求めるの？などの疑問が湧くかと思います。その辺りを少し補足しておきます。

双対ノルムについての補足

説明が前後しますが、後の説明のためにまずここで双対ノルムについて補足を書きます。そのための準備としてヘルダーの不等式という不等式を紹介します。
という2つの整数に、という関係があれば、

が成り立つというのがヘルダーの不等式です。実は、の双対ノルムはであり、その逆も成り立ちます。そのため、l2ノルムの双対ノルムはそのままl2ノルムですし、l1ノルムの双対ノルムはlノルムです。
このことはヘルダーの不等式から示すことができます。まずヘルダーの不等式より、任意のに対してが成り立ちます。lpノルムの双対ノルムはという条件下でのの最大値です。ここでをうまくとることにより、等式を成立させることができます。つまり、

が成立します。

何故、ノルム関数の共役関数が指示関数なのか？

まず、の時を考えます。式（8）の定義よりなわけですが、双対ノルムの定義式（10）と前述の前提より、かつとなるが存在することになります。これを満たすをとした場合、そのに定数をかけたも式（8）の最適解の候補として入ってきます。これはもちろん、とした場合も同様です。これらのことから、

となるため、のときはとなります。
一方で、の時は、式（12）のヘルダーの不等式よりなので、

です。前提より、なので、左辺（つまりは）が0以上になることはありません。ここで、という取り方もできるので、式（15）の左辺を0にすることができます。よって、のときはとなります。

結局のところproimal operatorはどうやって求めればよいか？

2つ方法があると思います。1つ目は式（2）から直接求める方法です。例えば、前回記したようの方法でl1ノルム正則化項のproximal operator式（5）を求めるというものです。
もう1つの方法としては、式（7）のMoreau decompositionの性質を使う方法です。特に正則化項がノルムなら、共役関数のproximal operatorはという超球への射影なので簡単に求められます。例えば、l2ノルムの双対ノルムはl2ノルムなので、

であり、求めたいproximal operatorは、

となります。

なお、式（2）にはステップ幅のが付いています。が付いていても基本的には上式のをに置き換えればよいだけですが、を含めたMoreau decompositionの一般形を書いておきましょう。

という形式になります。したがって、がノルム関数の場合は、

になります。ちなみに、多くの文献ではとしているので、として理論を展開している点に注意してください。

様々な正則化項のproximal operator

ここまでノルム関数が正則化項の時のproximal operator導出方法を見てきました。ここで、いくつかの正則化項に対するproximal operatorを記載しておきます。

l2ノルムの2乗正則化：

リッジ回帰等でよく使われる正則化項です。普通にproximal operatorの式を微分して0とおけば求められます。2乗しているので、以下のように式（16）とは異なる結果が得られます。

Elastic Net：

l1ノルム正則化とl2ノルム2乗正則化の混合です。l1ノルムは2つ変数が同程度の寄与度だったとしてもどちらかしか選ばれない（一方が0になってしまう）ケースが多いですが、Elastic Netではこの点を改善できるようです。このproximal operatorはl1ノルム正則化の時と同様、変数ごとに分離して求めることができます。

グループl1ノルム（グループに重複がない場合）：

変数がいくつかのグループに分けられるとき、それぞれのグループのlpノルムを取って、更にその和を取るノルムです（ ベクトルのうち、グループに所属する変数を並べたのが）。これを正則化項に用いるとグループごとにスパースになる、ならないが決まります。このグループl1ノルムを正則化項として使った回帰をgroup lassoと言います。
これはグループに重複がなく完全に分離ができるため、グループごとにlpノルムのproximal operatorを求めればよいです。p=2の時のグループのproximal operatorは下記の通りです。（式（16）とほぼ同じです。）

ちなみに、グループに重複がある場合も提案されていますが、ここでは割愛します（まだ追いきれてない）。

トレースノルム ：

最後は行列関数に関する正則化項です。トレースノルムはのトレースを正則化項とするものなので、の特異値をとすると、です。そのため、このトレースノルムを正則化項として使うと特異値が0になりやすく、結果低ランクな解が得られるようになります。
それではproximal operatorを見ていきます。の特異値分解をとすると、proximal operatorは以下のようになります。

特異値分解して特異値に対してソフトしきい値作要素を施した後に再構築すればよいわけです。

近接勾配法を実装

ということで式（1）を解く近接勾配法を実装しました。

近接勾配法を行う関数

ということで、まずは近接勾配法そのものの実装です。

def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        if (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        result = result_new
    return x_new, result_new

引数について説明します。
grad_fはにあたります。つまり、grad_fはxを引数として、その勾配を返す関数です。
proxはproximal operatorです。proxも関数で、入力変数vとステップ幅のgammaの2引数をとり、proximal operatorの計算結果を返します。
objectiveは目的関数です。objectiveはxを引数として、目的関数の計算結果を返す関数です。
init_xは初期値、tolは収束の許容誤差です。
やっていることは、式（6）の更新式をそのまま素直に実装しているだけです。

l1ノルム正則化項のproximal operator

次にproximal operatorの実装です。

def prox_norm1(v, gamma, lam):
    return soft_thresh(v, lam * gamma)

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

prox_norm1がl1ノルムのproximal operatorを計算する関数です。補助としてsoft_threshという関数を作っています。

lasso回帰問題を解く関数

以下は、というlasso回帰問題を解く関数です。なお、l1ノルム正則化項に対する近接勾配法はISTA（Iterative Shrinkage Thresholding Algorithm）と呼ばれるため、関数名をISTAとしています。

def ISTA(W, y, lam):
    objective = lambda x:np.sum(np.power(y - np.dot(W, x), 2))/2.0 + lam * np.sum(np.abs(x)) #(1)
    grad_f = lambda x: np.dot(W.T, np.dot(W, x) - y)  #(2)
    (u, l, v) = np.linalg.svd(W)   #(3)
    gamma = 1/max(l.real*l.real)   #(4)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam)  #(5)
    x_init = randn(W.shape[1])  #(6)
    (x_hat, result) = proximal_gradient(grad_f, prox, gamma, objective, x_init, 1e-5)  #(7)
    return x_hat, result

(1)ではラムダ式を使ってobjectiveをxを引数としてとり、上記の目的関数を計算する関数として定義しています。
また、であり、(2)はこの計算を行うgrad_fをxに引数とる関数として定義しています。
(3)、(4)はを定めています。最大特異値の逆数より小さい値をとして定めると収束するようですが、その辺りはまだ理解ができておりません。とりあえず、文献の通りに実装します。
(5)がproximal operatorです。prox_norm1関数はv,gamma,lamの3引数を受け取りますが、proximal_gradient関数に渡すproxはvとgammaの2引数のみ取る関数なので、ラムダ式を使って再定義しています。
(6)はランダムに初期値を決めて、(7)でproximal_gradient関数にて問題を解きます。

proximal operatorを別なものにする

上記ISTA関数のproximal operatorとobjectiveの部分だけ変えれば他の正則化にも対応可能です。例としてグループl1ノルム正則化（group lasso）の場合を示します。
まずは、グループl1ノルムのproximal operatorの実装です。

def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b)) * b

groupsは1～G（Gはグループ数）を要素に持つリストで、各変数がどのグループに属しているかを表します。後はグループごとにblock_soft_thresh関数にて、l2ノルムのproximal operatorを計算しています。
そして、上記のISTA関数内の(5)でproxとして、prox_group_normを指定すればgroup lassoになります。

    prox = lambda v, gamma:prox_group_norm(v, groups, gamma, lam) #(5)'

もし、Elastic Netを使いたかったらここにElastic Netのproximal operatorを持ってくればよいし、他のノルムの場合も同様です（本当はobjectiveも変更する必要がありますが上記では省略しています）。
また、を別なものにしたかったらを別途定義して、grad_fとすればよいです。
行列関数の場合も同様で、grad_fとproximal operator (あと、objective)さえ定義してあげたら機能します。そのため、proximal operatorとしてトレースノルムを指定してあげると、低ランク解を得ることができます。

まとめ

今回は近接勾配法とproximal operatorについて見てきました。上記の実装については次の機会に応用事例で実験したいと思います（ちなみに、疑似データで実験したところそれっぽく動いていました）。
また、今回は取り上げませんでしたが、この界隈ではISTAの高速版FISTAや拡張ラグランジュ法、その改良であるAlternating Direction Method of Multipliers (ADMM) など様々なアルゴリズムが提案されており、非常に奥が深いです。

応用例（2018/6/15追記）

近接勾配法の応用例の記事をいくつか書いたのでリンクを貼っておきます。

• スパースコーディング、辞書学習、超解像
• L1ノルム正則化項によるスパースロジスティック回帰
• トレースノルム正則化項付きロジスティック回帰による行列データの分類

参考文献

最後にこの記事を書くにあたって参考にした文献です（基本的に前回とほぼ同じです）。

• スパースモデリングのための凸最適化—近接勾配法による高速アルゴリズム：今回の件を調べるきっかけとなった解説論文
• Proximal Algorithms：共役関数、双対ノルム、Moreau decompositionについて非常に参考になった
• Optimization with Sparsity-Inducing Penalties：同上
• Convex Optimization：UCバークレーの講義資料？上２つの内容を補間するのに役にたった

また日本語の書籍で、proximal operator界隈が詳しく書かれているのは以下2冊かなと思います（他にもあったら教えてください）。

スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)

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

*1:もうひとつ、「単純」という条件も付きます。

l1ノルム正則化フーリエ変換を複素数のまま解く

はじめに

前回の記事では、離散フーリエ変換を係数をLasso回帰として求めスパースな解を得る方法を書きましたが、この際に複素数を取り扱いたくないため実数のみの問題に変換して、ソルバに解かせるということを行いました。今回は、複素数のまま問題を解くことを考え、複素数の最適解を得る方針で考えます。この方針で考察することで、以下のことが見えてきたので記事に残そうと思った次第です。

• 前回の記事では「この問題は解析的には解けず」などと書いているが、実はが離散フーリエの基底ベクトルを並べた行列であれば解析的に解けること
• 結局は、離散フーリエ変換してしきい値未満の要素を0にしていることと（ほぼ）等価であること
• 複素数としてのl1ノルムを考えると、別な解が得られること

ちなみに、自分は「この方法だとなくうまくいきそうだな」と思ったことはまず試してみて、後から理屈を考える（特にうまくいかなかったとき）というスタンスなので、理論的、数学的な厳密性はあまり考慮していません。したがって、本来理論的にはあり得ない式展開、考え方をしている場合もあるかもしれませんがご了承ください。

定式化

前回と同様、正則化項をつけた回帰問題として定式化します。

は元の時間軸の信号で、、が求めたい離散フーリエ係数であり、です。ここでの各要素は複素数です。は実信号なので実数ベクトルですが、虚部がすべて0の複素数と見なしてもよいでしょう。また、行列の各列ベクトルは

とします。前回の定義とほぼ一緒ですが、各要素を要素数で割っています。こうすると、はユニタリ行列になります。つまり、です(は随伴行列、は複素共役を表します)。
後は正則化項です。普通に考えると

（ただし、はの絶対値を表し ）としたいところですが、実はこの定義だと前回と同じ結果にはなりません。このことは後ほど説明するとして、今は

とします（はの実部、はの虚部を表す）。式（2）と式（3）の大きな違いは、実部と虚部が関連し合っているか、完全に独立なものとして扱うか、です。そのため、式（3）を正則化項として用いることは、実部と虚部を分けて解を求める前回の考え方と一致します。

解析解の導出と実装

当初は「スパースモデリングのための凸最適化」の解説論文を読んで近接勾配法で解こうなんて考えていましたが、解説を読んでいくと実は、がユニタリ行列であるため解析的に解が得られることがわかりました。
ということでその導出です。まず、解説に従うと、式（1）は以下のように式を変形できます。

このように式（1）の問題は各要素独立に最適解を選ぶ問題に変形されます。しかも、実部と虚部も分離しているので、こちらについても独立して最適解を選ぶことができます。ようするに、以下の最適化問題を各要素の実部、虚部それぞれ解けばよいわけです。

この解は次のような式で得られます。

解説論文と異なり式（1）および式（4）の第1項にがかかっていないため、最適解もにかけていることに注意してください。ここで式（5）のはソフトしきい値作用素と呼ばれる作用素です。
このソフトしきい値作用素を用いると、今回得たい式（1）の最適解は下記のように表せます。

さて、ここで通常の離散フーリエ変換の係数はで得られることを思い出すと、式（6）は通常の離散フーリエ変換の結果の実部、虚部それぞれに対してソフトしきい値作用素を作用させていることになります。
これは式（5）により離散フーリエ係数の絶対値が未満であれば0にすることに対応します。l1ノルム正則化項の離散フーリエ変換なんて言っていますが、紐解いていくと離散フーリエの結果に対ししきい値未満ならば0にするという割りとよくやる操作とほぼ同じになるわけですね（しきい値を超えている要素にはを加える点は異なりますが）。
ということで、以下は上式によるスパース離散フーリエ変換を行う関数の実装です。

def sparse_dft_with_thresh(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (soft_thresh((np.dot(W.T.conj(), y)).real, lam/2)
            +1j * soft_thresh((np.dot(W.T.conj(), y)).imag, lam/2))

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

sparse_dft_with_threshがスパースな離散フーリエ変換を行う関数で、soft_threshがソフトしきい値作用素を施す関数です。soft_threshを解説と合わせるために、引数で渡す時にlam/2をしています。

複素数l1ノルム正則化

ここで1つ疑問が残ります。正則化として式（2）のノルムを使うとどうなるか？です。これまでの議論から、式（2）のノルムを使った場合、式（1）は次のように変形できます。

この式を見ると要素同士は独立して最適解を選ぶことができますが、実部と虚部は独立させることができません（なので）。
この最適解を得るために、まず以下のような式を考えます。

この式（8）は近接作用素（proximal operator）と呼ばれるもので、近接勾配法などで重要な役割を持つ作用素です。今、式（7）について複素数を実部と虚部を並べた2次元ベクトルとして考えると（つまり、 ）、式（7）と式（8）の問題は（第1項のを除いて）等価になります（）。
ということで、式（8）を解くことを考えますが、実はこれはgroup lassoにおけるproximal operatorを求めることと等価です。group lassoのproximal operatorの導出やその基礎となる劣微分、凸解析、双対理論などは「ptimization with sparsity-inducing penalties」や「Proximal Algorithms」にまとまっています。正直、まだ理解できていない部分が多いですが、これらの資料に記載されているgroup lassoにおけるproximal operatorを使うと式（7）の解は下記のようになります（たぶん）。

式（6）では、実部、虚部それぞれでしきい値未満だと0になりましたが、こちらは複素数の絶対値がしきい値未満だと0になります。つまり、実部と虚部は運命共同体で、しきい値未満となる要素は実部も虚部一緒に0になり、その逆もまたしかりです。ということで実装です。

def group_lasso_prox(b, lam):
    x_hat = np.zeros(b.shape[0], "complex")
    x_hat[np.abs(b) >= lam] = (b[np.abs(b) >= lam] -
                                lam *  b[np.abs(b) >= lam]/np.abs(b)[np.abs(b) >= lam])
    return x_hat

def sparse_dft_with_complex_l1_norm(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (group_lasso_prox(np.dot(W.T.conj(), y), lam/2))

一応、実験した結果も載せておきます。実験には前回同様sin波とcos波とホワイトノイズを重ね合わせた疑似データを用います。

def main():
    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T))

    x1 = sparse_dft_with_thresh(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x1.real, "r")
    plt.hold(True)
    plt.plot(x1.imag, "b")

    x2 = sparse_dft_with_complex_l1_norm(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x2.real, "r")
    plt.hold(True)
    plt.plot(x2.imag, "b")

    plt.show()

実部と虚部を独立の変数として扱った正則化の結果

複素数l1ノルムによる正則化の結果

こうやって見ると、前者の方がスパース性という意味だと高いです。実部、虚部を運命共同体にする理由がない限り、前者のように独立して扱った方がよいように思えます。

最後に

ここまで見て感じたかと思いますが、結局のところ単なるしきい値での枝刈りに帰着するので、実装まで落としこんでしまうと正直あまり面白いものではありませんし、応用としてもあまり意味のあるものでもありません（これやるくらいなら、普通にFFTしてしきい値未満を0にする方が断然早いです）。
ただ、ここまで導出する過程で、ソフトしきい値作用素やproximal operator、劣微分、凸解析などなど最適化理論の色々なトピックスに触れることができました（まだ全然理解できていないですが）。
冒頭に書いた思いついたことはまず試してみるというスタンスは、試すために知らないトピックスを調べるきっかけになるというメリットがあります。やってみてダメなら何がダメか更に深追いするモチベーションにもなります。後はこういう試行錯誤は何より楽しいですね。ということで、今後は近接勾配法や凸解析の理論などの知識を深めていければと思います。