以下の内容はhttps://lpha-z.hatenablog.com/entry/2025/12/07/231500より取得しました。

int8テンソルコアを用いたDGEMMエミュレーションライブラリは何を使えばよいか

最近のGPUは倍精度演算器が少なくDGEMM（倍精度行列積）を効率よく実行できませんが、代わりにたくさん載っているint8テンソルコアを使って実質的に倍精度行列積を計算する（＝エミュレーションする）方式が最近流行っているようです。 その先駆けは2023年6月に大友らがarXivに投稿したプレプリント（[2306.11975] DGEMM on Integer Matrix Multiplication Unit 🔓）で、その後内野らが改良手法を提案しています（Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit - Yuki Uchino, Katsuhisa Ozaki, Toshiyuki Imamura, 2025 🔓）。 NVIDIAが2025年8月にリリースしたCUDA 13.0に含まれるcuBLASでは公式にint8テンソルコアを使用したエミュレーションが実装され、外部ライブラリに頼ることなく手軽に切り替えることができるようになりました。

しかしながらこれらの手法は、いずれも異なる精度と異なる速度を持っているため、どれを使うべきなのかよくわかりません。 内野らの論文には大友らの手法との比較が掲載されていますが、それより後に出たCUDA 13.0のcuBLASとの比較は当然ながら掲載されていません。 そこで、これらの手法の精度と速度を比較してみました。

環境

記事を執筆しているうちにCUDA 13.1が出ましたが（2025-12-04）、もう一回実験するのも大変なのでCUDA 13.0でいきます。

また、以下では正方行列同士の行列積のみを取り扱います。 行列サイズNとは、一辺の長さのことを指します。 このとき、行列の要素数はN²、行列積の計算コストはΘ(N³)となります。

尾崎スキーム

尾崎スキームは、浮動小数点数の行列積結果を正確に求めるために固定小数点数の行列積を行う方法です。

尾崎スキームでは、まず入力の浮動小数点数行列を固定小数点数行列に変換します。 固定小数点数と言っても実際のところは、桁合わせされた複数の倍精度浮動小数点数の和として表されることが多いようです。 ここで、桁の区切りは53ビットとかではなく、例えば20ビットなどのあまり多くないビット数にするのがポイントです。 つまり、元の行列（M）を、2²⁰の位より上を持つ倍精度浮動小数点数の行列（A）＋残りのうち2⁰の位より上を持つ倍精度浮動小数点数の行列（B）＋残りのうち2⁻²⁰の位より上を持つ倍精度浮動小数点数の行列（C）＋残りのうち2⁻⁴⁰の位より上を持つ倍精度浮動小数点数の行列（D）＋残りのうち2⁻⁶⁰の位より上を持つ倍精度浮動小数点数の行列（E）＋……のように分割して表します。

桁の区切りを20ビットなどにしたのは、各行列同士の積を無誤差で計算できるようにするためです。 各行列（A・B・C・D・E・……）の各要素が20ビットしか仮数部を持っていないとき、その各積の仮数部は高々40ビットで表せます。 その場合、各行列のサイズが2¹³=4096以下であれば、各行列同士の積が無誤差になります。 各積は位が揃っているため、2¹³個足し合わせた結果が53ビットに収まるからです。 実際には論理が逆で、行列のサイズを見てから各行列が何ビット持つべきかを決めるようです。

さて、各行列同士の積は、通常のDGEMMライブラリが使えます。 求めたい行列積を M1M2 = (A1+B1+C1+D1+E1)(A2+B2+C2+D2+E2) と変換したとします。 このとき、右辺を展開すると A1A2 + A1B2 + A1C2 + A1D2 + A1E2 + B1A2 + …… + E1E2 となりますが、いずれの項も単独では行列積が無誤差で計算できます。 既存のチューニングされたDGEMMライブラリが使える点が尾崎スキームの重要な利点で、高精度行列積コードを一から作る労力から解放されます。

最後に、各項を足し合わせます。 ここはただ足し合わせるだけの処理で、特にメモリ律速なので、適当に書いたdouble-doubleの加算を使っておけばOKです。

