2016-05-05

ベルマンフォード法とダイカストラ法について

メモ。
AOJのこの問題で、入力数が大きいパターンでタイムアウトになっていた。

有向グラフが与えられる。このとき、頂点rから各頂点までの最短経路を求めよ。

ダイカストラではなくベルマンフォードで解いていたのが原因。この問題では、Vが最大100000、Eが最大500000と大きいので、O(VE)であるベルマンフォードだと制限時間内(3s)に処理が終わらない。追記：ベルマンフォードを適切に実装したら問題なかったようだ。コードも差し替えておく。
ダイカストラ： O(E logV)
ベルマンフォード： O(VE)

教訓：負の辺がない場合はダイカストラの利用を考える

ベルマンフォードでの解答（抜粋）

/// <summary>
/// find the shortest path from start to all destination
/// </summary>
public class BellmanFord
{
    public List<edge> Edge = new List<edge>();
    public int V;   // number of vertex

    public class edge
    {
        public int from, to, cost;
        public edge(int _from, int _to, int _cost)
        {
            from = _from;
            to = _to;
            cost = _cost;
        }
    }

    public BellmanFord(int v)
    {
        V = v;
    }

    /// <summary>
    ///  return shortestPath[V] represents distance from startIndex
    /// </summary>
    public int[] GetShortestPath(int startIndex)
    {
        int[] shortestPath = new int[V];
        var INF = Int32.MaxValue;

        for (int i = 0; i < V; i++) shortestPath[i] = INF;
        shortestPath[startIndex] = 0;
        while (true)
        {
            bool update = false;
            foreach (edge e in Edge)
            {
                if (shortestPath[e.from] != INF && shortestPath[e.to] > shortestPath[e.from] + e.cost)
                {
                    shortestPath[e.to] = shortestPath[e.from] + e.cost;
                    update = true;
                }
            }
            if (!update) break;
        }

        return shortestPath;
    }

    /// <summary>
    ///  return true if it has negative close loop
    /// </summary>
    public bool HasNegativeLoop()
    {
        int[] d = new int[V];
        for (int i = 0; i < V; i++)
        {
            foreach (edge e in Edge)
            {
                if (d[e.to] > d[e.from] + e.cost)
                {
                    d[e.to] = d[e.from] + e.cost;
                    if (i == V - 1) return true;
                }
            }
        }
        return false;
    }
}

ダイカストラでの解答（抜粋）

/// <summary>
/// Get min cost between two points
/// </summary>
public class Dijkstra
{
    private class Edge
    {
        public int T { get; set; }
        public int C { get; set; }
    }

    private int maxIndex = -1;
    private const int INF = Int32.MaxValue;

    private Dictionary<int, List<Edge>> _edgeArray = new Dictionary<int, List<Edge>>();

    /// <summary>
    /// Dijkstra, get min cost between two points
    /// should not contain negatie cost path
    /// </summary>
    /// <param name="size">max index of vertice</param>
    public Dijkstra(int size)
    {
        maxIndex = size + 1;
    }

    // Add a path with its cost
    public void AddPath(int s, int t, int cost)
    {
        if (!_edgeArray.ContainsKey(s)) _edgeArray.Add(s, new List<Edge>());
        _edgeArray[s].Add(new Edge() { T = t, C = cost });
    }

    //Get the min cost from s
    //return Int32.MaxValue if no path
    public int[] GetMinCost(int s)
    {
        int[] cost = new int[maxIndex];
        for (int i = 0; i < cost.Length; i++) cost[i] = INF;
        cost[s] = 0;

        var priorityQueue = new PriorityQueue<ComparablePair<int, int>>(_edgeArray.Sum( e => e.Value.Count()) + 1);
        priorityQueue.Push(new ComparablePair<int, int>(0, s));

        while (priorityQueue.Count() > 0)
        {
            var costDestinationPair = priorityQueue.Pop();
            if (cost[costDestinationPair.Item2] < costDestinationPair.Item1) continue;
            if (!_edgeArray.ContainsKey(costDestinationPair.Item2)) continue;

            foreach (var edge in _edgeArray[costDestinationPair.Item2])
            {
                var newCostToi = costDestinationPair.Item1 + edge.C;
                if (newCostToi < cost[edge.T])
                {
                    cost[edge.T] = newCostToi;
                    priorityQueue.Push(new ComparablePair<int, int>(newCostToi, edge.T));
                }
            }
        }

        return cost;
    }
}

PriorityQueueを自作する必要があるのと、Dijkstra内を単純な配列にするとメモリが足りなくなるのが面倒だった。
こういう問題だとC++を使ったほうが楽そうだ。

2016-05-03

