Три решателя для линейной регрессии: решение наименьших квадратов, градиентный спуск и стохастический градиентный спуск. Градиентный спуск и стохастический градиентный спуск — отличные решения для многих задач машинного обучения. Это упражнение предназначено для того, чтобы понять, как эти два метода помогают решать коэффициенты в линейной регрессии, и сравнить результат с решением наименьших квадратов.

Данные, которые я использую здесь, представляют собой общедоступный набор данных от Kaggle-https://www.kaggle.com/datasets/karthickveerakumar/startup-logistic-regression.

Первое чтение в наборе данных:

import pandas as pd
df = pd.read_csv('/Users/keru/Downloads/50_Startups.csv')
df.head()

Есть 4 переменные — расходы на НИОКР, администрирование, расходы на маркетинг, состояние. Целевой переменной является прибыль.

Чтобы упростить эту задачу, я использую Прибыль как зависимую переменную, Расходы на НИОКР как независимую переменную и провожу простую линейную регрессию:

y = a + b*x.

Быстрая проверка линейности между прибылью и затратами на НИОКР:

Линейная регрессия

Модель между прибылью и расходами на НИОКР представляет собой линейную регрессию.

где у — прибыль, а х — затраты на НИОКР.

Необходимо оценить два параметра — b0 и b1.

Качество этой подгонки измеряется квадратом остаточной суммы, который представляет собой сумму квадрата ошибки между истинным y и прогнозируемым y.

Цель состоит в том, чтобы найти подобранную линию, которая может минимизировать этот член ошибки.

Мы могли бы построить RSS как ответ на b0 и b1

def rss(b0, b1, df):
    rss = 0
    for i in range(0, df.shape[0]):
        rss += (df.Profit.values[i] - (b0 * df['R&D Spend'].values[i] + b1)) ** 2
    return rss
grid_b0 = np.linspace(0.5, 1.2, 50)
grid_b1 = np.linspace(45000, 55000, 50)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
grid_b0, grid_b1 = np.meshgrid(grid_b0, grid_b1)
zs = np.array([rss(x,y, df) for x,y in zip(np.ravel(grid_b0), np.ravel(grid_b1))])
Z = zs.reshape(grid_b0.shape)
ax.plot_surface(grid_b0, grid_b1, Z, rstride=1, cstride=1, color='red', alpha = 0.5)
ax.set_xlabel('b1');ax.set_ylabel('b0');ax.set_zlabel('RSS')              
plt.show()

Оси x и y — это b1 и b0, а ось z — это RSS.

Из этого графика видно, что если оба b1 и b0 спустятся с холма до конца, они достигнут минимума RSS=0. Направление спуска можно определить по производной RSS по b0 и производной RSS по b1:

Этот минимум RSS находится у подножия холма, поскольку производная RSS, соответствующая b1 и b0, равна 0.

Наименьший квадрат (LS)

Это решение близкой формы, где мы устанавливаем производную RSS по b1 и производную RSS по b0 равными 0.

# least square solution 
y = df.Profit.values
x = df['R&D Spend'].values
n = df.shape[0]
# normalization y and x 
y_normalizer = np.sqrt(sum(y**2))
y_norm = y/y_normalizer
x_normalizer = np.sqrt(sum(x**2))
x_norm = x/x_normalizer
# least square solution 
b1 = (sum(y_norm*x_norm) - sum(y_norm)*sum(x_norm)/n)/(sum([t**2 for t in x_norm]) - sum(x_norm)*sum(x_norm)/n)
print(b1)
b0 = (sum(y_norm) - b1*sum(x_norm))/n
print(b0)

b1 = 0,622 и b0 = 0,058

Обратите внимание, что здесь я нормализовал y и x, используя длину их векторов. Масштабирование признаков важно для градиентного спуска. Я нормализовал x и y здесь, чтобы получить b0 и b1 для решения наименьших квадратов, чтобы я мог сравнить их с градиентным спуском.

Градиентный спуск

Алгоритм градиентного спуска:

Это вспомогательная функция для обновления b0 и b1, а также для вывода обоих градиентов — g(b0) и g(b1).

# gradient descent 
def gradient_descent(b0, b1, stepsize): 
    # calculate gradient descent 
    g_b0 = -2*sum(y_norm-(b0+b1*x_norm))
    g_b1 = -2*sum([(y_i-(b0+b1*x_i))*x_i for x_i, y_i in zip(x_norm, y_norm)])
    # update both b0 and b1
    b1_new = b1 - stepsize*g_b1
    b0_new = b0 - stepsize*g_b0
    return b0_new, b1_new, g_b0, g_b1

Теперь я начинаю итерации:

gradient_vector_length = []
# inital stepsize, b0, b1, gradient b0 and gradient b1 at iteration 0 
stepsize = 0.01
b0=0; b1=0
g_b0 = 1; g_b1 = 1
# while not converge (tolerance is 0.0001)
while np.sqrt(g_b0**2 + g_b1**2) > 0.0001:
    # get b0 and b1 at t+1
    b0, b1, g_b0, g_b1 = gradient_descent(b0, b1, stepsize)
    # collect gradient_vector_length for plot
    gradient_vector_length.append(np.sqrt(g_b0**2 + g_b1**2))
    
plt.plot(gradient_vector_length) 
plt.title('Gradient Vector Length')

Первый размер шага, который я пробовал, равен 0,1, и это размер вектора градиента:

Он никогда не падает! И в конце концов он взорвался. Размер этого шага слишком велик.

Попробуйте размер шага = 0,01

Он сходился довольно быстро с шагом 0,01.