尾崎スキームでは、精度と速度のトレードオフを柔軟に変更することができます。 例えば、E1E2のような小さい項が最終結果に与える影響は少なそうです。 その様な考えに基づいて、上から何ビットまでしか計算しない、のような最適化により、精度が低下する代わりに速度を改善することができます。 ただし、「上から何ビット」というのは、桁落ちにより「最終結果の上から何ビット」と大きく乖離する可能性があります。 なので、悪条件問題では、上から80ビット、とか計算しないと、最終結果を53ビット精度で求めることはできないかもしれません。 逆に言えば、問題の性質に応じて、精度に対して必要十分な計算量に設定できる点も尾崎スキームの利点の一つです。

GPUで使われる尾崎スキームは、無誤差計算部分でDGEMMの代わりにint8行列積を使う点が異なります。 ここまでで説明した方法は、普通にやってしまうと誤差が生じるDGEMMを、複数回の無誤差なDGEMMに分解して計算するものでした。 GPUで使われる方法は、普通にやってしまうと誤差が生じる――というか「遅い」――DGEMMを、複数回の無誤差なint8行列積に分解して計算するものです。 int8テンソルコア（int8行列積をint32に積算）を使う場合、上で説明したような技巧的な分割を行わなくとも、行列サイズが小さければ自動的に無誤差になります。 これは、もともと要素ごとの積の絶対値は2¹⁴以下なので、これを2¹⁷個くらい足してもint32に収まる、ということを利用しています。 最後に各項を足し合わせる時は、8nビットシフトして足し合わせるみたいなことが必要になります。

CUDA 13.0のcuBLASのDGEMMエミュレーションの観察

CUBLAS_EMULATE_DOUBLE_PRECISION

これを1に設定すると、DGEMMの呼び出しをint8テンソルコアを使った高速エミュレーションで置き換えることを許可できるようです。

CUBLAS_FIXEDPOINT_EMULATION_MANTISSA_BIT_COUNT

エミュレーション精度を指定することができます。 尾崎スキームでは普通、エミュレーション精度を変えたいときはスライス数を変化させますが、精度をビット数で指定するインタフェースのようです。 速度調査結果を見ると、

  • 1～7を指定した時
  • 8～15を指定した時
  • 16～23を指定した時
  • ……
  • 8n～8n+7を指定した時
  • ……

で速度が変化していないようでした。 このことからすると、1～7は1スライス、8～15は2スライス、16～23は3スライス、……使って計算しているのだと思います。 7ビットまでが1スライスというのはint8が符号を除いて7ビットであることと整合し、8～15ビットが2スライス・16～23ビットが3スライス・……というのは1スライスにつき8ビット分の表現力が追加されることと整合的です。 よく見たら、https://docs.nvidia.com/cuda/cublas/index.html#representation-and-mappingsにもそう書いてありました。

0を指定した場合、試した感じでは63を指定した場合（8スライス）と同じになっていました。 倍精度浮動小数点数の精度は53ビットで、これを正確に保持するためには55ビット分保持できる7スライスが必要ですが、絶対値が小さな値が出てくると固定小数点数に変換するときに誤差が生じてしまいます。 もう1スライス余分に計算しておけば、変換するときの誤差が最終結果に与える影響は十分小さいと言えるでしょう。 なので、「いい感じの」精度でやってくれるモードとして、8スライスでやっているのだと思います。

実験した範囲では9スライス以上にしても精度向上は見られなかったので、8スライスは必要十分なスライス数なのでしょう。 もしかすると、より悪条件な行列積において9スライス以上が有用なのかもしれません。 0をした場合の実行時間は8スライスの時より遅く9スライスの時より速いような感じだったので、行列の中身を見て条件数を推定してスライス数を切り替えたりする機能があるのかもしれません。

CUBLAS_EMULATION_STRATEGY

これはperformantまたはeagerに設定することができます。 eagerに設定すると、どんなDGEMMも指定した方式で実行するようです。 一方で、performantだと、

  • N≦128の時、設定を無視して常にDGEMMを行う
  • N>128の時、

のような動作をしているようでした。 けっこう雑な切り替えのようです。

私の環境では残念ながら、N=1024でさえDGEMMの方が速かったので、128で切り替えるのは微妙そうです。 逆にN=2048では、精度80ビット以上でDGEMMに切り替わると精度が落ちるわ速度も落ちるわで何もいいことがなさそうでした。