Topcoder Open 2016 Marathon Match(Round 1)の参加日記

C# プログラミングコンテストマラソンマッチ Topcoder

結果は55位/136人と振るわなかったが、今後のために記録を残しておく。

問題概要
https://community.topcoder.com/longcontest/?module=ViewProblemStatement&rd=16702&pm=10729

平面上に、木の幹がランダムに点として存在する
それぞれの幹からは、ランダムに根が生えている
根の終端からまた根もランダムに生えている
このとき、すべての木の幹が別のグループになるよう、平面を切りたい
切断線は無限大の長さになる（＝ある点からある点まで、とは切れない）
なるべく切られてしまう根を減らして、切断線を求めよ

実装内容
雑な貪欲。すべての幹の組み合わせを求め、幹間の距離が短い順に

幹間に線を引き、30等分くらいで切断開始点を決める
それぞれの開始点から、数度刻みでぐるっと一周で切断面を作ってみる
一番、切ってしまった根の長さ合計が少なかったものを採用

を繰り返しただけ。
あとは、すでに切られている幹間は無視したり、根が多すぎるときは単純化して計算量を減らしたりした。

感想と反省
開始3日くらいで実装・提出してから、終了までほぼ改善できずに終わってしまった。
7日目くらいで焼きなましをやってみようとしたのだが、

体調などのせいで時間が取れない（そして睡眠不足の悪循環）
データ構造が良くなくほぼ作り直し

のせいで日の目を見ずにタイムアップ。
まあ3日で何とか黄色を維持できるコードが書けたし、独自Visualizerも悪くない感じにできたので、上達はしているのかもしれない。

教訓

時間がない時こそ、実装を急がず戦略を立てて行動する
プロトタイプといえど、データ構造はなるべく柔軟性を持たせる

2016-03-30

k平均法（k-meansクラスタリング）実装

C# 機械学習

Courseraの機械学習コースでk平均法を学んだので、理解を深めるために実装してみた。
わざわざ日記にするほど難しくはないが、メモということで。ソースコードはここ。

メインは次のGetKMeans_withCost()で、tryNumの数だけk平均法を試して、コスト(最寄りのセントロイドまでの距離の合計)が最小になった結果を返している。返り値は各セントロイドの位置と最終的なコストのタプル。多次元対応。

入力はクラスタ数kとMatrix型の変数x。xは行が１サンプル、列が特徴を表す。

アルゴリズムは単純で、　　

xからk個ランダムに選び出し、この位置をセントロイドの初期値にする
それぞれのxに、最も近いセントロイドを割り当てる
セントロイドの位置を、割り当てられているxの平均値にアップデート

を繰り返しているだけ。セントロイドの位置が変わらなくなるまでループする。

static Tuple<int[], double> GetKMeans_withCost(int tryNum, int K, Matrix x)
{
    var r = new Random();

    int[] ret = null;
    var centroids = new Matrix[K];

    var minCost = double.MaxValue;
    while (tryNum > 0)
    {
        Initialize(r, x, centroids);

        var J = 0.0;
        int[] c = null;
        while (true)
        {
            var cTuple = AssignC(centroids, x);
            c = cTuple.Item1;

            var prevJ = J;
            J = cTuple.Item2;
            if (DoubleUtil.Equals(prevJ, J)) break;

            UpdateCentroids(centroids, x, c);
        }

        if (J < minCost || ret == null)
        {
            ret = c;
            minCost = J;
        }

        tryNum--;
    }

    return new Tuple<int[], double>(ret, minCost);
}

こちらはkの値をいろいろ試してみるコード。kの候補値(Ks[])ごとの結果を返す。
うまい具合にエルボー（これ以上増やしてもあまりコストが下がらないようなkの値）があれば、それを使うのが良いらしい。

public static double[] FindK(int tryNum, int[] Ks, Matrix x)
{
    var ret = new double[Ks.Length];

    for (int i = 0; i < Ks.Length; i++)
    {
        var result = GetKMeans_withCost(tryNum, Ks[i], x);
        ret[i] = result.Item2;
    }

    return ret;
}

2016-03-25

ニューラルネットワーク実装

C# 機械学習

Courseraの機械学習コースでニューラルネットワークを学んだので、習作としてC#で実装してみた。

多層パーセプトロン対応。Classification専用（シグモイド関数をベタ書きしてるので）。
ソースコードはここにおいてある。

実装の概略
ネットワーク各層の関数を、行列型の配列Thetas[]で保持している。

