Python 从头开始实现随机梯度下降。实施是否正确？
Stochastic Gradient Descent implementation in Python from scratch. is the implementation correct?
我知道这看起来与之前就同一主题提出的许多问题相似。我调查了他们中的大多数，但他们并没有完全回答我的问题。我的问题是我的梯度没有收敛到最优值，它甚至在非常低的 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
我不认为这是一个数据 issue.I 我使用的是相当简单的线性函数。
我认为是方程式，但我也仔细检查了它们，它们对我来说似乎没问题。
您的实施中有几处需要更正（大多数是出于效率原因）。当然，您可以通过简单地定义 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 对应于你的损失函数。但是，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 列表也是合乎逻辑的，因为这是对应于 W:
的列表
Costs = [J]
由于您想执行随机梯度下降，因此无需随机选择要从数据集中提取的样本。您有两个选择：要么在每个样本处更新权重，要么计算 J w.r.t 的梯度。你的体重。后者实现起来更简单一些，并且通常比前者收敛得更优雅。但是，由于您选择了前者，因此我将使用它。请注意，即使在这个版本中，您也不必随机选择样本，但我将使用与您相同的方法，因为这也应该有效。关于您的抽样，我认为最好确保您不要两次采用相同的索引。因此，您可能需要这样定义 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
运行 它 returns 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（一旦先前的修复已被made)，但不会停止，可能是因为 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)
通过更正此问题，您的代码可以正常工作。此错误也会出现在程序的开头，但只要完成两次以上的迭代就不会影响它。
我知道这看起来与之前就同一主题提出的许多问题相似。我调查了他们中的大多数，但他们并没有完全回答我的问题。我的问题是我的梯度没有收敛到最优值，它甚至在非常低的 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
我不认为这是一个数据 issue.I 我使用的是相当简单的线性函数。
我认为是方程式，但我也仔细检查了它们，它们对我来说似乎没问题。
您的实施中有几处需要更正（大多数是出于效率原因）。当然，您可以通过简单地定义 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 对应于你的损失函数。但是，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 列表也是合乎逻辑的，因为这是对应于 W:
Costs = [J]
由于您想执行随机梯度下降，因此无需随机选择要从数据集中提取的样本。您有两个选择：要么在每个样本处更新权重，要么计算 J w.r.t 的梯度。你的体重。后者实现起来更简单一些，并且通常比前者收敛得更优雅。但是，由于您选择了前者，因此我将使用它。请注意，即使在这个版本中，您也不必随机选择样本，但我将使用与您相同的方法，因为这也应该有效。关于您的抽样，我认为最好确保您不要两次采用相同的索引。因此，您可能需要这样定义 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
运行 它 returns 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（一旦先前的修复已被made)，但不会停止，可能是因为 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)
通过更正此问题，您的代码可以正常工作。此错误也会出现在程序的开头，但只要完成两次以上的迭代就不会影响它。