はじめに
はじめまして、サイカのソフトウェアエンジニアの松嶋です。
サイカでは数理・統計処理を Python、Julia、R で実装しています。 今回、保守の観点から、数値計算の数値ズレについてお話ししようと思います。
皆さんは、数値計算の結果が、パッケージの更新前後や実行環境の違いで微妙に変わったことはないでしょうか? numpy のバージョンを上げるとテストが失敗するようになったり、ローカル環境ではテストが通るのにCI環境ではテストが失敗したり。
この記事では、Python を使って以下の数値変化をご紹介したいと思います。
- 入力変数の順序
- 数値解析ライブラリのバグ修正
- 加算アルゴリズム
尚、ここで扱うのは、情報落ちなどの浮動小数点演算で発生する計算誤差そのものの話ではないことにご注意ください。
入力変数の順序
最初に入力変数の順序による数値変化をご紹介します。
Python の dict() オブジェクトは、Python 3.5 までは順序情報を持ちませんが、 Python 3.7 から挿入順序が保存されるようになりました(What’s New In Python 3.7)。
線形最小二乗法を例に変数順序による数値変化を見てみましょう。
サンプルコード
以下は、範囲制約付き線形最小二乗法のサンプルコードです。
各変数のデータと範囲制約を辞書で管理するようにしています。この辞書から行列と範囲制約の配列を作り、予め定義したベクトルとともに入力パラメータとして scipy.optimize.lsq_linear() に渡し、線形最小二乗法問題を解きます。
sample_dict.py
import pandas from scipy.optimize import lsq_linear def main(): xs = { "x1": { "data": [8, -6], "bounds": (1, 19) }, "x2": { "data": [-4, 2], "bounds": (-14, 82) }, "x3": { "data": [-9, -4], "bounds": (-51, -20) } } A = [[x["data"][0] for x in xs.values()], [x["data"][1] for x in xs.values()]] b = [-1, 9] df = pandas.DataFrame(A, columns=xs.keys()) bounds = ([xs[c]["bounds"][0] for c in df.columns], [xs[c]["bounds"][1] for c in df.columns]) result = lsq_linear(df.values, b, bounds=bounds, method="bvls") merged = [(v, r) for v, r in zip(df.columns, result.x)] for i in sorted(merged): print(i) if __name__ == "__main__": main()
このコードは Python 3.7 以降の環境では実行結果は一定ですが、Python 3.5 環境で実行すると実行結果が変動します。
実行環境
| version | |
|---|---|
| python | 3.5.10 |
| numpy | 1.13.3 |
| scipy | 0.19.1 |
| pandas | 0.20.3 |
実行結果
$ python sample_dict.py ('x1', 19.0) ('x2', 70.899999999999991) ('x3', -20.0) $ python sample_dict.py ('x1', 19.0) ('x2', 70.900000000000006) ('x3', -20.0) $ python sample_dict.py ('x1', 19.0) ('x2', 70.899999999999991) ('x3', -20.0)
Python 3.5 環境でも dict を OrderedDict に書き換えると、実行結果は変わらなくなります。
sample_dict.py
... from collections import OrderedDict def main(): ... xs = OrderedDict([ ("x1", { "data": [8, -6], "bounds": (1, 19) }), ("x2", { "data": [-4, 2], "bounds": (-14, 82) }), ("x3", { "data": [-9, -4], "bounds": (-51, -20) }) ]) ...
実行結果
$ python sample_dict.py ('x1', 19.0) ('x2', 70.899999999999991) ('x3', -20.0)
入力パラメータの行列は、辞書の各変数データから組み立てられるため、変数の順序が不定の場合は常に同じ行列になりません。そのため、使うデータは同じですが並べ方により行列計算の結果が変わっていました。
このような事象を防ぐためにも、使用するアルゴリズムが入力データをどのように処理するかを把握しておくことは大事だと思います。
数値解析ライブラリのバグ修正
続いて、数値解析ライブラリのアップデートによる数値の変化をご紹介します。
先ほどの線形最小二乗法のサンプルコードで使用した scipy.optimize.lsq_linear の method bvls は、バージョン 1.6.0 でバグが修正されました。
サンプルコード
以下は dict の変数順序の時と同様の、範囲制約付きの線形最小二乗法のサンプルコードです。
入力は行列 A とベクトル b、範囲制約 bounds を scipy.optimize.lsq_linear() に渡し、線形最小二乗法問題を解きます。
sample_bvls.py
import numpy from scipy.optimize import lsq_linear def main(): seed = 0 numpy.random.seed(seed) n, m = 3, 2 A = numpy.random.randint(-10, 10, size=[m, n]) b = numpy.random.randint(-10, 10, size=[m]) bounds = ((1, -14, -51), (19, 82, -20)) result = lsq_linear(A, b, bounds=bounds, method="bvls") merged = [(v, r) for v, r in zip(["x1", "x2", "x3"], result.x)] for i in sorted(merged): print(i) if __name__ == "__main__": main()
実行環境
| version | |
|---|---|
| python | 3.10.12 |
| numpy | 1.23.5 |
| scipy | 1.9.3 |
| pandas | 1.5.1 |
バグ修正後(scipy 1.6.0 以降)の実行結果
$ python sample_bvls.py ('x1', 14.73584905660377) ('x2', -14.0) ('x3', -20.0)
次に、バグ修正前の scipy 1.5.4 の bvls のコードを抜き出し、モジュールを置き換えることでバグ修正前の挙動を再現してみます。
bvls.py
import importlib import sys import numpy as np from numpy.linalg import lstsq, norm from scipy.optimize import OptimizeResult, lsq_linear from scipy.optimize._lsq import bvls, common # https://github.com/scipy/scipy/blob/19acfed431060aafaa963f7e530c95e70cd4b85c/scipy/optimize/_lsq/bvls.py#L17 def bvls_v1_5_4(A, b, x_lsq, lb, ub, tol, max_iter, verbose): m, n = A.shape x = x_lsq.copy() on_bound = np.zeros(n) mask = x < lb x[mask] = lb[mask] on_bound[mask] = -1 mask = x > ub x[mask] = ub[mask] on_bound[mask] = 1 free_set = on_bound == 0 active_set = ~free_set free_set, = np.nonzero(free_set) r = A.dot(x) - b cost = 0.5 * np.dot(r, r) initial_cost = cost g = A.T.dot(r) cost_change = None step_norm = None iteration = 0 if verbose == 2: common.print_header_linear() # This is the initialization loop. The requirement is that the # least-squares solution on free variables is feasible before BVLS starts. # One possible initialization is to set all variables to lower or upper # bounds, but many iterations may be required from this state later on. # The implemented ad-hoc procedure which intuitively should give a better # initial state: find the least-squares solution on current free variables, # if its feasible then stop, otherwise, set violating variables to # corresponding bounds and continue on the reduced set of free variables. while free_set.size > 0: if verbose == 2: optimality = bvls.compute_kkt_optimality(g, on_bound) common.print_iteration_linear(iteration, cost, cost_change, step_norm, optimality) iteration += 1 x_free_old = x[free_set].copy() A_free = A[:, free_set] b_free = b - A.dot(x * active_set) z = lstsq(A_free, b_free, rcond=-1)[0] lbv = z < lb[free_set] ubv = z > ub[free_set] v = lbv | ubv if np.any(lbv): ind = free_set[lbv] x[ind] = lb[ind] active_set[ind] = True on_bound[ind] = -1 if np.any(ubv): ind = free_set[ubv] x[ind] = ub[ind] active_set[ind] = True on_bound[ind] = 1 ind = free_set[~v] x[ind] = z[~v] r = A.dot(x) - b cost_new = 0.5 * np.dot(r, r) cost_change = cost - cost_new cost = cost_new g = A.T.dot(r) step_norm = norm(x[free_set] - x_free_old) if np.any(v): free_set = free_set[~v] else: break if max_iter is None: max_iter = n max_iter += iteration termination_status = None # Main BVLS loop. optimality = bvls.compute_kkt_optimality(g, on_bound) for iteration in range(iteration, max_iter): if verbose == 2: common.print_iteration_linear(iteration, cost, cost_change, step_norm, optimality) if optimality < tol: termination_status = 1 if termination_status is not None: break move_to_free = np.argmax(g * on_bound) on_bound[move_to_free] = 0 free_set = on_bound == 0 active_set = ~free_set free_set, = np.nonzero(free_set) x_free = x[free_set] x_free_old = x_free.copy() lb_free = lb[free_set] ub_free = ub[free_set] A_free = A[:, free_set] b_free = b - A.dot(x * active_set) z = lstsq(A_free, b_free, rcond=-1)[0] lbv, = np.nonzero(z < lb_free) ubv, = np.nonzero(z > ub_free) v = np.hstack((lbv, ubv)) if v.size > 0: alphas = np.hstack(( lb_free[lbv] - x_free[lbv], ub_free[ubv] - x_free[ubv])) / (z[v] - x_free[v]) i = np.argmin(alphas) i_free = v[i] alpha = alphas[i] x_free *= 1 - alpha x_free += alpha * z if i < lbv.size: on_bound[free_set[i_free]] = -1 else: on_bound[free_set[i_free]] = 1 else: x_free = z x[free_set] = x_free step_norm = norm(x_free - x_free_old) r = A.dot(x) - b cost_new = 0.5 * np.dot(r, r) cost_change = cost - cost_new if cost_change < tol * cost: termination_status = 2 cost = cost_new g = A.T.dot(r) optimality = bvls.compute_kkt_optimality(g, on_bound) if termination_status is None: termination_status = 0 return OptimizeResult( x=x, fun=r, cost=cost, optimality=optimality, active_mask=on_bound, nit=iteration + 1, status=termination_status, initial_cost=initial_cost) bvls.bvls = bvls_v1_5_4 importlib.reload(sys.modules['scipy.optimize._lsq.lsq_linear'])
sample_bvls.py
... def main(): ... import bvls result = lsq_linear(df.values, b, bounds=bounds, method="bvls") ... ...
実行結果
$ python sample_bvls.py ('x1', 6.576271186440669) ('x2', -14.0) ('x3', -20.0)
バグの修正前後の実行結果を比較すると、x1 の値は修正前の 6.576271186440669 から修正後は 14.73584905660377 に変わりました。
結構数値が変わりましたが、一方でこのバグ修正により bvls を使った lsq_linear の結果が全て変わる訳ではなく、扱う入力データによっては数値変化に気づきにくいかもしれません。
数値解析ライブラリのバグ修正は、場合によっては大きく数値が変化する可能性があるので、保守する場合は注意が必要です。
加算アルゴリズム
最後に加算アルゴリズムによる数値変化についてご紹介します。
Pandas や Numpy、Python の sum() は加算アルゴリズムが同じではありません。また Pandas でも resample() と DataFrame の sum() ではアルゴリズムが違います。
実行環境
| version | |
|---|---|
| python | 3.11.9 |
| numpy | 2.1.1 |
| pandas | 2.2.2 |
Pandas の sum() の比較
Pandas 1.2.0 では、 resample の sum() の加算アルゴリズムが Kahan summation になりました。一方、 DataFrame の sum() は、 Numpy ndarray.sum() と同じ計算方法になっています。
sample_sum1.py
import numpy import pandas df = pandas.DataFrame( data={"x": [0.1]*14}, index=pandas.date_range('2018-01-01', '2018-01-14', freq="D") ) a = df.sum().values[0] b = df.resample('W').sum().sum().values[0] numpy.testing.assert_equal(a, b)
実行結果
$ python sample_sum1.py AssertionError: Items are not equal: ACTUAL: np.float64(1.4000000000000004) DESIRED: np.float64(1.4000000000000001)
Numpy と Python の sum() の比較
また、 Numpy ndarray.sum() と Python の Built-in の sum() でも、加算アルゴリズムの違いにより、同じ結果にならないことがあります。
sample_sum2.py
import numpy arr = numpy.array([0.1]*14) a = arr.sum() b = sum(arr) numpy.testing.assert_equal(a, b)
実行結果
$ python sample_sum2.py Items are not equal: ACTUAL: np.float64(1.4000000000000004) DESIRED: np.float64(1.4000000000000001)
数値の加算処理は Python に限らず他の言語でも広く行われていますが、言語やライブラリによって採用している加算アルゴリズムが違うことがあり、サービス間でのやりとり等は注意した方がいいかもしれません。
まとめ
3つのケースをご紹介しましたが、数値ズレには様々な原因があります。 今回扱わなかったものとしては、CPU アーキテクチャやコンパイラ、OpenBLAS などの線形計算ライブラリなどは数値ズレの原因として特に注意した方がよいでしょう。
数値計算は、用途によって求められる数値の精度は様々ですが、数値が変化したときに原因を把握しておくことは大事だと思います。 数値の変化がアルゴリズム改善によるものなのか、パッケージの更新によるものなのか、数値変化の原因を把握しておくことは、数値計算の再現性の面でも重要になります。
また、今回ご紹介したような数値変化の検知については、数値の誤差を許容したテストでは見逃しやすく、気付くのも遅れやすいです。 扱うアルゴリズムや数値計算の特性を把握し微細な数値変化の最終的な影響が分かると、それが無視できるものなのか判断しやすくなり、保守もやり易くなると思います。
数値計算の細かい話しになりましたが、少しでもお役に立てれば幸いです。
最後まで読んでいただきありがとうございました。