public class NeuralNetwork
{
    public Matrix[] Thetas { get; set; }
    public int NumLayers { get; set; }
    public int[] NumNeurons { get; set; }

これに、（RandomizeThetas()で初期化したあと）関数Learn()に学習用データとパラメタを渡せば最急降下法で学習を開始する。

public void Learn(Matrix[] input, Matrix[] output, double alpha, double lambda, int maxItr)
{
    while (maxItr > 0)
    {
        var grad = GetGrad(Thetas, input, output, input.Length, lambda);
        UpdateThetas(input.GetLength(0), grad, alpha, lambda, input, output);
        maxItr--;
    }
}

上記のGetGrad()が、各ノードの偏微分を返す関数で、それを使ってUpdateThetas()内でThetas[]を更新する。

GetGrad()の実装も教科書通り。

public Matrix[] GetGrad(Matrix[] thetas, Matrix[] input, Matrix[] output, int m, double lambda)
{
    var L_deltas = InitializeLDeltas();

    for (int idx = 0; idx < input.GetLength(0); idx++)
    {
        Matrix[] z = null;
        var values = ForwardProp(input[idx], thetas, ref z);
        BackProp(L_deltas, output[idx], values, z);
    }

    var grad = GetGrad(thetas, L_deltas, input.Length, lambda);
    return grad;
}

ForwardProp()は名前の通りフォワードプロパゲーションを、BackProp()はバックプロパゲーションをする関数。
ForwardProp()は各層の中間値を返すので、それを使ってBackProp()で後ろの層から変化量もとめ、L_deltasを更新している。
全入力のループが終わったら、GetGrad()がL_deltasを使って偏微分を求める。

偏微分を求める実装がややこしいのだが、コースで推奨されているとおりに、近似値を求めて比較する方法でデバッグした。
${\theta}$ -eと ${\theta}$ +e（eは小さい数）のそれぞれでコストを求め、その傾きを使うというもの。

Matrix GetNumericGradient(NeuralNetwork nn, Matrix[] thetas, int idx, Matrix[] input, Matrix[] output)
{
    var t = thetas.ToArray();

    const double e = 1e-4;
    var ret = new Matrix(t[idx].RowNum, t[idx].ColNum);
    var diffMat = new Matrix(t[idx].RowNum, t[idx].ColNum);

    for (int row = 0; row < diffMat.RowNum; row++)
    {
        for (int col = 0; col < diffMat.ColNum; col++)
        {
            var orgValue = t[idx][row, col];

            t[idx][row, col] = orgValue - e;
            var loss1 = nn.J(thetas, input, output);

            thetas[idx][row, col] = orgValue + e;
            var loss2 = nn.J(thetas, input, output);

            ret[row, col] = (loss2 - loss1) / (2.0 * e);

            thetas[idx][row, col] = orgValue;
        }
    }

    return ret;
}

var numgrads = GetNumericGradients(nn, nn.Thetas, x, y);
var grads = nn.GetGrad(new Matrix[] { theta1, theta2 }, x, y, x.Length, lambda);

2通りの結果がほぼ同じなので、実装が合っていることを確認できる。

使ってみる

こちらのサイトにある手書き数字の学習と認識に挑戦してみる。

static void NeuralNetworkTest_Example()
{
    Console.WriteLine("**Neural Network test (test data from: http://archive.ics.uci.edu/ml/machine-learning-databases/pendigits/)**");

    var nn = new NeuralNetwork(3, new int[] { 16, 16, 10 });
    nn.RandomizeThetas();

    //train
    var trainingDataTuple = GetNNData(TraininData.Data);
    var inputs = trainingDataTuple.Item1;
    var outputs = trainingDataTuple.Item2;

    //test
    var testDataTuple = GetNNData(TestData.Data);
    var test_inputs = testDataTuple.Item1;
    var test_outputs = testDataTuple.Item2;

    var alphas = new double[] { 5, 7 };
    var lambdas = new double[] { 0.03, 0.1, 0.3, 1.0 };

    var parameters = nn.FindParameters(alphas, lambdas, inputs, outputs, test_inputs, test_outputs, 300);
            
    Console.WriteLine("Chosen params - alpha:{0}, lambda:{1}", parameters.Item1, parameters.Item2);

    nn.Learn(inputs, outputs, parameters.Item1, parameters.Item2, 300);

    var result = nn.GetResult(test_inputs, test_outputs);
    var err = result.Item1;
    var total = result.Item2;

    Console.WriteLine(string.Format("Total:{0}, Error:{1}, Success Rate: {2}", total, err, (double)(total - err) / (double)total));            
}

トレーニング用とテスト用のデータがすでに用意されていたので、そのまま利用した。
好ましいパラメタを選ぶために、本来ならクロスバリデーション用として別のデータを用意すべきだが、今回はテスト用を流用した。

3層ネットワークで、イテレーションは適当に300回。
合計5回の学習で1時間以上かかってしまった。

f:id:yambe2002:20160324232921p:plain

認識率は約90%だった。何も工夫しないとこんなものだろうか。または、イテレーション回数が少なくて収束しきっていないのかもしれない。
このあたりは、別にちゃんとしたライブラリを使っていろいろ試してみようと思う。

一から自分で組んでみてかなり勉強になったが、本格的に実装するならライブラリを使わないと難しいと改めて実感。（数学的に最適化されてないので実行がかなり遅いし、デバッグも難しい）

2016-03-03

HackerRank Week of Code - 19 参加日記

プログラミングコンテスト C#

HackerRankのWeek of Code - 19に参加したのでそのメモ。
結果は346位/3176人だった。

Fix the Cycles

4つのノードA,B,C,Dをもつ有向グラフを考える。D→Aのエッジをa、A→Bのエッジをb、B→Cのエッジをc、C→Dのエッジをd、A→Cのエッジをe、B→Dのエッジをfとする。
a~fの重みが与えられたとき、負の閉路をなくすためには、最小でいくら重みを足せばよいか。ただし、1つのエッジにのみ重みを足すことができる。
不可能な場合は-1を返こと。
制約：a~fの重みの値は-20から+20

全エッジに対し、足す重み全て(0~80)をチェックするだけでよい。
ちなみに、すべてのループがエッジaを通るため、答えが-1になることはない。

static bool Success(int a, int b, int c, int d, int e, int f, int val, int idx)
{
    switch (idx)
    {
        case 0: a += val; break;
        case 1: b += val; break;
        case 2: c += val; break;
        case 3: d += val; break;
        case 4: e += val; break;
        case 5: f += val; break;
    }

    return (b + f + a >= 0) && (b + c + d + a >= 0) && (e + d + a >= 0);
}

public static int Solve(int a, int b, int c, int d, int e, int f)
{
    for (int i = 0; i <= 80; i++)
    {
        for (int j = 0; j < 6; j++)
        {
            if (Success(a, b, c, d, e, f, i, j)) return i;
        }
    }

    return -1;  //should not be here since all loops have edge 'a'
}

Two Robots

無限個のキャンディーが入ったM個のコンテナが、１列に配置されている。そして2つのロボットがある。ロボットはコンテナ間を移動して、キャンディーを運ぶことができる。
キャンディー運ぶ順番（あるコンテナから別のコンテナへ）が指定されたとき、トータルで最小のロボット移動距離を求めよ（隣り合うコンテナ間の距離を距離1とする）。

キャンディーを運ぶ順番は次のように指定される。
N　#運ぶ回数
Ma1 Mb1　#コンテナMa1からMb1に運ぶ。必ず違う値
・・・
Man　Mbn

制約： ${ 1 \leq M \leq 1000, \, 1 \leq N \leq 1000 }$

次のDPを設定し、手順を後ろからループしていく。
DP[現在のロボットAの位置][現在のロボットBの位置] ＝以降の移動距離の最小（Aの位置＜Bの位置）

単純に実装するとメモリと実行時間が足りなくなるので、(1)DP配列の再利用(2)あり得ない条件は無視して枝刈り、を行った。

public class Solver
{
    public int n = 0;
    public int m = 0;
    public Tuple<int, int>[] queries;

