C#で行列と列ベクトルの加算を実装した話

概要
ブロードキャストとは
ブロードキャストの使いどころ
Math.NETにブロードキャストが見当たらない
自力で実装
実行速度
まとめ
付録
- C#のベンチマークコード
- Pythonのベンチマークコード

概要

C# + Math.NETでニューラルネットワークを実装しようと思ったけど, NumPyでいうブロードキャストに該当する機能が存在しない
↓
ブロードキャストは行列と列ベクトルを加算するときなどに必要
↓
自分で実装しよう

ブロードキャストとは

NumPyに実装されている便利な機能。線形代数においてMxNの行列AとM次元列ベクトルBについて, A+Bという演算は未定義。でもNumPyではエラーにならない。なぜか？NumPyではブロードキャストによってベクトルBをN個並べてMxNの行列に変換した後に行列Aに加算される*1。だから正しく計算できる。

ブロードキャストの使いどころ

例えばNNでは, $WX + B$ という計算を頻繁に行う。ここで $W$ はNNの重み, $X$ はNNへの入力をバッチの要素ごとに各列に並べたもの, そして $B$ はNNのバイアスベクトルを各列に並べたもの。しかし, 実際のところ $B$ の各列は同じ列ベクトルであるため, 実装上はわざわざバイアスベクトルを各列に並べるのはメモリの無駄になる。どうせならバイアスベクトルを $\boldsymbol{b}$ として, $WX + \boldsymbol{b}$ と書けた方が楽。このようなときにNumPyのブロードキャストが役に立つ。

Math.NETにブロードキャストが見当たらない

C#用(厳密には.NET用)の数学ライブラリとしてMath.NETがある。Math.NETは内部でIntel MKLやOpenBLASなどの数値計算ライブラリを用いているため, 行列積の計算などが非常に高速。しかし, ざっと調べたところ, Math.NETにはNumPyのブロードキャストに当たる機能が存在しなかった(自分が見逃しているだけかもしれないけど)。

自力で実装

以下, 行列とベクトルの各要素の型はfloat(32bit浮動小数点数)であることを前提に進める。
AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)というメソッドを実装することにした。このメソッドは行列leftと列ベクトルrightを加算した結果を行列sumに格納する。色々と試行錯誤したのだがそれを書くには余白が足りなすぎるので割愛。結論から示すと以下のような実装になった(例外処理は省略)。

static void AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)
{
    var leftArray = left.AsColumnMajorArray();
    var rightArray = right.AsArray();
    var sumArray = sum.AsColumnMajorArray();
    var (rowNum, colNum) = (sum.RowCount, sum.ColumnCount);

    unsafe
    {
        const int UNROLL_STEP_NUM = 8;
        fixed (float* leftPtr = leftArray)
        fixed (float* rightPtr = rightArray)
        fixed (float* sumPtr = sumArray)
        {
            for (var j = 0; j < colNum; j++)
            {
                var l = leftPtr + j * rowNum;
                var s = sumPtr + j * rowNum;
                for (var i = 0; i < rowNum; i += UNROLL_STEP_NUM)
                {
                    var leftVec = Avx.LoadVector256(l + i);
                    var rightVec = Avx.LoadVector256(rightPtr + i);
                    var sumVec = Avx.Add(leftVec, rightVec);
                    Avx.Store(s + i, sumVec);
                }
            }

            // 端数処理
            for (var j = 0; j < colNum; j++)
            {
                var l = leftPtr + j * rowNum;
                var s = sumPtr + j * rowNum;
                for (var i = UNROLL_STEP_NUM * (rowNum / UNROLL_STEP_NUM); i < rowNum; i++)
                    s[i] = l[i] + rightPtr[i];
            }
        }
    }
}

C#ではAVX組み込み関数を用いることができるので, それを用いて8要素同時に加算を行っている。列ベクトルの次元(rowNum)を8で割った余りは単純に2重ループを回して計算している(//端数処理以降)。

実行速度

測定方法

1024x1024の行列に1024次元の列ベクトルを加算する処理を1000回行う。行列とベクトルの各要素は0.0~1.0の一様乱数で初期化する。
また、比較対象としてNumPyでも同様の処理を行う。
ソースコードは"付録"の項を参照。

マシンスペック

CPU	Intel(R) Core(TM) i7-10510U
メモリ	DDR3 16GB

NumPyの情報

バージョン: 1.23.3
NumPyの設定によって結果が変わると思うので, numpy.__config__.show()の結果を載せておく。

>>> np.__config__.show()
openblas64__info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__info']
libraries = ['openblas64__info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
blas_ilp64_opt_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__info']
libraries = ['openblas64__info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
openblas64__lapack_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__lapack_info']
libraries = ['openblas64__lapack_info']
language = f77
lapack_ilp64_opt_info:
library_dirs = ['D:\\a\\numpy\\numpy\\build\\openblas64__lapack_info']
libraries = ['openblas64__lapack_info']
language = f77
define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None), ('HAVE_LAPACKE', None)]
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3
found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

※(2022/09/13 追記) NumPyをAVX2まで対応させたものに変更しました。以下が新たなnumpy.__config__.show()の結果です。