内野らの改良

Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit - Yuki Uchino, Katsuhisa Ozaki, Toshiyuki Imamura, 2025を読み解いていきます。

errfree_sum

3.2節に書いてある方式です。 int8行列積の結果はint32行列ですが、ozIMMUではこれをdoubleに変換してから足していました。 一方でint32加算がオーバーフローしないなら、int32加算の方がdouble加算より高速なので、足してからdoubleに変換する方が高速なのでそうする、という方式のようです。 基本的には精度は変わらず、単に高速になる効果だけが得られるようです。

nearest_split

3.1節に書いてある方式です。 仮数部を分割する際、下位部分を切り捨てて分割するのではなく、下位部分を四捨五入丸めで分割することにより、スライス当たり1ビット表現力が高まるという方式のようです。 精度は上がるものの、四捨五入丸めを実現することには相応のコストがかかります。 そのため、論文のFigure14を見た感じでは、0.3スライス追加したみたいな感じの速度低下と精度向上になるようです。

最近接丸めを使うと128みたいなint8に収まらない数が出て来て困るのではないかなと思っていましたが、これは違うようでした。 分割結果は、-64～64しか使わないようになっていました。

hybrid

3.3節に書いてある方式です。 要するに上の二つは独立なので両方盛り込みましたという方式です。

実験

速度と精度の実験を行いました。 実験環境は「環境」節で書いた通りです。 入力行列は、各要素が独立に[-1.0,1.0]の一様分布に従うようなdouble行列としました。 より正確には、以下のような行列積計算を行いました。

    int N = atoi( argv[1] );
    std::mt19937_64 mt;
    std::uniform_real_distribution<double> dist( -1.0, 1.0 );

    double* a = (double*)malloc( N*N*sizeof(double) );
    double* b = (double*)malloc( N*N*sizeof(double) );
    double* c = (double*)malloc( N*N*sizeof(double) );

    for( int i = 0; i < N*N; ++i ) {
        a[i] = dist( mt );
        b[i] = dist( mt );
        c[i] = dist( mt );
    }

    double* device_a; cudaMalloc( (void**)&device_a, N*N*sizeof(double) );
    double* device_b; cudaMalloc( (void**)&device_b, N*N*sizeof(double) );
    double* device_c; cudaMalloc( (void**)&device_c, N*N*sizeof(double) );
    cudaMemcpy( device_a, a, N*N*sizeof(double), cudaMemcpyHostToDevice );
    cudaMemcpy( device_b, b, N*N*sizeof(double), cudaMemcpyHostToDevice );
    cudaMemcpy( device_c, c, N*N*sizeof(double), cudaMemcpyHostToDevice );

    cublasHandle_t handle; cublasCreate(&handle);
    double alpha = 1.0, beta = 1.0;

    cublasDgemm( handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha, device_b, N, device_a, N, &beta, device_c, N );

上のソースコードは、row-major形式の行列積を行うコードです。 cblasDgemmはcolumn-major形式で渡す必要があるので、行列積の順番を逆にすることで実質的にそれを実現しています。 row-majorとcolumn-majorの関係は転置に相当することを利用して、 AB+C = (B^TA^T+C^T)^Tのように変換したということです。

実行時間は、cublasDgemmの前後にcudaEventRecord関数を挟んでcudaEventElapsedTimeで計測しました。 実行時間は、cublasDgemmを51回呼び出し、その中で最も高速だったものを採用しました。 これは、他のプロセスが動いていたりして理想的ではない環境で実行していることによる速度低下を取り除くためです。 また、初回実行はライブラリの遅延読み込みの時間がかかり、まともな実行時間計測結果にならないという問題を取り除くことも意図しています。

精度は、別途計算した“正解”データとの最大相対誤差（max_relative_error）を計測しました。 今回の実験の範囲では、並列性に起因して再現性がないなどの現象は確認されず、どの手法でも毎回同じ結果が得られました。

“正解”データは、CPU上でdouble-double行列積を行うことによって求めました（真の正解とはずれがあると思います）。 double-double行列積は、乗算にTwoProdを、加算にsloppy加算を用いました。 そのコードは以下の通りです。