    public int Solve()
    {
        int[][,] dp_all = new int[][,] { new int[m + 2, m + 2], new int[m + 2, m + 2] };

        int[] minIdx = new int[m + 1];
        for (int idx = 0; idx < m + 1; idx++)
        {
            minIdx[idx] = Int32.MaxValue;
        }

        for (int idx = 0; idx < queries.Length; idx++)
        {
            var val = queries[idx].Item2;
            if (minIdx[val] == Int32.MaxValue || idx < minIdx[val])
            {
                minIdx[val] = idx;
            }
        }

        int dpIdx = 0;
        for (int idx = queries.Length - 1; idx >= 0; idx--)
        {
            var q = queries[idx];

            var pos1 = idx == 0 ? 0 : queries[idx-1].Item2;

            for (int pos2 = 0; pos2 < m + 1; pos2++)
            {
                if (pos2 != 0 && idx < minIdx[pos2]) continue;  //cannot be

                var dp = dp_all[dpIdx];
                var dp_next = dp_all[dpIdx == 0 ? 1 : 0];

                dp[Math.Min(pos1, pos2), Math.Max(pos1, pos2)] = Math.Abs(q.Item1 - q.Item2);

                var movePos1 = (pos1 == 0 ? 0 : Math.Abs(pos1 - q.Item1))
                    + ((idx == queries.Length - 1) ? 0 : dp_next[Math.Min(q.Item2, pos2), Math.Max(q.Item2, pos2)]);

                var movePos2 = (pos2 == 0 ? 0 : Math.Abs(pos2 - q.Item1))
                    + ((idx == queries.Length - 1) ? 0 : dp_next[Math.Min(q.Item2, pos1), Math.Max(q.Item2, pos1)]);

                dp[Math.Min(pos1, pos2), Math.Max(pos1, pos2)] += Math.Min(movePos1, movePos2);
            }

            dpIdx = dpIdx == 0 ? 1 : 0;
        }

        return dp_all[dpIdx == 0 ? 1 : 0][0, 0];
    }
}

Scalar Products

次の数列を考える。
${ a_0 = 0, a_1 = C, a_{i+2} = (a_{i+1} + a_i) \% M,\,where\,0 \leq i \leq 2n}$
このとき、ベクトル ${v_1, v_2, ..., v_n}$ を次のように規定する。
${v_1 = (a_2, a_3), v_2 = (a_4, a_5), ..., v_n = (a_{2n}, a_{2n+1})}$
Sを ${v_i, v_j \, (v\leq i, j \leq n)}$ のスカラ積の集合としたとき、Sの剰余の個数をMでModしたものを求めよ。
制約： ${1 \leq C \leq 10^9, 1 \leq M \leq 10^9, 1 \leq n \leq 3 \times 10^5}$

${a_0 = 0,\, a_1 = C }$ で開始するフィボナッチ数列の変種と考えればよい。

フィボナッチ数列の特徴
${ a_x \times a_y + a_{x+1} \times a_{y+1} = a_{x+y+1} }$
を利用すると、この数列に対しては
${ a_x \times a_y + a_{x+1} \times a_{y+1} = C \times a_{x+y+1} }$
が導き出せる。

よって、 ${ C \times a_{x+y+1} }$ の部分を見ていけばよい。

public class Solver
{
    public long Solve(long c, long m, long n)
    {
        c %= m;
        var resultSet = new HashSet<long>();

        var prev = c;
        long prev2 = 0;
        for (int i = 2; i <= 4 * n - 1; i++)
        {
            var val = (prev2 + prev) % m;

            if (i % 2 == 1 && i >= 7) resultSet.Add((val * c) % m);

            prev2 = prev;
            prev = val;
        }

        return resultSet.Count() % m;
    }
}

これは解けなかったのが悔しかった。
「フィボナッチ数列」で検索してみるべきだった。

Coolguy and Two Subsequences

一次元配列A[]（長さN、Indexは1から始まる）が与えられたとき、次の疑似コードで得られるansを答えよ。
//f(a, b)はA[a, b]の最小値を返す
ans = 0
for a -> [1, N]
    for b -> [a, N]
        for c -> [b + 1, N]
            for d -> [c, N]
                ans = ans + min(f(a, b), f(c, d))

ans = ans % (10^9 + 1)
制約： ${ 1 \leq N \leq 2 \times 10^5, \, 1 \leq A_i \leq 10^9 }$

まった手が出なかった。
Editorialは理解できるのだが、最後が「これをSegmentTreeを使って解けばよい」で終わっていて困る・・・。このレベルなら自明であるべきなのだろうか。

仕方ないので、SegmentTreeの集中特訓をしてから、後でもう一度解いてみることにする。

Editorialの内容

配列Aにあるすべての数字xについて、 ${ min( f(a, \, b), \, f(c,\,d) ) = x }$ になる回数を求めればよい。
これを ${ cnt(x) }$ とすると、最終的な解答は ${ sum( x \times cnt(x) ) }$ になる。

ここで、
　 ${g(x) = min( f(a, \, b), \, f(c,\,d) ) \geq x になる回数}$
と定義すると、
${cnt(x) = cnt(x) - cnt(x + 1)}$ となる。

${g(x)}$ を求めるためには、別の配列Bを用意するとよい。
BはAと同じ長さで、 ${a_i \geq x}$ なら ${b_i = 1}$ 、それ以外なら ${b_i = 0}$ のもの。

すると、 ${g(x)}$ は
　 ${x1 \leq y1}$ かつ
　 ${y1 + 1 \leq x2}$ かつ ${ x2 \leq y2 }$
になる ${(x1,\,y1,\,x2,\,y2)}$ の組み合わせのうち、
${ B[x1:y1]}$ と ${ B[x2:y2] }$ の範囲内すべてが1になるものの数と等しくなる。

${g(x)}$ はSegmentTreeを使えば解くことができる。

2016-02-13

線形回帰とロジスティック回帰の実装

C# 機械学習

Courseraの機械学習コースで多変量線形回帰（最小二乗法）とロジスティック回帰を学んだので、理解を深めるためにC#で実装してみた。
特にロジスティック回帰は、コースで最急降下法の実装をやらなかった（Octave組み込みのfminuc()を使う）ため、練習として組んでみたかった。

行列クラス
まずは行列クラスを作成。
topcoder_template/MyLibrary_Matrix.cs at master · yambe2002/topcoder_template · GitHub

演算子のオーバーロードと、平均・標準偏差や合計を求める関数などもあり、ずいぶん大きくなってしまった。

逆行列を求める関数Inverse()は掃き出し法を使っているので、一般化逆行列を返すわけではない。
NormalizeFeature()は特徴量を標準化させるもので、行列クラスに入れるのは微妙だったかもしれない。

線形回帰

最急降下法（Gradient Descent）

各特徴量ごと、コストをその特徴量(j)にかかる係数 ${\theta_{j}}$ で偏微分し、コストが下がる方向にその係数 ${\theta_{j}}$ を変化させていく。返り値は ${\theta}$ のベクトルとコスト履歴のタプル。

public static Tuple<Matrix, double[]> GetResult_GradientDescent(Matrix x, Matrix y, double alpha, double lambda, int maxItr, bool logCosts=true)
{
    var theta = new Matrix(x.ColNum, 1);
    var costs = new List<double>();
    while (maxItr > 0)
    {
        var h = x * theta;

        if (logCosts)
        {
            var j = ((1.0 / (2.0 * (double)x.RowNum)) * ((h - y) ^ 2)).Sum();
            costs.Add(j);
        }
        var diff = (alpha / (double)x.RowNum) * (x.Transpose() * (h - y)) - ((alpha * lambda) / (double)x.RowNum) * theta;
        diff[0, 0] += ((alpha * lambda) / (double)x.RowNum) * theta[0, 0];
        theta = theta - diff;

        maxItr--;
    }

    return new Tuple<Matrix,double[]>(theta, costs.ToArray());
}

行列xがトレーニングデータ（データ数 ${ \times }$ 特徴数）、行列yが解のベクトル、alphaが降下パラメータ、lambdaが標準化係数。

正規方程式（Normal Equation）

一般解として計算する方法。 ${\theta}$ を求める係数のベクトルとすると、解は
${ \theta = (X ^{T} X) ^{-1} X^{T} Xy }$ 　になる。
なぜこれで求まるかは、線形回帰の Normal Equation（正規方程式）について - Qiitaが参考になった。
とても簡単で結果も厳密だが、データサイズが大きくなると計算に時間がかかる。

public static Tuple<Matrix, double[]> GetResult_NormalEquation(Matrix x, Matrix y, double lambda, bool logCosts = true)
{
    var xT = x.Transpose();
    var l = Matrix.Eye(x.ColNum);
    l[0, 0] = 0;

    var ret = (xT * x + lambda * l).Inverse() * xT * y;

    var costs = new List<double>();
    if (logCosts)
    {
        var h = x * ret;
        var j = ((1.0 / (2.0 * (double)x.RowNum)) * ((h - y) ^ 2)).Sum();
        costs.Add(j);
    }

    return new Tuple<Matrix,double[]>(ret, costs.ToArray());
}

ロジスティック回帰

線形回帰の最急降下法と似ているが、 ${h = g (\theta ^T x)}$ になる。
${g(z)}$ はシグモイド関数
${ g(z) = \frac{1}{ 1 + e^{-z} } }$
で、解が必ず0~1になる。
コストは ${ y = 0 }$ の場合と ${ y = 1 }$ の場合で場合わけして求める。

public static Tuple<Matrix, double[]> GetResult(Matrix x, Matrix y, double alpha, double lambda, int maxItr, bool logCosts=true)
{
    var theta = new Matrix(x.ColNum, 1);
    var costs = new List<double>();
    while (maxItr > 0)
    {
        var h = GetSigmoid(x * theta);
        var diff = h - y;

        if (logCosts)
        {
            var wktheta = theta.Clone();
            wktheta[0, 0] = 0;
            var j = ( 1.0 / (double)x.RowNum ) *
                    ( ( -1.0 * (y.Transpose() * h.Log())) - (1.0 - y).Transpose() * (1.0 - h).Log())
                    + ( lambda / ( 2.0 * (double)x.RowNum) ) * (wktheta^2).Sum();
            costs.Add(j[0, 0]);
        }

        var grad = ((alpha / (double)x.RowNum) * (diff.Transpose() * x).Transpose()) + ((alpha * lambda) / (double)x.RowNum) * theta;
        grad[0, 0] -= ((alpha * lambda) / (double)x.RowNum) * theta[0, 0];

        theta -= grad;

        maxItr--;
    }

    return new Tuple<Matrix, double[]>(theta, costs.ToArray());
}

static Matrix GetSigmoid(Matrix z)
{
    var a = Math.E ^ z;
    var b = 1.0 / a;
    var c = 1.0 + b;
    var d = 1.0 / c;

    var ret = 1.0 / ( 1.0 + (1.0 / (Math.E ^ ( z ))) );
    return ret;
}

検証
f:id:yambe2002:20160303041844p:plain
コースと同じデータを使って検証した。
Octaveの結果とほぼ一致。いずれも、それらしく動いているようだ。

メモ：当初、あるテストケースで最急降下法がうまく収束しなくて悩んでいたが、標準化を忘れていたのが原因だった。

2016-01-28

二分探索の練習

C# プログラミングコンテスト

二分探索の練習。
問題は「プログラミングコンテスト」「二分探索」あたりで適当にググったもの。

- UpperBoundとLowerBound
CのSTD関数をC#で実装。

public static int LowerBound(int[] ar, int val)
{
    var lb = -1;
    var ub = ar.Length;

    while (ub - lb > 1)
    {
        var mid = (lb + ub) / 2;
        if (ar[mid] >= val)
        {
            ub = mid;
        }
        else
        {
            lb = mid;
        }
    }

    return ub;
}

public static int UpperBound(int[] ar, int val)
{
    var lb = -1;
    var ub = ar.Length;

    while (ub - lb > 1)
    {
        var mid = (lb + ub) / 2;
        if (ar[mid] <= val)
        {
            lb = mid;
        }
        else
        {
            ub = mid;
        }
    }

    return lb + 1;
}

- Topcoder SRM676 Div1Easy(WaterTank)
https://community.topcoder.com/stat?c=problem_statement&pm=14019

Cリットル入る水がめと、配列tとxが与えられる。水がめには、t[i]の時間、x[i]リットルの水が注ぎこまれる。水がめが溢れないようになるべく小さい穴をあけたい。
この穴の単位時間あたりの水の排出量を求める。
配列の大きさは1～50、tとxは1～1,000,000、Cは1～10^9

LowerBoundと同じように二分探索する。
最初、1sごとにループするコードを書いてタイムアウトになってしまった。水の増え方が線形なので、t[i]単位でまとめて計算して問題ない。

public class WaterTank {
    double EPS = 1e-8;

