Pandas DataFrame (Series) における sum の違い

お前は何を言っているんだというタイトルだと思いますが、 DataFrame (Series) に対する sum に対するこの結果を見てみましょう。

import math
import pandas as pd

df = pd.DataFrame({"a": [1e-7] * int(1e7)})

df["a"].sum()
1.0000000000000173

sum(df["a"])
0.99999999975017

math.fsum(df["a"])
1.0

いかがでしたでしょうか。みたいな薄っぺらい商材ブログの文体とは異なり、ストイックな追求をここではしていく。

前提知識

そもそも単純だと思われるsum関数において、いくつかの加算アルゴリズムが存在しているという事実に注意が必要だ（と偉そうに書いていてこの記事を書くために調べて初めてきちんと理解した）。具体的には Kahan summation, pairwise summation が今回の題材である。

いきなり結論

• pandas Series.sum() は内部で numpy.sum() を呼び、それは pairwise summation を使う。
  • 一方で GroupBy sum() や Rolling sum() method は、kahan summation を使うため、上記と結果が一致しないことがある
• build-in sum関数(CPython)は単純ループ加算（呼び名不明）を使う
• math.fsum は kahan summation を使う

調査環境

>>> import sys
>>> sys.version
'3.9.7 (default, Dec  5 2021, 20:50:26) \n[GCC 9.3.0]'
>>> import numpy as np
>>> import pandas as pd
>>> np.__version__
'1.21.4'
>>> pd.__version__
'1.3.5'

個々の呼び出し元のソースコードの確認

pandas Series の sum method

（注意: 長文です）

breakpoint() の step-in で雑にコードをたどると
pandas/core/generic.py(10708)sum()
pandas/core/generic.py(10437)sum()
pandas/core/generic.py(104390)_min_count_stat_function()
pandas/core/series.py(4359)_reduce()
pandas/core/nanops.py(83)_f()
pandas/core/nanops.py(395)new_func()
pandas/core/nanops.py(546)nansum() -> (591) ndarray.sum()
numpy/core/_methods.py(46)_sum()

# numpy/core/_methods.py(48)
return umr_sum(a, axis, dtype, out, keepdims, initial, where)

ここで debug 上のstepで _methods.py(48)_sum()->1.0000000000000173 と値が返ってきてしまった。

さて、 umr_sum() とは事実上 numpy.core.umath.add.reduce を呼び出していることが同ファイルから伺える。

# numpy/core/_methods.py (10, 20)
from numpy.core import umath as um
...
umr_sum = um.add.reduce

numpy/core/umath.py を見ると __all__ で "add" が定義されているので、 from . import _multiarray_umath を遡ろうとするが、 numpy/core/_multiarray_umath というファイルは存在せず numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so のバイナリファイルがあった（調査環境は Ubuntu 20.04 on WSL2）。

一方で numpy.ufunc の method だということもわかる。

Pdb) p umr_sum
<built-in method reduce of numpy.ufunc object at 0x7fe743c56a90>

さて、行き詰ってしまったので numpy.core.umath.add を探してみると generate_umath.py(293) 付近になにやら怪しいものを発見

$ find -type f -exec grep -nH 'numpy.core.umath.add' {} \;
./numpy/core/code_generators/generate_umath.py:293:          docstrings.get('numpy.core.umath.add'),

C言語のコード生成で作っている

