Recap: Go * SIMDで高速化するベクトル検索 ~ ルーフラインモデルでSIMDが効く境界を探れ! ~

Table of Contents

はじめに

Go Conference 2026 のショートワークショップ Go * SIMDで高速化するベクトル検索 ~ ルーフラインモデルでSIMDが効く境界を探れ! ~ に参加してきた。オーガナイザはponさんだった。

資料は以下の通り:

しかし、前提知識がほぼない状態で参加したため分からないこと多く、本メモでは自己理解した時に詰まった部分や考えたことをメモとして残す。基本的な説明は上記スライドやリポジトリを参照すると良いと思う。

前提と計測環境

  • 探索は全探索(10 万件すべてと内積)
  • クエリは 1 本ずつ
  • 1 コアのみ。並列化は付録扱い
  • ワークスペースで割
  • データはメモリ上の配列、384 次元 float32

つまり、1 回の検索ごとに、この 153.6 MB(10 万[件] * 384 [次元] * 4 [byte]) を全部メモリから読むことになる。これらの演算をワークショップ内で改善させていくことが目的だった。

SIMD と Go のパッケージ

SIMD

SIMD(Single Instruction, Multiple Data)は、1 つの命令で複数のデータをまとめて処理する CPU の機能。Multiple Instructions, Multiple Data だと MIMD になる。詳しくはパタヘネを読んでください。

普通に書いたループは値を 1 個ずつ処理する。a[i]b[i] を 1 個取り出して 1 回掛けて足す、を繰り返すので、1 命令で float32 を 1 個しか扱わない。これをスカラ処理と呼ぶ。

対して CPU の中には普通の変数より大きなベクトルレジスタがあり、256 [bit] のレジスタには float32 (32 [bit]) が 8 個入る。8 個詰めて掛け算命令を 1 回実行すると、8 個分の掛け算が同時に走る。この詰めた 1 個ぶんの区画をレーンと呼ぶ。

つまり 384 次元の内積なら、スカラで 384 回かかる掛け算が 48 回で済むので、理論上は 8 倍になるはず(実際はそう上手くはいかないが…w)。

scalar-vs-simd.png

https://github.com/po3rin/gocon2026-simd-search/より引用

Go の SIMD パッケージ

Go 1.27 では GOEXPERIMENT=simd を付けてビルドすると SIMD のパッケージが使える。アセンブリも cgo も書かずにベクトル命令を直接扱えて、**メソッド呼び出しがほぼそのまま 1 つの CPU 命令になるのが嬉しいところだと思う。

そして、パッケージはアーキテクチャ固有用の (archsimd) とポータブル用の(simd) の2つある。#73787によると、ハードが持つ SIMD 演算は本質的に複雑で非ポータブルであるが、Go 自身はシンプルでポータブルな思想があるため、このような棲み分けになったらしい。

simd/archsimdsimd
型名Float32x8(レーン数が入る)Float32s(入らない)
レーン幅型で固定実行時に CPU が決める (4 / 8 / 16)
CPU 機能チェック必須不要
使える命令そのアーキの命令ほぼ全部全アーキ共通に持てる演算だけ
対象アーキamd64 / arm64 / wasm で API が別1 ソースで全アーキ

関連する Issue は以下の通り:

  • #78979(open): AMD64 の archsimd をデフォルト有効にする提案。通れば GOEXPERIMENT が要らなくなる
  • #79781(open): ARM64 の SVE 対応。SVE はレジスタ幅が実装依存なので、ポータブル API の設計とも絡みそう
  • #79413(open): 標準ライブラリ crypto の手書きアセンブリを Go の SIMD に置き換える。これが本来の狙いだと思う
  • #80857(closed / completed): min/max の畳み込みループをコンパイラが自動ベクトル化する

ただし、ワークショップでは、ポータブルな simd だと Stage 2 以降で必要な命令が足りないため archsimd で進めていく。

性能の上限を判定する

flop (op)

FLOP は FLoating-point OPeration の略で、浮動小数点演算 1 回を意味する。例えば、x := a + b は 1 [flop]、z := a * b + c なら 2 [flop] 。浮動小数点数ではない場合は単に OP を用いる。

内積のループで数えると、1 要素あたり掛け算 1 回と足し算 1 回で 2 flop。1 回の検索なら以下になる:

2 [flop] * 384 [次元] * 100,000 [件] = 76.8M [flop]

この 76.8M を実行時間で割れば GFLOP/s が出る。ベンチではこの割り算をしている。

なお FMA(Fused Multiply-Add、a * b + c を 1 命令で行う命令)は1 命令だが 2 flopと数えていた。

算術(演算)強度

算術強度(Arithmetic Intensity, AI)は flop / byte、つまり「メモリから 1 バイト運ぶごとに何回計算するか」を意味する指標。他の文献をあたると「演算強度」という単語が用いられていることが多かったが、このメモではワークショップ資料にあわせて算術強度という言葉を用いる。

この単位系の数値を持つ理由は、CPU に計算する能力(flop/s)とデータを運ぶ能力(byte/s)という独立した 2 つの指標があり、どんな処理も必ずどちらかが先に枯れるためだと思う。今回の内積なら 2 [flop] / 4 [byte] = 0.5 [flop/byte]。「4 バイトも運んできて、計算はたった 2 回」で、かなり低い側になる。

ルーフラインモデル

高速化の手を打つ前に 「いま何が理由で詰まっているか」と「追加でどの程度伸びる余地があるか」 を明確にするためのモデル。

roofline-concept.png()

https://github.com/po3rin/gocon2026-simd-search/より引用

Claude Code くん曰く

初出は Roofline: An Insightful Visual Performance Model for Floating-Point Programs and Multicore Architectures(Williams et al. 2008)。翌年に CACM 52(4), pp.65-76 へ「Roofline: An Insightful Visual Performance Model for Multicore Architectures」として公刊されており(タイトルから “Floating-Point Programs and” が落ちている)、引用されるのは通常こちら。

