python+numpy利用梯度下降算法实现逻辑回归（对比sklearn官方demo）

概要描述

在逻辑回归中，可以使用梯度下降算法来求解模型（即确定sigmoid函数中的各系数和截距）。
本文用python实现梯度下降算法，并使用numpy做矩阵计算求解逻辑回归模型，最后与sklearn官方代码结果验证结果正确性。

详细说明

python实现梯度下降算法

# !/usr/bin/python
# -*- coding: utf-8 -*-
"""
PROJECT_NAME = Datawhale
Author : sciengineer
Email : 821072960@qq.com
Time = 2021/7/22 18:09
"""


import numpy as np

# General a toy dataset:s it's just a straight line with some Gaussian noise:
xmin, xmax = -5, 5
n_samples = 100
np.random.seed(0)
X = np.random.normal(size=n_samples)
y = (X > 0).astype(np.int64)
X[X > 0] *= 4
X += .3 * np.random.normal(size=n_samples)
X = X[:, np.newaxis]


# learning rate: lr
lr = 1


# gradient descent algorithm
def grad_desc(x_ndarr, y_ndarr):


    """
    :param x_ndarr: ndarray of features (for example, 422 samples with 3 fetures, then the shape of x_ndarr
    is (422,3)
    :param y_ndarr: ndarray of label (for example, 422 samples with 1 label, then the shape of y_ndarr
    is (422,)
    :return: a list contains interception and coeffients of the logistic regression model
    """
    x0 = np.ones(x_ndarr.shape[0])
    # in order to use matrix multiplication, add this column x0
    x_ndarr = np.insert(x_ndarr, 0, x0, axis=1)
    theta_arr = np.zeros(x_ndarr.shape[1])

    # to find a reasonable number to initialize the j_loss_last, I
    # calculate the j_loss when i == 0 in the for loop (69), and
    # set j_loss_last the  same order of magnitude with j_loss.
    j_loss_last = 1e2
    delta = 1e-10

    for i in range(10 ** 10):

        # z = theta0 * x0 + theta1 * x1 +theta2 * x2+ ... +thetan * xn
        # z = np.dot(x_ndarr, theta_arr)
        z = np.dot(x_ndarr, theta_arr)
        # z = np.dot( theta_arr, x_ndarr.T)
        #
        y_hat = 1/(1+np.exp(-z))
        # use the clip to bound the y_hat in (0,1) , and to avoid warning "divide by zero encountered in log" or
        # "invalid value encountered in log". Otherwise this will ruin the optimization.
        y_hat = np.clip(y_hat,delta,1-delta)
        # j_loss = np.dot((y_hat - y_ndarr).T, (y_hat - y_ndarr))
        j_loss = -np.dot(y_ndarr , np.log(y_hat))  - np.dot((1-y_ndarr),np.log(1-y_hat))
        delta_j_loss = j_loss_last - j_loss
        rate = abs(delta_j_loss / j_loss)
        # partial derivative of function j_loss with respect to variable theta_arr
        pd_j2theta_arr = np.dot(y_hat - y_ndarr, x_ndarr)

        # theta_arr updates each interation
        theta_arr = theta_arr - lr * 0.01*pd_j2theta_arr
        j_loss_last = j_loss

        # I choose the rate as the condition of convergence

        if rate < 5 * 1e-10:
            break
    return theta_arr



theta_arr = grad_desc(X, y)
# The coeffients: theta1, theta2,..., thetan
print('Coeffients: \n', theta_arr[1:])
# The interception: theta0
print('Interception: \n', theta_arr[0])

运行结果：

Coeffients: 
 [6.90478134]
Interception: 
 -1.6481480918181262

sklearn官方demo

只选取这个官方demo中的前面部分，主要是为了验证手写的逻辑回归算法实现的正确性。

# !/usr/bin/python
# -*- coding: utf-8 -*-
# Code source: Gael Varoquaux
# License: BSD 3 clause

import numpy as np
import matplotlib.pyplot as plt

from sklearn import linear_model
from scipy.special import expit

# General a toy dataset:s it's just a straight line with some Gaussian noise:
xmin, xmax = -5, 5
n_samples = 100
np.random.seed(0)
X = np.random.normal(size=n_samples)
y = (X > 0).astype(float)
X[X > 0] *= 4
X += .3 * np.random.normal(size=n_samples)

X = X[:, np.newaxis]

# Fit the classifier
clf = linear_model.LogisticRegression(C=1e10)
clf.fit(X, y)

print('Coeffients: \n', clf.coef_)
print('Interception: \n',clf.intercept_)

运行结果：

Coeffients: 
 [[6.90879439]]
Interception: 
 [-1.64913083]

可以看到，使用python实现的逻辑回归与sklearn库的运行结果一致。当然，这个实现肯定有许多需要优化的地方，欢迎讨论。