3
$\begingroup$

Objective

Seeking for help, advise why the gradient descent implementation does not work below.

Background

Working on the task below to implement the logistic regression.

enter image description here

Gradient descent

Derived the gradient descent as in the picture.

Typo fixed as in the red in the picture. The cross entropy log loss is $- \left [ylog(z) + (1-y)log(1-z) \right ]$

enter image description here

Implemented the code, however it says incorrect.

import numpy as np  


def LogitRegression(arr):
  # code goes here
  x = arr[0]
  y = arr[1]
  a = arr[2]
  b = arr[3]

  z = 1.0 / (1.0 + np.exp(-a * x - b))
  lr = 1.0

  a = a - x * (z-y)
  a = np.round(a, decimals=3)
  b = b - (z-y)
  b = np.round(b, decimals=3)
  
  return ", ".join([str(a), str(b)])

# keep this function call here 
print(LogitRegression(input()))

enter image description here

If I reverted the sign of the gradient update, it works. However, not sure why.

  # a = a - x * (z-y)
  a = a + x * (z-y)
  a = np.round(a, decimals=3)

  # b = b - (z-y)
  b = b + (z-y)     
  b = np.round(b, decimals=3)

enter image description here

$\endgroup$

2 Answers 2

2
$\begingroup$

You are missing a minus sign before your binary cross entropy loss function. The loss function you currently have becomes more negative (positive) if the predictions are worse (better), therefore if you minimize this loss function the model will change its weights in the wrong direction and start performing worse. To make the model perform better you either maximize the loss function you currently have (i.e. use gradient ascent instead of gradient descent, as you have in your second example), or you add a minus sign so that a decrease in the loss is linked to a better prediction.

$\endgroup$
1
  • $\begingroup$ Thanks a lot but please help understand. The cross entropy loss function is $-1 * [ ylog(z) + (1-y)log(1-z) ]$ (apology that there was a typo as in $-1 * [ ylog(z) - (1-y)log(1-z) ]$ beflre). Please explain if this -1 at the top is not the one? and $\frac{dL}{dz}=\frac{z−y}{z(1−z)}$ for the backpropagation from the loss function where z is the sigmoid(-ax-b)? – $\endgroup$
    – mon
    Dec 7, 2021 at 23:26
1
$\begingroup$

I think your implementation is correct and the answer provided is just wrong.

Just for reference, the below figure represents the theory / math we are using here to implement Logistic Regression with Gradient Descent:

enter image description here

Here, we have the learnable parameter vector $\theta = [b,\;a]^T$ and $m=1$ (since a singe data point), with $X=[1,\; x]$, where $1$ corresponds to the intercept (bias) term.

Just making your implementation a little modular and increasing the number of epochs to 10 (instead of 1):

def update_params(a, b, x, y, z, lr):
    a = a + lr * x * (z-y)
    a = np.round(a, decimals=3)
    b = b + lr * (z-y)
    b = np.round(b, decimals=3)
    return a, b          

def LogitRegression(arr):
  x, y, a, b = arr
  lr = 1.0
  num_epochs = 10
  #losses, preds = [], []
  for _ in range(num_epochs):    
      z = 1.0 / (1.0 + np.exp(-a * x - b))    
      bce = -y*np.log(z) -(1-y)*np.log(1-z)
      #losses.append(bce)
      #preds.append(z)
      print(bce, y, z, a, b)
      a, b = update_params(a, b, x, y, z, lr)            
  
  return ", ".join([str(a), str(b)])

LogitRegression([1,1,1,1])
# 0.12692801104297263 1 0.8807970779778823 1 1
# 0.10135698320837692 1 0.9036104015620354 1.119 1.119 # values after 1 epoch
# 0.08437500133718023 1 0.9190865327845347 1.215 1.215
# 0.0721998635352405 1 0.9303449352007099 1.296 1.296
# 0.06305834631954188 1 0.9388886913913739 1.366 1.366
# 0.05601486564909184 1 0.9455250799418752 1.427 1.427
# 0.05042252914331105 1 0.9508275872468411 1.481 1.481
# 0.04582166273506799 1 0.9552122969502131 1.53 1.53
# 0.041959389233941616 1 0.958908721799535 1.575 1.575
# 0.03871910934525996 1 0.962020893877162 1.616 1.616

If you plot the BCE loss and the predicted y (i.e., z) over iterations, you get the following figure (as expected, BCE loss is monotonically decreasing and z is getting closer to ground truth y with increasing iterations, leading to convergence):

