01行列積はなぜGPU向きか
行列積 C = A × B の定義はシンプルだ。出力の各要素は、Aの1行とBの1列の内積で決まる。
C[i][j] = Σ_k A[i][k] × B[k][j]
ここで重要なのは、出力 C[i][j] が互いに独立に計算できる点だ。C[0][0] を計算するのに C[0][1] の結果は要らない。つまり、出力の全要素を同時並列に計算できる。M×N個の出力要素があれば、理屈の上ではM×N個の計算を一斉に走らせられる——これがGPUの大量の演算ユニットと噛み合う。
GPUは「1スレッド＝1出力要素」という割り当てで、数千〜数万のスレッドを同時に動かす。CPUが要素を1つずつ順に計算する間に、GPUは全要素をまとめて片付ける。行列が大きいほど、この並列性の利得が大きくなる。
02naive実装 — 1スレッド1要素
まず素朴な実装から。各スレッドが自分の担当する C[row][col] を1つ計算する。WGSL（WebGPU Shading Language）のコンピュートシェーダーで書くとこうなる。
WGSL — naive matmul
// 各スレッドが C[row][col] を1つ計算
@compute @workgroup_size(16, 16)
fn main(@builtin(global_invocation_id) gid : vec3<u32>) {
let col = gid.x; let row = gid.y;
if (row >= D.M || col >= D.N) { return; } // 範囲外は何もしない
var sum : f32 = 0.0;
for (var k : u32 = 0u; k < D.K; k = k + 1u) {
sum = sum + A[row * D.K + k] * B[k * D.N + col]; // 内積
}
C[row * D.N + col] = sum;
}
global_invocation_id は、このスレッドがグリッド全体で何番目かを表す座標だ。gid.x を列、gid.y を行に対応させ、各スレッドが自分の (row, col) を担当する。@workgroup_size(16, 16) は、スレッドを16×16＝256個ずつの塊（ワークグループ）にまとめる指定だ。
行列はGPUには1次元配列として渡すので、2次元の [row][col] は row * 幅 + col で添字を計算する（row-major）。境界チェック（row >= D.M）は、行列サイズがワークグループの倍数でないときに、はみ出したスレッドが配列外を触らないための保険だ。
naiveの弱点
この実装は正しく動くが、同じデータを何度もグローバルメモリから読む。たとえばCの同じ行の要素は、すべてAの同じ行を読む。隣のスレッドが直前に読んだ値を、また読み直す。GPUのグローバルメモリは広いが遅い。この「読み直し」の無駄が、大きな行列で効いてくる。
03タイル化 — 共有メモリで再利用する
そこでタイル化（tiling）の出番だ。アイデアはこうだ。行列を小さなタイル（ブロック）に区切り、ワークグループ内で共有できる高速なメモリ（共有メモリ）にタイルを一度だけ読み込み、その中で何度も再利用する。
タイル化：AとBのタイルを共有メモリに読み込み、その中で内積を計算。K方向にタイルをずらしながら累積する。
WGSLでは var<workgroup> が共有メモリにあたる。ワークグループ内の全スレッドがアクセスでき、グローバルメモリよりずっと高速だ。コードの核心はこうなる。
WGSL — tiled matmul（TILE = タイルの一辺）
// 共有メモリ：A・Bのタイルを格納
var<workgroup> Asub : array<array<f32, TILE>, TILE>;
var<workgroup> Bsub : array<array<f32, TILE>, TILE>;
fn main(...) {
var sum : f32 = 0.0;
let numTiles = (D.K + TILE - 1u) / TILE;
for (
var t = 0u; t < numTiles; t++) {
// ① 各スレッドがタイルの1要素を共有メモリへロード
Asub[lid.y][lid.x] = A[ ... ];
Bsub[lid.y][lid.x] = B[ ... ];
workgroupBarrier();
// ② 全スレッドのロード完了を待つ
// ③ 共有メモリ上のタイルで内積（高速）
for (
var k = 0u; k < TILE; k++) {
sum = sum + Asub[lid.y][k] * Bsub[k][lid.x];
}
workgroupBarrier();
// ④ 次のタイルロード前に全員待つ
}
C[row * D.N + col] = sum;
}
ポイントは workgroupBarrier() だ。これはワークグループ内の全スレッドの同期点で、「全員がここに到達するまで待つ」。①で各スレッドが共有メモリにタイルを書き込むが、②のバリアで「全員が書き終わった」ことを保証してから、③で読む。これがないと、まだ書かれていない領域を読んでしまう。④のバリアは、次のタイルで共有メモリを上書きする前に、全員が③を読み終えるのを待つ。
なぜ速くなるか
naiveでは、Cの1要素の計算でAとBをK回ずつグローバルメモリから読む。タイル化では、タイルを一度だけ共有メモリに載せ、ワークグループ内の全スレッドがそれを共有して使い回す。グローバルメモリへのアクセス回数が、おおよそタイルの一辺の分だけ減る。遅いメモリを触る回数を減らすこと——これがGPU最適化の王道だ。
04なぜタイルサイズ32なのか
ところが実測すると、タイルサイズ次第で結果が一変した。ベンチマークでは、タイル16では naive に負けることすらあり、タイル32で一気に最速になった（詳細は結果編）。同じタイル化なのに、なぜここまで違うのか。
| タイルサイズ | ワークグループのスレッド数 | 傾向 |
|---|
| 8×8 | 64スレッド | 遅い（naive以下のことも） |
| 16×16 | 256スレッド | 中途半端 |
| 32×32 | 1024スレッド | 最速（両機種で） |
理由は共有メモリの再利用率にある。タイルが小さいと、共有メモリに載せた値を使い回す回数が少なく、ロードと同期（バリア）のコストばかりがかさむ。タイル8では、その手間が naive の素朴さにすら負ける。
タイルを大きくすると、一度ロードしたタイル内で行う計算（内積）の回数が増え、ロード1回あたりの「元が取れる」度合いが上がる。32×32（1024スレッド/ワークグループ）は、今回の2機種のモバイルGPUがちょうど活かせる上限近くで、ここで再利用率と並列度のバランスが最適になった。
教訓
「並列化すれば速い」「タイル化すれば速い」は素朴すぎる。最適なパラメータはGPUのアーキテクチャ（共有メモリ量、ワークグループ上限、メモリ帯域）に依存し、机上では決まらない。実測して初めて分かる——これがGPUプログラミングの現実だ。本エンジンは実測結果を受けて、デフォルトのタイルサイズを32に設定している。
05転置を避けた共分散行列
ここから金融計算だ。ポートフォリオ最適化の入力となる共分散行列は、資産のリターン時系列 X（n期間×m資産）から作る。
Σ = (1/(n-1)) X'ᵀ X' （X' は列ごと中心化）
XᵀX は転置を含む積だ。素直にやるならXを転置してから行列積にかけるが、それは余分なメモリコピーを生む。だが、よく見ると (XᵀX)[i][j] は「Xの列iと列jの内積」だ。つまり転置せずに、列同士を直接かけ合わせればいい。
WGSL — 共分散行列（転置を省略）
fn main(@builtin(global_invocation_id) gid : vec3<u32>) {
let i = gid.y; // 資産 i
let j = gid.x; // 資産 j
if (i >= D.m || j >= D.m) { return; }
let mi = MEAN[i]; let mj = MEAN[j]; // 列平均（事前計算）
var sum : f32 = 0.0;
for (var t = 0u; t < D.n; t++) { // 期間方向に走査
let xi = X[t * D.m + i] - mi; // 中心化した列i
let xj = X[t * D.m + j] - mj; // 中心化した列j
sum = sum + xi * xj;
}
S[i * D.m + j] = sum / f32(D.n - D.ddof); // 不偏分散は n-1
}
各スレッドが共分散行列の1要素 Σ[i][j] を担当し、期間方向（t）に走査して、中心化した列iと列jの積を足し込む。中心化（平均を引く操作）も式に組み込んでいるので、Xを書き換える必要がない。列平均だけは事前にCPUで計算してGPUに渡している（資産数mは小さく、ここはGPU化の旨味が薄いため）。
06ポートフォリオ並列評価 — 行列演算×独立並列
共分散行列Σと期待リターンμが揃えば、配分w（各資産の比率）に対するポートフォリオのリスクとリターンが計算できる。
リスク σ = √(wᵀΣw) 期待リターン μp = wᵀμ
効率的フロンティアを描くため、数万通りの配分wを一気に評価する。ここで設計が面白いことになる。各スレッドが1つの配分wを担当する——これは行列積（1スレッド1要素）ではなく、独立並列（各試行が独立、モンテカルロと同じ構造）だ。スレッド内部では wᵀΣw という行列・ベクトル演算を行う。行列演算と独立並列が合流する。
WGSL — ポートフォリオ並列評価（各スレッド＝1配分）
fn main(@builtin(global_invocation_id) gid : vec3<u32>) {
let p = gid.x; // p番目の配分を担当
if (p >= D.N) { return; }
let base = p * D.m; // 配分pの先頭位置
// 期待リターン wᵀμ
var mp : f32 = 0.0;
for (var i=0u; i<D.m; i++){ mp = mp + W[base+i]*MU[i]; }
// 分散 wᵀΣw（二重ループ）
var varp : f32 = 0.0;
for (var i=0u; i<D.m; i++){
var row : f32 = 0.0;
for (var j=0u; j<D.m; j++){ row = row + W[base+j]*SIG[i*D.m+j]; }
varp = varp + W[base+i]*row;
}
let sd = sqrt(max(0.0, varp));
RISK[p] = sd; RET[p] = mp;
if (sd > 1e-9) { SHARPE[p] = (mp - D.rf)/sd; } else { SHARPE[p] = 0.0; }
}
1スレッドが wᵀΣw（m×mの二重ループ）と wᵀμ を計算し、リスク・リターン・シャープレシオを出力する。N個の配分を別々のスレッドが同時に処理するので、5万通りの配分でも一瞬だ。各スレッドの計算量は資産数mの二乗 O(m²) なので、資産数が増えるほどGPUの並列性が効く（結果編で実証）。
設計の妙
逆行列を使う解析解を避け、力任せに数万配分を並列評価することで、効率的フロンティアを散布図の外縁として描く。逆行列はGPUの並列化に向かない（逐次性が強い）が、この「総当たり評価」はGPUの独立並列と完全に噛み合う。アルゴリズムをGPUの得意な形に寄せる、という発想だ。
07エンジンへの組み込み
これらのシェーダーは、共通エンジン gwaw-compute.js にメソッドとして組み込んでいる。利用側は数行で呼べる。
JavaScript — 呼び出し例
const gc = new GwawCompute();
await gc.detect(); // WebGPU/CPU を自動判定
// 行列積（デフォルトで tiled-32 を使用）
const { C } = await gc.matmul(A, B, M, K, N);
// 共分散行列
const { Sigma } = await gc.covariance(X, n, m);
// N通りの配分を一括評価
const { risk, ret, sharpe } = await gc.evalPortfolios(Sigma, mu, W, N, m);
内部では、WebGPUが使えれば各シェーダーを実行し、使えなければ同じ計算のCPU版に自動フォールバックする。どちらでも結果が一致することは、開発中つねにCPU版を基準に検証した。GPUの計算結果は「速いが正しいか不安」になりがちなので、CPUという答え合わせの基準を常に持っておくのは、GPUプログラミングで欠かせない習慣だと思う。
以上が、ベンチマークであの数字を出した実装の中身だ。要点をまとめると——行列積は1スレッド1要素の独立性でGPUに向く。タイル化は共有メモリの再利用で速くなるが、サイズ選びが命。共分散は転置を省いて列の内積で計算。ポートフォリオ評価は逆行列を避け、数万配分の独立並列に持ち込む。いずれも「GPUの得意な形に計算を寄せる」という一貫した発想でできている。
▶ このページで実際に試す — 行列積ベンチ
あなたのデバイスで CPU / WebGPU naive / WebGPU tiled(32) を実行します。記事の数字を追体験してください。
WebGPU
CPU
※ 計測は単発のため、初回はシェーダーのコンパイル時間を含み遅めに出ることがあります。2回目以降が実力値に近くなります。