    int[] t;
    int[] x;
    double C;

    bool NoOverflow(double rate)
    {
        var tank = 0.0;
        for (int n = 0; n < t.Length; n++)
        {
            //これだと間に合わない！
            //線形に増えるので1sごとにループする必要はない
            //var water = (double)x[n];
            //for (int sec = 0; sec < t[n]; sec++)
            //{
            //    tank += water;
            //    tank -= rate;
            //    if (tank > C) return false;
            //    if (tank < 0) tank = 0.0;
            //}

            tank += (double)x[n] * t[n];    //キャストを忘れない！（忘れてもサンプルテストは通る）
            tank -= rate * t[n];
            if (tank > C) return false;
            if (tank < 0) tank = 0.0;
        }

        return true;
    }

    public double minOutputRate(int[] pt, int[] px, int pC)
    {
        t = pt;
        x = px;
        C = (double)pC;

        double lb = 0.0, ub = 1000001.0;

        while (ub - lb > EPS)
        {
            var mid = (lb + ub) / 2;
            if (NoOverflow(mid))
                ub = mid;
            else
                lb = mid;
        }

        return ub;
    }
}

- Topcoder SRM456 Div2Hard(CutSticks)
https://community.topcoder.com/stat?c=problem_statement&pm=10563&rd=13909

いろいろな長さの棒sticks[]がある。そして最大でC回、棒を取り出して2つに割ることができる（割ってできた棒をさらに割ることもできる）。
最終的にできた棒を、増加しない順にならべて、K番目の棒を取り出す。このとき、その取り出した棒の最大長を求める。
なお、最低でもK本の棒ができるように、割っていく必要がある。
sticks[]の長さ：1～50
sticks[]の各要素：1～1000000000
C：1～1000000000
K:1～sticks[]の長さ+C

「ある長さにできるかどうか」を判定しながら二分探索する。
終了条件を上限と下限の差にすると無限ループになってしまうので、十分大きな固定回数のループにする。

public class CutSticks {
    double EPS = 1e-10;

