一、简介
  在深度学习领域，我们往往采用梯度下降(或上升)法来优化训练函数模型，梯度下降法尤其是在解决凸优化问题上表现极佳。模型优化涉及到反向传播过程，反向传播过程需要先推导出梯度计算公式然后利用机器进行代数运算。这篇博文的工作是详细推导了逻辑回归反向传播梯度计算公式(什么是梯度？简单来讲就是成本函数对所有未知参数所组成向量的导数，梯度是成本函数变化最快的方向)，最后使用标准梯度下降算法来优化逻辑回归并附上逻辑回归优化实现代码。




  逻辑回归问题描述如Figure1所示，它是在多元线性回归的基础上对回归结果使用sigmoid函数进行激活计算得到条件概率$p(y=0|x)或p(y=1|x)$，用来解决二分类问题。这里多说一点如果你想实现多分类可以使用softmax，soft就体现在每一类别都有非零置信度(条件概率)，而不是像hard那样只有一个类别有1置信度其余置信度均为0。以一个样本作为输入的logistic regression forward propagation数学形式化表达为:




二、反向传播梯度计算公式推导
  接下来是logistic regression back propagation，反向传播梯度计算要求1.熟知前向传播代数计算过程，2.复合函数偏导数的链式法则[3]，3.函数矩阵或函数向量微分理论[4]。在反向传播推导梯度时可以(1)先计算一个输入样本的梯度，然后根据变量与变量的偏导数有相同shape(形式)来推广计算多个输入样本的梯度也就是实现梯度的向量化或叫矩阵化[1][2]，之所以要将梯度向量化是为了使用相关矩阵理论来简化代数运算。(2)直接对向量化后的多输入样本根据链式法则和向量求导理论推导梯度即可。接下来我将使用这两种方法来详细的推导梯度计算公式。

  首先，按方法(1)推导梯度。由于逻辑回归parameters就只有$w$和$b$，故而反向传播过程只需计算$\frac {∂L} {∂w}\equiv dw$和$\frac {∂L} {∂b}\equiv db$，从而获得梯度$grad(L)=(dw,db)$。过程如下所示。





  其次，按方法(2)计算得到梯度$grad(J)=(dw,db)$，推导过程如下所示。


三、逻辑回归标准梯度下降优化算法


  逻辑回归优化实现，详见我的github仓库：https://github.com/yiguanxian/implement-logistic-regression 。关于优化算法除了标准梯度下降法还有很多据此法改进的深度学习优化算法，我后面会写一篇关于深度学习中的优化算法的文章，待续…

  由于博主水平有限，博文中难免有错误和不当之处，欢迎小伙伴们给与批评和指正，也欢迎小伙伴们补充其他更简单的方法，可以直接留言评论。希望我的这篇博文能对你有所帮助，祝好！
参考：

[1]https://mooc.study.163.com/learn/2001281002?tid=2001392029#/learn/content?type=detail&id=2001702011

[2]https://mooc.study.163.com/learn/2001281002?tid=2001392029#/learn/content?type=detail&id=2001702012

[3] 欧阳光中，朱学炎等.数学分析下册[M].北京：高等教育出版社，2007：168.

[4]吴孟达，李兵等.高等工程数学[M].北京：科学出版社，2004:71-75.

1. $\operatorname {softmax}$ 函数的定义
对于任意一个 $n$ 维向量 $X = \begin {pmatrix} \vdots \\ {x}_{i} \\ \vdots \end {pmatrix}$ ，定义
$\operatorname {softmax} (X) = \begin {pmatrix} \vdots \\ \dfrac {e ^ {x_i}} { \sum \limits _{k} {e ^ {x_k} } } \\ \vdots \end {pmatrix}$
2. $\operatorname {softmax}$ 的梯度
记 $n$ 维向量 $\widehat {Y} = \begin {pmatrix} \vdots \\ \hat {y}_{i} \\ \vdots \end {pmatrix} = \operatorname {softmax} (X)$ ，则
$\dfrac {\partial } {\partial {x_j}} {\hat {y}_{i} } = \begin {cases} \dfrac {e ^ {x_j} \sum \limits _{k} {e ^ {x_k} } - e ^ {x_j} e ^ {x_j} }  { { \left ( \sum \limits _{k} {e ^ {x_k} } \right ) } ^2 }, & i = j  \\ \\ - \dfrac { e ^ {x_i} e ^ {x_j} }  { { \left ( \sum \limits _{k} {e ^ {x_k} } \right ) } ^2 }, & i \neq j \end {cases} = \begin {cases} \hat {y}_{j}  \left ( 1 - \hat {y}_{j} \right ), & i = j \\ - \hat {y}_{i} \hat {y}_{j} , & i \neq j \end {cases}$
3. $\text {CrossEntropyLoss}$ 的定义
设 $n$ 维向量 $Y = \begin {pmatrix} \vdots \\ {y}_{i} \\ \vdots \end {pmatrix}$ 是一个样本的标签，$n$ 维向量 $\widehat {Y} = \begin {pmatrix} \vdots \\ \hat {y}_{i} \\ \vdots \end {pmatrix}$ 是该样本的预测值，
且 $\sum \limits _{i} { \widehat y_{i} } = \sum \limits _{i} {y_{i}} = 1$ ，