とのこと。

こちらのルーフラインモデルでは、性能の上限が以下で決まる:

上限 [flop/s] = min(演算ピーク, 算術強度 [flop/byte] * メモリ帯域 [byte/s])

横軸に算術強度、縦軸に性能を取ると、左が右上がりの斜線(メモリ律速)、右が水平線(演算律速)の折れ線になる。屋根の形をしているのが名前の由来で、どのコードもこの線より上には行けない

そして min() の形が、そのまま打つべき手を教えてくれる:

点の位置詰まっているもの打つ手ワークショップでの該当部分
斜線に張り付いているデータ転送運ぶ量を減らす / 運んだデータを使い回す(算術強度を上げる)Stage 1b → 2
水平線に張り付いている演算実装効率を上げる(SIMD、アキュムレータ分割)Stage 2
どちらからも遠い実装の無駄(依存連鎖など)まず素直に実装効率を上げるStage 0 → 1a

また、折れ曲がる点をリッジと呼ぶ。2 本の線が交わる点なので等号を解くだけで出る:

算術強度 * メモリ帯域 = 演算ピーク
→ リッジ = 演算ピーク / メモリ帯域

単位を見ると割り算になる理由が分かる。[flop/s] / [byte/s] で秒が約分され、算術強度と同じ flop/byte になる。コードの要求(算術強度)とマシンの供給(リッジ)が同じ単位に揃うので、そのまま大小比較できるのが嬉しいポイントだと思う。

演算ピークの算出方法

演算ピークの出し方は理論値と実測値の 2 つある。理論値は掛け算 4 つで出る:

理論ピーク = FMA ユニット数 [op/cycle] * レーン数 [lane/op] * FMA 1 個の flop [flop/lane] * クロック [cycle/s]

EPYC 7763 なら 2 * 8 * 2 * 3.5GHz ≒ 110 GFLOP/s になる。一方、実測値はメモリに触らずレジスタ上で FMA だけを回し続けるベンチで出す。教材の make roofline-ceiling がそれで、同じ 4 コア Codespace での教材の実測値は以下だった:

演算ピーク  25.59 GFLOP/s
メモリ帯域  20.80 GB/s(読むだけ)
→ リッジ    25.59 / 20.80 ≒ 1.23 [flop/byte]

ここで、実測が理論の 1/4 にしかならないことが分かる。原因はハードの限界ではなく、Go のコンパイラがアキュムレータをレジスタに置き続けられず、FMA のたびにスタックと往復させるためであると考えられる。これを register spill と呼ぶ。

具体的なアセンブリを見ると以下のようになっており、FMA 1 個ごとに load と store が付くので、load/store ポートが先に飽和することが分かる。

VMOVDQU 0x1b8(SP), Y2    ; スタックからレジスタへ load
VFMADD213PS Y1, Y0, Y2   ; FMA
VMOVDQU Y2, 0x1b8(SP)    ; レジスタからスタックへ store

この問題は既知の issue golang/go#76969 として報告されているが、closed as not planned 扱いであった。

この全探索の上限は 10 GFLOP/s 前後

道具が揃ったので、実際に当てはめてみる。

リッジ   = 演算ピーク / メモリ帯域 = 25.59 / 20.80 ≒ 1.23 [flop/byte]
算術強度 = 0.5 [flop/byte]  <  リッジ 1.23          → メモリ律速
上限     = 算術強度 * メモリ帯域 = 0.5 * 20.8 ≒ 10 [GFLOP/s]

コードを 1 行も書かないうちに「この全探索はどう頑張っても 10 GFLOP/s 前後が天井」という見積りが立つ。 以降の Stage では、この数字と実測を比べていく。

ワークショップ

ここから、ワークショップの内容に入る。基本的にはリポジトリに参考コードがあるので、そちらを参照すると良いと思う。

Stage 0: スカラ基準

最適化は基準値から始まる:

func Dot(a, b []float32) float32 {
    var sum float32
    for i := range a {
        sum += a[i] * b[i]
    }
    return sum
}

func (ix *Index) Search(q []float32, k int) []Result {
    t := newTopK(k)
    for id := 0; id < ix.N; id++ {
        t.push(id, vec.Dot(q, ix.Vec(id)))
    }
    return t.results()
}
BenchmarkSearchNaive-4    68    34877996 ns/op    4403.92 MB/s    0.5000 AI(flop/byte)    2.202 GFLOP/s    153.6 MB/query

34.878 ms、2.202 GFLOP/s。見積りの 10 GFLOP/s に対して 22 % しか出ていない。 帯域に直すと 4.40 GB/s で、これも 20.8 GB/s の 2 割しかない。つまりメモリ帯域に届く前に詰まっている

ここでの課題は2つ。

  • 1.浮動小数点演算を 1 つずつしか実行できないこと
  • 2.sum の演算が逐次処理で、前の演算が終わるまで次を実行できないこと

次の Stage でこれらの課題解決に取り組む。

archsimd の読み方

1. 型が「データの形」を表す

型名そのものが「何 bit 幅のレジスタに、どの型を何個詰めるか」を意味する。型を選ぶことが、使う命令幅とレーン数を選ぶことになる。

var a archsimd.Float32x8    // float32 を 8 レーン  = 256 [bit] (AVX2)
var b archsimd.Float32x16   // float32 を 16 レーン = 512 [bit] (AVX-512)
var c archsimd.Uint64x4     // uint64 を 4 レーン

2. メソッドが 1 つの CPU 命令に対応する

各メソッドはベクトル命令にほぼ 1 対 1 で変換されるので、メソッド名から出てくる機械語の見当がつく。

va := archsimd.LoadFloat32x8(xs)   // スライス -> レジスタ (ロード)
va = va.MulAdd(vb, acc)            // 積和     -> VFMADD
xo := vc.Xor(vd)                   // XOR      -> VPXOR
po := xo.OnesCount()               // popcount -> VPOPCNTQ
va.Store(xs)                       // レジスタ -> スライス (ストア)

