累積和

累積和（るいせきわ、英: prefix sum, cumulative sum, scan）とは、計算機科学において、数列 ${\displaystyle x_{0},x_{1},x_{2},\dots }$ に対して、その先頭部分の総和を求めることによって得られる数列 ${\displaystyle y_{0},y_{1},y_{2},\dots }$ の事。

${\displaystyle y_{0}=x_{0}}$

${\displaystyle y_{1}=x_{0}+x_{1}}$

${\displaystyle y_{2}=x_{0}+x_{1}+x_{2}}$

${\displaystyle \dots }$

例えば、自然数の累積和は三角数になる。

自然数	1	2	3	4	5	6	...
自然数の累積和	1	3	6	10	15	21	...

累積和は、単に足し算だけで無く、二項演算子 ${\displaystyle \oplus }$ に一般化することが可能であり、そのため幅広い応用が可能である。これにより、関数型プログラミング言語では、scanと呼ばれる基本的な処理となっている。なお、途中の計算過程を記録する必要が無く、最終結果だけが必要な場合はfoldと呼ばれる。^[1][2]

いくつかのアルゴリズムにおいて有用な基本処理であり、カウントソートなどのアルゴリズムで利用されている。二項演算子 ${\displaystyle \oplus }$ に対して引き算 ${\displaystyle \ominus }$ が存在する場合、事前に累積和を求めておくと、m番目からn番目までの総和が「n番目の累積和 ${\displaystyle \ominus }$ (m-1)番目の累積和」により、高速に求めることができる。2次元の場合はこれをsummed-area tableと呼ぶ。^[3][4]

累積和は、逐次計算においては、単に前の結果と計算するだけで簡単に求まるが、並列計算の分野でも広く研究されており、foldやscanの二項演算子が ${\displaystyle (a\oplus b)\oplus c=a\oplus (b\oplus c)}$ という結合法則を満たすと並列化することが可能であり、並列アルゴリズムの有用な基本処理になっている。^[5][6][7]

高階関数のscan

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

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

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

inclusiveとexclusive

プログラミング言語やライブラリにおける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の実装の一覧。

プログラミング言語	inclusive scan	exclusive scan
APL	+\
C++	std::inclusive_scan	std::exclusive_scan
Clojure	reductions init無し	reductions init有り
CUDA	thrust::inclusive_scan cub::DeviceScan::InclusiveScan	thrust::exclusive_scan cub::DeviceScan::ExclusiveScan
F#		scan
Haskell	scanl1	scanl
HPF	sum_prefix	sum_prefix(..., exclusive=.true.)
Java	Gatherers.scan Arrays.parallelPrefix
Kotlin		scan
MPI	MPI_Scan	MPI_Exscan
Python	itertools.accumulate 引数initialがNoneの場合	itertools.accumulate 引数initialがNoneでない場合
Rust	scan
Scala		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)}$ とする。

アルゴリズム1：短いスパン、高い並列性

並列性の高い16入力の並列累積和の回路表現

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]

アルゴリズム2：ワーク効率が高い方法

ワーク効率が高い16入力の並列累積和の回路表現

ワーク効率の良い並列累積和は、以下の手順で計算できる。^[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つのプロセッサが複数のインデックスを処理するように割り当てることで実現される。

より一般的なfoldやscanに適用する方法

本項では、要素 ${\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;
}