    bool CanMake(double len, double[] sticks, int C, int K)
    {
        var cnt = 0;
        foreach (var stick in sticks)
        {
            if (stick - len > EPS)
            {
                var maxDevideNum = (int)(stick / len) - 1;
                var devideNum = Math.Min(C, maxDevideNum);

                cnt += (devideNum + 1);
                C -= devideNum;

                if (cnt >= K) return true;
            }
        }

        return false;
    }

    public double maxKth(int[] sticks, int C, int K) {
        double lb = 0;
        double ub = 1000000001;

        var dsticks = sticks.Select(a => (double)a).OrderByDescending(b => b).ToArray();

        //while(ub - lb > EPS)  //これだと (new int[] { 1000000000 }, 1000000000, 1) のケースで終了しない
        for (int i = 0; i < 500; i++)
        {
            var mid = (lb + ub) / 2;

            if (CanMake(mid, dsticks, C, K))
                lb = mid;
            else
                ub = mid;
        }

        return lb;
    }
}

- Topcoder SRM636 Div2Hard(ChocolateDividingHard)
https://community.topcoder.com/stat?c=problem_statement&pm=13388&rd=16079

チョコレートの板String[] chocolateがある。
chocolateの各char要素は、そのセルのクオリティを表している('0'～'9')。
縦に４つ、横に４つに割ったとき、１つの塊のクオリティの合計の最小値を最大にしたい。この値を答えよ。
chocolate[]の縦横の要素数は4～75

最大のクオリティを求める二分探索を行う。
その中で縦（または横）全パターンのループを回し、パターンごとにクオリティを満足する最小の横Indexを二分探索で求めていく。
横に4分割できればOK。

二分探索の中で二分探索をしていて、はまるとデバッグが大変になる問題だと思う。

public class ChocolateDividingHard {
    int[,] sums;
    int rowNum;
    int colNum;
    int sum;