3. 使う前にその CPU が対応しているか確かめる

これが一番怖い部分だと思う。未対応の CPU でメソッドを呼ぶと、不正命令 (SIGILL) でプロセスごと落ちる。recover できる panic にもならないので、実行時に機能フラグでガードする必要がある。

var hasSIMD = archsimd.X86.AVX2() && archsimd.X86.FMA()

x86 では MulAdd が使う FMA 命令が AVX2 に含まれておらず別の拡張として提供されているので、ここでは 2 つ確認している。arm64 では FMA 相当の命令が Neon 自体に含まれるため、この区別はない。

さらに厄介なのが、この機能フラグ自体のバグが複数報告されていることだった。正しくチェックを書いたつもりでも、検出側が間違っていれば意味がない。

  • #78772(open): X86.AVX512() が darwin/amd64 で常に false になる
  • #79437(open): internal/cpuHasGFNIHasAVX512F の内側に誤って入れ子になっている
  • #81008(closed / completed): GODEBUG=simd=+ が CLMUL の有無を反転して報告していた

そして重要なのは、SIMD を使えばすぐ速くなるわけではないということ。効くかどうかは「何が性能を抑えているか」で決まる。それを見極めるのが次章のルーフラインモデルである。

Stage 1a: AVX2 で 8 個まとめる

まず Stage 0 の課題 1(浮動小数点演算を 1 つずつしか実行できない)に手を入れる。Float32x8 でスライスから 8 要素をベクトルレジスタにロードし、MulAdd(FMA)でアキュムレータに足し込む形に置き換える。

func Dot(a, b []float32) float32 {
    var acc archsimd.Float32x8
    for len(a) >= 8 {
        va := archsimd.LoadFloat32x8(a)
        vb := archsimd.LoadFloat32x8(b)
        acc = va.MulAdd(vb, acc)
        a, b = a[8:], b[8:]
    }

    var buf [8]float32
    acc.Store(buf[:])
    archsimd.ClearAVXUpperBits()
    sum := buf[0] + buf[1] + buf[2] + buf[3] + buf[4] + buf[5] + buf[6] + buf[7]
    for i := range a {
        sum += a[i] * b[i]
    }
    return sum
}
BenchmarkSearchNaive-4    238    10174585 ns/op    15096.44 MB/s    0.5000 AI(flop/byte)    7.548 GFLOP/s    153.6 MB/query

34.878 ms から 10.175 ms で 3.43 倍、7.548 GFLOP/s まで来た。ただし理論上は 8 倍のはずである。届かないのは Stage 0 の課題が丸ごと残っているからで、acc = va.MulAdd(vb, acc) は前の acc を読んで新しい acc を書くため、鎖が 8 本束になっただけで、鎖であることは変わっていない

Stage 1b: アキュムレータを 2 本にする

アキュムレータを 2 本に分けると、互いの結果を待たなくなる:

func Dot(a, b []float32) float32 {
    var acc, acc2 archsimd.Float32x8
    for len(a) >= 8*2 {
        va := archsimd.LoadFloat32x8(a)
        vb := archsimd.LoadFloat32x8(b)
        acc = va.MulAdd(vb, acc)

        va = archsimd.LoadFloat32x8(a[8:])
        vb = archsimd.LoadFloat32x8(b[8:])
        acc2 = va.MulAdd(vb, acc2)
        a, b = a[16:], b[16:]
    }

    var buf [8]float32
    acc.Add(acc2).Store(buf[:])
    archsimd.ClearAVXUpperBits()
    sum := buf[0] + buf[1] + buf[2] + buf[3] + buf[4] + buf[5] + buf[6] + buf[7]
    for i := range a {
        sum += a[i] * b[i]
    }
    return sum
}
BenchmarkSearchNaive-4    237    9825801 ns/op    15632.31 MB/s    0.5000 AI(flop/byte)    7.816 GFLOP/s    153.6 MB/query

10.175 ms から 9.826 ms で、3.5 % しか速くなっていない。

正直に書くと期待外れだった。教材は内積単体で 348 ns から 55.2 ns(6.3 倍)になると書いていて、アキュムレータ分割は効くはずなのだ。理由は、全探索の粒度ではもう別の上限に当たっているからだと考えられる。帯域が 15.6 GB/s まで来ていて DRAM 側が効き始めている。ここが Stage 1 の到達点ということだろう。

そのため、次はルーフトップモデルにおける点を上に動かせないなら右に動かす方針をとる。つまり算術強度を上げるのだが、分数なので方法は 2 つしかない。分母(運ぶ量)を小さくするか、分子(計算回数)を大きくするかだ。今回は前者で進める。

Stage 2: int8 量子化でバイトを削る

各要素を float32(4 byte)から int8(1 byte)に表し直す。これを量子化と呼ぶ:

算術強度 = 2 [op] / 1 [byte] = 2.0 [flop/byte]   (4 倍・リッジ 1.23 の右へ)
転送量   = 153.6 MB → 38.4 MB/query              (1/4)

使うのは対称スカラ量子化だ。比率を変えずに値を縮めるイメージで、ベクトルごとに絶対値の最大値 maxAbs を測り、scale = maxAbs / 127 として各要素を round(v / scale) で -127〜127 に丸める。以下で示す実装を見ればイメージがつきやすいだろう。

また、こちらの量子化は非可逆圧縮なので、外れ値が 1 つあると他の要素の分解能が落ちる。しかし、ベクトル検索で知りたいのは「上位 N 件」という大小関係なので、大きな問題にはならないと考えられる。丸め誤差はおおよそ独立でゼロ平均なので、384 次元を足し合わせる間に絶対誤差は $\sqrt{384}$ 程度で増える一方、信号は 384 で増える。相対誤差が薄まる方向に働くはずだと思う。また、float から int の演算になるので、評価指標も FLOP/s から(単なる)op/s になる。

