I am currently implementing logistic regression from scratch and I'm comparing my model with SKLearn's logistic regression. Since this is just an exercise, I decided to use toy data, specifically using the Palmer's Penguins dataset. Comparing my implementation's coefficients with SKLearn's I get this: (the dict is my implementation's coefficients, while the normal list is SKLearn's)

{'Adelie': array([ 2.71869553e+01, -5.10530642e+00,  8.95099131e+00,  5.24749815e-02,
        6.27858209e-03]), 'Chinstrap': array([-2.92432584e+01,  2.77557691e+00,  1.09108267e-01, -1.99461744e-01,
       -1.52392419e-02]), 'Gentoo': array([  0.68028231,  -0.15432404, -14.95178004,   0.67591593,
[[ 2.71960599e+01 -5.10686666e+00  8.95378324e+00  5.24706034e-02
 [-3.05419279e-02  2.36553040e+00 -2.89090975e-01 -2.96850695e-01
 [-2.37512840e-01 -1.37551011e+00 -1.51342172e+01  8.53889447e-01

As you can see, the One vs Rest coefficients for "Adelie" are pretty much the same for both models, but the other classes are different. I'm wondering why this is the case.

I'm aware that I have to use the paramater multi_class='ovr' as SKLearn's default parameter for multi_class is biased towards something (I forgot what bias it has). I also used penalty='none' for the SKLearn model, as I had a regularization parameter of 0. As for the optimization, I used SKLearns minimize function with method='BFGS' as SKLearn's logistic regression uses 'LBFGS' as it's default method for gradient descent. Looking at the user guide for SKLearn's Logistic Regression, it looks like I implemented the same equations they are using for logistic regression, so I don't think there are any errors here.

My question is why are the coefficients different for the 'Chinstrap' and 'Gentoo' classes? Using google, I found out that, for logistic regression, a log-loss cost function is convex, so wouldn't this mean that the coefficients for all classes for both models should be the same if one of them ended up being the same? Why would the coefficients for one class be the same but be different for other classes? I want to say that this is due to the errors in my implementation or because there was one parameter that I overlooked when using SKLearn's model, but I wasn't able to find any errors. Below is my code. Thanks for reading this!

import numpy as np
from palmerpenguins import load_penguins
from sklearn.linear_model import LogisticRegression
from scipy.special import expit
from scipy import optimize as op
from sklearn.utils import shuffle

def logisticRegression(X, y, classes):
    m, n = X.shape  # m = number of examples, n = number of features, not including bias
    theta = {}
    for currClass in classes:
        yCurr = (y == currClass)
        yCurr = yCurr.astype(int)
        theta_ini = np.zeros(n + 1)  # theta_0 is the bias variable
        theta_ini[0] = 1  # convention is that theta_0 = 1
        X_bias = np.insert(X, 0, np.ones(X.shape[0]), axis=1)  # adding column of ones for the bias variable
        # gradient descent
        currTheta = op.minimize(costFunction, theta_ini, (X_bias, yCurr), method='BFGS').x
        theta[currClass] = currTheta
    return theta

def costFunction(theta, X, y, lambda_=0.0):
    m = X.shape[0]
    yp = expit(np.matmul(X, theta))
    eps = 1e-5
    leftOp = np.sum(y * np.log(yp + eps))
    rightOp = np.sum((1 - y) * np.log(1 - yp + eps))
    # reg = np.sum(np.square(theta[1:])) * lambda_ / 2
    cost = -1*(leftOp + rightOp) / m
    return cost

def main():
    df = load_penguins()
    # Data cleaning:
    # There are columns with NaN, so we drop those
    df = df.dropna()
    # Dataset has island, sex, year, all of which are not needed so we drop these variables
    features = df.drop(columns=['island', 'sex', 'year'])
    target = features.drop(columns=['bill_length_mm', 'bill_depth_mm', 'flipper_length_mm', 'body_mass_g'])
    features = features.drop(columns='species')

    # Convert features and target to numpy arrays
    X = features.to_numpy()
    y = target.to_numpy()
    y = y.reshape(y.shape[0])
    species = np.unique(y)

    theta = logisticRegression(X, y, species)
    model =  LogisticRegression(max_iter=10000, penalty='none', multi_class='ovr').fit(X,y)
    skTheta = np.insert(model.coef_, 0, model.intercept_, axis=1)

if __name__ == "__main__":
  • $\begingroup$ You're computing $\sigma(x) = \frac{1}{1+\exp(-x)}$ and then $\log \sigma(x)$ or $\log(1 - \sigma(x))$. If you work through the algebra, you can simplify the expression greatly. This will solve the $\log 0$ problem & the bias that arises from the $\log(\sigma(x) + \epsilon)$ kludge because the simplification means won't lose numerical precision from round-tripping: $\log \exp (x)$. $\endgroup$
    – Sycorax
    Apr 30, 2023 at 23:29

1 Answer 1

  • Scikit-learn is open-source so you can check the source code yourself for the differences.
  • If you look at the implementations of statistical algorithms in high-quality software (R, scikit-learn) vs naive implementations, you will notice that they are much more complicated. They do a lot of clever tricks to improve the vanilla algorithms.
  • One obvious culprit in your code may be using hardcoded + 1e-5 to avoid problems with logarithms. The constant is arbitrary and may alone be a reason for the difference, but there are many other differences if you look at the scikit-learn code.
  • $\begingroup$ Yeah, I did reckon that Scikit-learn would have different algorithms but I was just curious why it was the same for one class but different for the others. In my eyes, if it was the same for one, then wouldn't the algorithm be the same? But yeah, I should go look at the source code more carefully. As for the hardcoded + 1e-5, I have that on there to prevent getting log(0), as my code would frequently throw out errors of trying to do something with infinity. $\endgroup$
    – Jed
    Sep 11, 2022 at 23:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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