enter image description here

Now, if you change your update_params() to the following:

def update_params(a, b, x, y, z, lr):
    a = a + lr * x * (z-y)
    a = np.round(a, decimals=3)
    b = b + lr * (z-y)
    b = np.round(b, decimals=3)
    return a, b

and call LogitRegression() with the same set of inputs:

LogitRegression([1,1,1,1])
# 0.12692801104297263 1 0.8807970779778823 1 1
# 0.15845663982299638 1 0.8534599691639768 0.881 0.881 # provided in the answer
# 0.2073277757451888 1 0.8127532055353431 0.734 0.734
# 0.28883714051459425 1 0.7491341990786479 0.547 0.547
# 0.4403300268044629 1 0.6438239068707556 0.296 0.296
# 0.7549461015956136 1 0.4700359482354282 -0.06 -0.06
# 1.4479476778575628 1 0.2350521962362353 -0.59 -0.59
# 2.774416770021533 1 0.0623858513799944 -1.355 -1.355
# 4.596141947283801 1 0.010090691161759239 -2.293 -2.293
# 6.56740642634977 1 0.0014054377957286094 -3.283 -3.283

and you will end up with the following figure if you plot (clearly this is wrong, since the loss function increases with every epoch and z goes further away from ground-truth y, leading to divergence):

enter image description here

Also, the above implementation can easily be extended to multi-dimensional data containing many data points like the following:

def VanillaLogisticRegression(x, y): # LR without regularization
    m, n = x.shape
    w = np.zeros((n+1, 1))
    X = np.hstack((np.ones(m)[:,None],x)) # include the feature corresponding to the bias term
    num_epochs = 1000 # number of epochs to run gradient descent, tune this hyperparametrer
    lr = 0.5 # learning rate, tune this hyperparameter
    losses = []
    for _ in range(num_epochs):
        y_hat = 1. / (1. + np.exp(-np.dot(X, w))) # predicted y by the LR model
        J = np.mean(-y*np.log2(y_hat) - (1-y)*np.log2(1-y_hat)) # the binary cross entropy loss function
        grad_J = np.mean((y_hat - y)*X, axis=0) # the gradient of the loss function
        w -= lr * grad_J[:, None] # the gradient descent step, update the parameter vector w
        losses.append(J)
        # test corretness of the implementation
        # loss J should monotonically decrease & y_hat should be closer to y, with increasing iterations
        # print(J)            
    return w

m, n = 1000, 5 # 1000 rows, 5 columns
# randomly generate dataset, note that y can have values as 0 and 1 only
x, y = np.random.random(m*n).reshape(m,n), np.random.randint(0,2,m).reshape(-1,1)
w = VanillaLogisticRegression(x, y)
w # learnt parameters
# array([[-0.0749518 ],
#   [ 0.28592107],
#   [ 0.15202566],
#   [-0.15020757],
#   [ 0.08147078],
#   [-0.18823631]])

If you plot the loss function value over iterations, you will get a plot like the following one, showing how it converges.

enter image description here

Finally, let's compare the above implementation with sklearn's implementation, which uses a more advanced optimization algorithm lbfgs by default, hence likely to converge much faster, but if our implementation is correct both of then should converge to the same global minima, since the loss function is convex (note that sklearn by default uses regularization, in order to have almost no regularization, we need to have the value of the input hyper-parameter $C$ very high):

from sklearn.linear_model import LogisticRegression
clf = LogisticRegression(random_state=0, C=10**12).fit(x, y)
print(clf.coef_, clf.intercept_)
# [[ 0.28633262  0.15256914 -0.14975667  0.08192404 -0.18780851]] [-0.07612282]

Compare the parameter values obtained from the above implementation and the one obtained with sklearn's implementation: they are almost equal.

Also, let's compare the predicted probabilities obtained using these two different implementations of LR (one from scratch, another one from sklearn's library function), as can be seen from the following scatterplot, they are almost identical:

pred_probs = 1 / (1 + np.exp(-X@w))
plt.scatter(pred_probs, clf.predict_proba(x)[:,1])
plt.grid()
plt.xlabel('pred prob', size=20)
plt.ylabel('pred prob (sklearn)', size=20)
plt.show()

enter image description here

Finally, let's compute the accuracies obtained, they are identical too:

print(sum((pred_probs > 0.5) == y) / len(y)) 
# [0.527]
clf.score(x, y)   
# 0.527

This also shows the correctness of the implementation.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.