func Dot(a, b []int8) int32 {
    // アキュムレータは Int32x8 で固定。DotProductPairs の戻り型が Int32x8 なので他に選べない
    var acc0, acc1 archsimd.Int32x8
    for len(a) >= 16*2 {
        acc0 = acc0.Add(archsimd.LoadInt8x16(a).ExtendToInt16().
            DotProductPairs(archsimd.LoadInt8x16(b).ExtendToInt16()))
        acc1 = acc1.Add(archsimd.LoadInt8x16(a[16:]).ExtendToInt16().
            DotProductPairs(archsimd.LoadInt8x16(b[16:]).ExtendToInt16()))
        a, b = a[32:], b[32:]
    }
    if len(a) >= 16 {
        acc0 = acc0.Add(archsimd.LoadInt8x16(a).ExtendToInt16().
            DotProductPairs(archsimd.LoadInt8x16(b).ExtendToInt16()))
        a, b = a[16:], b[16:]
    }

    var buf [8]int32
    acc0.Add(acc1).Store(buf[:])
    archsimd.ClearAVXUpperBits()
    sum := buf[0] + buf[1] + buf[2] + buf[3] + buf[4] + buf[5] + buf[6] + buf[7]
    for i := range a {
        sum += int32(a[i]) * int32(b[i])
    }
    return sum
}

スコアの復元には scale を掛け戻す:

t.push(id, qScale*ix.Scales[id]*float32(Dot(q8, ix.Code8(id))))

これは量子化の定義をそのまま代入した形になっている。

$$ \sum_i d_i q_i \approx \sum_i \bigl(\mathrm{code8}_i \cdot \mathrm{dScale}\bigr)\bigl(\mathrm{q8}_i \cdot \mathrm{qScale}\bigr) = \mathrm{dScale} \cdot \mathrm{qScale} \cdot \sum_i \mathrm{code8}_i \cdot \mathrm{q8}_i $$

ここで $d_i$ は DB ベクトルの、$q_i$ はクエリの元の float32 の要素で、$\mathrm{code8}_i$ と $\mathrm{q8}_i$ がそれぞれを量子化した int8、$\mathrm{dScale}$ と $\mathrm{qScale}$ がその scale である。

Claude Code くん曰く、↓らしい。

$\mathrm{dScale}$ と $\mathrm{qScale}$ が総和の外に出られるのは scale が「要素ごと」ではなく「ベクトルごとに 1 つ」だから成り立つ。要素ごとに違う scale を持っていたら整数の積和命令は使えない。教材が対称(ゼロ点なし)のスカラ量子化を選んでいるのはこの理由だと考えられる。

BenchmarkDotInt8Naive-4     6374116        371.9 ns/op
BenchmarkDotInt8SIMD-4     75458914         31.38 ns/op
BenchmarkSearchInt8-4           604    3965408 ns/op    9683.75 MB/s    2.000 AI(flop/byte)    19.37 Gop/s    38.40 MB/query

内積単体は 371.9 ns から 31.38 ns で 11.85 倍。float32 の SIMD 内積より速い。全探索は 9.826 ms から 3.965 ms で 2.48 倍になった。

精度に関しては、上位10件の正しさを確認するための Recall@10 を採用しており、 Recall@10 = 0.948 であった。これは、float32 の上位 10 件を正解としたとき、平均 9.5 件が一致するという意味であって問題なさそうだ。

一方で、運ぶ量を 1/4 にしたのに 4 倍にはなっていない。これは、リッジを越えたので、律速がメモリから int8 内積の命令スループットに移ったためと考えられる。実際に、教材では「時間を決めているのがメモリから内積の計算に変わった」と書かれている。

Stage 3 では、次なる高速化として、int8 でも 10 万件 * 384 byte = 38.4 MB で、L3 キャッシュに全て乗らないという課題に取り組む。

Stage 3: バイナリ量子化でキャッシュに乗せる

各要素を符号 1 bit だけにする。正なら 1、そうでなければ 0。

表現1 ベクトル10 万件キャッシュに
float321536 byte153.6 MB乗らない
int8384 byte38.4 MB乗らない
1 bit48 byte4.8 MB乗る

とはいえ、本当に正負の符号だけで検索できるのだろうか?調べてみると、答えは、精度が落ちるが Yes らしい。なぜなら、検索のベクトル演算において「意味が近い = 向きが近い」であり、各軸について正側か負側かを記録する符号は、向きの粗い記述になっていると考えられるためとのこと。

ここで嬉しいのが、ベクトルの内積計算の代わりにビットの一致のみを確認すれば良くなったことである。したがって、元々 SIMD で高速化していた部分はハミング距離を確認すれば良くなり、ハミング距離は XOR をして各ビットを数え上げるだけで良くなるので高速化が期待できる。

ちなみに、この考えの起源は Charikar の Similarity estimation techniques from rounding algorithms(STOC 2002, pp.380-388)で、後に SimHash という通称で広まった論文とのこと。

Claude Code くん曰く↓らしい。

2 本のベクトル $u$, $v$ のなす角を $\theta$ とすると、ランダムな超平面の法線 $r$ に対して $\Pr[\operatorname{sign}(r \cdot u) \neq \operatorname{sign}(r \cdot v)] = \theta / \pi$ が成り立つ、という結果である。ビットが食い違う割合がそのまま角度に比例するので、ハミング距離を数えるだけでコサイン類似度を推定できる。

func Quantize(v []float32, out []uint64) {
    for i := range out {
        out[i] = 0
    }
    for i, x := range v {
        if x > 0 {
            out[i/64] |= 1 << (i % 64)
        }
    }
}

func Hamming(a, b []uint64) int {
    var d int
    for i := range a {
        d += bits.OnesCount64(a[i] ^ b[i])
    }
    return d
}

func (ix *Index) Code(id int) []uint64 {
    return ix.Codes[id*ix.Words : (id+1)*ix.Words]
}

