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

Question

ある25個のデータでz軸に観測結果、x, yで観測値の位置三次元で表現されるものを三次元のガウス関数の積分によってfittingをおこないました。

このプログラムの最小二乗法を行っているoptimeze.leastsqの部分での所要時間が10秒ほどかかっています。ここで、同様の処理を1800回ほど連続して行うため、総処理時間が5時間ほどになってしまいます。

leastsqの引数を変更してみたりしたのですが、あまり変化が見らなかったです（中身を正確に理解しているとは言えないので、方法が悪い可能性はあり）。

よって処理速度の向上を目標としているのですが、良い方法が見つからないため、詳しい方がおりましたらご指摘していただければ幸いです。


修正後のプログラム
```python
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

Accepted Answer

下記でinteg_gaussian()の速度を改善できそうです。
- func_gauss は、毎回作り直す必要はないので、関数の外で作成
- func_gauss は ufunc になっている(はずな)ので、forは使わず、引数にarrayを渡して計算
- 二重のforループはブロードキャストを使った演算により削除

一番上の、sympyの部分を関数の外に出す変更だけでも、かなり変わります。

```python
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
```

関連した質問