则定义

$\text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) = - \sum \limits _{i} {y_{i} \ln \widehat y_{i}}$
4. $\text {CrossEntropyLoss with softmax}$ 的梯度
在多分类任务的机器学习中，通常取 $\widehat {Y} = \operatorname {softmax} \left ( X \right )$ ，此时
$\dfrac {\partial} {\partial x_{j}} {\left ( y_{i} \ln \widehat y_{i} \right )} = y_{i} \cdot \dfrac {1} {\widehat y_{i} }\cdot  \dfrac {\partial} {\partial x_{j}} {\widehat y_{i}}$

$= \begin {cases} y_{j} \cdot \dfrac {1} {\widehat y_{j} }\cdot  \widehat y_{j}  \left ( 1 - \widehat y_{j}   \right ), & i = j \\ - y_{i} \cdot \dfrac {1} {\widehat y_{i} }\cdot \widehat y_{i}  \widehat y_{j}, & i \neq j \end {cases}$

$= \begin {cases} y_{j}  \left ( 1 - \widehat y_{j}   \right ), & i = j \\ - y_{i} \widehat y_{j} , & i \neq j \end {cases}$
因此

$\dfrac {\partial} {\partial x_{j}} {\text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) }$

$= - \sum \limits _{i}  { \dfrac {\partial} {\partial x_{j}} {\left ( y_{i} \ln \widehat y_{i} \right )} }$
$= - \left [ y_{j}  \left ( 1 - \widehat y_{j}   \right ) + \sum \limits _{\underset {i \neq j} {i} } \left (- y_{i} \widehat y_{j} \right ) \right ]$
$= - \left [ y_{j}  \left ( 1 - \widehat y_{j}   \right ) - \widehat y_{j} \sum \limits _{\underset {i \neq j} {i} } {y_{i}} \right ]$
$= - \left [ y_{j}  \left ( 1 - \widehat y_{j}   \right ) - \widehat y_{j} \left (1 - y_{j} \right ) \right ]$

$= - \left ( y_{j} - \widehat y_{j} \right )$

$= \widehat y_{j} - y_{j}$
因此

$\dfrac {\partial} {\partial X} {\text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) } = \widehat {Y} - Y$
5. 论文  Distilling the Knowledge in a Neural Network 中公式 (2) 的推导
在论文 Distilling the Knowledge in a Neural Network 中，相当于取

$X = \dfrac {1} {T} \cdot Z$

$q_{i} =\widehat y_{i} = { \left ( \operatorname {softmax} (X) \right ) } _{i}$

$p_{i} = y_{i}$

则

$C = { \text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) }$

于是

$\dfrac {\partial C} {\partial z_{i}} = \dfrac {\partial} {\partial z_{i}} {\text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) }$
$= \dfrac {\partial x_{i}} {\partial z_{i}}  \cdot \dfrac {\partial} {\partial x_{i}} {\text {CrossEntropyLoss} \left (Y, \widehat {Y} \right ) }$
$= \dfrac {1} {T} \cdot \left ( \widehat y_{i} - y_{i} \right )$
$= \dfrac {1} {T} \cdot \left ( q_{i} - p_{i} \right )$
参考文献
Classification and Loss Evaluation - Softmax and Cross Entropy Loss

一. 改为向量化公式

上节中我们求得梯度下降的计算公式如下：



对上式子进行向量化处理。

（1）首先将第0项和其他项进行统一：（X0恒等于1）



进一步化成矩阵的乘法：





（2）最后转置得梯度为：



二. 修改梯度的函数实现方式：

（1）原本计算的代码：

