Scipy Minimize andf Least Square results on parameter covariance matrix

Question

I experience a problem using scipy.optimize.minimize to estimate parameter covariance matrix. Here is a small exercise snippet:

import numpy as np

# Dataset
N=20
rng = np.random.default_rng(2022)
ti = 10.0 * rng.random(N)
ti = np.sort(ti)
sigma_e = 1.
e = rng.normal(0, sigma_e, ti.shape)
param_true = np.array([3.5, 1.0])
yi =param_true[1] + param_true[0]*ti +e

#Least SQ estmation
def test(params, X, y):
    X = jnp.c_[ X, np.ones(len(X)) ] 
    residuals = jnp.dot(X, params) - y
    return residuals

res_lsq=scipy.optimize.leastsq(test, jnp.array([0.,0.]), args=(ti,yi), 
Dfun=None, full_output=True, col_deriv=0, ftol=1.49012e-08, 
xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

print(res_lsq[0]) # parameters
print(res_lsq[1]) # covariance mtx

I get

[3.51045968 0.9103981 ]

[[ 0.00432748 -0.02341142]
 [-0.02341142  0.17665437]]

that I have validated by an other method. Now using minimize L-BFGS-B

# Minimize 
def lik(params, X, y):
    X = np.c_[ X, np.ones(len(X)) ] 
    residuals = np.dot(X, params) - y
    return np.mean(residuals ** 2)

lik_model = scipy.optimize.minimize(lik, jnp.array([0.,0.]), args=(ti,yi), method='L-BFGS-B',
                     options={'gtol': 1e-6,'disp': False})

print(lik_model.x)
print(lik_model.hess_inv.todense())

I get

[3.51045967 0.91039815]

[[ 0.04327373 -0.23411416]
 [-0.23411416  1.76654345]]

As you may probably notice, the parameters are very close to the leastsq method but the Inverse Hessian coefficients are 10 times bigger.

Does someone can explain me this feature? or is it a bug? or is there conditionning factor of the Hessian that is internally used that I have missed?

Thanks

user_11635654 · Answer

I found how to proceed and it may useful for other people, and at first using sigma_e=1 is masking part of the problem. So, in the above exercise you can set sigma_e = 3.0, and proceed to the replacement

def lik(params, X, y):
    X = np.c_[ X, np.ones(len(X)) ] 
    residuals = np.dot(X, params) - y
    return 0.5*np.sum((residuals/sigma_e) ** 2)

to be used with

lik_model = minimize(lik, jnp.array([0.,0.]), args=(ti,yi), method='BFGS', options={'gtol': 1e-6,'disp': False})

And

def test(params, X, y):
    X = jnp.c_[ X, np.ones(len(X)) ] 
    residuals = jnp.dot(X, params) - y
    return residuals/sigma_e

if you use

res_lsq=scipy.optimize.leastsq(test, jnp.array([0.,0.]), args=(ti,yi), 
Dfun=None, full_output=True, col_deriv=0, ftol=1.49012e-08, 
xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

Then, res_lsq[1] and lik_model.hess_inv will gives the same results up to numerical uncertainty. (in case you use L-BFGS-B then you get the inverse hessian by lik_model.hess_inv.todense(). I must confess that I have made 2 mistakes before 1) not scaling by the residuals by sigma_e, 2) the use of np.mean() instead of np.sum(). But, I have not yet figure out the need of the 0.5 factor in the minimize case, although it is an indication of the normal distribution factor, but I do not manage to get it in the code.

Ok, so far so good now I hope you will get the right covariance matrix in your own use-cases.

Scipy Minimize andf Least Square results on parameter covariance matrix

1 Answers