【强基固本】基础算法：使用numpy实现逻辑回归随机梯度下降（附代码）

深度学习算法工程师面试，记录一道较为基础的笔试题：

输入：目标向量Y(N*1），矩阵X(N*K)；输出：使用随机梯度下降求得的逻辑回归系数W(K+1)。

分析：该问题需要先列出逻辑回归的函数解析式，再选择损失函数，最后算出损失函数关于更新参数的导数，即可开始随机梯度下降。

作者：Algernon

地址：https://www.zhihu.com/people/thisiszhou

01

逻辑回归解析式

其中

02

Loss函数

由于逻辑回归的输出值在[0, 1]之间，并且逻辑回归虽然名字为回归，但实际是分类任务，所以损失函数使用交叉熵。其中交叉熵函数解析式为：

y为label，y‘为预测值，N为y的长度。

03

关于更新参数的导数

更新参数有  ，b，Loss的解析式为：

为了便于求导，我们换一种写法，其中sum函数相当于给后面的列向量乘一个都为1的行向量S:

其中  为矩阵乘法，CE和  都是逐元素函数，则：

其中

（其中  为逐元素乘法，在最后的计算中，逐元素乘法优先级比矩阵乘法优先级高）

04

参数更新（梯度下降）

其中  为学习率。

05

代码实现

先封装一下涉及到的函数：

import numpy as np
from random import shuffle

# sigmoid函数
def sigmoid(x):
    return 1 / (1 + np.exp(-x))


# 交叉熵函数
def ce(y, y_label):
    return -y_label*np.log(y) - (1-y_label)*np.log(1-y)


# 交叉熵函数导数
def ce_grad(y, y_label):
    return -y_label/y + (1-y_label)/(1-y)

初始化参数以及输入：

# 确定随机性唯一
np.random.seed(32)  
n = 6
k = 7
# 随机初始化输入与参数
X = np.random.rand(n,k) 
y_label = np.random.randint(0,2,size=(n,1)).astype(np.float32)
w = np.random.rand(k,1)
b = np.random.rand(1,1)

前向函数：

def forward(X, w, b):
    y1 = np.dot(X, w) + b
    y2 = sigmoid(y1)
    y3 = ce(y2,y_label)
    loss = sum(y3)
    return y1, y2, y3, loss

反向求导：

def gradients(y1, y2, y3, X, y_label):
    grad1 = np.ones(len(y3))
    grad2 = ce_grad(y2, y_label)
    grad3 = sigmoid(y1)*(1-sigmoid(y1))
    grad4_w = X
    grad4_b = 1
    return (
        np.dot(grad1, grad2*grad3*grad4_w),
        np.dot(grad1, grad2*grad3*grad4_b)
    )

grad_w数值：

array([2.34286961, 2.97101168, 1.98692618, 1.81275096, 2.52826215,
       2.42595535, 1.9706045 ])

使用tensorflow进行验证：

import tensorflow as tf

def ce(y, y_label):
    return -y_label*tf.log(y) - (1-y_label)*tf.log(1-y)
    
X = tf.Variable(X)
w = tf.Variable(w)
b = tf.Variable(b)

y1 = tf.matmul(X, w) + b
y2 = tf.sigmoid(y1)
y3 = ce(y2, y_label)
loss = tf.reduce_sum(y3)
grad = tf.gradients(loss, w)
with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    ret = sess.run(grad)

grad数值：

array([[2.34286961],
        [2.97101168],
        [1.98692618],
        [1.81275096],
        [2.52826215],
        [2.42595535],
        [1.9706045 ]])

可见我们自己计算的导数没有问题。

使用随机梯度下降进行参数更新。随机梯度下降，一般会随机选择batch，这里为了简便，直接将所有向量进行BP：

训练过程：

yita = 1e-2
train_num = 10000
for i in range(train_num):
    y1, y2, y3, loss = forward(X, w, b)
    g_w, g_b = gradients(y1, y2, y3, X, y_label)
    w -= yita*g_w.reshape([-1, 1])
    b -= yita*g_b
    if i % 1000 == 0:
        print("loss:", loss)

输出：

loss: [11.6081676]
loss: [1.18844796]
loss: [0.71728752]
loss: [0.49936237]
loss: [0.37872785]
loss: [0.30340733]
loss: [0.25233963]
loss: [0.21561081]
loss: [0.18800623]
loss: [0.16654284]

之后看看forward输出：

>>> forward(X, w, b)[1]
array([[0.01485668],
       [0.00538101],
       [0.01436137],
       [0.01684294],
       [0.0247247 ],
       [0.93002105]])

y_label：

>>> y_label
array([[0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [1.]], dtype=float32)

06

总结

BP过程相对基础，但确实不是很简单。例如，此题loss关于w的导数，是典型的标量关于向量求导。关于向量、矩阵的求导，推荐以下文章：

https://zhuanlan.zhihu.com/p/24709748

本文目的在于学术交流，并不代表本公众号赞同其观点或对其内容真实性负责，版权归原作者所有，如有侵权请告知删除。

直播预告

左划查看更多

历史文章推荐

分享、点赞、在看，给个三连击呗！