b0 и b1 после покрытия:

b0 = 0.058
b1 = 0.622

Стохастический градиентный спуск (SGD)

В SGD градиент вычисляется с использованием только одной выборки вместо всего набора данных при градиентном спуске.

Алгоритм очень похож на градиентный спуск, за исключением того, что мы перебираем каждую выборку в наборе данных (b0, b1) каждый раз обновляются:

Одно важное замечание здесь заключается в том, что данные должны быть перемешаны в первую очередь. Поскольку обновление производится каждой выборкой, здесь важна случайность выборки.

import sklearn
df_shuffled=sklearn.utils.shuffle(df)
y = df_shuffled.Profit.values
x = df_shuffled['R&D Spend'].values
# normalize 
y_normalizer = np.sqrt(sum(y**2))
y_norm = y/y_normalizer
x_normalizer = np.sqrt(sum(x**2))
x_norm = x/x_normalizer

Я попробовал 5 разных размеров шага от 0,5 до 0,001.

stepsize = [0.5, 0.1, 0.05, 0.01, 0.001]
# save gradient vector length for each stepsize
gradient_vector_length_dict = dict()
b1_dict = dict()
b0_dict = dict()
# save number of pass it went over the dataframe
number_of_pass_dict = dict()

Теперь для каждого размера шага проведите SGD, сохраните коэффициенты (b0, b1) и размер вектора градиента. Обратите внимание, что первые 200 сгенерированных коэффициентов здесь отбрасываются, а в качестве предсказанных коэффициентов используется среднее значение остальных. Это сделано для минимизации шума, поскольку последние сгенерированные коэффициенты могут быть либо очень хорошими, либо плохими, в зависимости от конкретного образца.

for stepsize_i in stepsize:
    
    gradient_vector_length = []
    b0_sgd = []
    b1_sgd = []
    
    # inital stepsize, b0, b1, gradient b0 / b1 at iteration 0 
    stepsize = stepsize_i
    b0=0; b1=0
    g_b0 = 1; g_b1 = 1
    
    # while not converge (tolerance is 0.0001)
    satisfied = False
    tol = 0.0001
    number_of_pass = 0
    while satisfied == False:
        # if go over the entire dataframe, then add one
        number_of_pass +=1
        for i in range(df.shape[0]):
            xi = x_norm[i]
            yi = y_norm[i]    
            # get b0 and b1 at t+1
            b0, b1, g_b0, g_b1 = stochastic_gradient_descent(b0, b1, stepsize, yi, xi)
            b0_sgd.append(b0)
            b1_sgd.append(b1)
            # collect gradient_vector_length for plot
            gradient_vector_length.append(np.sqrt(g_b0**2 + g_b1**2))
            # check is gradient vector size is smaller than tolerance 
            satisfied = np.sqrt(g_b0**2 + g_b1**2) < tol
            if satisfied:
                print('gradient vector size: ', np.sqrt(g_b0**2 + g_b1**2))
                break
                
    graident_vector_length_dict[str(stepsize)] = gradient_vector_length
    # discard first 200 generated b0 and b1 
    b1_dict[str(stepsize)] = np.mean(b1_sgd[200:])
    b0_dict[str(stepsize)] = np.mean(b0_sgd[200:])
    number_of_pass_dict[str(stepsize)] = number_of_pass

Затем постройте размер вектора градиента для каждого размера шага:

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1,5, figsize=(20,5))
for i, (key, value) in enumerate(graident_vector_length_dict.items()):
    axes[i].plot(value)
    axes[i].set(title=key.upper(), xlabel='Gradient Vector Length')
plt.show()

Размер шага 0,5 сходился намного лучше, чем другие 4 размера шага здесь, которые имеют наименьшее колебание ближе к концу.

По количеству проходов видно, что:

{'0.5': 10, '0.1': 5, '0.05': 29, '0.01': 30, '0.001': 20}

шаг 0,5 прошел через весь набор данных 10 раз.

Два коэффициента (b0 и b1) при размере шага 0,5:

b0: 0.065
b1: 0.565

Сравните результаты по 3 методам

С очень низким допуском — 0,0001, градиентный спуск очень хорошо соответствует наименьшим квадратам, в то время как между SGD и двумя другими есть небольшие различия.

Решение методом наименьших квадратов всегда будет находить наилучшую подогнанную линию, поскольку градиент = 0. Однако для многих задач машинного обучения размер объекта может быть очень большим, тогда решение градиента = 0 может оказаться неосуществимым. Например, если N ‹ D, то мы не можем решить градиент = 0. Даже если мы можем это решить, градиентный спуск в вычислительном отношении более эффективен, когда размер объекта действительно высок.

Однако градиентный спуск может быть очень медленным, когда набор данных большой, потому что весь набор данных используется для вычисления RSS для обновления коэффициентов. Представьте, что для вычисления RSS для одной выборки требуется 1 миллисекунда, при 1 миллиарде выборок потребуется ~ 11 дней, чтобы сделать один шаг для градиентного спуска, тогда как для SGD потребуется всего 1 миллисекунда.

SGD отлично подходит для скорости, но он чувствителен к параметрам, таким как допуск и размер шага в этом посте. Выбор правильного размера шага для SGD потребует много испытаний.

Вот несколько советов:

  1. Попробуйте сначала найти большие и маленькие размеры шагов, а промежуточные — экспоненциально разнесенные между ними.
  2. Уменьшите размер шага по мере увеличения итерации:
stepsize = stepsize / i, where i is i_th iteration

На этом все, надеюсь, вам понравилось :)