C#で、ディープラーニング　～誤差逆伝播法～

#C#で誤差逆伝播法

以前にもC#でディープラーニングするプログラムを記事にしました。
そのときは、誤差逆伝播法を使いませんでした。
今回は、前回のプログラムを改造して、誤差逆伝播法によるディープラーニングをしようと思います。
一応、プログラム実行時に誤差逆伝播をするかしないかということと、学習回数を指定できるようにしました。

なお、今回参考にした本は、「ゼロから作る　Deep Learning(オライリー)」です。
理論から説明があるので、Pythonではなくても、様々な言語でプログラムが書けます。

##行列計算用のクラスのプログラム

まずは、行列用のクラスです。

```CSharp:Matrix.cs
using System.Linq;

namespace Matrix
{
    class Mat
    {
        //プロパティR,Cは、行列の行と列の数を示す

        private int r = 0;
        public int R
        {
            get { return r; }
        }
        private int c = 0;
        public int C
        {
            get { return c; }
        }
        private bool err = false;
        public bool Err
        {
            get { return err; }
        }
        private double[][] matrix_data;
        public double[][] Matrix_data    //行列データを取得、設定
        {
            get {
                double[][] a = new double[2][];
                a[0] = new double[] { 0, 0 };
                a[1] = new double[] { 0, 0 };
                if (err) return a;
                else return matrix_data;
            }
            set
            {
                matrix_data = value;
            }
        }
        public double[][] Zero_matrix    //同じ形の、全ての要素が0の行列を示す
        {
            get
            {
                double[][] zm = new double[this.r][];
                for (int i = 0; i < this.r; i++)
                {
                    zm[i] = new double[this.c];
                    for (int j = 0; j < this.c; j++)
                    {
                        zm[i][j] = 0;
                    }
                }
                return zm;
            }
        }
        public double[][] T      //転置行列を取得
        {
            get
            {
                double[][] t = new double[this.c][];
                for (int i = 0; i < this.c; i++)
                {
                    t[i] = new double[this.r];
                    for (int j = 0; j < this.r; j++)
                    {
                        t[i][j] = matrix_data[j][i];
                    }
                }
                return t;
            }
        }
        public Mat(params double[][] vs)
        {
            int len = vs[0].Length;

            for (int i = 0; i < vs.Length; i++)
            {
                if (i != 0 && len != vs[i].Length)
                {
                    err = true;
                }
            }
            if (!err)
            {
                r = vs.Length;
                c = vs[0].Length;
                matrix_data = vs;
            }
        }

        private double[][] sigmoid_out;   //シグモイド関数実行時の結果を一時保存する
        public double[][] sigmoid()   //シグモイド関数
        {
            double[][] sig = new double[1][];
            sig[0] = new double[this.c];

            for(int i = 0; i < this.c; i++)
            {
                sig[0][i] = 1 / (1 + System.Math.Exp(-this.matrix_data[0][i]));
            }

            sigmoid_out = sig;

            return sig;
        }
        public double[][] sigmoid_backword(double[][] dout)  //誤差逆伝播時のシグモイド関数の微分
        {
            double[][] sig = new double[1][];
            sig[0] = new double[this.c];

            for (int i = 0; i < this.c; i++)
            {
                sig[0][i] = dout[0][i] * (1.0 - this.sigmoid_out[0][i]) * this.sigmoid_out[0][i];
            }

            return sig;
        }
        public double[][] softmax()   //ソフトマックス関数
        {
            double[][] sm = new double[1][];
            sm[0] = new double[this.c];

            double m = this.matrix_data[0].Max();

            double[] exp_a = new double[this.c];
            for (int i = 0; i < this.c; i++)
            {
                exp_a[i] = System.Math.Exp(this.matrix_data[0][i] - m);
            }

            double sum = 0.0;
            for (int i = 0; i < this.c; i++)
            {
                sum = sum + exp_a[i];
            }

            for (int i = 0; i < this.c; i++)
            {
                sm[0][i] = exp_a[i] / sum;
            }

            return sm;
        }
        public double cross_etp_err(Mat t)   //公差エントロピ誤差
        {
            double delta = 0.0000001;
            double sum = 0.0;
            for (int i = 0; i < this.c; i++)
            {
                sum = sum + t.matrix_data[0][i] * System.Math.Log(this.matrix_data[0][i] + delta);
            }

            return -sum;
        }
        public double[][] numerical_gradient(System.Func<double> loss)  //勾配を求める
        {
            double h = 0.0001;
            double[][] grad = new double[this.r][];
            double tmp_val = 0.0;
            double fxh1 = 0.0;
            double fxh2 = 0.0;

            for(int i = 0; i < this.r; i++)
            {
                grad[i] = new double[this.c];
                for(int j = 0; j < this.c; j++)
                {
                    tmp_val = this.matrix_data[i][j];
                    this.matrix_data[i][j] = tmp_val + h;
                    fxh1 = loss();

                    this.matrix_data[i][j] = tmp_val - h;
                    fxh2 = loss();

                    grad[i][j] = (fxh1 - fxh2) / (2 * h);
                    this.matrix_data[i][j] = tmp_val;
                }
            }

            return grad;
        }
        //以下　演算子オーバーロード
        public static double[][] operator +(Mat p1, Mat p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2.C && p1.R == p2.R)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] + p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator +(double[][] p1, Mat p2)
        {
            double[][] d = new double[p1.Length][];

            if (p1[0].Length == p2.C && p1.Length == p2.R)
            {
                for (int i = 0; i < p1.Length; i++)
                {
                    d[i] = new double[p1[0].Length];
                    for (int j = 0; j < p1[0].Length; j++)
                    {
                        d[i][j] = p1[i][j] + p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.Length; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator +(Mat p1, double[][] p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2[0].Length && p1.R == p2.Length)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] + p2[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator +(double p1, Mat p2)
        {
            double[][] d = new double[p2.R][];
            for (int i = 0; i < p2.R; i++)
            {
                d[i] = new double[p2.C];
                for (int j = 0; j < p2.C; j++)
                {
                    d[i][j] = p2.Matrix_data[i][j] + p1;
                }
            }

            return d;
        }
        public static double[][] operator +(Mat p1, double p2)
        {
            double[][] d = new double[p1.R][];
            for (int i = 0; i < p1.R; i++)
            {
                d[i] = new double[p1.C];
                for (int j = 0; j < p1.C; j++)
                {
                    d[i][j] = p1.Matrix_data[i][j] + p2;
                }
            }

            return d;
        }
        public static double[][] operator -(Mat p1, Mat p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2.C && p1.R == p2.R)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] - p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator -(double[][] p1, Mat p2)
        {
            double[][] d = new double[p1.Length][];

            if (p1[0].Length == p2.C && p1.Length == p2.R)
            {
                for (int i = 0; i < p1.Length; i++)
                {
                    d[i] = new double[p1[0].Length];
                    for (int j = 0; j < p1[0].Length; j++)
                    {
                        d[i][j] = p1[i][j] - p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.Length; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator -(Mat p1, double[][] p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2[0].Length && p1.R == p2.Length)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] - p2[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }

        public static double[][] operator -(double p1, Mat p2)
        {
            double[][] d = new double[p2.R][];
            for (int i = 0; i < p2.R; i++)
            {
                d[i] = new double[p2.C];
                for (int j = 0; j < p2.C; j++)
                {
                    d[i][j] = p1 - p2.Matrix_data[i][j];
                }
            }

            return d;
        }
        public static double[][] operator -(Mat p1, double p2)
        {
            double[][] d = new double[p1.R][];
            for (int i = 0; i < p1.R; i++)
            {
                d[i] = new double[p1.C];
                for (int j = 0; j < p1.C; j++)
                {
                    d[i][j] = p1.Matrix_data[i][j] - p2;
                }
            }

            return d;
        }
        public static double[][] operator *(Mat p1, Mat p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2.C && p1.R == p2.R)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] * p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator *(double[][] p1, Mat p2)
        {
            double[][] d = new double[p1.Length][];

            if (p1[0].Length == p2.C && p1.Length == p2.R)
            {
                for (int i = 0; i < p1.Length; i++)
                {
                    d[i] = new double[p1[0].Length];
                    for (int j = 0; j < p1[0].Length; j++)
                    {
                        d[i][j] = p1[i][j] * p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.Length; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator *(Mat p1, double[][] p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2[0].Length && p1.R == p2.Length)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] * p2[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }

        public static double[][] operator *(double p1, Mat p2)
        {
            double[][] d = new double[p2.R][];
            for (int i = 0; i < p2.R; i++)
            {
                d[i] = new double[p2.C];
                for (int j = 0; j < p2.C; j++)
                {
                    d[i][j] = p1 * p2.Matrix_data[i][j];
                }
            }

            return d;
        }
        public static double[][] operator *(Mat p1, double p2)
        {
            double[][] d = new double[p1.R][];
            for (int i = 0; i < p1.R; i++)
            {
                d[i] = new double[p1.C];
                for (int j = 0; j < p1.C; j++)
                {
                    d[i][j] = p1.Matrix_data[i][j] * p2;
                }
            }

            return d;
        }
        public static double[][] operator /(Mat p1, Mat p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2.C && p1.R == p2.R)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] / p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator /(double[][] p1, Mat p2)
        {
            double[][] d = new double[p1.Length][];

            if (p1[0].Length == p2.C && p1.Length == p2.R)
            {
                for (int i = 0; i < p1.Length; i++)
                {
                    d[i] = new double[p1[0].Length];
                    for (int j = 0; j < p1[0].Length; j++)
                    {
                        d[i][j] = p1[i][j] / p2.Matrix_data[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.Length; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double[][] operator /(Mat p1, double[][] p2)
        {
            double[][] d = new double[p1.R][];

            if (p1.C == p2[0].Length && p1.R == p2.Length)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p1.C];
                    for (int j = 0; j < p1.C; j++)
                    {
                        d[i][j] = p1.Matrix_data[i][j] / p2[i][j];
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }

        public static double[][] operator /(double p1, Mat p2)
        {
            double[][] d = new double[p2.R][];
            for (int i = 0; i < p2.R; i++)
            {
                d[i] = new double[p2.C];
                for (int j = 0; j < p2.C; j++)
                {
                    d[i][j] = p1 / p2.Matrix_data[i][j];
                }
            }

            return d;
        }
        public static double[][] operator /(Mat p1, double p2)
        {
            double[][] d = new double[p1.R][];
            for (int i = 0; i < p1.R; i++)
            {
                d[i] = new double[p1.C];
                for (int j = 0; j < p1.C; j++)
                {
                    d[i][j] = p1.Matrix_data[i][j] / p2;
                }
            }

            return d;
        }
        public static double[][] dot(Mat p1, Mat p2)    //行列のドット積
        {
            double[][] d = new double[p1.R][];
            double temp = 0;

            if (p1.C == p2.R)
            {
                for (int i = 0; i < p1.R; i++)
                {
                    d[i] = new double[p2.C];
                    for (int j = 0; j < p2.C; j++)
                    {
                        for(int a = 0; a < p1.C; a++)
                        {
                            temp = temp + p1.Matrix_data[i][a] * p2.Matrix_data[a][j];
                        }
                        d[i][j] = temp;
                        temp = 0.0;
                    }
                }
            }
            else
            {
                for (int k = 0; k < p1.R; k++)
                {
                    d[k] = new double[2] { 0, 0 };
                }
            }

            return d;
        }
        public static double sum(Mat mat)     //1次元のMatクラスの、全ての和
        {
            double s = 0.0;

            for (int i = 0; i < mat.c; i++)
            {
                s = mat.Matrix_data[0][i] + s;
            }

            return s;
        }
    }
}
```