    //row1 < row2, col1 < col2
    //(row1 + 1, col1 + 1) to (row2, col2)
    int GetArea(int row1, int row2, int col1, int col2)
    {
        var ret = sums[row2, col2] -
            (col1 >= 0 ? sums[row2, col1] : 0) -
            (row1 >= 0 ? sums[row1, col2] : 0) +
            (col1 >= 0 && row1 >= 0 ? sums[row1, col1] : 0);

        return ret;
    }

    bool CanDevide(int col1, int col2, int size, int[] rows)
    {
        for (int rowIdx = 0; rowIdx < rows.Length; rowIdx++)
        {
            var row1 = rowIdx > 0 ? rows[rowIdx - 1] : -1;
            var row2 = rows[rowIdx];
            if (GetArea(row1, row2, col1, col2) < size) return false;
        }
        return true;
    }

    int GetMinColIdx(int startIdx, int size, int[] rows)
    {
        //binary search (l, u]
        int l = startIdx - 1;   //ここを l = startIdx にすると、startIdxが答えの時に間違った解を返す！
        int u = colNum;

        while (l < u - 1)
        {
            var mid = (l + u) / 2;
            if (CanDevide(startIdx - 1, mid, size, rows))
                u = mid;
            else
                l = mid;
        }

        return (u == colNum) ? -1 : u;
    }

