Example of scipy.optimize on least square linear regression

least_square.py

import numpy as np
import scipy as sp
import scipy.optimize as opt

# Generate linear regression
a = 3
b = 2
n = 1000
np.random.seed(42)
X = np.random.normal(loc=0, scale=4, size=(n,1))
y = a + b*X + np.random.randn(n,1)

# regression function
def regression(X, y):
	# define the loss function
	n, d = X.shape
	X = np.hstack([np.ones((n,1)), X])
	def loss(params):
		residual = y - (X @ params.reshape(-1,1))
		loss = np.square(residual).mean()
		gradient = (-2 * residual * X).mean(axis=0)
		return loss, gradient
	# find initial params
	params0 = np.random.randn(d+1)
	return opt.minimize(loss, params0, jac=True, method="BFGS")

# Run
reg = regression(X, y)
print(reg)