累積和

関数型プログラミング言語の観点では、累積和は加算に限らず任意の二項演算子へと一般化できる。この一般化によって得られる高階関数はscanと呼ばれ、foldと密接に関連している。scanとfoldはどちらも与えられた二項演算子を同じ数列に適用するが、両者には違いがある。scanは演算の途中結果を含む全ての中間値の列を返すのに対し、foldは最終結果のみを返す。

例えば、階乗数列は、自然数列に対して加算の代わりに乗算を用いたscanを行うことで生成できる。

入力	1	2	3	4	5	6	...
累積の積	1	2	6	24	120	720	...

プログラミング言語やライブラリにおけるscanの実装には、inclusiveとexclusiveの2種類が存在する。

入力	1	2	3	4	5	6	...
inclusive	1	3	6	10	15	21	...
exclusive	0	1	3	6	10	15	...

inclusive scanでは、出力 ${\displaystyle y_{i}}$ を計算する際に入力 ${\displaystyle x_{i}}$ を含めるのに対し（すなわち、${\textstyle y_{i}=\bigoplus _{j=0}^{i}x_{j}}$）、exclusive scanでは ${\displaystyle x_{i}}$ を含めない（すなわち、${\textstyle y_{i}=\bigoplus _{j=0}^{i-1}x_{j}}$）。後者の場合、${\displaystyle y_{0}}$ を未定義のままとするか、scanの初期値として特別な値 ${\displaystyle x_{-1}}$ を受け取る実装が一般的である。

inclusive scanとexclusive scanは相互に変換可能である。inclusive scanをexclusive scanに変換するには、scanで得られた配列を右に1つシフトし、左端に単位元（identity value）を挿入すればよい。逆に、exclusive scanをinclusive scanに変換するには、scanで得られた配列を左に1つシフトし、右端に「scanの最後の要素と入力配列の最後の要素の和」を挿入すればよい。^[8]

以下はプログラミング言語でのscanの実装の一覧。

二項演算子 ${\displaystyle \oplus }$ が ${\displaystyle (a\oplus b)\oplus c=a\oplus (b\oplus c)}$ という結合法則を満たす場合は並列化が可能である。

この並列化手法はGPUでも利用可能である。NVIDIAのGPU Gems 3のFigure 39-7によると、2007年当時、要素数nが1,000の場合はCPUの方が高速だが、要素数nが1,000,000ある場合はGPU(CUDA)の方が高速である。^[8]

並列計算において、累積和を求めるためのアルゴリズムは多数存在する。NVIDIAでは2013年にNVIDIA CUDAのCUB 1.0.1で実装され^[9]、2016年に論文が書かれたdecoupled look-back法を使用していて、二項演算子が単純な足し算の場合は、要素数nが十分大きければ、メモリ帯域がボトルネックとなっていて、要素数nのメモリコピーと同じ速度で動く^[10]。この手法はIntel GPU用のIntel oneAPI DPC++でも使用されている^[11]。

以下、2007年に出版されたNVIDIAのGPU Gems 3のChapter 39で紹介されているアルゴリズムを説明する^[8]。これらはより効率が良いものが発見されているので現在はNVIDIAは使用していない。1つ目のアルゴリズムは、より短いスパン（計算の依存関係の深さ）を持ち、高い並列性を実現できるが、ワーク効率（計算量の総和）が低い。2つ目のアルゴリズムは、ワーク効率が高いものの、スパンが2倍となり並列性が低下する。以下に、それぞれのアルゴリズムについて説明する。二項演算子 ${\displaystyle \oplus }$ の計算量は ${\displaystyle O(1)}$ とする。

HillisとSteeleは、以下の並列累積和アルゴリズムを提案している。^[12]

for i <- 0 to floor(log2(n)) do
    for j <- 0 to n - 1 do in parallel
        if j < 2i then
            xi+1
j <- xi
j
        else
            xi+1
j <- xi
j ${\displaystyle \oplus }$ xi
j - 2i

ここで、${\displaystyle x_{j}^{i}}$ は、ステップ ${\displaystyle i}$ における配列 ${\displaystyle x}$ の ${\displaystyle j}$ 番目の要素の値を表す。

このアルゴリズムを単一プロセッサで実行した場合、計算量は ${\displaystyle O(n\log n)}$ となる。しかし、少なくとも ${\displaystyle n}$ 個のプロセッサを用いて内側のループを並列実行できる環境であれば、外側のループの繰り返し回数に等しい ${\displaystyle O(\log n)}$ 時間で計算を完了できる。^[13]

ワーク効率の良い並列累積和は、以下の手順で計算できる。^[5][14][15]

1. 隣接する要素の和を計算する（ペアの先頭要素のインデックスが偶数であるものを対象とする）。
   • 例：${\displaystyle z_{0}=x_{0}\oplus x_{1},z_{1}=x_{2}\oplus x_{3},\dots }$
2. ステップ1で得た数列 ${\displaystyle z_{0},z_{1},z_{2},\dots }$ に対して、再帰的に累積和を計算する。
   • 結果として、新たな数列 ${\displaystyle w_{0},w_{1},w_{2},\dots }$ を得る。