##ディープラーニングのメインプログラム

```Csharp:Main.cs
namespace Matrix
{
    class Program
    {
        static Mat X = new Mat(
                new double[] { 6.3, 10.4, 11.1, 16.4 }
                ),
            W1 = new Mat(
                new double[] { 4.0, 2.0, 6.0, 1.0, 5.0 },
                new double[] { 1.0, 3.0, 7.0, 2.0, 1.0 },
                new double[] { 4.0, 2.0, 6.0, 1.0, 5.0 },
                new double[] { 1.0, 3.0, 7.0, 2.0, 1.0 }
                ),
            B1 = new Mat(
                new double[] { 2.0, 6.0, 4.0, 4.0, 1.0 }
                ),
            W2 = new Mat(
                new double[] { 3.0, 5.0 },
                new double[] { 6.0, 3.0 },
                new double[] { 1.0, 5.0 },
                new double[] { 3.0, 5.0 },
                new double[] { 6.0, 3.0 }
                ),
            B2 = new Mat(
                new double[] { 3.0, 6.0 }
                ),
            W3 = new Mat(
                new double[] { 7.0, 5.0, 7.0, 6.0 },
                new double[] { 4.0, 2.0, 4.0, 2.0 }
                ),
            B3 = new Mat(
                new double[] { 4.0, 6.0, 5.0, 4.0 }
                ),
            T = new Mat(
                new double[] { 0, 0, 1, 0 }
                );
        static Mat A1 = new Mat(
                new double[] { 0.0, 0.0, 0.0, 0.0, 0.0 }
                ),
                Z1 = new Mat(
                new double[] { 0.0, 0.0, 0.0, 0.0, 0.0 }
                ),
                A2 = new Mat(
                new double[] { 0.0, 0.0 }
                ),
                Z2 = new Mat(
                new double[] { 0.0, 0.0 }
                ),
                A3 = new Mat(
                new double[] { 0.0, 0.0, 0.0, 0.0 }
                ),
                Y = new Mat(
                new double[] { 0.0, 0.0, 0.0, 0.0 }
                );

        static double loss()     //交差エントロピー誤差を求める
        {
            A1.Matrix_data = Mat.dot(X, W1) + B1;
            Z1.Matrix_data = A1.sigmoid();
            A2.Matrix_data = Mat.dot(Z1, W2) + B2;
            Z2.Matrix_data = A2.sigmoid();
            A3.Matrix_data = Mat.dot(Z2, W3) + B3;
            Y.Matrix_data = A3.softmax();

            return Y.cross_etp_err(T);
        }
        static Mat loss_forward()     //誤差逆伝播時に使用する　交差エントロピー誤差とソフトマックス関数の微分
        {
            Mat dx = new Mat(
                new double[] { 0.0, 0.0, 0.0, 0.0 }
                );

            dx.Matrix_data = Y - T;
            dx.Matrix_data = dx / 2;

            return dx;
        }
        static void Main(string[] args)
        {
            double rate = 0.1;     //学習率

            Mat W1_grad = new Mat(W1.Zero_matrix),
                B1_grad = new Mat(B1.Zero_matrix),
                W2_grad = new Mat(W2.Zero_matrix),
                B2_grad = new Mat(B2.Zero_matrix),
                W3_grad = new Mat(W3.Zero_matrix),
                B3_grad = new Mat(B3.Zero_matrix);

            System.Console.WriteLine("誤差逆伝播：しない　0　、する　1");
            var deeplearning = int.Parse(System.Console.ReadLine());
            System.Console.WriteLine("学習回数");
            var learn = int.Parse(System.Console.ReadLine());

            for (int i = 0; i < learn; i++)
            {
                if (deeplearning == 0)            //誤差逆伝播法を使用しない場合
                {
                    W1_grad.Matrix_data = W1.numerical_gradient(loss);
                    B1_grad.Matrix_data = B1.numerical_gradient(loss);
                    W2_grad.Matrix_data = W2.numerical_gradient(loss);
                    B2_grad.Matrix_data = B2.numerical_gradient(loss);
                    W3_grad.Matrix_data = W3.numerical_gradient(loss);
                    B3_grad.Matrix_data = B3.numerical_gradient(loss);

                    W1.Matrix_data = W1 - (rate * W1_grad);
                    B1.Matrix_data = B1 - (rate * B1_grad);
                    W2.Matrix_data = W2 - (rate * W2_grad);
                    B2.Matrix_data = B2 - (rate * B2_grad);
                    W3.Matrix_data = W3 - (rate * W3_grad);
                    B3.Matrix_data = B3 - (rate * B3_grad);
                }
                else                     //誤差逆伝播法を使用する場合
                {
                    loss();
                    Mat affine3 = loss_forward();
                    Mat affine3_dx = new Mat(Mat.dot(affine3, new Mat(W3.T)));   //癖のあるプログラムなので、Matクラスの設計を考え直す必要あり。
                    W3_grad.Matrix_data = Mat.dot(new Mat(Z2.T), affine3);
                    B3_grad.Matrix_data = affine3.Matrix_data;

                    Mat sigmoid_2 = new Mat(A2.sigmoid_backword(affine3_dx.Matrix_data));

                    Mat affine2 = sigmoid_2;
                    Mat affine2_dx = new Mat(Mat.dot(affine2, new Mat(W2.T)));
                    W2_grad.Matrix_data = Mat.dot(new Mat(Z1.T), affine2);
                    B2_grad.Matrix_data = affine2.Matrix_data;

                    Mat sigmoid_1 = new Mat(A1.sigmoid_backword(affine2_dx.Matrix_data));

                    Mat affine1 = sigmoid_1;
                    Mat affine1_dx = new Mat(Mat.dot(affine1, new Mat(W1.T)));
                    W1_grad.Matrix_data = Mat.dot(new Mat(X.T), affine1);
                    B1_grad.Matrix_data = affine1.Matrix_data;

                    W1.Matrix_data = W1 - (rate * W1_grad);
                    B1.Matrix_data = B1 - (rate * B1_grad);
                    W2.Matrix_data = W2 - (rate * W2_grad);
                    B2.Matrix_data = B2 - (rate * B2_grad);
                    W3.Matrix_data = W3 - (rate * W3_grad);
                    B3.Matrix_data = B3 - (rate * B3_grad);
                }

                loss();            //パラメータ更新後、最終的な出力を得、表示。
                System.Console.WriteLine(i.ToString() + ":" + Y.Matrix_data[0][0].ToString() + "==" + Y.Matrix_data[0][2].ToString());
            }
        }
    }
}
```