Реализация стохастического градиентного спуска в Python с нуля. правильная реализация?
Я знаю, что это могло бы походить на многие вопросы, заданные ранее по той же теме. Я опросил большинство из них, но они не совсем ответили на мой вопрос. Моя проблема в том, что мой градиент не сходится к оптимуму, он довольно расходится и колеблется даже при очень низких значениях альфа.
Моя функция генерации данных ниже
X = [[float(np.random.randn(1)) for i in range(0,100)] for j in range(0,5)]
X = np.array(X).transpose()
Y = [float(0) for i in range(0,100)]
Y = 2*X[:,0] + 3*X[:,1] + 1*X[:,2] + 4*X[:,3] + 1*X[:,4] + 5
fig, ax = plt.subplots(1,5)
fig.set_size_inches(20,5)
k = 0
for j in range(0,5):
sns.scatterplot(X[:,k],Y,ax=ax[j])
k += 1
Моя реализация SGD приведена ниже
def multilinreg(X,Y,epsilon = 0.000001,alpha = 0.01,K = 20):
Xnot = [[1] for i in range(0,len(X))]
Xnot = np.array(Xnot)
X = np.append(Xnot,X, axis = 1)
vars = X.shape[1]
W = []
W = [np.random.normal(1) for i in range(vars)]
W = np.array(W)
J = 0
for i in range(len(X)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + X[i,j] * W[j]
J = J + (0.5/(len(X)))*((Y[i]-Yunit)**2)
err = 1
iter = 0
Weights = []
Weights.append(W)
Costs = []
while err > epsilon:
index = [np.random.randint(len(Y)) for i in range(K)]
Xsample, Ysample = X[index,:], Y[index]
m =len(Xsample)
Ypredsample = []
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + X[i,j] * W[j]
Ypredsample.append(Yunit)
Ypredsample = np.array(Ypredsample)
for i in range(len(Xsample)):
for j in range(vars):
gradJunit = (-1)*(Xsample[i,j]*(Ysample[i] - Ypredsample[i]))
W[j] = W[j] - alpha*gradJunit
Jnew = 0
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + Xsample[i,j]*W[j]
Jnew = Jnew + (0.5/(len(Xsample)))*((Ysample[i]-Yunit)**2)
Weights.append(W)
err = abs(float(Jnew - J))
J = Jnew
Costs.append(J)
iter += 1
if iter % 1000 == 0:
print(iter)
print(J)
Costs = np.array(Costs)
Ypred = []
for i in range(len(X)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + X[i,j] * W[j]
Ypred.append(Yunit)
Ypred = np.array(Ypred)
return Ypred, iter, Costs, W
Гиперпараметры указаны ниже
epsilon = 1*(10)**(-20)
alpha = 0.0000001
K = 50
Я не думаю, что это проблема данных. Я использую довольно простую линейную функцию.
Я думаю, что это уравнения, но я их дважды проверил, и они мне кажутся хорошими.
1 ответ
В вашей реализации необходимо исправить несколько вещей (большинство из них по соображениям эффективности). Конечно, вы выиграете время, просто определивw = np.array([5, 2, 3, 1, 4, 1]), но это не отвечает на вопрос, почему ваша реализация SGD не работает.
Прежде всего, вы определяете X при выполнении:
X = [[float(np.random.randn(1)) for i in range(0,100)] for j in range(0,5)]
X = np.array(X).transpose()
Более быстрый способ выполнить эту операцию - просто выполнить:
X = np.random.randn(100, 5)
Затем вы определяете Y:
Y = [float(0) for i in range(0,100)]
Y = 2*X[:,0] + 3*X[:,1] + 1*X[:,2] + 4*X[:,3] + 1*X[:,4] + 5
Первая инициализация Y = [float(0) for i in range(0,100)] бесполезно, так как вы мгновенно отменяете Yсо второй строкой. Более сжатый способ записи этой строки также мог бы быть таким:
Y = X @ np.array([2, 3, 1, 4, 1]) + 5
Теперь по поводу вашей реализации SGD. Линии:
Xnot = [[1] for i in range(0,len(X))]
Xnot = np.array(Xnot)
X = np.append(Xnot,X, axis = 1)
можно более эффективно переписать как:
X = np.hstack((np.ones(len(X)).reshape(-1, 1), X))
Аналогично линии
W = []
W = [np.random.normal(1) for i in range(vars)]
W = np.array(W)
можно переписать с помощью numpyфункции. Обратите внимание, чем первая строкаW = [] бесполезно, поскольку вы отменяете W сразу после этого, не используя его. np.random.normal может напрямую сгенерировать более 1 выборки, используя sizeаргумент ключевого слова. Кроме того, обратите внимание, что при использованииnp.random.normal(1), вы выбираете из нормального распределения со средним значением 1 и стандартным 1, в то время как вы, вероятно, хотите выбрать из нормального распределения со средним значением 0 и стандартным 1. Следовательно, вы должны определить:
W = np.random.normal(size=vars)
Yunit прогноз, который вы делаете, используя W. По определению, вы можете вычислить это, выполнив следующие действия:
Yunit = X @ W
что позволяет избежать вложенных forпетли. Как вы вычисляетеJхотя странно. Если я не ошибаюсь,Jсоответствует вашей функции потерь. Однако формула дляJ, предполагая, что потеря MSE равна J = 0.5 * sum from k=1 to len(X) of (y_k - w*x_k) ** 2. Следовательно, эти два вложенныхfor циклы можно переписать как:
Yunit = X @ W
J = 0.5 * np.sum((Y - Yunit) ** 2)
В качестве дополнительного замечания: именование err таким образом может меня ввести в заблуждение, так как error- это обычно стоимость, в то время как она обозначает прогресс, достигнутый здесь на каждом этапе. Линии:
Weights = []
Weights.append(W)
можно переписать как:
Weights = [W]
Было бы логично также добавить J на ваш Costs list, так как это тот, который соответствует W:
Costs = [J]
Поскольку вы хотите выполнить стохастический градиентный спуск, нет необходимости выбирать произвольно, какие образцы вы хотите взять из своего набора данных. У вас есть два варианта: либо вы обновляете свои веса для каждого образца, либо вы можете вычислить градиентJпротив ваших весов. Последний вариант немного проще в реализации и обычно сходится более изящно, чем первый. Однако, поскольку вы выбрали первое, я буду работать именно с ним. Обратите внимание, что даже в этой версии вам не нужно выбирать образцы случайным образом, но я буду использовать тот же метод, что и вы, поскольку он также должен работать. Что касается вашей выборки, я думаю, что лучше убедиться, что вы не берете один и тот же индекс дважды. Следовательно, вы можете определитьindex как это:
index = np.random.choice(np.arange(len(Y)), size=K, replace=False)
m бесполезно, так как всегда будет равно Kв этом случае. Вы должны сохранить его, если вы выполняете выборку, не гарантируя, что у вас нет одного и того же индекса дважды. Если вы хотите выполнить выборку без проверки, что вы дважды выбрали один и тот же индекс, просто введитеreplace=True в choice функция.
Еще раз, вы можете использовать умножение матриц для вычисления Yunitболее эффективно. Следовательно, вы можете заменить:
Ypredsample = []
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + X[i,j] * W[j]
Ypredsample.append(Yunit)
по:
Ypredsample = X @ W
Точно так же вы можете вычислить обновление своего веса, используя numpyфункции. Следовательно, вы можете заменить:
for i in range(len(Xsample)):
for j in range(vars):
gradJunit = (-1)*(Xsample[i,j]*(Ysample[i] - Ypredsample[i]))
W[j] = W[j] - alpha*gradJunit
по:
W -= alpha * np.sum((Ypredsample - Ysample).reshape(-1, 1) * Xsample, axis=0)
Как и раньше, вычисление затрат может быть выполнено с помощью умножения матриц. Однако обратите внимание, что вы должны вычислитьJпо всему набору данных. Следовательно, вам следует заменить:
Jnew = 0
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + Xsample[i,j]*W[j]
Jnew = Jnew + (0.5/(len(Xsample)))*((Ysample[i]-Yunit)**2)
по:
Jnew = 0.5 * np.sum((Y - X @ W) ** 2)
Наконец, вы можете использовать умножение матриц, чтобы делать свои прогнозы. Следовательно, ваш окончательный код должен выглядеть так:
import numpy as np
X = np.random.randn(100, 5)
Y = X @ np.array([2, 3, 1, 4, 1]) + 5
def multilinreg(X, Y, epsilon=0.00001, alpha=0.01, K=20):
X = np.hstack((np.ones(len(X)).reshape(-1, 1), X))
vars = X.shape[1]
W = np.random.normal(size=vars)
Yunit = X @ W
J = 0.5 * np.sum((Y - Yunit) ** 2)
err = 1
Weights = [W]
Costs = [J]
iter = 0
while err > epsilon:
index = np.random.choice(np.arange(len(Y)), size=K, replace=False)
Xsample, Ysample = X[index], Y[index]
Ypredsample = Xsample @ W
W -= alpha * np.sum((Ypredsample - Ysample).reshape(-1,1) * Xsample, axis=0)
Jnew = 0.5 * np.sum((Y - X @ W) ** 2)
Weights.append(Jnew)
err = abs(Jnew - J)
J = Jnew
Costs.append(J)
iter += 1
if iter % 10 == 0:
print(iter)
print(J)
Costs = np.array(Costs)
Ypred = X @ W
return Ypred, iter, Costs, W
Запуск возвращается W=array([4.99956786, 2.00023614, 3.00000213, 1.00034205, 3.99963732, 1.00063196]) за 61 итерацию с финальной стоимостью 3.05e-05.
Теперь, когда мы знаем, что этот код верен, мы можем использовать его, чтобы определить, где у вас произошла ошибка. В этом фрагменте кода:
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + X[i,j] * W[j]
Ypredsample.append(Yunit)
Ypredsample = np.array(Ypredsample)
ты используешь X[i, j] вместо того Xsample[i, j], что не имеет смысла. Плюс, если вы печатаетеW вместе с J а также iter в вашем цикле вы можете видеть, что программа находит правильный W довольно быстро (после внесения предыдущего исправления), но не останавливается, вероятно, потому что Jвычисляется неправильно. Ошибка в том, что эта строка:
Jnew = Jnew + (0.5/(len(Xsample)))*((Ysample[i]-Yunit)**2)
с неправильным отступом. В самом деле, это не должно быть частьюfor j in range(vars) цикл, но он должен быть частью for i in range(len(Xsample)) только цикл, например:
Jnew = 0
for i in range(len(Xsample)):
Yunit = 0
for j in range(vars):
Yunit = Yunit + Xsample[i,j]*W[j]
Jnew = Jnew + (0.5/(len(Xsample)))*((Ysample[i]-Yunit)**2)
Исправив это, ваш код будет работать правильно. Эта ошибка также присутствует в начале вашей программы, но не влияет на нее, если выполняется более двух итераций.