>>> np.__config__.show()
blas_armpl_info:
NOT AVAILABLE
blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
NOT AVAILABLE
accelerate_info:
NOT AVAILABLE
atlas_3_10_blas_threads_info:
NOT AVAILABLE
atlas_3_10_blas_info:
NOT AVAILABLE
atlas_blas_threads_info:
NOT AVAILABLE
atlas_blas_info:
NOT AVAILABLE
blas_info:
NOT AVAILABLE
blas_src_info:
NOT AVAILABLE
blas_opt_info:
NOT AVAILABLE
lapack_armpl_info:
NOT AVAILABLE
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
NOT AVAILABLE
openblas_clapack_info:
NOT AVAILABLE
flame_info:
NOT AVAILABLE
atlas_3_10_threads_info:
NOT AVAILABLE
atlas_3_10_info:
NOT AVAILABLE
atlas_threads_info:
NOT AVAILABLE
atlas_info:
NOT AVAILABLE
lapack_info:
NOT AVAILABLE
lapack_src_info:
NOT AVAILABLE
lapack_opt_info:
NOT AVAILABLE
numpy_linalg_lapack_lite:
language = c
define_macros = [('HAVE_BLAS_ILP64', None), ('BLAS_SYMBOL_SUFFIX', '64_')]
Supported SIMD extensions in this NumPy install:
baseline = SSE,SSE2,SSE3,SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
found =
not found = XOP,FMA4,AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

測定結果

平均実行速度[ms]

マイコード 0.473

NumPy(SSE3まで対応)のadd関数 1.233

※(2022/09/13 追記) この実験で用いたNumPyはSSE3までにしか対応していないため, NumPyにとってかなり不利な比較になっていました。後日, AVXに対応したものを用意して再実験した結果を掲載する予定です。

※(2022/09/13 追記) 新たな実験結果を追加しました。NumPyをAVXに対応させても, 特に大きな速度差はありませんでした。

	平均実行速度[ms]
マイコード	0.473
NumPyのadd関数	1.248

NumPyよりも2倍以上高速に動作した。

まとめ

C#を用いて行列と列ベクトルの加算を実装した。実は行列と行ベクトルの加算も実装したが, 内容はほとんど同じなので割愛。自分は高速化に関して素人なので, もっとよい実装方法があったらコメントで指摘してくださるとうれしいです。

付録

C#のベンチマークコード

using System;
using System.Diagnostics;
using System.Runtime.Intrinsics.X86;

using MathNet.Numerics.LinearAlgebra.Single;

class Program
{
    const int DIM = 1024;
    const int SAMPLE_NUM = 1000;

    public static void Main()
    {
        Benchmark(AddColumnVector);
    }

    static void Benchmark(Action<DenseMatrix, DenseVector, DenseMatrix> f)
    {
        var result = new DenseMatrix(DIM, DIM);

        var sw = new Stopwatch();
        for (var i = 0; i < SAMPLE_NUM; i++)
        {
            var left = InitMatrixWithRand(DIM, DIM);
            var right = DenseVector.OfArray(InitMatrixWithRand(DIM, 1).AsColumnMajorArray());
            sw.Start();
            f(left, right, result);
            sw.Stop();
        }
        Console.WriteLine($"Ellapsed = {(float)sw.ElapsedMilliseconds / SAMPLE_NUM}[ms]");
    }

    static DenseMatrix InitMatrixWithRand(int rowNum, int colNum)
    {
        var ret = new DenseMatrix(rowNum, colNum);
        var array = ret.AsColumnMajorArray();
        for (var i = 0; i < array.Length; i++)
            array[i] = Random.Shared.NextSingle();
        return ret;
    }

    static void AddColumnVector(DenseMatrix left, DenseVector right, DenseMatrix sum)
    {
        var leftArray = left.AsColumnMajorArray();
        var rightArray = right.AsArray();
        var sumArray = sum.AsColumnMajorArray();
        var (rowNum, colNum) = (sum.RowCount, sum.ColumnCount);

        unsafe
        {
            const int UNROLL_STEP_NUM = 8;
            fixed (float* leftPtr = leftArray)
            fixed (float* rightPtr = rightArray)
            fixed (float* sumPtr = sumArray)
            {
                for (var j = 0; j < colNum; j++)
                {
                    var l = leftPtr + j * rowNum;
                    var s = sumPtr + j * rowNum;
                    for (var i = 0; i < rowNum; i += UNROLL_STEP_NUM)
                    {
                        var leftVec = Avx.LoadVector256(l + i);
                        var rightVec = Avx.LoadVector256(rightPtr + i);
                        var sumVec = Avx.Add(leftVec, rightVec);
                        Avx.Store(s + i, sumVec);
                    }
                }

                // 端数処理
                for (var j = 0; j < colNum; j++)
                {
                    var l = leftPtr + j * rowNum;
                    var s = sumPtr + j * rowNum;
                    for (var i = UNROLL_STEP_NUM * (rowNum / UNROLL_STEP_NUM); i < rowNum; i++)
                        s[i] = l[i] + rightPtr[i];
                }
            }
        }
    }
}

Pythonのベンチマークコード

import time
import numpy as np

DIM = 1024
SAMPLE_NUM = 1000

ellapsed = 0
s = np.zeros((DIM, DIM))
for _ in range(SAMPLE_NUM):
    left = np.random.rand(DIM, DIM)
    right = np.random.rand(DIM, 1)
    start = time.perf_counter()
    np.add(left, right, s)
    end = time.perf_counter()
    ellapsed += end - start
    
print(f"Ellapsed = {(ellapsed / SAMPLE_NUM) * 1000}[ms]")

*1:もちろん, わざわざMxN行列用のメモリを確保してベクトルBをそこに書き込んでいるわけではなく, 実装上は普通に行列Aの各列にベクトルBを加算している