    bool CanDevide(int size, int[] rows)
    {
        var col = -1;
        for (int i = 0; i < 4; i++)
        {
            col = GetMinColIdx(col + 1, size, rows);
            if (col == -1) return false;
        }
        return true;
    }

    bool CanDevide(int size)
    {
        for (int row1 = 0; row1 < rowNum; row1++)
            for (int row2 = row1 + 1; row2 < rowNum; row2++)
                for (int row3 = row2 + 1; row3 < rowNum; row3++)
                    if (CanDevide(size, new int[] { row1, row2, row3, rowNum - 1 }))
                        return true;
        return false;
    }

    int findBest()
    {
        //binary search [l, u)
        int l = 0;
        int u = sum + 1;

        while (u > l + 1)
        {
            var mid = (l + u) / 2;

            if (CanDevide(mid))
                l = mid;
            else
                u = mid;
        }

        return l;
    }

    public int findBest(string[] chocolate) {
        rowNum = chocolate.Length;
        colNum = chocolate[0].Length;

        sums = new int[rowNum, colNum];
        for (int row = 0; row < rowNum; row++)
        {
            for (int col = 0; col < colNum; col++)
            {
                sums[row, col] = 
                    (row > 0 ? sums[row - 1, col] : 0) +
                    (col > 0 ? sums[row, col - 1] : 0) +
                    ((int)(chocolate[row][col] - '0')) -
                    (row > 0 && col > 0 ? sums[row - 1, col - 1] : 0);
                sum += (int)(chocolate[row][col] - '0');
            }
        }

        return findBest();
    }
}