How to use scipy least_squares to get the estimation of unknow variables

Question

I am a newbie in using scipy.optimize. I have the following function calls func. I have x and y values given as a list and need to get the estimated value of a, b and c. I could use curve_fit to get the estimation of a, b and c. However, I want to explore the possibilities of using least_squares. When I run the following code, I get the below error. It'd be great if anyone could point me to the right direction.

import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import least_squares

np.random.seed(0)

x = np.random.randint(0, 100, 100) # sample dataset for independent variables
y = np.random.randint(0, 100, 100) # sample dataset for dependent variables  

def func(x,a,b,c):
    return a*x**2 + b*x + c

def result(list_x, list_y):
      popt = curve_fit(func, list_x, list_y)
    
sol = least_squares(result,x, args=(y,),method='lm',jac='2-point',max_nfev=2000)   

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be 
safely coerced to any supported types according to the casting rule ''safe''

Answer 1

The following code uses the least_squares() routine for optimization. The most important change in comparison to your code is ensuring that func() returns a vector of residuals. I also compared the solution to the linear algebra result to ensure correctness.

import numpy as np
from scipy.optimize import curve_fit
from scipy.optimize import least_squares

np.random.seed(0)

x = np.random.randint(0, 100, 100) # sample dataset for independent variables
y = np.random.randint(0, 100, 100) # sample dataset for dependent variables

def func(theta, x, y):
    # Return residual = fit-observed
    return (theta[0]*x**2 + theta[1]*x + theta[2]) - y

# Initial parameter guess
theta0 = np.array([0.5, -0.1, 0.3])

# Compute solution providing initial guess theta0, x input, and y input
sol = least_squares(func, theta0, args=(x,y))
print(sol.x)

#------------------- OPTIONAL -------------------#
# Compare to linear algebra solution
temp = x.reshape((100,1))
X = np.hstack( (temp**2, temp, np.ones((100,1))) )
OLS = np.linalg.lstsq(X, y.reshape((100,1)), rcond=None)
print(OLS[0])