scipy - 如何使用 SciPy 最小化重用 Jacobian 和 Inverse Hessian

Question

我目前正在使用 SciPy.optimize 中的最小化 (BFGS) 函数来校准我的模型。校准后，参数会不时受到轻微扰动，我想重新使用之前校准的 Jacobian 和 IHessian 作为我的初始 ste(一种对 jac 和 hessian 的初始猜测，然后切换到这些的常规计算)

知道怎么做吗？

Answer 1

你必须稍微修改 scipy 才能做到这一点。 BFGS 的 Hessian 近似存储在变量 Hk 中，首先在文件 _optimize.py 中函数 _minimize_bfgs() 的第 1308 行定义:https://github.com/scipy/scipy/blob/65d8361519e3f1166d81dceef9f1d6375b4ed808/scipy/optimize/_optimize.py#L1308

你将不得不修改 scipy 的代码返回或存储这个矩阵，并在下次将它作为参数再次传递给函数。

---

或者，您可以改用 scipy 的非线性求根，因为最小化函数与找到点 x 相同，其中函数的梯度 g(x) 为零。需要注意的是，从算法的角度来看，这远不够稳健。寻根器可能会收敛到鞍点或最大值等。

用于优化的 BFGS 算法与用于寻找非线性函数根的 Broyden 方法非常相似。 Broyden 是 BFGS 中的 B，BFGS 基本上是应用于方程 g(x)=0 的 Broyden 求根方法的修改版本。

从代码的角度来看，与 BFGS 相比，scipy 可以更好地访问 Broyden 非线性求根器的内部结构。 Broyden Jacobian 矩阵(与 BFGS Hessian 类似)存储为“Jacobian”类型的对象，如果您覆盖“setup()”方法不执行任何操作，则可以重新使用该对象，因此它不会重设。例如:

import numpy as np
from scipy.optimize import nonlin, rosen, rosen_der


x0 = np.array([0.0, 0.0])

counter = [1]
def callback(x, Fx):
    print('k=', counter[0], '   xk=', x, '   f(xk)=', Fx)
    counter[0] += 1

print('first solve:')
jac = nonlin.BroydenFirst()
sol1 = nonlin.nonlin_solve(rosen_der, x0, jacobian=jac, callback=callback, line_search='wolfe')

print('')
print('Jacobian after first solve:')
print(jac.todense())


print('')
print('second solve:')
jac.setup = lambda x,y,z: None # Prevent Jacobian from resetting. Try commenting this line out!
counter = [1]
sol2 = nonlin.nonlin_solve(rosen_der, x0, jacobian=jac, callback=callback, line_search='wolfe')

print('')
print('Jacobian after second solve:')
print(jac.todense())

这会产生以下输出:

first solve:
k= 1    xk= [-0.5  0. ]    f(xk)= [-53. -50.]
k= 2    xk= [-0.32283933  0.16713271]    f(xk)= [ 5.47792447 12.58149553]
k= 3    xk= [-0.28061925  0.08236319]    f(xk)= [-2.1553475  0.7232059]
k= 4    xk= [-0.21641127  0.02839649]    f(xk)= [-4.02884252 -3.68746961]
k= 5    xk= [-0.18109429  0.02229409]    f(xk)= [-3.12286088 -2.10021066]
k= 6    xk= [-0.19320298  0.02730297]    f(xk)= [-3.16110526 -2.00488442]
k= 7    xk= [-0.1500595   0.01281555]    f(xk)= [-2.88248831 -1.94046136]
k= 8    xk= [-0.14063472  0.00987684]    f(xk)= [-2.83825517 -1.98025685]
k= 9    xk= [-0.10559036  0.00122285]    f(xk)= [-2.63043649 -1.98529385]
k= 10    xk= [-0.07978769 -0.00399875]    f(xk)= [-2.49036939 -2.07296403]
k= 11    xk= [-0.04033837 -0.00840453]    f(xk)= [-2.24254189 -2.00634254]
k= 12    xk= [-0.01281811 -0.01000148]    f(xk)= [-2.07775868 -2.03315754]
k= 13    xk= [ 0.04612425 -0.00598848]    f(xk)= [-1.75801514 -1.62318494]
k= 14    xk= [ 0.03215952 -0.00653348]    f(xk)= [-1.83833134 -1.51354257]
k= 15    xk= [0.08036797 0.0021266 ]    f(xk)= [-1.6999891  -0.86648289]
k= 16    xk= [-0.04367326 -0.01098396]    f(xk)= [-2.31254872 -2.57826183]
k= 17    xk= [-0.07820406 -0.00365169]    f(xk)= [-2.46195331 -1.95351227]
k= 18    xk= [-0.02940764 -0.00337684]    f(xk)= [-2.10871    -0.84832937]
k= 19    xk= [-0.0314977  -0.00316138]    f(xk)= [-2.11532548 -0.83069659]
k= 20    xk= [-0.00053675 -0.00472907]    f(xk)= [-2.00208887 -0.94587179]
k= 21    xk= [ 0.01086898 -0.00493684]    f(xk)= [-1.95628507 -1.01099561]
k= 22    xk= [ 0.03866879 -0.00408101]    f(xk)= [-1.83641115 -1.11525685]
k= 23    xk= [ 0.06232922 -0.00227391]    f(xk)= [-1.72179137 -1.23176741]
k= 24    xk= [0.09084425 0.0019575 ]    f(xk)= [-1.5895594 -1.2590346]
k= 25    xk= [0.11524109 0.00670698]    f(xk)= [-1.46650132 -1.3147068 ]
k= 26    xk= [0.17879539 0.0246649 ]    f(xk)= [-1.12012011 -1.46057767]
k= 27    xk= [0.17505925 0.02929594]    f(xk)= [-1.5553635  -0.26996006]
k= 28    xk= [0.18893989 0.03342207]    f(xk)= [-1.4500936  -0.45524167]
k= 29    xk= [0.22546628 0.04688035]    f(xk)= [-1.19240743 -0.79093866]
k= 30    xk= [0.23877529 0.05259369]    f(xk)= [-1.10029946 -0.88399008]
k= 31    xk= [0.25428437 0.06010509]    f(xk)= [-1.02807924 -0.91109027]
k= 32    xk= [0.27116407 0.06871439]    f(xk)= [-0.9353492  -0.96311185]
k= 33    xk= [0.28827306 0.07809215]    f(xk)= [-0.84584567 -1.00184213]
k= 34    xk= [0.30453673 0.08760267]    f(xk)= [-0.76480491 -1.0279903 ]
k= 35    xk= [0.32072413 0.09764762]    f(xk)= [-0.68934851 -1.04326926]
k= 36    xk= [0.33639066 0.10790526]    f(xk)= [-0.6203388  -1.05068295]
k= 37    xk= [0.35201809 0.11866338]    f(xk)= [-0.55625337 -1.05067108]
k= 38    xk= [0.36724683 0.12964277]    f(xk)= [-0.49759833 -1.04549305]
k= 39    xk= [0.38246707 0.14110658]    f(xk)= [-0.44343805 -1.03489669]
k= 40    xk= [0.39728775 0.15273218]    f(xk)= [-0.39410312 -1.02107526]
k= 41    xk= [0.41215838 0.16485939]    f(xk)= [-0.34886931 -1.00302937]
k= 42    xk= [0.42650931 0.17699212]    f(xk)= [-0.30793999 -0.9836144 ]
k= 43    xk= [0.44109504 0.18976381]    f(xk)= [-0.27072811 -0.9602033 ]
k= 44    xk= [0.45477201 0.20212697]    f(xk)= [-0.23719208 -0.93812269]
k= 45    xk= [0.4681153  0.21458789]    f(xk)= [-0.21291526 -0.9088083 ]
k= 46    xk= [0.47908973 0.22505828]    f(xk)= [-0.1854592  -0.89373795]
k= 47    xk= [0.50398499 0.2497432 ]    f(xk)= [-0.13370942 -0.85153389]
k= 48    xk= [0.51439976 0.26044589]    f(xk)= [-0.11498655 -0.83224564]
k= 49    xk= [0.53852539 0.28603837]    f(xk)= [-0.06750659 -0.79424541]
k= 50    xk= [0.55785196 0.30748713]    f(xk)= [-0.05606785 -0.74233693]
k= 51    xk= [0.56092826 0.31093463]    f(xk)= [-0.0466508  -0.74117561]
k= 52    xk= [0.61261365 0.37369717]    f(xk)= [-0.38311202 -0.31966368]
k= 53    xk= [0.52465657 0.27398909]    f(xk)= [-0.68302112 -0.25508661]
k= 54    xk= [0.74386124 0.52529992]    f(xk)= [ 7.82778407 -5.60592561]
k= 55    xk= [0.54293777 0.30014958]    f(xk)= [-2.07995492  1.07363175]
k= 56    xk= [0.58511821 0.34990237]    f(xk)= [-2.59425851  1.50781064]
k= 57    xk= [0.35868093 0.08830817]    f(xk)= [ 4.50558888 -8.06876881]
k= 58    xk= [0.47983856 0.24586376]    f(xk)= [-4.0381075   3.12374291]
k= 59    xk= [0.3775732 0.1412489]    f(xk)= [-1.04660917 -0.26252449]
k= 60    xk= [0.36511619 0.13121294]    f(xk)= [-0.96352471 -0.41937734]
k= 61    xk= [ 0.1114212  -0.05351408]    f(xk)= [  1.16118684 -13.18575168]
k= 62    xk= [0.25437135 0.02356053]    f(xk)= [ 2.69511008 -8.22885014]
k= 63    xk= [0.16777908 0.02828527]    f(xk)= [-1.67353225  0.0270904 ]
k= 64    xk= [0.21546576 0.04636992]    f(xk)= [-1.56427885 -0.01111462]
k= 65    xk= [0.24867847 0.0590896 ]    f(xk)= [-1.22895955 -0.55027587]
k= 66    xk= [0.26241736 0.06557769]    f(xk)= [-1.13033017 -0.65703565]
k= 67    xk= [0.27594487 0.07265633]    f(xk)= [-1.06297502 -0.69784814]
k= 68    xk= [0.29626637 0.08383401]    f(xk)= [-0.94058127 -0.78794969]
k= 69    xk= [0.31324088 0.09384403]    f(xk)= [-0.83777351 -0.85516414]
k= 70    xk= [0.3305509  0.10477546]    f(xk)= [-0.74543613 -0.89768638]
k= 71    xk= [0.34571472 0.11489762]    f(xk)= [-0.66954478 -0.9242097 ]
k= 72    xk= [0.36107381 0.12568167]    f(xk)= [-0.60009902 -0.9385247 ]
k= 73    xk= [0.37526559 0.13609304]    f(xk)= [-0.53928208 -0.94624546]
k= 74    xk= [0.38966873 0.14710916]    f(xk)= [-0.48301011 -0.94651221]
k= 75    xk= [0.40322511 0.15787327]    f(xk)= [-0.43270865 -0.94344462]
k= 76    xk= [0.41704718 0.16925377]    f(xk)= [-0.3860966  -0.93491704]
k= 77    xk= [0.4300686  0.18033373]    f(xk)= [-0.34418922 -0.92505424]
k= 78    xk= [0.44351101 0.1921499 ]    f(xk)= [-0.30541382 -0.9104218 ]
k= 79    xk= [0.45596    0.20341705]    f(xk)= [-0.27054816 -0.89649514]
k= 80    xk= [0.46923896 0.2158001 ]    f(xk)= [-0.23845909 -0.87701904]
k= 81    xk= [0.48085841 0.22691624]    f(xk)= [-0.2095579  -0.86171446]
k= 82    xk= [0.49450418 0.24035192]    f(xk)= [-0.18369471 -0.83649135]
k= 83    xk= [0.50433095 0.25022809]    f(xk)= [-0.15987461 -0.82432328]
k= 84    xk= [0.52923466 0.27606143]    f(xk)= [-0.08884812 -0.80558078]
k= 85    xk= [0.54817807 0.29670844]    f(xk)= [-0.07244057 -0.75815082]
k= 86    xk= [0.55196258 0.30088075]    f(xk)= [-0.06107966 -0.7563875 ]
k= 87    xk= [0.58409242 0.3383646 ]    f(xk)= [-0.17778262 -0.55987076]
k= 88    xk= [0.57415074 0.32741847]    f(xk)= [-0.33941812 -0.44612013]
k= 89    xk= [0.5759845  0.32975938]    f(xk)= [-0.38752971 -0.39975147]
k= 90    xk= [0.54936923 0.30264256]    f(xk)= [-1.08497164  0.16720094]
k= 91    xk= [0.59588566 0.35541025]    f(xk)= [-0.88701248  0.06610646]
k= 92    xk= [0.62557699 0.38945577]    f(xk)= [-0.27571094 -0.37815896]
k= 93    xk= [0.62659512 0.39086428]    f(xk)= [-0.30639535 -0.35143459]
k= 94    xk= [0.60750698 0.37014954]    f(xk)= [-1.04859775  0.21696188]
k= 95    xk= [0.63679781 0.40648354]    f(xk)= [-0.97401477  0.19441837]
k= 96    xk= [0.67603801 0.4552325 ]    f(xk)= [-0.16255821 -0.35897816]
k= 97    xk= [0.67756175 0.457579  ]    f(xk)= [-0.23538087 -0.30218326]
k= 98    xk= [0.66719264 0.44635928]    f(xk)= [-0.98940486  0.24265117]
k= 99    xk= [0.68116519 0.46516345]    f(xk)= [-0.95848113  0.23548731]
k= 100    xk= [0.72613353 0.52567098]    f(xk)= [-0.08332021 -0.31978466]
k= 101    xk= [0.72763167 0.52819241]    f(xk)= [-0.17934169 -0.25108513]
k= 102    xk= [0.72088263 0.52080286]    f(xk)= [-0.88438909  0.22621877]
k= 103    xk= [0.72930705 0.53299447]    f(xk)= [-0.86394012  0.22113746]
k= 104    xk= [0.77143759 0.59395855]    f(xk)= [-0.09998023 -0.2314799 ]
k= 105    xk= [0.772421   0.59561435]    f(xk)= [-0.14005778 -0.20396923]
k= 106    xk= [0.7641847  0.58495593]    f(xk)= [-0.77048152  0.19553579]
k= 107    xk= [0.77477743 0.60123328]    f(xk)= [-0.74585811  0.19064376]
k= 108    xk= [0.81159937 0.65777517]    f(xk)= [-0.07866333 -0.1836731 ]
k= 109    xk= [0.81222326 0.65888019]    f(xk)= [-0.10705352 -0.16528704]
k= 110    xk= [0.80623427 0.65086575]    f(xk)= [-0.66231087  0.17040917]
k= 111    xk= [0.81361762 0.66280639]    f(xk)= [-0.64378248  0.16655104]
k= 112    xk= [0.84654107 0.71601533]    f(xk)= [-0.09817672 -0.12329062]
k= 113    xk= [0.85700521 0.73465556]    f(xk)= [-0.35373609  0.03952514]
k= 114    xk= [0.83758317 0.70213926]    f(xk)= [-0.52374073  0.1187387 ]
k= 115    xk= [0.90388038 0.81436146]    f(xk)= [ 0.76163845 -0.52765703]
k= 116    xk= [0.86767736 0.75453523]    f(xk)= [-0.84467781  0.33424436]
k= 117    xk= [0.90054052 0.81082101]    f(xk)= [-0.1440905  -0.03044198]
k= 118    xk= [0.90932277 0.82694243]    f(xk)= [-0.20846799  0.01490863]
k= 119    xk= [0.87018732 0.75667299]    f(xk)= [-0.06714585 -0.1105966 ]
k= 120    xk= [0.82135459 0.67121157]    f(xk)= [ 0.76362535 -0.68235826]
k= 121    xk= [0.86037593 0.74128628]    f(xk)= [-0.63700442  0.20790696]
k= 122    xk= [0.83873931 0.70465215]    f(xk)= [-0.71455427  0.23370366]
k= 123    xk= [1.03860171 1.04308697]    f(xk)= [14.86961426 -7.12131062]
k= 124    xk= [0.84811163 0.7215493 ]    f(xk)= [-1.06909728  0.45119093]
k= 125    xk= [0.86097865 0.74383034]    f(xk)= [-1.15489858  0.50922045]
k= 126    xk= [0.66388223 0.40558191]    f(xk)= [ 8.66399552 -7.03154158]
k= 127    xk= [0.83144109 0.69757602]    f(xk)= [-2.42627201  1.25634529]
k= 128    xk= [0.78711785 0.62687425]    f(xk)= [-2.73036387  1.46394824]
k= 129    xk= [1.21594817 1.31516318]    f(xk)= [ 79.89010983 -32.67335542]
k= 130    xk= [0.80183352 0.64738549]    f(xk)= [-1.82311642  0.88970056]
k= 131    xk= [0.81130929 0.66135898]    f(xk)= [-1.39515741  0.62724291]
k= 132    xk= [0.84379242 0.71145948]    f(xk)= [-0.13482093 -0.10523573]
k= 133    xk= [0.8488667  0.72108011]    f(xk)= [-0.47388572  0.1010872 ]
k= 134    xk= [0.83939571 0.70529745]    f(xk)= [-0.56036629  0.14245826]
k= 135    xk= [0.90713661 0.81941977]    f(xk)= [ 1.07593817 -0.69541067]
k= 136    xk= [0.86501385 0.74991553]    f(xk)= [-0.84661343  0.33331324]
k= 137    xk= [0.89854797 0.80733425]    f(xk)= [-0.18342061 -0.01084162]
k= 138    xk= [0.91010896 0.82821754]    f(xk)= [-0.15037671 -0.01615486]
k= 139    xk= [0.92061081 0.84725285]    f(xk)= [-0.05883515 -0.05428094]
k= 140    xk= [0.93145142 0.8677041 ]    f(xk)= [-0.17522727  0.02046812]
k= 141    xk= [0.90995165 0.82823094]    f(xk)= [-0.25978644  0.04378789]
k= 142    xk= [0.9791429  0.95553256]    f(xk)= [ 1.20698881 -0.63765106]
k= 143    xk= [0.92368943 0.85431748]    f(xk)= [-0.56470048  0.22306163]
k= 144    xk= [0.9419506 0.8881915]    f(xk)= [-0.46295033  0.18411344]
k= 145    xk= [0.9577381  0.91746347]    f(xk)= [-0.16160146  0.04023942]
k= 146    xk= [0.96701862 0.93506488]    f(xk)= [-0.04270618 -0.01202489]
k= 147    xk= [0.97141166 0.94372517]    f(xk)= [-0.0900311   0.01691066]
k= 148    xk= [0.96094308 0.92354637]    f(xk)= [-0.12991058  0.026951  ]
k= 149    xk= [0.99570993 0.99061839]    f(xk)= [ 0.3179641  -0.16397559]
k= 150    xk= [0.97199417 0.94517136]    f(xk)= [-0.21102069  0.07973764]
k= 151    xk= [0.98787082 0.97590466]    f(xk)= [-0.03054315  0.00318098]
k= 152    xk= [0.99093975 0.98195051]    f(xk)= [-0.01372648 -0.00221709]
k= 153    xk= [0.99398549 0.98801939]    f(xk)= [-0.016895    0.00244771]
k= 154    xk= [0.95542339 0.91169344]    f(xk)= [ 0.3466759  -0.22808167]
k= 155    xk= [0.99166933 0.98353597]    f(xk)= [-0.06740292  0.02558392]
k= 156    xk= [0.9849434  0.97039446]    f(xk)= [-0.14080314  0.05619102]
k= 157    xk= [0.99794854 0.99583192]    f(xk)= [ 0.02359149 -0.01387567]
k= 158    xk= [0.99621505 0.99248275]    f(xk)= [-0.02284649  0.00766731]
k= 159    xk= [0.99710996 0.99424387]    f(xk)= [-0.0119997   0.00311882]
k= 160    xk= [0.99816951 0.99634977]    f(xk)= [-0.00661416  0.0014793 ]
k= 161    xk= [0.99951767 0.99903483]    f(xk)= [-0.00066891 -0.00014795]
k= 162    xk= [0.99970257 0.99940648]    f(xk)= [-0.001095    0.00025015]
k= 163    xk= [0.99900072 0.99800724]    f(xk)= [-0.00391391  0.00095864]
k= 164    xk= [0.99997795 0.99995574]    f(xk)= [ 1.74913071e-05 -3.08012179e-05]
k= 165    xk= [0.99997761 0.99995513]    f(xk)= [-8.70938135e-06 -1.80379225e-05]
k= 166    xk= [0.99997777 0.99995565]    f(xk)= [-9.03605856e-05  2.29477018e-05]
k= 167    xk= [0.99997745 0.99995502]    f(xk)= [-9.16328188e-05  2.32707922e-05]
k= 168    xk= [1. 1.]    f(xk)= [ 2.02059708e-07 -1.00387831e-07]

Jacobian after first solve:
[[ 116.48889476  -56.35457991]
 [-225.91025031  112.73024426]]

second solve:
k= 1    xk= [0.56260435 1.12745333]    f(xk)= [-183.36781519  162.18593492]
k= 2    xk= [1. 1.]    f(xk)= [ 9.83445546e-08 -5.00097297e-08]

Jacobian after second solve:
[[ 380.39702492 -133.25513499]
 [-329.18184311  142.82269653]]

您可以看到第一次需要 168 次迭代，第二次重复使用第一次的雅可比行列式时只需要 2 次迭代。

我发现使用“Wolfe 条件”线性搜索并使用“Broyden 的第一种方法”(而不是 Broyden 的第二种方法)往往会使求解器在这里更加稳健。