# This dictionary describes all the ufunc implementations, generating
# all the function names and their corresponding ufunc signatures.  TD is
# an object which expands a list of character codes into an array of
# TypeDescriptions.
defdict = {
'add':
    Ufunc(2, 1, Zero,
          docstrings.get('numpy.core.umath.add'),
          'PyUFunc_AdditionTypeResolver',
          TD(notimes_or_obj, simd=[('avx2', ints)], dispatch=[('loops_arithm_fp', 'fdFD')]),
          [TypeDescription('M', FullTypeDescr, 'Mm', 'M'),
           TypeDescription('m', FullTypeDescr, 'mm', 'm'),
           TypeDescription('M', FullTypeDescr, 'mM', 'M'),
          ],
          TD(O, f='PyNumber_Add'),
          ),
# 略
# L1206
if __name__ == "__main__":
    filename = __file__
    code = make_code(defdict, filename)
    with open('__umath_generated.c', 'w') as fid:
        fid.write(code)

make_code() でc言語のコードを作成する DSL 的パートであろう。注目するのは PyNumber_Add と dispatch=[('loops_arithm_fp', 'fdFD')] の 'loops_arithm_fp だ。

loops_arithm_fp をファイル名で調べると numpy/core/src/umath/loops_utils.h.src というファイルが見つかる。そ、そしてついに L78 に @TYPE@_pairwise_sum(char *a, npy_intp n, npy_intp stride) という目的っぽいC言語の関数を見つける！！！（ここまで2.5時間）

---

以下ボツとなった調査プロセス（供養）。探索的調査方法が誰かの参考になるため残しておく。

さて、ここで行き詰まるのも勿体ないので `numpy/core/include/numpy/*.h` あたりを読むと何やら `PyArray_Sum()` なるものを `numpy/core/include/numpy/__multiarray_api.h` にて発見できた。とはいえローカルファイル上ではソースが存在しないため、ここから github のソースコードを探すとにする

numpy/core/src/multiarray/calculation.c#L521 にて PyArray_Sum() が定義され、その中で PyArray_GenericReduceFunction() を呼び出している。それは numpy/core/src/multiarray/number.c#L228 で定義され、その中で

 meth = PyObject_GetAttrString(op, "reduce");
    if (meth && PyCallable_Check(meth)) {
        ret = PyObject_Call(meth, args, kwds);
    }

とある。ここで op = numpy.add と推測されることから、結局は numpy.add.reduce を呼び出しているのであろうか。コードリーディング力が低すぎて循環してしまった。。。（分かる方コメントを頂けると助かります）。

ということで stackoverflow に助けを借りにいくと、今回の調べたいことそのもの の回答を見つけた。ずばり numpy/core/src/umath/loops_utils.h.src#L78 にて @TYPE@_pairwise_sum(...) で定義されている。

---

この @TYPE@_pairwise_sum(...) を読むと、次の特徴がわかる

1. 配列要素数が8未満の場合、単純な足し算ループ
2. 配列要素数が PW_BLOCKSIZE (=128) 以下の場合、　r[0] から r[7] までの8つの配列を用意し、r[0] には i = 0, 8, 16, ..., r[1] には i = 1, 9, 17, ..., [...], r[7] には i = 7, 15, 23, ... が詰め込まれ、最後に res = ((r[0] + r[1]) + (r[2] + r[3])) + ((r[4] + r[5]) + (r[6] + r[7])) の順番で足し合わされる。
3. 配列要素が 128 超の場合、　前半と後半の約半分に（前半部分が8の倍数になるように）分け、再帰的に足し合わせる。

Groupby.Sum

pandas/core/groupby/groupby.py(1838)sum()
pandas/core/groupby/groupby.py(1352)_agg_general()
pandas/core/groupby/generic.py(1056)_cython_agg_general()
pandas/core/internals/managers.py(530)get_numeric_data()
pandas/core/internals/managers.py(539)_combine()
pandas/core/groupby/generic.py(1066)_cython_agg_general()
pandas/core/internals/managers.py(1213)grouped_reduce()
pandas/core/internals/blocks.py(374)apply()
pandas/core/groupby/generic.py(1066)array_func()
pandas/core/groupby/ops.py(991)_cython_operation()
pandas/core/groupby/ops.py(631)cython_operation()
pandas/core/groupby/ops.py(488)_cython_op_ndim_compat()
pandas/core/groupby/ops.py(525)_call_cython_op()
(580) の func の call で結果が返っている

(Pdb) p func
<cyfunction group_add at 0x7f2f9a657ac0>

$ find -type f -exec grep -nH 'group_add' {} \;
# 略
/pandas/_libs/groupby.pyx:482:def group_add

/pandas/_libs/groupby.pyx を確認 すると

For floats, use Kahan summation to reduce floating-point

とあるので、 kahan summation が使われていることがわかる。実装はこの部分を読めば良い（解説略）

Rolling.sum()

pandas/core/window/rolling.py(1670)sum()
pandas/core/window/rolling.py(1221)sum()
pandas/core/window/rolling.py(482)_apply()
pandas/core/window/rolling.py(411)_apply_blockwise()
pandas/core/window/rolling.py(392)_apply_series()
pandas/core/window/rolling.py(514)homogeneous_func()
pandas/core/window/rolling.py(520)calc()
pandas/core/window/rolling.py(531)calc()

(Pdb) p func
<built-in function roll_sum>

$ find -type f -exec grep -nH 'roll_sum' {} \;
# 略
pandas/_libs/window/aggregations.pyx:118:def roll_sum

pandas/_libs/window/aggregations.pyx を確認 すると """ remove a value from the sum calc using Kahan summation """ と記述があることが確認できた。

math.fsum

CPython 上の math.fsum を見てみると（少し拡張された） kahan summation であることがわかる。

// L1474-1485
for (i = j = 0; j < n; j++) {       /* for y in partials */
    y = p[j];
    if (fabs(x) < fabs(y)) {
        t = x; x = y; y = t;
    }
    hi = x + y;
    yr = hi - x;
    lo = y - yr;
    if (lo != 0.0)
        p[i++] = lo;
    x = hi;
}

sum 関数

bltinmodule.c の builtin_sum_impl(...) の中の L2564 あたりの for 文。単純なループでの順番での足し上げ。

    for(;;) {
        item = PyIter_Next(iter);
        if (item == NULL) {
            /* error, or end-of-sequence */
            if (PyErr_Occurred()) {
                Py_DECREF(result);
                result = NULL;
            }
            break;
        }
        temp = PyNumber_Add(result, item);
        Py_DECREF(result);
        Py_DECREF(item);
        result = temp;
        if (result == NULL)
            break;
    }
    Py_DECREF(iter);
    return result;

PyNumber_Add は、内部で BINARY_OP1(v, w, NB_SLOT(nb_add), "+") を call する。以降の解説は Python処理系入門 〜1 + 1 で学ぶ処理系解読の基礎〜を参考にすると良い。

再現

pairwise summation by numpy

python, C言語でも結果が再現できない。原因がわかる方、コメントもらえると大変うれしいです。

python code (再現できず)

def my_pairwise(lst):
    PW_BLOCKSIZE = 128
    n = len(lst)

    if n < 8:
        out = 0.0
        for i in range(n):
            out += lst[i]
        return out
    elif n <= PW_BLOCKSIZE:
        r = [0.0] * 8
        r[0] = lst[0]
        r[1] = lst[1]
        r[2] = lst[2]
        r[3] = lst[3]
        r[4] = lst[4]
        r[5] = lst[5]
        r[6] = lst[6]
        r[7] = lst[7]
        for i in range(8, n - (n % 8), 8):
            r[0] += lst[0 + i]
            r[1] += lst[1 + i]
            r[2] += lst[2 + i]
            r[3] += lst[3 + i]
            r[4] += lst[4 + i]
            r[5] += lst[5 + i]
            r[6] += lst[6 + i]
            r[7] += lst[7 + i]
        res = ((r[0] + r[1]) + (r[2] + r[3])) +\
              ((r[4] + r[5]) + (r[6] + r[7]))
        
        # rest of 8-blocks
        for i in  range(n - (n % 8), n, 1):
            res += lst[i]
        
        return res
    else:
        n2 = n // 2
        n2 -= n2 % 8
        return my_pairwise(lst[:n2]) + my_pairwise(lst[n2:])

out = my_pairwise([1e-7] * int(1e7))
"%0.20f" % out
# '0.99999999999999988898'
# 再現できない

c言語 (こちらも再現できず)

https://paiza.io/projects/9ELMbLsizKFydjQTTMgaug?language=c

#include <stdio.h>

#define PW_BLOCKSIZE 128
#define IN_SIZE 10000000

double my_pairwise_sum(double *a, int n, int stride);


int main(void){
    double in_array[IN_SIZE];
    for (int i = 0; i < IN_SIZE; i++) {
        in_array[i] = 1.0 / IN_SIZE;
    }
    double out = my_pairwise_sum(in_array, IN_SIZE, 0);
    printf("%.20f\n", out);
}

double my_pairwise_sum(double *a, int n, int stride)
{
    if (n < 8) {
        int i;
        double res = 0.;
        for (i = 0; i < n; i++) {
            res += *(a + i * stride);
        }
        return res;
    }
    else if (n <= PW_BLOCKSIZE) {
        int i;
        double r[8], res;

        /*
         * sum a block with 8 accumulators
         * 8 times unroll reduces blocksize to 16 and allows vectorization with
         * avx without changing summation ordering
         */
        r[0] = *(a + 0 * stride);
        r[1] = *(a + 1 * stride);
        r[2] = *(a + 2 * stride);
        r[3] = *(a + 3 * stride);
        r[4] = *(a + 4 * stride);
        r[5] = *(a + 5 * stride);
        r[6] = *(a + 6 * stride);
        r[7] = *(a + 7 * stride);

        for (i = 8; i < n - (n % 8); i += 8) {
            /* small blocksizes seems to mess with hardware prefetch */
            r[0] += *(a + (i + 0) * stride);
            r[1] += *(a + (i + 1) * stride);
            r[2] += *(a + (i + 2) * stride);
            r[3] += *(a + (i + 3) * stride);
            r[4] += *(a + (i + 4) * stride);
            r[5] += *(a + (i + 5) * stride);
            r[6] += *(a + (i + 6) * stride);
            r[7] += *(a + (i + 7) * stride);
        }

        /* accumulate now to avoid stack spills for single peel loop */
        res = ((r[0] + r[1]) + (r[2] + r[3])) +
              ((r[4] + r[5]) + (r[6] + r[7]));

        /* do non multiple of 8 rest */
        for (; i < n; i++) {
            res += *(a + i * stride);
        }
        return res;
    }
    else {
        /* divide by two but avoid non-multiples of unroll factor */
        int n2 = n / 2;

        n2 -= n2 % 8;
        return my_pairwise_sum(a, n2, stride) +
               my_pairwise_sum(a + n2 * stride, n - n2, stride);
    }
}
// 0.99999999999999988898 

fsum

python code

def my_fsum(lst):
    res = 0.0 
    c = 0.0
 
    for l in lst:
        y = l - c
        t = res + y
        c = (t - res) - y
        res = t
 
    return res        

1.0
# 一致

basic sum

python code

def my_sum(lst):
    out = 0.0
    for l in lst:
        out += l
    return out

l =  [1/1e7] * int(1e7)
my_sum(l)
# 0.99999999975017

その他おまけ

# https://github.com/python/cpython/blob/c602c1be439e295fed9ebab47e895ef1d9df28be/Lib/statistics.py#L153
>>> import statistics
>>> statistics._sum([1e50, 1, -1e50] * 1000)
(<class 'float'>, Fraction(1000, 1), 3000)

import numpy as np
from math import fsum
from statistics import _sum

for i in range(1, 21):
    lst = [0.1] * i
    print(f"{sum(lst)=:.18f}, {np.sum(lst)=:.18f}, {fsum(lst)=:.18f}, {_sum(lst)[1]=}")

sum(lst)=0.100000000000000006, np.sum(lst)=0.100000000000000006, fsum(lst)=0.100000000000000006, _sum(lst)[1]=Fraction(3602879701896397, 36028797018963968)
sum(lst)=0.200000000000000011, np.sum(lst)=0.200000000000000011, fsum(lst)=0.200000000000000011, _sum(lst)[1]=Fraction(3602879701896397, 18014398509481984)
sum(lst)=0.300000000000000044, np.sum(lst)=0.300000000000000044, fsum(lst)=0.300000000000000044, _sum(lst)[1]=Fraction(10808639105689191, 36028797018963968)
sum(lst)=0.400000000000000022, np.sum(lst)=0.400000000000000022, fsum(lst)=0.400000000000000022, _sum(lst)[1]=Fraction(3602879701896397, 9007199254740992)
sum(lst)=0.500000000000000000, np.sum(lst)=0.500000000000000000, fsum(lst)=0.500000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 36028797018963968)
sum(lst)=0.599999999999999978, np.sum(lst)=0.599999999999999978, fsum(lst)=0.600000000000000089, _sum(lst)[1]=Fraction(10808639105689191, 18014398509481984)
sum(lst)=0.699999999999999956, np.sum(lst)=0.699999999999999956, fsum(lst)=0.700000000000000067, _sum(lst)[1]=Fraction(25220157913274779, 36028797018963968)
sum(lst)=0.799999999999999933, np.sum(lst)=0.800000000000000044, fsum(lst)=0.800000000000000044, _sum(lst)[1]=Fraction(3602879701896397, 4503599627370496)
sum(lst)=0.899999999999999911, np.sum(lst)=0.900000000000000022, fsum(lst)=0.900000000000000022, _sum(lst)[1]=Fraction(32425917317067573, 36028797018963968)
sum(lst)=0.999999999999999889, np.sum(lst)=1.000000000000000000, fsum(lst)=1.000000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 18014398509481984)
sum(lst)=1.099999999999999867, np.sum(lst)=1.100000000000000089, fsum(lst)=1.100000000000000089, _sum(lst)[1]=Fraction(39631676720860367, 36028797018963968)
sum(lst)=1.199999999999999956, np.sum(lst)=1.200000000000000178, fsum(lst)=1.200000000000000178, _sum(lst)[1]=Fraction(10808639105689191, 9007199254740992)
sum(lst)=1.300000000000000044, np.sum(lst)=1.300000000000000266, fsum(lst)=1.300000000000000044, _sum(lst)[1]=Fraction(46837436124653161, 36028797018963968)
sum(lst)=1.400000000000000133, np.sum(lst)=1.400000000000000355, fsum(lst)=1.400000000000000133, _sum(lst)[1]=Fraction(25220157913274779, 18014398509481984)
sum(lst)=1.500000000000000222, np.sum(lst)=1.500000000000000444, fsum(lst)=1.500000000000000000, _sum(lst)[1]=Fraction(54043195528445955, 36028797018963968)
sum(lst)=1.600000000000000311, np.sum(lst)=1.600000000000000089, fsum(lst)=1.600000000000000089, _sum(lst)[1]=Fraction(3602879701896397, 2251799813685248)
sum(lst)=1.700000000000000400, np.sum(lst)=1.700000000000000178, fsum(lst)=1.700000000000000178, _sum(lst)[1]=Fraction(61248954932238749, 36028797018963968)
sum(lst)=1.800000000000000488, np.sum(lst)=1.800000000000000266, fsum(lst)=1.800000000000000044, _sum(lst)[1]=Fraction(32425917317067573, 18014398509481984)
sum(lst)=1.900000000000000577, np.sum(lst)=1.900000000000000355, fsum(lst)=1.900000000000000133, _sum(lst)[1]=Fraction(68454714336031543, 36028797018963968)
sum(lst)=2.000000000000000444, np.sum(lst)=2.000000000000000444, fsum(lst)=2.000000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 9007199254740992)

公式ドキュメントでの言及

15. 浮動小数点演算、その問題と制限
を読むと良い。

x = 0.1
x.as_integer_ratio()
# (3602879701896397, 36028797018963968)

x.hex()
# '0x1.999999999999ap-4'

最後に

後半が疲れからかグダグダになった感は否めない。

参考

https://qiita.com/jo7ueb/items/25ede2bd48c4a2b322e2
https://qiita.com/yucatio/items/c03def631bc5eaaa9887
https://github.com/numpy/numpy/pull/10635/
https://yigarashi.hatenablog.com/entry/cpython
https://note.nkmk.me/python-float-hex/