bits.OnesCount64 はスカラの POPCNT 命令 1 個にコンパイルされ、1 命令で 64 次元ぶん数える。

BenchmarkSearchNaive-4              63    37578911 ns/op    4087.40 MB/s    0.5000 AI    2.044 GFLOP/s    153.6 MB/query
BenchmarkSearchBinary-4           2318      899242 ns/op    5337.83 MB/s                                   4.800 MB/query

同じ実行内の比較で 37.579 ms から 0.899 ms、41.8 倍である。

しかし、やはり Recall が81%程度低下している。

Recall@10: binary=0.180

Stage 4: float32 SIMD で rerank して精度を戻す

L3キャッシュでハンドリングするため、次のようなStage3 と SIMD のハイブリッドアプローチをとる。

  1. 1 bit のハミング距離で 10 万件すべてを比べ、上位 100 件(返したい 10 件の 10 倍)に絞る
  2. その 100 件だけ float32 のベクトルを読み、Stage 1 の SIMD 内積で採点し直して上位 10 件を返す。

2 段目が安い理由が肝で、読む float32 は 100 件 * 1536 byte = 154 KB しかなくキャッシュに乗る。10 万件全部に float32 内積をすると 9.4 ms かかるが、100 件なら 1/1000 である。

BenchmarkSearchBinary-4              2318     899242 ns/op
BenchmarkSearchBinaryRerank-4        2530     969063 ns/op

0.899 ms から 0.969 ms で、増加は 0.07 ms だけだ。また、Recall@10 = 0.868 となっており、良さそうである。

ただし、この 0.868 は 2 万件で測った数字だった。ベンチと同じ 10 万件だとどうなるかは教材の数字を自分で確かめたに分けて書いた。

教材の数字を自分で確かめた

ここまでは教材の説明を追ってきたが、調べていくうちに教材が載せている数字そのものが気になった。以下はワークショップ後に自分で測り直した部分である。

Recall 0.868 は 2 万件での数字だった

ここで引っかかった。教材のまとめ表は速度と Recall を同じ行に並べているが、速度は 10 万件のベンチ、Recall は 2 万件のテストで測られている(脚注に明記はある)。factor は 10 で固定なので、件数が増えたら候補 100 件で足りるのだろうか。

気になったので、教材のコードをそのまま使って N と factor を振ってみた。

// 教材の clusteredIndex / SearchNaive / SearchBinaryRerank をそのまま使う
func TestRecallScaleSweep(t *testing.T) {
	const dim, nq, k = 384, 50, 10
	for _, n := range []int{20_000, 100_000} {
		ix, sample := clusteredIndex(n, dim, 128, 7, 8)

		qs := make([][]float32, nq)
		for i := range qs {
			qs[i] = sample()
		}
		exact := make([]map[int]bool, nq)
		for i := range qs {
			exact[i] = idSet(ix.SearchNaive(qs[i], k)) // float32 全探索の上位 k を正解とする
		}

		for _, f := range []int{10, 20, 50, 100} {
			var r float64
			for i := range qs {
				r += overlap(exact[i], ix.SearchBinaryRerank(qs[i], k, f))
			}
			t.Logf("N=%7d factor=%3d Recall@10=%.3f", n, f, r/float64(nq*k))
		}
	}
}

結果はこうなった。

Nbinary 単体factor=10factor=20factor=50factor=100
20,000(教材の Recall テスト)0.1800.8681.0001.0001.000
100,000(教材のベンチ)0.0880.4280.6140.9301.000

N=20,000 の 0.180 / 0.868 は教材が載せている数字そのものなので、条件は再現できている。その上で ベンチと同じ 10 万件にすると 0.868 が 0.428 に落ちた

なぜ落ちるのかも測ってみた。各クエリについて「float32 での真の上位 10 件」が「1bit のハミング距離順」で何位に沈んでいるかを数えると、次のようになる。

N = 20,000N = 100,000
中央値38 位137 位
90 パーセンタイル108 位463 位
最悪164 位777 位

クラスタ数は 128 で固定なので、N が 5 倍になるとクラスタあたりの件数が約 156 件から約 780 件に増える。1bit のハミング距離はクラスタを当てるところまでしかできず、その中の順位は当てられないので、同一クラスタの 780 件が団子になって真の上位 10 件が埋もれる。最悪順位 777 位がクラスタサイズとほぼ一致しているのが分かりやすい。

つまり factor はクラスタあたりの件数(≒ N)に比例して増やす必要がある。10 万件で 0.87 を出したいなら factor=50(候補 500 件)が必要で、そのとき 0.930 になる。

救いは、これが安い対処だということ。1bit の走査は factor に関係なく必ず全件見るので、増えるのは 2 段目で読む float32 だけだ。

factor候補数2 段目で読む float32索引サイズ
10100150 KB4.8 MB
50500750 KB4.8 MB

索引は 4.8 MB のままなので、Stage 3 で得た「キャッシュに乗る」という成果は壊れない。

ただし 1 点注意があって、教材の topK は挿入ソートで実装されていて、コメントに「k は小さい(〜100)前提」と書かれている。factor=50 にすると 1 段目の newTopK(k*factor) が 500 件になり、この前提を超える。10 万回の push のうち上位入りするたびに最大 500 要素のシフトが走るので、絞り込み自体が遅くなる可能性がある。ここは測っていないので要確認。効いているならヒープに変えるのが定石だと思う。

詰まったポイント

1. archsimd と simd のどちらを使うのか

Go 1.27 には SIMD のパッケージが 2 つある。アーキテクチャ固有の simd/archsimd と、ポータブルな simd だ。後者は型名にレーン数がなく(Float32s)、幅は実行時に決まり、CPU 機能チェックも要らない。前者は型で幅が固定で、機能チェックを怠ると SIGILL でプロセスが死ぬ。

最初は「ポータブルなほうが遅いのでは」と思ったが、機械語を見て納得した。ポータブル版は幅ごとに関数が複製される:

