読者です 読者をやめる 読者になる 読者になる

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

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の結果とほぼ一致。いずれも、それらしく動いているようだ。



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