UnityでPytorchライクの機械学習ライブラリを作る。3日目：行列積の高速化

# はじめに
　今回はドット積と単項演算を実装する予定でしたが、ドット積はライブラリの根幹をなす部分なので実装するにあたって少し速度を追求しようと思います。[こちら](https://qiita.com/nocturnality/items/92b3cf7d7fc5dab4a3a7)の記事を参考にして現時点で実行可能なものを試してみることにしました。具体的には以下の三つです。

* ループ交換
* ループアンローリング
* 並列化

　これらを実装してみて、速度を計測した上で最も早いものを使用することにします。

# 比較方法について
　比較の方法は15要素の配列$[50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 600, 700, 800, 900, 1000]$について、それぞれの計算方式でより短い時間で計算できた回数が多いものが優れているものとします。

# テスト用のプログラム
　時間の計測には以下の関数を用いてDot1とDot2のところに比較したい関数が入ります。

```cs:TensorTest.cs
using R = System.Double;
using System.Diagnostics;

namespace Rein.Tests
{
    public class TensorTest
    {
        public void ProfileDot(){
            Stopwatch sw1 = new Stopwatch();
            Stopwatch sw2 = new Stopwatch();
            // int[] shape = new int[2]{ , 300 };
            int loopNum = 10;
            int[] sizes = new int[]{ 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 600, 700, 800, 900, 1000 };

            foreach (int size in sizes){
                Console.WriteLine(size);
                sw1.Reset();
                sw2.Reset();
                for (int i=0; i < loopNum; i++){
                    Tensor left = new Tensor(
                        new int[2]{size, size}
                    );
                    Tensor right = new Tensor(
                        new int[2]{size, size}
                    );

                    sw1.Start();
                    Tensor out1 = new Rein.Functions.Dot1().Forward(left, right);
                    sw1.Stop();

                    sw2.Start();
                    Tensor out2 = new Rein.Functions.Dot2().Forward(left, right);
                    sw2.Stop();

                    System.Threading.Thread.Sleep(100);
                }

                Console.WriteLine($"ループ交換前：{sw1.ElapsedMilliseconds}/{loopNum}ms");
                Console.WriteLine($"ループ交換後：{sw2.ElapsedMilliseconds}/{loopNum}ms");
            }
        }
    }
}
```

# ループ交換の有無
　行列として表している配列はメモリでは以下のような順番で格納されて居ます。縦方向に移動してアクセスすると移動幅が大きくなり時間がかかります。そこで横方向へのアクセスを増やすことで速度を上げる手法がループ交換です。

![行列.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179010/c87130a5-6e75-fad7-2dc3-3fdcacc04360.jpeg)

## 高速化なしの実装
　まず最初に高速化を施さない場合の実装を書いておきます。

```cs:Dot.cs
namespace Rein.Functions{
    public class Dot: BinaryFunction{

        public Dot():base("Dot"){
            
        }

        protected override Tensor BinaryForward(Tensor left, Tensor right)
        {
            int N = left.Size / left.Shape.Last();
            int M= right.Shape[1];
            int L = right.Shape[0];
            R[] data = new R[N * M];

            for (int i=0; i < N; i++){
                for (int j=0; j < M; j++){
                    for (int k=0; k < L; k++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }
            }

            List<int> shape = new List<int>(left.Shape);
            shape[shape.Count - 1] = right.Shape[1];
            return new Tensor(data, shape);
        }

        protected override void BinaryBackward()
        {
            
        }
    }
}
```

## ループ交換を行った場合の実装
　ループ内部のjとkを交換して実装します。

```cs:Dot.cs
{
    public class Dot_optim1: BinaryFunction{
        public Dot_optim1():base("Dot"){

        }

        protected override Tensor BinaryForward(Tensor left, Tensor right)
        {
            int N = left.Size / left.Shape.Last();
            int M= right.Shape[1];
            int L = right.Shape[0];
            R[] data = new R[N * M];

            for (int i=0; i < N; i++){
                for (int k=0; k < L; k++){
                    for (int j=0; j < M; j++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }
            }

            List<int> shape = new List<int>(left.Shape);
            shape[shape.Count - 1] = right.Shape[1];
            return new Tensor(data, shape);
        }

        protected override void BinaryBackward()
        {
            throw new System.NotImplementedException();
        }
    }
}
```


## 計測
　以下が計測の結果です。ほぼ全ての結果でループ交換ありの場合の方が高速でした。そこで以降はループ交換ありで実装していきます。

|サイズ|ループ交換なし|ループ交換あり|
|---|---|---|
|50|14ms|12ms|
|100|92ms|86ms|
|150|388ms|417ms|
|200|660ms|496ms|
|250|987ms|755ms|
|300|2066ms|1611ms|
|350|2755ms|2144ms|
|400|4631ms|3651ms|
|450|7748ms|5386ms|
|500|10486ms|7460ms|
|600|19710ms|13072ms|
|700|48664ms|19344ms|
|800|56212ms|31548ms|
|900|116427ms|41787ms|
|1000|127298ms|58175ms|

![ループ交換による速度変化 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179010/40a3f9e8-d8ff-6ad7-b8bf-050179326e7f.png)


# ループアンローリング
　ループアンローリングはfor文の刻み幅を大きくすることで条件判断の回数を少なくする方法です。基本的にはアンローリングする数を増やすほど速度は早くなりますが、その分コードが長くなってしまいます。
　例えばループで13回処理をしたい時に刻み幅を2にして処理をすると、以下のようになります。

```cs
int i;
for(i = 0;i < (13 - 1);i+=2){
    DoSomething(i);
    DoSomething(i + 1);
}

// 余った文の処理を行う
for (;i < 13;i++){
    DoSomething(i);
}
```

　このようにfor文一つで済んでいた物が複数必要になるので実装は面倒になります。今回は2つのループでアンローリングを行い、それぞれの刻み幅は4にしています。


## ループアンローリングの実装
　以下が実装です。

```cs:Dot.cs
    public class Dot_optim6: BinaryFunction{
        public Dot_optim6():base("Dot"){

        }

        protected override Tensor BinaryForward(Tensor left, Tensor right)
        {
            int N = left.Size / left.Shape.Last();
            int M= right.Shape[1];
            int L = right.Shape[0];
            R[] data = new R[N * M];
            int k, k1, k2, k3, j;

            for (int i=0; i < N; i++){
                for (k=0; k < L - 3; k+=4){
                    k1 = k + 1;
                    k2 = k + 2;
                    k3 = k + 3;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3];

                        data[i*M + j] += left.Data[i*L + k1] * right.Data[k1*M + j];
                        data[i*M + j+1] += left.Data[i*L + k1] * right.Data[k1*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k1] * right.Data[k1*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k1] * right.Data[k1*M + j+3];

                        data[i*M + j] += left.Data[i*L + k2] * right.Data[k2*M + j];
                        data[i*M + j+1] += left.Data[i*L + k2] * right.Data[k2*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k2] * right.Data[k2*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k2] * right.Data[k2*M + j+3];

                        data[i*M + j] += left.Data[i*L + k3] * right.Data[k3*M + j];
                        data[i*M + j+1] += left.Data[i*L + k3] * right.Data[k3*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k3] * right.Data[k3*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k3] * right.Data[k3*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }

                for (; k < L; k++){
                    k1 = k + 1;
                    k2 = k + 2;
                    k3 = k + 3;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*L + j] += left.Data[i*L + k] * right.Data[k*L + j];
                    }
                }
            }

            List<int> shape = new List<int>(left.Shape);
            shape[shape.Count - 1] = right.Shape[1];
            return new Tensor(data, shape);
        }

        protected override void BinaryBackward()
        {
            throw new System.NotImplementedException();
        }
    }
```

## 計測
　今回はループ交換のみとループ交換プラス+ループアンローリングを比較します。大体１〜２割程度の改善がみられたので今後はループ交換+アンローリングを暫定的に最善の手法とします。

|サイズ|交換|交換+アンローリング|
|---|---|---|
|50|16ms|14ms|
|100|106ms|62ms|
|150|266ms|163ms|
|200|531ms|393ms|
|250|891ms|717ms|
|300|1591ms|1361ms|
|350|2284ms|1935ms|
|400|3752ms|3252ms|
|450|5328ms|4579ms|
|500|7560ms|6697ms|
|600|13445ms|11000ms|
|700|21393ms|18280ms|
|800|31945ms|26829ms|
|900|46691ms|37558ms|
|1000|61885ms|50993ms|

![アンローリング.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179010/4cb328c8-6eb9-fa08-e7fc-d14715d37a7a.png)


# 並列化
　並列化は真面目に実装すると手間なので、C#のSystemライブラリで行える手軽な方法を試そうと思います。
　for文などの複数回繰り返し行われる処理を並列化したい時にはSystem.Threading.TasksのParallelを使用します。例として以下のようなfor文を並列化する場合

```cs:hoge.cs
for(int i = 0; i < N; i++){
    DoSomething(i);
}
```

 for文を丸ごとParallel.Forに置き換えて

```cs:hoge.cs
Parallel.For(0, N, i =>{
    DoSomething(i);
});
```

　このように簡単に導入できます。

## 並列化した実装
　実際にループ交換+アンローリングを並列化したものがこちらです。

```cs:Dot.cs
public class Dot_optim7: BinaryFunction{
        public Dot_optim7():base("Dot"){

        }

        protected override Tensor BinaryForward(Tensor left, Tensor right)
        {
            int N = left.Size / left.Shape.Last();
            int M= right.Shape[1];
            int L = right.Shape[0];
            R[] data = new R[N * M];

            Parallel.For(0, N, i => {
                int k;
                for (k=0; k < L - 3; k+=4){
                    int k1 = k + 1;
                    int k2 = k + 1;
                    int k3 = k + 1;
                    int j;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3];

                        data[i*M + j] += left.Data[i*L + k1] * right.Data[k1*M + j];
                        data[i*M + j+1] += left.Data[i*L + k1] * right.Data[k1*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k1] * right.Data[k1*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k1] * right.Data[k1*M + j+3];

                        data[i*M + j] += left.Data[i*L + k2] * right.Data[k2*M + j];
                        data[i*M + j+1] += left.Data[i*L + k2] * right.Data[k2*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k2] * right.Data[k2*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k2] * right.Data[k2*M + j+3];

                        data[i*M + j] += left.Data[i*L + k3] * right.Data[k3*M + j];
                        data[i*M + j+1] += left.Data[i*L + k3] * right.Data[k3*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k3] * right.Data[k3*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k3] * right.Data[k3*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }

                for (; k < L; k++){
                    int k1 = k + 1;
                    int k2 = k + 1;
                    int k3 = k + 1;
                    int j;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }
            });

            List<int> shape = new List<int>(left.Shape);
            shape[shape.Count - 1] = right.Shape[1];
            return new Tensor(data, shape);
        }

        protected override void BinaryBackward()
        {
            throw new System.NotImplementedException();
        }
    }
```

## 測定
　ループ交換+アンローリングとループ交換+アンローリング+並列化を比較していきます。結果としは200~800は並列化した方が早くなりました。ちょっと意味が分からないです。何度か試してみたのですが同様の結果が出ました。おそらく並列化する際にThreadの立ち上げに時間がかかるのでサイズ50の時がサイズ100の時より遅くなったと思われます。
　正直私のPCスペックでは性能があまり上がりませんし速度も安定しない上に、私の用途ではサイズ100以下の物がほとんどなので基本的には並列化なしを採用しようと思います。

|サイズ|交換+アンローリング|+並列化|
|---|---|---|
|50|17ms|134ms|
|100|95ms|107ms|
|150|204ms|257ms|
|200|460ms|407ms|
|250|850ms|663ms|
|300|1322ms|964ms|
|350|2227ms|1566ms|
|400|3213ms|2211ms|
|450|4959ms|3397ms|
|500|6824ms|4911ms|
|600|11920ms|9986ms|
|700|19104ms|14949ms|
|800|27855ms|26552ms|
|900|39310ms|43938ms|
|1000|50518ms|60170ms|

![並列化.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179010/87ec4363-5ad3-6a5d-9b39-0b58ffbe8747.png)


# (2021/2/21追加)代入回数を減らす
　上の三つを実装してみて気付いたのですが、ループアンローリングの際に同じインデックスに代入するような計算を一度に足し合わせることでより早く計算できるのではないかと考えました。
　実際に試してみます

## 実装(ループ交換+アンローリング+代入回数を減らす)

```cs:Dot.cs
public class Dot_optim8: BinaryFunction{
        public Dot_optim8():base("Dot"){

        }

        protected override Tensor BinaryForward(Tensor left, Tensor right)
        {
            int N = left.Size / left.Shape.Last();
            int M= right.Shape[1];
            int L = right.Shape[0];
            R[] data = new R[N * M];
            int k, k1, k2, k3, j;

            for (int i=0; i < N; i++){
                for (k=0; k < L - 3; k+=4){
                    k1 = k + 1;
                    k2 = k + 2;
                    k3 = k + 3;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j]
                            + left.Data[i*L + k1] * right.Data[k1*M + j]
                            + left.Data[i*L + k2] * right.Data[k2*M + j]
                            + left.Data[i*L + k3] * right.Data[k3*M + j];
                        
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1]
                            + left.Data[i*L + k1] * right.Data[k1*M + j+1]
                            + left.Data[i*L + k2] * right.Data[k2*M + j+1]
                            + left.Data[i*L + k3] * right.Data[k3*M + j+1];

                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2]
                            + left.Data[i*L + k1] * right.Data[k1*M + j+2]
                            + left.Data[i*L + k2] * right.Data[k2*M + j+2]
                            + left.Data[i*L + k3] * right.Data[k3*M + j+2];

                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3]
                            + left.Data[i*L + k1] * right.Data[k1*M + j+3]
                            + left.Data[i*L + k2] * right.Data[k2*M + j+3]
                            + left.Data[i*L + k3] * right.Data[k3*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                    }
                }

                for (; k < L; k++){
                    k1 = k + 1;
                    k2 = k + 2;
                    k3 = k + 3;
                    for (j=0; j < M - 3; j+=4){
                        data[i*M + j] += left.Data[i*L + k] * right.Data[k*M + j];
                        data[i*M + j+1] += left.Data[i*L + k] * right.Data[k*M + j+1];
                        data[i*M + j+2] += left.Data[i*L + k] * right.Data[k*M + j+2];
                        data[i*M + j+3] += left.Data[i*L + k] * right.Data[k*M + j+3];
                    }
                    
                    for(; j < M; j++){
                        data[i*L + j] += left.Data[i*L + k] * right.Data[k*L + j];
                    }
                }
            }

            List<int> shape = new List<int>(left.Shape);
            shape[shape.Count - 1] = right.Shape[1];
            return new Tensor(data, shape);
        }
    }
```

## 測定
　左が(交換＋アンローリング)、右が(交換＋アンローリング＋代入減らし)です。代入減らしを行うことで３割程度の高速化ができるようなのでこれを採用します。

|サイズ|交換+アンローリング|+代入減らし|
|---|---|---|
|50|13ms|12ms|
|100|74ms|52ms|
|150|233ms|142ms|
|200|502ms|247ms|
|250|747ms|444ms|
|300|1482ms|1033ms|
|350|2031ms|1357ms|
|400|2862ms|1957ms|
|450|4387ms|3068ms|
|500|6215ms|4330ms|
|600|11164ms|7555ms|
|700|17970ms|12607ms|
|800|27703ms|19008ms|
|900|39243ms|27557ms|
|1000|54042ms|36383ms|


![代入減らし.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179010/30602e2b-c69a-2bc1-037f-98e2a4e4d86d.png)


# 結論
　基本はループ交換+アンローリング＋代入減らしを使用して必要になったら、並列化を使用するのがベスト。

# 終わりに
　ということで今回は、行列積の高速化を行いました。今回高速化をするにあたって色々調べてみたのですが突き詰めるほど低階層というかCPUのキャッシュなどを意識した実装が必要になってくるので終わりが無いなと思い、今回は言語レベルやメモリレベルの改良に留めておきました。GPUを使用するという方法もありますが、私がGPU持ってないので実装してもしょうがないということでそちらも後回しにしました。
　とにかく動けば問題ないのです。
　次回は、今回の結果を踏まえてドット演算の勾配処理の実装を行います。