最小二乗法をscipy.optimize.leastsqを用いて行ったプログラムの処理速度を上げたい

Question

ある25個のデータでz軸に観測結果、x, yで観測値の位置三次元で表現されるものを三次元のガウス関数の積分によってfittingをおこないました。

このプログラムの最小二乗法を行っているoptimeze.leastsqの部分での所要時間が10秒ほどかかっています。ここで、同様の処理を1800回ほど連続して行うため、総処理時間が5時間ほどになってしまいます。

leastsqの引数を変更してみたりしたのですが、あまり変化が見らなかったです（中身を正確に理解しているとは言えないので、方法が悪い可能性はあり）。

よって処理速度の向上を目標としているのですが、良い方法が見つからないため、詳しい方がおりましたらご指摘していただければ幸いです。

修正後のプログラム

import time
import numpy as np
from scipy import optimize
from sympy import *

a = Symbol('a')
f = 1/sqrt(2)*exp(-(a**2/2))
F = integrate(f, a)
func_gauss = lambdify(a, F)

#二次元ガウスの積分
def integ_gaussian(param, x_range, y_range):
    height = param[0]
    sigma = param[1]
    cen_x = param[2]
    cen_y = param[3]
    constant = param[4]
    #x_list, y_listに範囲内の不定積分を格納する
    #積分する回数を減らすことができる
    x_list = [func_gauss((x-cen_x)/sigma) for x in x_range]
    y_list = [func_gauss((y-cen_y)/sigma) for y in y_range]
    #np.diffは前後の差を求める
    #二次元ガウスの積分値が求まる
    diff_x = np.diff(x_list)
    diff_y = np.diff(y_list)
    #np.meshgridで格子座標を生成し、x, yの二次元ガウスの積分値の掛け算を網羅する
    mesh_x, mesh_y = np.meshgrid(diff_x, diff_y)
    np_integ = (height * mesh_x * mesh_y + constant).ravel()
    return np_integ

#実測値 ― 計算値から最小二乗法を行う
def residuals_2D(param, z, x_range, y_range):
    return (z - integ_gaussian(param, x_range, y_range))