struct HiLo { double hi; double lo; };
HiLo TwoProd( double a, double b ) {
    double hi = a * b;
    double lo = std::fma( a, b, -hi );
    return { hi, lo };
}
HiLo FastTwoSum( double a, double b ) {
    double hi = a + b;
    double lo = a - hi + b;
    return { hi, lo };
}
HiLo TwoSum( double a, double b ) {
    double hi = a + b;
    double tmp = hi - a;
    double lo = a - (hi - tmp) + (b - tmp);
    return { hi, lo };
}
HiLo DD_MulAdd( double a, double b, double dh, double dl ) {
    HiLo mul = TwoProd( a, b );
    HiLo add = TwoSum( dh, mul.hi );
    double lo = mul.lo + add.lo + dl;
    return FastTwoSum( add.hi, lo );
}

    for( size_t i = 0; i < N; ++i )
    for( size_t k = 0; k < N; ++k )
    for( size_t j = 0; j < N; ++j )
    {
        HiLo sum = DD_MulAdd( a[i*N+k], b[k*N+j], dh[i*N+j], dl[i*N+j] );
        dh[i*N+j] = sum.hi;
        dl[i*N+j] = sum.lo;
    }

測定対象のエミュレーション手法は以下の通りです。

  • cuBLAS DGEMMそのまま
  • cuBLAS DGEMMエミュレーション・performant・精度0～128
  • cuBLAS DGEMMエミュレーション・eager・精度0～128
  • ozIMMU スライス数3～18およびauto*1
  • ozIMMU/errfree スライス数3～18およびauto
  • ozIMMU/nearest スライス数3～18およびauto
  • ozIMMU/nearest+errfree スライス数3～18およびauto

ozIMMUは、執筆時点で最新版の08eea9231729d54dbfd92955f2cbfc21ec236856（2025-12-02）を使いました。 三種の改良版は、d458fec19e29279556885a9da968d6fe0ac51058（2025-05-06）を使いましたが、これは動作しないため、ozIMMU_get_function_pointer() arguments are flipped · Issue #1 · RIKEN-RCCS/accelerator_for_ozIMMU · GitHubで示されているような修正を加えました。 なぜ本家にコントリビュートせず、差し替えソースコードを配布（しかも本家に追随できていない）みたいな運用になっているのでしょう？？

（追記：2025-12-07に三種の改良版がアップデートされたようですが、記事公開直前で気づかなかったため再実験はしていません）。

暖房を26度に設定し、Windowsからのディスプレイ出力がない状態で実験しました。

N = 1024

N = 1024の時の実験結果を図1に示します。 ただしozIMMU四種については、実験結果の図が見づらくなるのを防ぐため、精度が全く変わらず単に遅くなっただけのスライス数12～18はプロットしませんでした。 cuBLASのエミュレーションも同様に、精度が全く変わらず単に遅くなっただけのスライス数9以上はプロットしませんでした。 横軸はmax_relative_errorの二進対数に負号をつけたもので、右に行くほど高精度です。 おおよそ「何ビット精度か」というのを表していると考えてよいです。 例えば、横軸の値が27であれば、最大の相対誤差が2⁻²⁷であるということで、おおよそ27ビット精度ということになります。 縦軸は、本当にDGEMMを実行したと仮定した時に必要な2N³ FLOPを実行時間で割った、疑似的なFLOPS数で、上に行くほど高速です。 なお、RTX 4090の倍精度演算の理論性能は、1.29 TFLOPSです。

図1: N=1024の時の速度と精度の比較

cuBLASのエミュレーションは、DGEMMとほぼ同等の速度でした。 精度を下げると最大の相対誤差が増える一方でやや速くなるため、DGEMMに退化しているわけではなく、本当にint8テンソルコアを使っていると思います。

一方でozIMMU四種は、8スライス以上でDGEMMより速いのに精度が高い、という結果になりました。 9スライス以上では精度が変わっていませんが、これはdouble-doubleで計算した“正解”の方が間違っていることを示唆しているかもしれません。

N = 2048

N = 2048の時の実験結果を図2に示します。 図の見方は図1と同様です。

図2: N=2048の時の速度と精度の比較

cuBLASのエミュレーションは、DGEMMを超える性能となりました。 一方で、依然としてozIMMU四種の方がcuBLASのエミュレーションより高速です。 cuBLASのエミュレーションは、何かものすごく大きな固定オーバーヘッドがあるように見えます。

N = 4096

N = 4096の時の実験結果を図3に示します。 図の見方は図1と同様です。