def dJ(theta, X_b, y):
    res = np.empty(len(theta))
    res[0] = np.sum(X_b.dot(theta) - y)
    for i in range(1, len(theta)):
        res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
    return res * 2 / len(X_b)

（2）根据新的公式修改为：

def dJ(theta, X_b, y):
    return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)

（3）修改完的完整的线性回归模型：

import numpy as np
from sklearn.metrics import r2_score


class LinearRegression:

    # 初始化LinearRegression模型
    def __init__(self):
        self.coef_ = None
        self.interception_ = None
        # theta为私有变量
        self._theta = None

    # 使用正规化方程0
    def fit_normal(self, X_train, y_train):
        assert X_train.shape[0] == y_train.shape[0], \
            "X_train和y_train大小必须一致"

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)

        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict(self, X_predict):
        assert self.interception_ is not None and self.coef_ is not None, \
            "必须先拟合才能进行预测！"
        assert X_predict.shape[1] == len(self.coef_), \
            "预测值的特征数必须和训练集个数一致！"
        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return X_b.dot(self._theta)

    # 根据测试数据集X_test和y_test确定当前模型的准确度
    def score(self, X_test, y_test):
        y_predict = self.predict(X_test)
        return r2_score(y_test, y_predict)

    def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """
        根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型
        :param X_train:
        :param y_train:
        :param eta:
        :param n_iters:
        :return:
        """
        assert X_train.shape[0] == y_train.shape[0], "训练集X和y必须大小一致"

        def J(theta, X_b, y):
            """
            求代价函数
            :param theta: 参数值
            :param X_b:
            :param y:
            :return: 代价函数值
            """
            try:
                return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            # res = np.empty(len(theta))
            # res[0] = np.sum(X_b.dot(theta) - y)
            # for i in range(1, len(theta)):
            #     res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])
            # return res * 2 / len(X_b)
            return X_b.T.dot(X_b.dot(theta) - y) * 2. / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
            theta = initial_theta
            i_iter = 0

            while i_iter < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient
                if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
                    break
                i_iter += 1
            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])  # 为原先的X添加一列,(())
        init_theta = np.zeros(X_b.shape[1])

        self._theta = gradient_descent(X_b, y_train, init_theta, eta, n_iters)
        self.interception_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

 

 
黎曼局部差分向量，SPD流形局部差分向量计算
 
前言
$\quad$在欧氏空间中，对于任意两个向量 $x,y \in R^k$，它们之间的差分向量为：$x-y$，该差分向量也可以等价为欧氏空间距离公式的平方在 $x$ 点的梯度 $\triangledown_{x}\delta^2(x,y)$ 。

$\quad$相应的，SPD流形上的黎曼局部差分向量可以定义为：$\mathop \zeta(X,Y) = {\delta}(X,Y)\frac{\triangledown_{X}\delta^2(X,Y)}
{\parallel \triangledown_{X}\delta^2(X,Y)\parallel}$, 其中，这里的 $\delta$ 为SPD流形上的距离度量 [link]
 
黎曼局部差分向量
$\quad$由于SPD流形的特殊性，我们不能通过直接相减来得到差分向量，对应SPD流形上四种常用的距离度量 [link]，相应的黎曼局部差分向量有如下的几种具体形式： 

 
1.使用AIRM度量，黎曼局部差分向量有如下的表达形式：
$\mathop \zeta(X,Y) = 2X^{1/2} \log(X^{-1/2}YX^{-1/2})X^{1/2}$
2.使用Stein散度，黎曼局部差分向量有如下的表达形式：
$\mathop \zeta(X,Y) = X(X+Y)^{-1}X-\frac{1}{2}X$
3.使用Jeffrey散度，黎曼局部差分向量有如下的表达形式：
$\mathop \zeta(X,Y) = \frac{1}{2}X(Y^{-1}-X^{-1}YX^{-1})X$
4.使用LEM度量，黎曼局部差分向量有如下的表达形式：
$\mathop \zeta(X,Y) = X(\log(X)-\log(Y))+(\log(X)-\log(Y))X$
 

 

示例代码：https://github.com/Kai-Xuan/SPD-OPERATIONS/tree/master/SPD-LDV/ [link]
如果这个内容对于您的研究工作有帮助，我们将非常感激您可以引用我们的文章：[1].
 

 参考:

1. Chen K X, Ren J Y, Wu X J, et al. Covariance Descriptors on a Gaussian Manifold and their Application to Image Set Classification[J]. Pattern Recognition, 2020, 107: 107463. [link]