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

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

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

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

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

線形回帰

• 最急降下法（Gradient Descent）

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

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がトレーニングデータ（データ数特徴数）、行列yが解のベクトル、alphaが降下パラメータ、lambdaが標準化係数。

• 正規方程式（Normal Equation）

一般解として計算する方法。を求める係数のベクトルとすると、解は
　になる。
なぜこれで求まるかは、線形回帰の 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());
}

ロジスティック回帰

線形回帰の最急降下法と似ているが、 になる。
はシグモイド関数

で、解が必ず0~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;
}

検証

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

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