从头开始在 Python 中实现随机梯度下降. 实现是否正确

Question

我知道这似乎与之前就同一主题提出的许多问题相似。我已经调查了他们中的大多数人，但他们并没有完全回答我的问题。我的问题是我的梯度没有收敛到最优值，而是在非常低的 alpha 值下发散和振荡。

我的数据生成功能如下

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

我不认为这是一个数据问题。我正在使用一个相当简单的线性函数。

我认为这是方程式，但我也仔细检查了它们，它们对我来说似乎很好。

Answer 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可以使用size关键字参数直接生成 1 个以上的样本。另外，请注意，使用np.random.normal(1)时，您是从均值为 1 和标准 1 的正态分布中采样，而您可能希望从均值为 0 和标准 1 的正态分布中采样。因此，您应该定义：

W = np.random.normal(size=vars)

Yunit是您使用W进行的预测。根据定义，您可以通过执行以下操作来计算它：

Yunit = X @ W

这避免了嵌套的for循环。不过，您计算J的方式很奇怪。如果我没记错的话，J对应于你的损失函数。然而，假设MSE损失为J = 0.5 * sum from k=1 to len(X) of (y_k - w*x_k) ** 2，则J的公式。因此，这两个嵌套的for循环可以重写为：

Yunit = X @ W
J = 0.5 * np.sum((Y - Yunit) ** 2)

作为附带评论：以这种方式命名err可能会误导我，因为error通常是成本，而它表示这里每一步取得的进展。台词：

Weights = []
Weights.append(W)

可以改写为：

Weights = [W]

将J添加到您的Costs列表中也是合乎逻辑的，因为这是对应于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

运行它会在 61 次迭代中返回W=array([4.99956786, 2.00023614, 3.00000213, 1.00034205, 3.99963732, 1.00063196])，最终成本为 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)

通过更正此问题，您的代码可以正常工作。此错误也出现在程序的开头，但只要完成两次以上的迭代，就不会影响它。