vec.DotPortable(SB)           ← 幅を見て振り分けるだけのスタブ
vec.DotPortable@simd0(SB)     ← エミュレーション版
vec.DotPortable@simd128(SB)   ← 128 bit 専用
vec.DotPortable@simd256(SB)   ← 256 bit 専用
vec.DotPortable@simd512(SB)   ← 512 bit 専用

スタブは simd.maxVectorSize を読んで該当版を CALL するだけなので、各複製の中では acc.Len() はコンパイル時定数になる。コストはグローバル 1 回読みと数回の比較と 1 回の CALL だけで、ループ本体は archsimd 版と同じ命令になる。

この多版化はコンパイラの midway パッケージが担っていて、根拠は次のとおり。

  • 複製名の @simd<N> を組み立てている箇所: cmd/compile/internal/midway/rewrite.gofmt.Sprintf("%s@simd%d", ...))。同じ関数がディスパッチ用の switch 文も生成している
  • どの幅で複製するかを決めている箇所: cmd/compile/internal/midway/midway.gorewriteSizes()amd64 は {0, 128, 256, 512}、arm64 と wasm は {0, 128} を返す。手元で観測した 4 つの複製がちょうどこれに対応する(0 がエミュレーション版)
  • ディスパッチの判定に使う変数: src/simd/midway_common.gomaxVectorSize
  • 複製される側のコード(本記事で測った関数): internal/vec/dot_portable.go

自分で確かめるなら次のコマンドでシンボルが見える。

GOEXPERIMENT=simd GOARCH=amd64 go build -gcflags=-S ./internal/vec 2>&1 | grep 'TEXT.*DotPortable'

なお rewriteSizes() の arm64 のところには // this will change for SVE and cannot just be a size-based choice. というコメントが付いていて、#79781 の SVE 対応がこの設計に影響することが示されている。

では何が違うのか。ポータブル版には必要な命令がない:

  • int8 → int16 の幅拡張(Int8x16.ExtendToInt16())がない
  • int8 のペア積和(Int16x16.DotProductPairs())がない
  • popcount(Uint64x4.OnesCount())がない

しかも罠がある。simd.Int8s には Mulあるが、これは int8 → int8 だ。内積の積は最大 127 * 127 = 16129 で int8 に入らないので、コンパイルは通るのに結果が壊れる。

これらは将来追加されるのか

ポータブル版の提案スレッド #78902 に、何を入れるかの基準が明記されている。

In this version, the supported vector methods are those in the intersection of the wasm SIMD API and the current amd64 SIMD API … The portable API will be expanded over time by various architecture-specific APIs with emulations to fill in the intersection.

「wasm と amd64 の共通部分」が基準ということなので、各アーキの archsimd を実際に引いて共通部分を調べてみた。

GOEXPERIMENT=simd GOARCH=arm64 go doc simd/archsimd.Int8x16
GOEXPERIMENT=simd GOOS=js GOARCH=wasm go doc simd/archsimd.Int8x16
操作amd64arm64wasm共通部分
幅拡張 ExtendLo8ToInt16(int8→int16・同幅)入る
popcount OnesCount(バイト単位)入る
ペア積和 DotProductPairsなしなし入らない

幅拡張と popcount は 3 アーキすべてにあるので、方針どおりなら追加される見込みがありそうだ。ただし形は変わる。

  • ExtendToInt16()Int8x16(128 bit)から Int16x16(256 bit)を返す、レジスタ幅が倍になる操作である。#78902 は「Within a given execution of a program all the vectors have the same length in bits」と書いているので、この形はポータブルには持ち込めない。入るとしたら ExtendLo / ExtendHi 系(同じ幅でレーン数が半分)になり、amd64 で 1 命令だったものが Lo と Hi の 2 回に分かれる
  • popcount も、64bit 単位(Uint64x4.OnesCount、AVX-512 VPOPCNTDQ)は arm64 に無い。バイト単位なら 3 つ揃うので、入るとしたらそちらで、64 bit ごとに畳むには自分で足す必要がある

一方 DotProductPairs は amd64 にしかない。 これは記事の後半で触れる「arm64 版の int8 内積が 3 段構成になる」話の裏返しで、arm64 にそもそも相当命令が無いことが API にそのまま出ている。Stage 2 の int8 内積がポータブルに書けないのは一時的な未整備ではなく、ハードウェアの差そのものが理由だった。方針には「emulation で埋めながら拡張する」ともあるので理屈上は入りえるが、1 命令とエミュレーションで性能特性が変わるので優先度は低そうに見える。

なお #78902 の本文は 2026-04 時点のもので、すでに現行 API と名前が食い違っているAddPairsGroupedConcatAddPairsGroupedToSimd()Float32sFromArch[T]()LoadInt8SliceLoadInt8s など)。「小さく始めて拡張する」が実際に動いている証拠ではある。

使い分けの結論

判断は「使いたい命令がポータブル側にあるか」の一点になる。float32 の四則と FMA で足りるなら simd、整数の幅拡張・ペア積和・popcount が要るなら archsimd だ。

そして切り替えは関数ごとに分ける必要すらなかった。 #78902 が ToArch() という逃げ道を用意していて、式の途中で archsimd の型に落とせる。現行 API にも func (x Float32s) ToArch() any が存在する。つまり 「ポータブルで書いて、足りない箇所だけ ToArch() で落とす」が公式に想定された使い方だった。

2. レーンを 1 つの float32 に畳む API がない

アキュムレータに足し込んだ後、8 レーンを 1 つの float32 にする段で手が止まった。ReduceSum を探したが見つからないのだ。

調べた結果、amd64 には Float32x8 にも Float32x16 にも ReduceSum がない。arm64 の Float32x4 にもない(ReduceMaxReduceMin だけある)。唯一あるのは arm64 の Int32x4.ReduceSum()(ADDV)で、そのため int8 の arm64 版だけ Store なしで書ける。

