I am using lmfit to look for the parameters that optimize the fit of a model to molecular spectra. The residual program invokes a Fortran code that computes the energy level and matches them with the available experimental data, providing the set of residuals that are the input to lmfit-minimize.
There are two possible cases: linear and non-linear molecules. I prepare one example for each of them. With linear molecules the output is what I would expect, reaching the same minimum than a previous code of mine that uses the Fortran version of Minuit for parameter optimization. This is an example, where that values printed are the Chi square values obtained each time the function is called.
import numpy as np
from lmfit import minimize, Parameters, fit_report
from subprocess import Popen, PIPE
#
%run ../../bin/residuals_U3.py
input_filename = "HNC_input_file"
BENT = ".F."
exp_data_file="exp_HNC_Danielle.dat"
N_val=60
LMAX=7
VMAX=7
EMINL=".F."
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=2390.0, vary=True)
U3H_FIT_param.add('P21',value=-37.0,vary=True)
U3H_FIT_param.add('P22',value=21.0, vary=True)
U3H_FIT_param.add('P23',value=-9.0, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT,
exp_data_file, N_val, LMAX, VMAX, EMINL))
31298390.0981
31298390.0981
31298390.0981
31298405.9262
31298393.7425
31298390.2979
31298399.7958
1239848240.71
....
5.7028931985
5.70223335575
5.68090132895
5.67973533085
5.68050720657
5.68089119394
5.68021995958
5.57489186835
5.57489701562
5.57489288067
5.57489188595
5.57489412245
5.57489158478
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 169
# data points = 19
# variables = 4
chi-square = 5.575
reduced chi-square = 0.372
Akaike info crit = -15.297
Bayesian info crit = -11.519
[[Variables]]
P11: 2375.84727 +/- 59.96417 (2.52%) (init= 2390)
P21: -37.5581932 +/- 1.065396 (2.84%) (init=-37)
P22: 19.2964371 +/- 0.426936 (2.21%) (init= 21)
P23: -9.63266677 +/- 0.261630 (2.72%) (init=-9)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P11, P23) = -1.000
C(P11, P21) = -1.000
C(P21, P23) = 1.000
C(P22, P23) = -1.000
C(P11, P22) = 1.000
C(P21, P22) = -1.000 code here
So far so good. The problem arises once I tried to fit data from a non-linear molecular species. Then, after very few evaluations, the program seems to converge to the initial values provided. If we call this example 1:
import numpy as np
from lmfit import minimize, Parameters, fit_report
from subprocess import Popen, PIPE
%run ../../bin/residuals_U3.py
input_filename = "NCNCS_input_file"
BENT = ".T."
exp_data_file="test_exp_NCNCS.dat"
N_val=70
LMAX=9
VMAX=6
EMINL=".F."
# Parameters far from the minimum I already know
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
33820.0608635
33820.0608635
33820.0608635
33942.2383744
34935.7037896
33810.3126988
33357.9527561
33303.3818229
34574.0192135
34526.2057525
33337.2494285
34458.8982496
34798.9915284
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 13
# data points = 70
# variables = 4
chi-square = 33303.382
reduced chi-square = 504.597
Akaike info crit = 439.544
Bayesian info crit = 448.538
[[Variables]]
P11: 199.999988 +/- 4.02e-05 (0.00%) (init= 200)
P21: -2.60000161 +/- 6.71e-07 (0.00%) (init=-2.6)
P22: 1.50000025 +/- 3.57e-07 (0.00%) (init= 1.5)
P23: -0.79999987 +/- 2.12e-07 (0.00%) (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P21, P23) = -0.675
C(P21, P22) = 0.506
C(P22, P23) = -0.430
C(P11, P21) = -0.420
C(P11, P23) = -0.317
C(P11, P22) = -0.216
If we repeat example 1 with a more verbose output that includes, apart from Chi^2 as previously, the parameter values and the residuals obtained we have:
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
.....
&INP1b P11=199.999988473 /
&INP2b P21=-2.60000161305, P22=1.50000025484, P23=-0.799999875067 /
RESIDUALS [ 0. -3.99268261 -11.73264434 -19.05089687 -21.8743998
-23.98786565 -28.05424445 -0.80354366 -4.29117203 -10.41169177
-16.11475578 -20.76739082 -24.72907688 -29.96025254 -1.19554552
-4.77165975 -10.04051152 -15.40288703 -20.18713017 -25.8307694
-31.91088782 -2.46101195 -5.74047469 -9.9850673 -15.24025329
-20.72837057 -26.68827531 -33.31621124 -3.66948178 -6.53021075
-10.89700164 -15.81327816 -21.81723225 -28.08793728 -35.76616676
-4.99405722 -7.65935238 -11.55129453 -16.5529505 -22.40531267
-29.36326533 -38.30271257 -7.21713077 -8.76345313 -12.77142684
-17.71182548 -24.01827081 -31.19676872 -40.10579542 -7.61666334
-9.92207061 -13.58164103 -18.52451949 -25.03610104 -33.39256778
-43.32754575 -8.67692476 -11.12324272 -15.00726898 -20.04660062
-27.02514438 -35.41059799 -46.51995156 -10.18456526 -12.0963503
-15.60152462 -21.17712212 -28.94048104 -37.94146875 -49.83299737]
34798.9915284
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 13
# data points = 70
# variables = 4
chi-square = 33303.382
reduced chi-square = 504.597
Akaike info crit = 439.544
Bayesian info crit = 448.538
[[Variables]]
P11: 199.999988 +/- 4.02e-05 (0.00%) (init= 200)
P21: -2.60000161 +/- 6.71e-07 (0.00%) (init=-2.6)
P22: 1.50000025 +/- 3.57e-07 (0.00%) (init= 1.5)
P23: -0.79999987 +/- 2.12e-07 (0.00%) (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P21, P23) = -0.675
C(P21, P22) = 0.506
C(P22, P23) = -0.430
C(P11, P21) = -0.420
C(P11, P23) = -0.317
C(P11, P22) = -0.216
There is clearly a better minimum as the example (example 2) shows:
# Parameters closer to the minimum
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=201.0, vary=True)
U3H_FIT_param.add('P21',value=-2.570,vary=True)
U3H_FIT_param.add('P22',value=1.4980, vary=True)
U3H_FIT_param.add('P23',value=-0.81, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
866.43578979
866.43578979
866.43578979
880.37396166
890.227624754
852.194331523
780.074380341
801.333830227
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 8
# data points = 70
# variables = 4
chi-square = 801.334
reduced chi-square = 12.141
Akaike info crit = 178.645
Bayesian info crit = 187.639
[[Variables]]
P11: 200.999990 +/- 7.47e-06 (0.00%) (init= 201)
P21: -2.57000013 +/- 1.08e-07 (0.00%) (init=-2.57)
P22: 1.49800011 +/- 5.88e-08 (0.00%) (init= 1.498)
P23: -0.80999997 +/- 8.11e-09 (0.00%) (init=-0.81)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P11, P21) = -0.386
C(P11, P22) = -0.379
C(P21, P23) = 0.263
C(P11, P23) = -0.125
However, if in this bent case I change the minimization method, e.g. to nelder, the results are as expected, even if I start the minimization quite far from the minimum as in example 1:
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL), method = "nelder")
33820.0608635
386363.244872
540126.50767
27839.1627683
21180.6786213
398690.128256
105124.802101
98763.3226762
332146.050613
3785.18423466
175555.611803
45570.4112121
15015.0822576
20437.4181311
76313.1360929
3192.21680729
...
355.980240836
332.571174789
332.571174789
355.980240836
350.020618839
358.718935527
332.571174789
332.571174789
332.571174789
332.571174789
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 280
# data points = 70
# variables = 4
chi-square = 332.571
reduced chi-square = 5.039
Akaike info crit = 117.085
Bayesian info crit = 126.079
[[Variables]]
P11: 201.428112 (init= 200)
P21: -2.58733894 (init=-2.6)
P22: 1.52918377 (init= 1.5)
P23: -0.81586079 (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
In this case the more verbose output for some of the iterations is the following
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL), method = "nelder")
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=210.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -26.16300891 -35.92631045 -18.8046449 7.84306164
29.85328814 49.35183801 5.44122598 -7.20663588 -7.51532458
5.28787982 22.79715815 41.85783891 59.23158371 14.12838386
9.2311941 12.82257834 24.05402553 39.05684174 54.70325573
69.54410104 24.02558592 24.14278484 30.06592218 40.05349953
53.44068403 67.14018094 80.79972754 33.95012425 37.30408452
44.53432899 54.88610982 67.23500121 79.31052971 91.39583861
44.20244335 49.60557888 58.08586594 68.72834808 79.93361024
91.24600502 101.35283675 53.69459696 61.88966887 70.56905834
80.97053725 91.89016079 102.16450603 111.61917172 64.64821697
73.47362861 82.86652184 93.10114386 103.25734084 112.65601867
120.69284205 75.35070819 84.29816695 94.28958114 104.26041597
113.8639484 122.29150128 129.26978921 85.42971752 95.71850413
105.5315729 115.31685183 124.24099756 132.05550135 137.78892407]
386363.244872
&INP1b P11=200.0 /
&INP2b P21=-2.73, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -6.51482403 -11.84722003 -17.94556936 -34.33607138
-57.67632401 -81.40734149 -1.63845707 -8.55295725 -17.17160577
-29.53548333 -47.97384287 -69.24130491 -94.24014416 -4.42106668
-13.37478867 -25.1140549 -40.25940819 -59.31786691 -81.45795179
-106.40406773 -8.47241409 -19.88180645 -33.26697111 -49.98215925
-70.28802414 -93.70400896 -119.18639263 -14.91219834 -26.42593577
-41.64563448 -60.30755635 -81.61793552 -105.65083382 -132.51508599
-20.93474386 -34.42278364 -50.99246305 -70.22964079 -92.3174835
-117.92944692 -146.19607146 -29.37468772 -42.38500645 -59.95420309
-80.42758498 -103.61555313 -130.08637485 -159.51597624 -35.8953851
-51.09068006 -69.33552966 -90.56701074 -114.87032877 -142.75145306
-173.6936809 -43.62317666 -59.85538811 -79.26959621 -101.53354992
-126.78686941 -155.31065489 -187.76051864 -52.18050662 -68.57148933
-88.53569045 -111.87500818 -138.66419476 -168.68410045 -203.08199158]
540126.50767
........................
&INP1b P11=201.428112052 /
&INP2b P21=-2.58733894638, P22=1.52918377018, P23=-0.815860796044 /
RESIDUALS [ 0. 2.70481025 2.5647914 -0.66779985 -1.84848475 -2.294670
-2.8511186 -0.06763408 2.30399185 2.4069086 0.60429218 -1.17025016
-1.7678735 -2.25363581 -0.97533799 1.89312414 2.19413958 1.4505511
-0.1054865 -1.2606611 -2.4558844 -1.59327497 0.99634674 1.79137313
1.42330317 0.81319402 -0.59261343 -1.9298181 -2.60320618 0.49113672
1.90950601 1.62687398 1.43875694 0.03109259 -2.37265779 -3.45057765
0.19965851 1.72320248 2.2898178 1.34013775 0.40188506 -2.25368524
-4.79999944 -0.29638363 1.49036912 2.26290378 1.88161992 0.70066741
-2.00055299 -4.2727996 -0.37954819 2.06048067 3.136022 2.91270888
0.68341221 -2.54241374 -4.11982953 -0.25818249 2.16914822 3.46019775
2.91502288 1.04311452 -2.62395441 -3.82618485 0.00951645 2.95514821
3.92300846 3.37200094 1.280833 -3.34939407]
332.571174789
&INP1b P11=201.428112052 /
&INP2b P21=-2.58733894638, P22=1.52918377018, P23=-0.815860796044 /
RESIDUALS [ 0. 2.70481025 2.5647914 -0.66779985 -1.84848475 -2.29467012
-2.8511186 -0.06763408 2.30399185 2.4069086 0.60429218 -1.17025016
-1.7678735 -2.25363581 -0.97533799 1.89312414 2.19413958 1.4505511
-0.1054865 -1.2606611 -2.4558844 -1.59327497 0.99634674 1.79137313
1.42330317 0.81319402 -0.59261343 -1.9298181 -2.60320618 0.49113672
1.90950601 1.62687398 1.43875694 0.03109259 -2.37265779 -3.45057765
0.19965851 1.72320248 2.2898178 1.34013775 0.40188506 -2.25368524
-4.79999944 -0.29638363 1.49036912 2.26290378 1.88161992 0.70066741
-2.00055299 -4.2727996 -0.37954819 2.06048067 3.136022 2.91270888
0.68341221 -2.54241374 -4.11982953 -0.25818249 2.16914822 3.46019775
2.91502288 1.04311452 -2.62395441 -3.82618485 0.00951645 2.95514821
3.92300846 3.37200094 1.280833 -3.34939407]
332.571174789
I would appreciate any explanation or guess of why the default (Levenberg-Marquardt method) may be failing or if there is a way of tuning it to overcome this problem.
