从头开始在Python中进行随机梯度下降实现实施正确吗？

问题描述

我知道这似乎与之前在同一主题上提出的许多问题相似。我已经调查了其中大多数，但他们并未完全回答我的问题。我的问题是我的梯度没有收敛到最优值，而是在非常低的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,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

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

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

解决方法

在您的实现中有几件事需要更正（大多数出于效率考虑）。当然，您只需定义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,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关键字参数直接生成多个样本。另外，请注意，在使用np.random.normal(1)时，您是从均值1和std 1的正态分布中采样的，而您可能想从均值0和std 1的正态分布中采样。因此，您应该定义：

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

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

    Yunit = X @ W

避免嵌套的for循环。尽管您计算J的方式很奇怪。如果我没记错的话，J对应于您的损失函数。但是，假设MSE损失为J，则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,1]) + 5

def multilinreg(X,Y,epsilon=0.00001,alpha=0.01,K=20):
    X = np.hstack((np.ones(len(X)).reshape(-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)),replace=False)
        Xsample,Ysample = X[index],Y[index]
        Ypredsample = Xsample @ W
        W -= alpha * np.sum((Ypredsample - Ysample).reshape(-1,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)

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