# Hand-crafted regression overshoots despite learning rate provided from hessian of the loss

I tried to do my own linear regression but I have seen that there are differences with LinearRegression from sklearn. Indeed, this is what gives my own linear regression:

>>> beta = estimate_beta(X_train.values, y_train.values)
>>> (X_test*beta).sum(axis = 1)
1142949     777.876173
696990     1010.574350
1437036     852.938428
1404619     439.875784
1151040    1002.271470
...
848794      396.552199
291118     3137.906208
128267      541.926482
459668      752.088114
1120216     424.453391
Length: 482738, dtype: float64


Which is very different from what sklearn returns :

array([[3.24519438],
[2.12164625],
[0.98814432],
...,
[1.13965629],
[2.09961495],
[1.84152796]])


And, of course, of the true values:

array([[2.],
[1.],
[1.],
...,
[1.],
[4.],
[2.]])


How come I have results so far from the real values with my own model?

Here is my model:

import numpy as np

def predict(x_i, beta):
"""assumes that the first element of each x_i is 1"""
return np.dot(x_i, beta)

def error(x_i, y_i, beta):
return y_i - predict(x_i, beta)

def squared_error(x_i, y_i, beta):
return error(x_i, y_i, beta) ** 2

"""The gradient (with respect to beta)
correspond to the ith squared error term"""
return [-2 * x_ij * error(x_i, y_i, beta)
for x_ij in x_i]

def estimate_beta(x, y):
beta_initial = [random.random() for x in x[0]]
return minimize_stochastic(squared_error,
x, y,
beta_initial,
0.01)

def in_random_order(data):
"""generator that returns the elements of data in random order"""
indexes = [i for i, _ in enumerate(data)] # create a list of indexes
random.shuffle(indexes)                   # suffle them
for i in indexes:
yield data[i]

def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
data = zip(x, y)
theta = theta_0                           # the initial guess
alpha = alpha_0                           # initial step size
min_theta, min_value = None, float("inf") # the minimum so far
iterations_with_no_improvement = 0

# if we ever go 100 iterations with no improvement, stop
while iterations_with_no_improvement <100:
value = sum(target_fn(x_i, y_i, theta) for x_i, y_i in data)

if value < min_value:
# if we've found a new minimum, remeber it
# and go back to the original step size
min_theta, min_value = theta, value
iterations_with_no_improvement = 0
alpha = alpha_0
else:
# otherwise we're not improving, so try shinking the step_size
iterations_with_no_improvement += 1
alpha *= 0.9

# and take a gradient step for each of the data points
for x_i, y_i in in_random_order(data):
theta = np.substract(theta, alpha * gradient_i)

return min_theta


I tried with a minimal and reproducible code:

>>>x = [[1, 49, 4, 2500],[1, 12, 2, 125], [1, 35, 2, 3790], [1, 60, 4, 4500], [1, 10, 4, 5000]]
>>>spendings_customer = [2000, 0, 3600, 0, 3500]
>>>import random
>>>random.seed(0)
>>>beta = estimate_beta(x, spendings_customer)


Which returns to me:

[0.8444218515250481,
0.7579544029403025,
0.420571580830845,
0.25891675029296335]


Which is strange ...

# With a learning rate from the hessian of the loss

The problem could have been because of my learning rate, which could have been too big. The best learning rate of a cost function is strictly less than 2/λ, where λ is the largest eigenvalue of the Hessian. I want to get this learning rate of my gradient descent algorithm. So I tried to compute it using this answer:

$$H_L(w) = X^T X$$

# Hessian is X.t*X
h = np.dot(X_test.T,X_test)
from numpy import linalg as LA
w, v = LA.eig(np.array(h))
max(w)


which returns:

(74119951381184.84+0j)


So I changed alpha_0 to something less that 2/λ But the predictions still look very different from the actual results:

1142949    1883.752457
696990     1015.531962
1437036    3342.723244
1404619     397.374124
1151040    3485.172794
...
848794     2366.912144
291118     2037.073368
128267     1853.395186
459668     4395.533697
1120216     786.328257
Length: 482738, dtype: float64