#第１引数はガウス関数のheightで体積を表す。観測データ中心の最も値が高いものの1/3としている
#第２引数はガウス関数のσでグラフの拡がりを表す。
#第３引数はガウス関数の中心のx座標
#第４引数はガウス関数の中心のy座標
#第５引数はガウス関数のコンスタントで、5×5の観測データの最外のデータ値の平均をとったもの
param_2D = [3221.33/3, 0.3, 0, 0, 307.30625

#観測データ
actual_data = [276.67, 336.11, 248.22, 294.89, 290.44, 
               285.33, 438.78, 472.89, 289.33, 416.78, 
               343.0, 2456.44, 3618.44, 308.78, 281.56,
               279.78, 456.56, 492.89, 273.22, 282.56,
               272.22, 444.0, 281.89, 269.67, 313.78]

#積分範囲を指定するのに使用
np_range = np.arange(-2.5, 2.6)

#optimised_param_2Dにfitting結果のパラメータが格納される
#optimised_param_2Dの所要時間を記録する。
time_1 = time.time()
optimised_param_2D = optimize.leastsq(residuals_2D, param_2D, args = (actual_data, np_range, np_range, ftol=1.49012e-1, xtol=1.49012e-1, gtol=1.49012e-1, maxfev=40))
time_2 =time.time()
print(time_2-time_1)

処理時間
time_2 - time_1 = 0.007936954498291016

fitting結果
height:1943.76
sigma:0.31
center_x:-0.41
center_y:0.01
constant:304.71

Comments

[jbpb0]

「scipy.optimize.leastsq」のパラメータ「maxfev」で、最適化で「residuals_2D」を最大何回実行するのかを指定できます

デフォルト値は
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html
に書いてあるように「200*(最適化変数の数+1)」なので、質問の条件だとN=5で200*(5+1)=1200回になりますから、それよりも小さい値を指定したら早く終わるかもしれません

ただし、あくまでも最大の回数なので、実際は1200回よりも少ない回数で終わってるかもしれません
一旦、(結果の精度は気にせずに)思いっきり少なく指定してみて、どれくらいの数値だと効果があるのかを確認してみてください

たとえば、10を指定して早くなったら次は100を指定して、それでもまだ早いなら200を指定して、みたいに徐々に数値を増やして、どこまで数値を増やすと効果が無くなるのかを確認する
また、各数値での「fit」の精度を確認して、精度的に十分になる数値も確認する
回数が少ないと、当然精度は低下しますが、それでも必要な精度を満たしてればいいわけです
必要以上の精度まで出そうとガンバって計算して、必要以上に時間かけても意味無いですから

[jbpb0]

上記のように直接回数を制限する方法の他に、精度を落として少ない回数で終わるようにする方法もあります
それを行うのが、「scipy.optimize.leastsq」のパラメータ「ftol」「xtol」「gtol」の指定です

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html
に書かれてるデフォルト値よりも大きな数値を指定したら、早く精度が達成できて早く終わるかもしれません

ただし、当然のことながら本当に必要な精度よりも落としてはダメです
本当に必要な精度ギリギリ達成したら止まるように設定するのがポイントです

こちらも、一旦は思いっきり大きく指定してみて早く終わるようにして、その時の精度を確認して、精度が悪すぎるなら指定する数値を小さくして確認、とやるといいと思います

[bsdfan]

単純にinteg_gaussianに時間がかかりすぎているのでは？
sympyを使ったり、forループを使ったりしている部分は改善できないでしょうか？

[python_wakaran]

jbpb0さん、bsdfanさんご指摘ありがとうございます。

まずjbpb0さんのおっしゃる2通りの方法を試して、最終的に
optimised_param_2D =
optimize.leastsq(residuals_2D, param_2D, args = (z, np_range, np_range), ftol=1.49012e-1, xtol=1.49012e-1, gtol=1.49012e-1, maxfev=40)としました。
今回はfittingのためのデータが少なかったからか、精度をだいぶ落としても問題なく作動しました。もう一桁大きくしてしまうと、さすがに値がおかしくなってしまいました。
またmaxfevはかなり小さくしても導かれる値としては問題なかったのですが、
RuntimeWarning: Number of calls to function has reached maxfev = 10.
などのエラーが発生してしまうため、このエラーが発生しない値で可能な限り小さくしました。

bsdfanさんに指摘されたsympyに関しては、不定積分も用いて積分を行いたいと考えているので、こちらは変えるのは難しいかなと考えています。しかし、積分に関しては以前scipy.integrateを用いて行っていたので、そちらと比較してみるのも一つの方法かなと考えることができました。後日実装してみようと思います。
また、ご指摘の通りforループに関しては内包表現に変更いたしました。私自身、なんとなく内包表現は避けてきたのですが（なんかめんどくさそうという理由）、調べてみると内包表現のほうがforループと比較した際に速度が速いことが学べました。

以上お二人のご指摘の結果、速度は7.15秒ほどに短縮されました。
同時にこれが現在のプログラムでの限界速度かなとも感じています。
ご指摘のほどありがとうございました。

[jbpb0]

・heightを掛ける
・constantを足す
はinteg_listの全要素に共通なので、一括で行えるような

・x_list[i+1] - x_list[i]
・y_list[j+1] - y_list[j]
はnumpy.diffで一括計算できるので、その後でnumpy.meshgridでそれぞれ2次元化して、(要素同士の)掛け算したら、
integ_list =…
の二重のforループは無くせるような

参考
https://deepage.net/features/numpy-diff.html
https://deepage.net/features/numpy-meshgrid.html
https://algorithm.joho.info/programming/python-numpy-multiplication/

[python_wakaran]

jbpd0さんのご指摘の通りnumpy.diffとnumpy.meshgridを使用する方法に変更しました。
丁寧なご指摘のおかげで、大変勉強になりました。ありがとうございました！

[jbpb0]

np_integに対して、「.ravel()」を2回やってますよ

[python_wakaran]

本当ですね。失礼しました。

Answer 1

下記でinteg_gaussian()の速度を改善できそうです。

• func_gauss は、毎回作り直す必要はないので、関数の外で作成
• func_gauss は ufunc になっている(はずな)ので、forは使わず、引数にarrayを渡して計算
• 二重のforループはブロードキャストを使った演算により削除

一番上の、sympyの部分を関数の外に出す変更だけでも、かなり変わります。

import time
import numpy as np
from scipy import optimize
#from scipy import integrate
from sympy import *


#func_gaussで不定積分を定義する
a = Symbol('a')
f = 1/sqrt(2)*exp(-(a**2/2))
F = integrate(f, a)
func_gauss = lambdify(a, F)

#二次元ガウスの積分
def integ_gaussian(param, x_range, y_range):
    height = param[0]
    sigma = param[1]
    cen_x = param[2]
    cen_y = param[3]
    constant = param[4]
 
    #x_list, y_listに範囲内の不定積分を格納する
    #積分する回数を減らすことができる
    x_list = func_gauss((x_range - cen_x) / sigma)
    y_list = func_gauss((y_range - cen_y) / sigma)
    #不定積分に積分したい範囲を代入して積分値を求める
    np_integ = (height * np.diff(x_list)[:, None] * np.diff(y_list) + constant).ravel()
    return np_integ

Comments

[python_wakaran]

回答ありがとうございます。
以上のプログラムを参考にして、jbpb0さんの指摘も参考にしたプログラムをに修正しました。
その結果、処理速度は0.0079と圧倒的に早くなりました。
bsdfanさんのおっしゃる通りsympy部分を外に出したことが、これほど高速化できた要因であると考えられます。丁寧なアンサーありがとうございました。