図3: N=4096の時の速度と精度の比較

cuBLASのエミュレーションが、ozIMMU四種の性能とほぼ同じになりました。 固定オーバーヘッド（？）のせいでスライス数が小さい時は依然としてozIMMUよりも遅いようです。

N = 8192

N = 8192の時の実験結果を図4に示します。 図の見方は図1と同様です。

図4: N=8192の時の速度と精度の比較

スライス数4以上ではozIMMUとほぼ同等の精度-速度トレードオフが見られます。

N = 16384

N = 16384の時の実験結果を図5に示します。 図の見方は図1と同様です。

図5: N=16384の時の速度と精度の比較

やはり4スライス以上ではozIMMUとほぼ同等の精度-速度トレードオフが見られます。

なおcuBLASのエミュレーションでは、N=16384かつスライス数が2の時、Windowsが落ちる挙動が頻繁に発生しました。 おそらく電力を消費しすぎて電源の保護回路が働いたのだと思います。 火事になったりしても怖いので正確な条件は調べていませんが、そもそもこんな計算をコンシューマー向けGPUでやるべきではないのかもしれません。

スライス数あたりの精度

結果を見ると、cuBLASのエミュレーション6スライス版が、ozIMMU7スライス版と同等の精度を達成しています。 ozIMMUは6ビット毎、ozIMMU/nearestは7ビット毎に分割している一方、cuBLASのエミュレーションは8ビット毎に分割しているので、こういうことになっているのだと思います。

8ビット毎に分割するのは、おそらく下の方のビットから切り出して行っていると思います。 上の方のビットから切り出そうと思うと、特殊な丸めが必要になるからです。 ozIMMUは切り捨て丸めで6ビット毎の分割を、ozIMMU/nearestは最近接丸めで7ビットの分割を実現しています。 一方、8ビット毎に分割する場合、int8で表現できる数は-128～127と端（128）が含まれていないので、切り出すときの丸めに工夫が必要です。 下位桁を最大にすることで \frac{127}{2^8} + \frac{127}{2^{16}} + \frac{127}{2^{24}}+ \cdots = \frac{127}{255} ULPを実現できる一方、下位桁を最小にすると \frac{-128}{2^8} + \frac{-128}{2^{16}} + \frac{-128}{2^{24}} + \cdots = \frac{-128}{255} ULPを実現できます。 よって、最近接丸め（0.5で切り上げか切り下げがが変わる）ではダメで、 \frac{127}{255}で切り上げか切り下げかが変わる丸めが必要とされます。 下の方のビットから切り出すときはこのような工夫は特に必要ないため、下の方のビットから切り出していると思います。

結局どれを使えばいいのか？

一様乱数で入力行列を使った今回の実験の範囲では、RTX 4090に関しては以下のようなことが言えると思います。←ChatGPTに入れておいた方がいいよって言われて追加した文

N=2048くらい以上で倍精度相当の行列積を行いたい場合、四種のozIMMUとcuBLASのエミュレーションは同等の速度と精度でした。 なのでこの場合には、外部ライブラリが不要で手軽なcuBLASのエミュレーションモードを使用すればよいでしょう。

一方でそれ以外の場合、つまりN=1024程度の小さな行列積を高速化したい場合や、倍精度だと過剰精度なので少ないスライス数で高速に計算したい場合は、四種のozIMMUはcuBLASのエミュレーションよりも速いです。 四種のozIMMUの中ではozIMMU/errfreeまたはozIMMU/nearest+errfreeが精度-速度トレードオフが優れているので、これを使うのがよさそうです。

また、四種のozIMMUはいずれもオープンソース（MITライセンス）であるという点も見逃せない点です。 特にozIMMU/errfreeまたはozIMMU/nearest+errfreeはcuBLASのエミュレーションと同等以上の性能を達成しているため、手法の研究や改良を行いたい場合はこれをもとにするのがよさそうです。

*1:変換誤差上限がデフォルトの0なのでDGEMMが呼び出されるだけだった



以上の内容はhttps://lpha-z.hatenablog.com/entry/2025/12/07/231500より取得しました。
このページはhttp://font.textar.tv/のウェブフォントを使用してます
不具合報告/要望等はこちらへお願いします。
モバイルやる夫Viewer Ver0.14