float32 では Store してスカラで足すしかない。ポータブル版で書く場合は buf の大きさに注意が必要で、[8]float32 と決め打ちすると AVX-512 機(16 レーン)で足りないため、[16]float32 で確保して Store(buf[:n]) に渡すStorelen(s) >= Len() を要求し、短いと panic する。StorePart なら min(len(s), Len()) 個書いて書いた個数を返すので、そちらでも書ける。

issue を探したが、水平加算そのものを扱った専用の issue は見つからなかった。 似た形の「ポータブル API に操作が足りない」という報告としては #80149(open・Mask8s のビット抽出を公開してほしい)があり、ポータブル側は今も操作を足しながら固めている段階に見える。

ただし、想定された解決策は提案スレッドに書かれていた。 #78902 の本文が ToArch() でアーキ固有の型に落とす方法を示していて、しかも挙げられている例が水平和そのものだった。

func sum(x simd.Float32s) float32 {
	switch a := x.ToArch().(type) {
	case archsimd.Float32x8:
		a = a.AddPairsGrouped(a)   // 現行では ConcatAddPairsGrouped
		a = a.AddPairsGrouped(a)
		return a.GetLo().GetElem(0) + a.GetHi().GetElem(0)
	case archsimd.Float32x16:
	...
	}

つまり「ポータブル API に無いから諦めて Store する」ではなく、そこだけ ToArch() で降りるのが本来の書き方だった。自分は Store してスカラで足す方を選んだが、ポータブル性を保ったまま水平加算だけアーキ固有にできたことになる。

3. ClearAVXUpperBits を自分で呼ぶ必要がある

コードに出てくる archsimd.ClearAVXUpperBits() の意味が最初は分からなかった。中身は VZEROUPPER という命令 1 つで、SIMD で使ったレジスタの上位 128 bit をゼロに掃除する。

これを置かないと、SIMD からスカラの計算に戻る境界で Intel 機が大きく遅くなる。教材の付録には、これ 1 命令の有無で内積単体が 167.4 ns と 23.4 ns(7.1 倍)変わった実測が載っている。YMM レジスタの上位 128 bit に値が残った状態(dirty)でレガシー SSE 命令を実行すると、命令ごとに偽の依存とマージ μop が挿入されて蓄積するからだ。

Go 1.27 のコンパイラはこれを自動挿入しない。 golang/go#80835 に「archsimd の intrinsics を使う関数にレガシー SSE 符号化が出力され、AVX-SSE 遷移ペナルティを起こす」として上がっているが、執筆時点で open のままなので、自分で呼ぶ必要がある。置く位置は「最後のベクトル演算のあと、最初のスカラ float 演算の前」だ。

SIMD まわりのコード生成品質については他にも open の報告があって、この領域がまだ固まっていないことが分かる。

  • #79984(open): simd の演算がループ不変式の巻き上げ(hoisting)の対象として扱われていない
  • #78138(open): VPSRLW の定数畳み込み規則が無く、NotEqual のコード生成も最適でない

4. int8 の内積なのにアキュムレータが Int32x8

ここが一番引っかかった。「int8 に削ったのに 32 bit で計算するなら、削った意味がないのでは」と思ったのだ。

まず型の流れを追うと、選択の余地がないことが分かる。LoadInt8x16Int8x16(16 * 8 bit = 128 bit)を読み、.ExtendToInt16()Int16x16(16 * 16 bit = 256 bit)に広げ、.DotProductPairs()Int32x8(8 * 32 bit = 256 bit)を返す。戻り型が Int32x8 なので、アキュムレータはそれ以外にできない(型が合わずコンパイルできない)。メソッド名の Pairs は「隣り合う 2 レーンの積を足して 1 つにまとめる」ことを表していて、だからレーン数が 16 から 8 に半減する。

そして幅ごとに「積が何個入るか」を見ると、なぜ int32 なのかがはっきりする:

範囲積(最大 16129)が何個入るか
int8-128 〜 1270 個(1 個すら入らない)
int16±327672 個(ペア和 32258 でぎりぎり)
int32±約 21.5 億約 133,000 個

384 次元でも最大 384 * 16129 = 6,193,536 で int32 の 0.3 % しか使わない。int64 にしない理由は、256 bit に入るレーン数が半分になってスループットが落ちるだけだからだ。

そして「削った意味がない」という懸念は的外れだった。算術強度を決めるのは DRAM から運ぶバイト数だけである。メモリ上の Codes8 は 10 万件 * 384 byte = 38.4 MB で転送に効くが、アキュムレータはレジスタの中に 2 本しか存在しない(合計 64 byte)。これを 32 bit にしても DRAM から運ぶ量は 1 バイトも増えない。「保存・転送は狭く、計算は広く」が量子化カーネルの定石なのだと理解した。

5. int8 から int16 を経由する理由

上の型の流れで ExtendToInt16 を挟むのが最初は冗長に見えた。int8 から直接 int32 にできないのか、と。

調べると、AVX2 には int8 の積和命令が存在しないInt8x16 のメソッドは Mul(int8 → int8、溢れる)と MulSign だけで、DotProduct* が 1 つもない。積和を持っているのは int16 のほうなので、int16 にするのは「積和命令に到達するため」だった。直接 int32 にする Int8x16.ExtendToInt32()存在するが、戻り型が Int32x16 で 512 bit になる。AVX-512 が必要で、AVX2 機では使えない。

ここに気づくと納得が早い。256 bit / 16 要素 = 16 bit/要素 なので、int8(128 bit)では積が入らず、int32(512 bit)では収まらない。「1 命令で 16 要素を処理する」と決めた時点で、中間表現は int16 以外にありえない。

なお arm64 は経路が逆だった。Int8x16.MulWidenLo(SMULL、int8 * int8 → int16 の幅拡張つき掛け算)がある代わりにペア積和がないので、「掛けながら広げて、そのあともう一度広げる」3 段構成になる。amd64 側に MulWiden* は 1 つもない。

6. Code8 が何なのか分からなかった

Index が同じ 10 万件を 3 つの表現で同時に持っていることに気づいていなかった:

func (ix *Index) Vec(id int) []float32 { return ix.Data[id*ix.Dim : (id+1)*ix.Dim] }
func (ix *Index) Code8(id int) []int8  { return ix.Codes8[id*ix.Dim : (id+1)*ix.Dim] }
func (ix *Index) Code(id int) []uint64 { return ix.Codes[id*ix.Words : (id+1)*ix.Words] }

この 3 つは完全に対応する兄弟で、Code8 は「id 番目の DB ベクトルの int8 版」を切り出すだけである。Codes8[][]int8 ではなく全ベクトルを 1 本につなげた平らな配列なので、切り出しにコピーもアロケートも発生しない。フラットにしている理由は、DRAM から連続して読めるようにするためだ。スライスのスライスだと各ベクトルがヒープ上に散らばり、順次読みにならないので帯域が出ない。

そして Vec ではなく Code8 と内積を取るのが Stage 2 の目的そのものだった。Vec(id) は 1 件 1536 byte、Code8(id) は 384 byte で、float32 の原本を読まずに済ませるために別表現を持っている。ちなみに Add のたびに Codes8Scales は破棄される。Data が増えたのに int8 表現が古いままだと id と実体がずれた結果を返すからだ。Codes(1 bit)のほうは Add の中でその場で 1 件追記されるので破棄されない。

結果

最終的に到達した点を並べる。すべて 4 コア Codespace(AMD EPYC 7763)の自分の実測である:

Stage表現1 クエリ転送量AI達成Recall@10
0 スカラfloat3234.878 ms153.6 MB0.52.202 GFLOP/sexact
1a SIMD(acc 1 本)float3210.175 ms153.6 MB0.57.548 GFLOP/sexact
1b SIMD(acc 2 本)float329.826 ms153.6 MB0.57.816 GFLOP/sexact
2 int8 量子化int83.965 ms38.40 MB2.019.37 Gop/s0.948
3 バイナリ量子化1 bit0.899 ms4.80 MB0.180
4 1 bit + rerank1 bit + float320.969 ms4.80 MB + 154 KB0.868
  • スカラ全探索から 38.8 倍(同一実行内の 37.579 ms → 0.969 ms)、精度は Recall@10 = 0.868
  • 当たった上限は順に 依存連鎖 → メモリ帯域 → int8 内積の命令スループット → 精度 の 4 つ
  • SIMD が効いたのは Stage 0→1(速さの主因)、Stage 1→2(int8 内積 11.85 倍)、Stage 3→4(精度の回復)
  • 未計測のまま残したもの: 自分の Codespace の演算ピークとメモリ帯域(make roofline-ceiling)、内積単体のベンチ(BenchmarkDot*)、クエリのバッチ化と goroutine 並列(教材の付録)

学び

SIMD もベクトル検索も関係なく持ち帰れそうだと思ったことを挙げる。

1. 天井は自分で測る。理論値で判断しない。

教材が屋根に使っていたのは理論ピーク 110 GFLOP/s ではなく実測の 25.6 だった。理論値で線を引くと、どう頑張っても届かない高さになってモデルが予測の役に立たなくなる。しかも実測が 1/4 に留まる原因は register spill というコンパイラ側の事情で、ハードの限界ではなかった。測るべきは「この言語・この処理系で到達できる上限」で、CPU のスペック表ではない。

2. 上限を越える手は 2 種類しかない。

実装効率を上げる(SIMD、アキュムレータ分割)か、算術強度を上げる(量子化、バッチ化)かだ。点が屋根のどこにいるかで、打つべき手が決まる。斜線に張り付いていたら実装をいくら速くしても無駄で、逆に屋根から遠ければ実装に無駄がある。

3. 「SIMD が効かない」ではなく「いまの算術強度では効かない」。

Stage 1b で 3.5 % しか伸びず一度は行き止まったが、int8 化して算術強度を上げた先で SIMD が再び効いた(内積単体で 11.85 倍)。手法の是非ではなく、適用する場所の問題だった。

4. 速度と精度は別の軸として管理する。

Stage 2 で量子化を入れた瞬間に「速くなったが間違った文書を返す」が起こりうる状態になった。速度だけ見ていると気づけない。指標を 2 本持って初めて Stage 3 の「46 倍だが Recall 0.18」が異常だと判断できる。

5. 「安い一次フィルタ + 高価な二次判定」は領域を越えて現れる構造。

1bit で粗く 100 件に絞り、float32 で正確に採点し直す 2 段構成は、自分が普段扱っている認証のリスクスコアリングと同じ形だった。安いシグナルで大量に粗く弾いて、残ったものだけ重い検査に回す。コストと精度のトレードオフがある領域に共通して現れるのだと思う。

6. 主張の粒度と計測の粒度を合わせる。

これは反省。自分は make roofline-ceiling を自分の Codespace で回さなかったので、この記事の屋根は教材の値の引き写しになっている。Stage 1a → 1b の効果も、全探索しか測っておらず内積単体を測っていない。そして教材の「Recall 0.868」に引っかかったのも同じ構図で、速度と精度を違う N で測っていたことが原因だった。 測る対象と主張する対象がずれていないかは、毎回確認する価値がある。

おわりに

また、こちらのワークショップは、「Go Far, Go Together」の通り、Go 単体から派生してベクトル検索の分野の素養についても調べるきっかけとなった。業務等ですぐに利用する領域ではないが、素養として学ぶきっかけになったので非常に良かったと思う。

本ワークショップのオーガナイザである pon さん、ならびに関係者のみなさま、本当に良い機会を提供してくださり、ありがとうございました!

参考

教材とセッション:

ルーフラインモデル:

Go の SIMD:

ベクトル検索と量子化:

Posted on 12, 2026