3. 最終的な累積和 ${\displaystyle y_{0},y_{1},y_{2},\dots }$ を、中間数列の値を用いて求める。
   • 具体的には、各 ${\displaystyle y_{i}}$ は、これまでに計算された中間数列の要素の和で表される。
   • 例：
     • ${\displaystyle y_{0}=x_{0}}$
     • ${\displaystyle y_{1}=z_{0}}$
     • ${\displaystyle y_{2}=z_{0}\oplus x_{2}}$
     • ${\displaystyle y_{3}=w_{1}}$
   • 最初の値 ${\displaystyle y_{0}}$ を決めた後、それ以降の各 ${\displaystyle y_{i}}$ は、数列 ${\displaystyle w}$ の半分の位置にある値をコピーするか、直前の値に数列 ${\displaystyle x}$ の一部の値を加えることで求める。

入力数列の長さを ${\displaystyle n}$ とすると、このアルゴリズムの再帰の深さは ${\displaystyle O(\log n)}$ となる。したがって、並列実行時の計算時間も ${\displaystyle O(\log n)}$ に抑えられる。

このアルゴリズムの総ステップ数は ${\displaystyle O(n)}$ であり、${\displaystyle O(n/\log n)}$ 個のプロセッサを持つ並列ランダムアクセス機械（PRAM）上で、非対称的な遅延なしに実装可能である。これは、プロセッサの数よりも要素数が多い段階では、1つのプロセッサが複数のインデックスを処理するように割り当てることで実現される。

本項では、要素 ${\displaystyle \oplus }$ 要素（Haskellの表記では型a -> a -> a^[16]）として書いているが、関数型プログラミング言語では、

• 初期状態（型b^[17]）
• f(状態, 要素) = 次の状態（型b -> a -> b^[17]）

と考えることにより、foldやscanで逐次処理を表現している。foldは最終状態（型b）で、scanは状態遷移列（型[b]）である。

ここで、結合法則を満たす

• 状態 ${\displaystyle \oplus }$ 状態 = 結合した状態（型b -> b -> b）

が存在すると、

1. まず、f(状態, 要素) は f(初期状態, 要素) でしか計算しないので、要素から f(初期状態, 要素) へのmap変換を行う。
2. その上で、状態 ${\displaystyle \oplus }$ 状態 -> 結合状態を使用すると、上記の並列アルゴリズムが適用できる。

これにより、逐次処理のfoldやscanで表現していたものに対して、並列アルゴリズムが使えるようになる。

具体例として、指数移動平均 ${\displaystyle s_{t}=(1-\alpha )s_{t-1}+\alpha x_{t},\ s_{0}=0}$ を考える。要素は ${\displaystyle x_{t}}$ だが、状態を (指数移動平均の最後の値, 指数移動平均の長さ) のタプルとすると下記変換式により並列計算ができる。

• 初期状態 = ${\displaystyle (0,0)}$
• f(初期状態 ${\displaystyle (0,0)}$, 要素 ${\displaystyle x_{t}}$) = ${\displaystyle (\alpha x_{t},1)}$
• 状態 ${\displaystyle (s_{t},l_{t})}$ ${\displaystyle \oplus }$ 状態 ${\displaystyle (s_{u},l_{u})}$ = 結合状態 ${\displaystyle ((1-\alpha )^{l_{u}}s_{t}+s_{u},l_{t}+l_{u})}$

これをNVIDIA CUDAのThrust(C++)で実装する。要素数nが十分大きく、二項演算子の計算量が小さい場合はinclusive_scanはメモリコピーの速度で動くが^[10]、要素から状態への並列map、状態の並列scan、状態から要素への並列mapが必要で、単純に実装すると3回分のメモリコピーが発生するが、Thrustではtransform_iteratorを使用すると1回分のメモリコピーにまとめることができるので、それを使用している。

#include <thrust/device_vector.h>
#include <thrust/iterator/transform_output_iterator.h>

struct State { float v; int len; };

struct toState { // 要素→状態
    const float alpha;
    __device__ State operator()(const float x) const { return State{alpha * x, 1}; }
};
struct toElement { // 状態→要素
    __device__ float operator()(const State &s) const { return s.v; }
};
struct plusStates { // 状態＋状態
    const float alpha;
    __device__ State operator()(const State &a, const State &b) const {
        return State{std::pow(1 - alpha, (float) b.len) * a.v + b.v, a.len + b.len};
    }
};

int main() {
    const float alpha = 0.2f;
    thrust::device_vector<float> input = {3, 1, 4, 1, 5, 9, 2, 6};
    thrust::device_vector<float> output(input.size());

    // 指数移動平均の計算
    thrust::inclusive_scan(
        thrust::make_transform_iterator(input.begin(), toState{alpha}),
        thrust::make_transform_iterator(input.end(), toState{alpha}),
        thrust::make_transform_output_iterator(output.begin(), toElement{}),
        plusStates{alpha}
    );

    // 計算結果の出力
    for (float x: output) {
        std::cout << x << " ";
    }
    std::cout << std::endl;

    return 0;
}