ML入门-Logistic优化求解与牛顿法

Posted: 2019-02-03 Categories: MachineLearning , LogisticRegression Tags: AI , MachineLearning , 人工智能 , 机器学习 , Logistic回归 Symbols: 15k Duration ≈ 17 mins.

（1）牛顿法和拟牛顿法；（2）Logistic 回归的优化求解。

ML入门-Logistic优化求解与牛顿法

前言：在机器学习的优化求解算法中，梯度下降法是无约束优化问题优化求解最常用的方法之一，还有一种求解方法就是牛顿法。相比梯度下降法，牛顿法的收敛速度更快，但同时，每次迭代需要的计算量也更大。

在牛顿法之前，需要先了解：泰勒公式。泰勒公式通俗地讲，就是当函数 $f$ 在点 $x$ 处的一阶导数、二阶导数……n 阶导数已知时，即可使用 n 阶泰勒展开来逼近函数 $f$ 在点 $x$ 的邻域的函数值，因此泰勒公式求的是一个点的邻域的近似函数值。

对应优化求解中，通常我们习惯于使用迭代法来求解，而迭代法的本质即：每次向正确的方向移动一小段，直到达到给定条件。这就与泰勒公式不谋而合——每次移动一小段之后的函数值可以使用移动前该点的泰勒展开来逼近，因此牛顿法或拟牛顿法都是在泰勒公式的基础上进行的。

1. 牛顿法

牛顿法（Newton - Raphson，牛顿 - 拉夫逊）是牛顿在 17 世纪提出的用于求解方程的根的方法。其求解思想如下：

假设点 $x^{\ast}$ 为函数 $f(x)$ 的根，则 $f(x) = 0$。
将函数 $f(x)$ 在点 $x_0$ 处进行一阶泰勒展开有：$f(x) \approx f(x_0) + (x - x_0) f'(x_0)$
假设点 $x$ 为 $x_0$ 邻域内一点，且 $x$ 为函数的根，则有：$f(x) \approx f(x_0) + (x - x_0) f'(x_0) = 0$
将上式变换即可得：$x = x_0 - \dfrac {f(x_0)} {f'(x_0)}$

上述牛顿法得到的结论，拓展到迭代的过程中，假设当前处在迭代第 t 轮，则可以得到下一轮 (t + 1) 时刻的解的表达式为：

$
x^{(t + 1)} = x^{(t)} - \dfrac {f(x^{(t)})} {f'(x^{(t)})}
$

这就是牛顿法优化求解的基本思想。下图展示了牛顿法求解方程 $f(x) = 0$ 的根的过程（图自 Wiki）：

牛顿法迭代过程

回到最优化问题中，通常会将问题转化成求极小值（误差、损失最小等），极小值对应了函数的导数为 0，因此需要适当调整牛顿法的目标，从求 $f(x) = 0$ 变为求 $f'(x) = 0$，因此原问题变为求 $f'(x)$ 的根。

令 $g(x) = f'(x)$，则关于 $x$ 的迭代条件变为：

$
x^{(t + 1)} = x^{(t)} - \dfrac {g(x^{(t)})} {g'(x^{(t)})} = x^{(t)} - \dfrac {f'(x^{(t)})} {f''(x^{(t)})}
$

在实际问题中，通常输入 $X$ 的维度都大于 1，因此将一阶导数替换为梯度（即 $f$ 分别对每个 $x_i$ 求偏导后组成向量）：

$
\nabla f(x_1, \cdots, x_D)
$

将二阶导数替换为海森（Hessian）矩阵 H：

$
H(X) = \left[
\begin{matrix}
\dfrac {\partial^2 f} {\partial^2 x^2_1} & \dfrac {\partial^2 f} {\partial x_1 \partial x_2} & \cdots & \dfrac {\partial^2 f} {\partial x_1 \partial x_D}
\\
\dfrac {\partial^2 f} {\partial x_2 \partial x_1} & \dfrac {\partial^2 f} {\partial^2 x^2_2} & \cdots & \dfrac {\partial^2 f} {\partial x_2 \partial x_D}
\\
\vdots & \vdots & \ddots & \vdots
\\
\dfrac {\partial^2 f} {\partial x_D \partial x_1} & \dfrac {\partial^2 f} {\partial x_D \partial x_2} & \cdots & \dfrac {\partial^2 f} {\partial^2 x^2_D}
\end{matrix}
\right]
$

Hessian 矩阵即：第 i 行第 j 列的元素为 $f$ 先对 $x_i$ 求偏导后再对 $x_j$ 求偏导。由于 $f$ 先对 $x_i$ 再对 $x_j$ 求偏导和 $f$ 先对 $x_j$ 再对 $x_i$ 求偏导相等，即 $\dfrac {\partial^2 f} {\partial x_i \partial x_j} = \dfrac {\partial^2 f} {\partial x_j \partial x_i}$，因此 Hessian 矩阵是对称的。

这样，牛顿法的迭代公式就变换为：

$$
x^{(t + 1)} = x^{(t)} - H^{-1} (X^{(t)}) \ \nabla f(X^{(t)})
$$

二阶导转换为 Hessian 矩阵 $H(X^{(t)})$ 后作为分母，使用逆运算 $H^{-1} (X^{(t)})$ 来表示。

总结牛顿法求解目标函数极值的迭代步骤如下：

从 $t = 0$ 开始，初始化 $X^{(0)}$ 为随机值
计算目标函数 $f(X)$ 在点 $X^{(t)}$ 的梯度：$g^{(t)} = \nabla f(X^{(t)})$，以及 Hessian 矩阵：$H^{(t)} = H(X^{(t)})$
计算移动方向：$d^{(t)} = (H^{(t)})^{-1} \ g^{(t)}$
根据迭代公式更新 $X$ 的值：$X^{(t + 1)} = X^{(t)} - d^{(t)}$
判断是否满足迭代终止条件（是否到达最大迭代次数，或相邻两次迭代的相对变化量或绝对变化量小于预设值，通常使用绝对量：$\dfrac {f(X^{(t + 1)}) - f(X^{(t)})} {f(X^{(t)})} \le \varepsilon$），若满足则循环计数，返回最佳参数 $X^{(t + 1)}$ 和目标函数极小值 $f(X^{(t + 1)})$，否则跳转到第 2 步

其中，第 3 步计算移动方向 $d^{(t)}$ 时，由于矩阵的逆求解困难，因此常用线性方程组计算：$H^{(t)} d^{(t)} = g^{(t)}$，当 $X$ 维度比较小时，可采用解析法求解 $d^{(t)}$，当 $X$ 维度比较高时，可采用梯度下降法或共而梯度下降法求解，因此对 $d^{(t)}$ 的求解又是一个迭代的计算过程。

对比梯度下降法中的移动方向：$d^{(t)} = - \eta g^{(t)}$，牛顿法：$d^{(t)} = - (H^{(t)})^{-1} \ g^{(t)}$，Hessian 矩阵相比学习率（步长）$\eta$ 包含的信息更多，因此牛顿法收敛速度更快，但从上述步骤也可明显看出牛顿法每次迭代的计算量都大幅增加。

由于梯度下降法仅使用了一阶导数，而牛顿法使用了二阶导数矩阵，因此梯度下降法是一阶最优化算法，而牛顿法是二阶最优化算法。

2. 拟牛顿法

牛顿法虽然收敛速度比梯度下降法更快，但在高维的情况下，计算目标函数二阶偏导数的复杂度很大，而且有时候目标函数的 Hessian 矩阵无法保持正定，不存在逆矩阵，此时牛顿法将不再能使用。

为此提出：拟牛顿法（Quasi-Newton Methods），拟牛顿法旨在：不用二阶偏导数，而构造出可以近似 Hessian 矩阵（或 Hessian 矩阵的逆矩阵）的正定对称矩阵，再逐步优化目标函数。不同的近似 Hessian 矩阵构造方法产生了不同的拟牛顿法：BFGS / L-BFGS。

扩展阅读：

谈谈常见的迭代优化方法

Mathematical optimization: finding minima of functions

重新考虑迭代条件，假设有某点 $x_0$，目标函数 $f(x)$ 在该点（已知点）进行二阶泰勒展开：

$
f(x) \approx f(x_0) + f'(x_0) (x - x_0) + \dfrac {1} {2} f''(x_0) (x - x_0)^2
$

当 $X$ 为向量 / 矩阵时，导数即为梯度，上式转换为：

$
f(X) \approx f(X_0) + \nabla f(X_0) (X - X_0) + \dfrac {1} {2} (X - X_0)^T \nabla^2 f(X_0) (X - X_0)
$

对上式取梯度运算（$X_0$ 为已知点，因此其与其导数均可视为常数项）得：

$
\nabla f(X) \approx \nabla f(X_0) + \nabla^2 f(X_0) (X - X_0)
$

由于函数在点 $X_0$ 处的泰勒展开可以近似 $X_0$ 邻域内的函数值，假设 $X_0$ 即为迭代 t 次后得到的 $X^{(t + 1)}$，则迭代前一轮的 $X^{(t)}$ 即为邻域内一点，可用二阶泰勒展开逼近，因此可得迭代关系：

$
\nabla f(X^{(t)}) \approx \nabla f(X^{(t + 1)}) + \nabla^2 f(X^{(t + 1)}) (X^{(t)} - X^{(t + 1)})
$

再使用 gradiant 和 Hessian 分别表示 $X$ 的一阶梯度和二阶梯度矩阵：$g^{(t)} = \nabla f(X^{(t)})$，$H^{(t)} = \nabla^2 f(X^{(t)})$，整理可得：

$
g^{(t + 1)} - g^{(t)} \approx H^{(t + 1)} (X^{(t + 1)} - X^{(t)})
$

进一步，引入记号：$s^{(t)} = X^{(t + 1)} - X^{(t)}$ 表示 $X$ 的变化量，$y^{(t)} = g^{(t + 1)} - g^{(t)}$ 表示梯度变化量，则可得简洁迭代关系：

$$
y^{(t)} \approx H^{(t + 1)} s^{(t)}
$$

由于牛顿法中 Hessian 矩阵的逆难以计算，因此在拟牛顿法中，令 $B$ 表示 $H$ 的近似，$D$ 表示 $H^{-1}$ 的近似，代入上式即可得到 拟牛顿法的条件< 为：

$$
y^{(t)} = B^{(t + 1)} s^{(t)}
$$

或

$$
s^{(t)} = D^{(t + 1)} y^{(t)}
$$

实际上，拟牛顿法的条件给出了 Hessian 矩阵的近似需要满足的条件。

2.1 BFGS

BFGS（Broyden, Fletcher, Glodfarb, Shanno）被认为是数值效果最好的拟牛顿法，且具有全局收敛性和超线性收敛速度。

BFGS 算法采用迭代法逼近 Hessian 矩阵：$B^{(t + 1)} = B^{(t)} + \Delta B^{(t)}$，初始值 $B^{(0)} = I$ 为单位矩阵，因此关键在于如何构造 $\Delta B^{(t)}$。

为保证矩阵 $B$ 的正定性，令 $\Delta B^{(t)} = \alpha u u^T + \beta v v^T$，代入上述拟牛顿法条件可得：

$
\begin{aligned}
y^{(t)} &= B^{(t + 1)} s^{(t)} = (B^{(t)} + \Delta B^{(t)}) s^{(t)} = B^{(t)} s^{(t)} + \Delta B^{(t)} s^{(t)}
\\
&= B^{(t)} s^{(t)} + \alpha u u^T s^{(t)} + \beta v v^T s^{(t)}
\\
&= B^{(t)} s^{(t)} + u (\alpha u^T s^{(t)}) + v (\beta v^T s^{(t)})
\end{aligned}
$

当 $u^T$ 或 $v^T$ 的维数与 $s^{(t)}$ 的维数一致（由于 $\Delta B^{(t)}$ 是构造的，因此可以构造为相同维数）时，$u^T s^{(t)}$ 以及 $v^T s^{(t)}$ 均为一个标量数值（向量的转置 x 向量 = 数值）。

令 $\alpha u^T s^{(t)} = 1$，$\beta v^T s^{(t)} = -1$，即 $\alpha = \dfrac {1} {u^T s^{(t)}}$，$\beta = - \dfrac {1} {v^T s^{(t)}}$，得到：

$
u - v = y^{(t)} - B^{(t)} s^{(t)}
$

不妨令 $u = y^{(t)}$，$v = B^{(t)} s^{(t)}$，代入 $\alpha$ 和 $\beta$ 的表达式得：

$
\alpha = \dfrac {1} {u^T s^{(t)}} = \dfrac {1} {(y^{(t)})^T s^{(t)}}
$

$
\beta = - \dfrac {1} {v^T s^{(t)}} = - \dfrac {1} {(B^{(t)} s^{(t)})^T s^{(t)}} = - \dfrac {1} {(s^{(t)})^T (B^{(t)})^T s^{(t)}}
$

代入 $\Delta B^{(t)}$ 的表达式得：

$$
\begin{aligned}
\Delta B^{(t)} &= \alpha u u^T + \beta v v^T
\\
&= \dfrac {y^{(t)} (y^{(t)})^T} {(y^{(t)})^T s^{(t)}} - \dfrac {B^{(t)} s^{(t)} (B^{(t)} s^{(t)})^T} {(s^{(t)})^T (B^{(t)})^T s^{(t)}}
\end{aligned}
$$

Sherman-Morrison 公式：若 $A$ 为非奇异方阵，$1 + v^T A^{-1} u \ne 0$，则有：

$
(A + uv^T)^{-1} = A^{-1} - \dfrac {A^{-1} u v^T A^{-1}} {1 + v^T A^{-1}}
$

由于牛顿法迭代过程需要计算 Hessian 矩阵的逆矩阵，因此根据 Sherman-Morrison 公式可得：

$$
(B^{(t + 1)})^{-1} = D^{(t + 1)} = \left( I - \dfrac {s^{(t)} (y^{(t)})^T} {(y^{(t)})^T s^{(t)}} \right) D^{(t)} \left( I - \dfrac {y^{(t)} (s^{(t)})^T} {(y^{(t)})^T s^{(t)}} \right) + \dfrac {s^{(t)} (s^{(t)})^T} {(y^{(t)})^T s^{(t)}}
$$

综上，对于 $\Delta B^{(t)}$ 和 $D^{(t + 1)}$，均可使用 $s^{(t)}$（$X$ 的变化量）和 $y^{(t)}$（函数梯度的变化量）表示，而初始化 $B^{(0)} = D^{(0)} = I$ 是一指的，$s^{(t)}$ 和 $y^{(t)}$ 在 $X$ 已知（优化求解的目的是寻找最优的 $X$，而 $X$ 本来就是已知的）时均是可计算的。

整理得 BFGS 更新参数的流程如下：

从 $t = 0$ 开始，初始化 $D^{(0)} = I$
计算移动方向：$d^{(t)} = D^{(t)} g^{(t)}$

先用 $B^{(t)} = B^{(t - 1)} + \Delta B^{(t - 1)}$ 迭代解出 $B^{(t)}$，再用 Sherman-Morrison 公式解出 $D^{(t)}$。
更新 $X$ 的值：$X^{(t + 1)} = X^{(t)} - d^{(t)}$
$s^{(t)} = d^{(t)}$
若 $||g^{(t + 1)}|| \le \varepsilon$，则迭代终止
计算：$y^{(t)} = g^{(t + 1)} - g^{(t)}$
$t = t + 1$，跳转第 2 步

2.2 L-BFGS

在 BFGS 中，每一轮迭代需要存储 Hessian 矩阵或其近似（$B^{(t)}$ 或 $D^{(t)}$），但当 $X$ 维数很高时，矩阵的维度也会很高，需要耗费大量存储空间。

L-BFGS（Limited memory BFGS）不直接存储 Hessian 矩阵或其近似（$B^{(t)}$ 或 $D^{(t)}$），而是存储迭代计算过程中的 $s^{(t)}$ 和 $y^{(t)}$ 来计算，从而减少参数存储所需空间。

在 BFGS 中，Hessian 矩阵的更新公式为：

$
(B^{(t + 1)})^{-1} = D^{(t + 1)} = \left( I - \dfrac {s^{(t)} (y^{(t)})^T} {(y^{(t)})^T s^{(t)}} \right) D^{(t)} \left( I - \dfrac {y^{(t)} (s^{(t)})^T} {(y^{(t)})^T s^{(t)}} \right) + \dfrac {s^{(t)} (s^{(t)})^T} {(y^{(t)})^T s^{(t)}}
$

令 $\rho^{(t)} = \dfrac {1} {(y^{(t)})^T s^{(t)}}$，$V^{(t)} = \left( I - \dfrac {y^{(t)} (s^{(t)})^T} {(y^{(t)})^T s^{(t)}} \right) = I - \rho^{(t)} y^{(t)} (s^{(t)})^T$

由于向量 $s^{(t)}$（X的变化量）和向量 $y^{(t)}$（梯度变化量）的维数相同（都与 $X$ 维度相等），因此 $s^{(t)} (y^{(t)})^T$ 和 $s^{(t)} (s^{(t)})^T$ 是矩阵，而 $(y^{(t)})^T s^{(t)}$ 是一个数值，因此 $\rho^{(t)}$ 是一个常数。

则有：

$
\begin{aligned}
(V^{(t)})^T &= \left( I - \rho^{(t)} y^{(t)} (s^{(t)})^T \right)^T
\\
&= I - \rho^{(t)} \left( y^{(t)} (s^{(t)})^T \right)^T
\\
&= I - \rho^{(t)} s^{(t)} (y^{(t)})^T
\\
&= I - \dfrac {s^{(t)} (y^{(t)})^T} {(y^{(t)})^T s^{(t)}}
\end{aligned}
$

因此原 Hessian 更新公式变为：

$$
D^{(t + 1)} = (V^{(t)})^T D^{(t)} V^{(t)} + \rho^{(t)} s^{(t)} (s^{(t)})^T
$$

将上述迭代更新公式展开得：

$
\begin{aligned}
D^{(t + 1)} &= (V^{(t)})^T D^{(t)} V^{(t)} + \rho^{(t)} s^{(t)} (s^{(t)})^T
\\
& \Downarrow
\\
D^{(1)} &= (V^{(0)})^T D^{(0)} V^{(0)} + \rho^{(0)} s^{(0)} (s^{(0)})^T
\\ \\
D^{(2)} &= (V^{(1)})^T D^{(1)} V^{(1)} + \rho^{(1)} s^{(1)} (s^{(1)})^T
\\
&= (V^{(1)})^T \left( (V^{(0)})^T D^{(0)} V^{(0)} + \rho^{(0)} s^{(0)} (s^{(0)})^T \right) V^{(1)} + \rho^{(1)} s^{(1)} (s^{(1)})^T
\\
&= (V^{(1)})^T (V^{(0)})^T D^{(0)} V^{(0)} V^{(1)} + (V^{(1)})^T \rho^{(0)} s^{(0)} (s^{(0)})^T V^{(1)} + \rho^{(1)} s^{(1)} (s^{(1)})^T
\\
& \ \ \vdots
\\
D^{(t + 1)} &= \left[ (V^{(t)})^T (V^{(t - 1)})^T \dots (V^{(0)})^T \right] D^{(0)} \left[ V^{(0)} V^{(1)} \dots V^{(t)} \right]
\\
& \ + \left[ (V^{(t)})^T (V^{(t - 1)})^T \dots (V^{(1)})^T \right] \left[ \rho^{(0)} s^{(0)} (s^{(0)})^T \right] \left[ V^{(1)} V^{(2)} \dots V^{(t)} \right]
\\
& \ + \left[ (V^{(t)})^T (V^{(t - 1)})^T \dots (V^{(2)})^T \right] \left[ \rho^{(1)} s^{(1)} (s^{(1)})^T \right] \left[ V^{(2)} V^{(3)} \dots V^{(t)} \right]
\\
& \ + \ \cdots \cdots
\\
& \ + (V^{(t)})^T \left[ \rho^{(t - 1)} s^{(t - 1)} (s^{(t - 1)})^T \right] V^{(t)}
\\
& \ + \rho^{(t)} s^{(t)} (s^{(t)})^T
\end{aligned}
$

从上述迭代过程可知，计算 $D^{(t + 1)}$ 需要用到 $\left \{ s^{(k)} y^{(k)} \right \}^t_{k = 0}$，若存储空间有限，仅能存储 m 组 $\left \{ s^{(k)} y^{(k)} \right \}$，当 $t > m$ 时，应当丢弃较早生成的 $\left \{ s^{(k)} y^{(k)} \right \}$。当然，由于丢弃了部分信息，此时计算的 $D^{(m + 1)}$ 是 $D^{(t + 1)}$ 的近似，也即 Hessian 矩阵的逆矩阵的进一步近似。

当 $t > m + 1$ 时，构造近似公式：

$
\begin{aligned}
D^{(t + 1)} &= \left[ (V^{(t)})^T (V^{(t - 1)})^T \dots (V^{(t - m + 1)})^T \right] D^{(0)} \left[ V^{(t - m + 1)} V^{(1)} \dots V^{(t)} \right]
\\
& \ + \left[ (V^{(t)})^T (V^{(t - 1)})^T \dots (V^{(t - m + 2)})^T \right] \left[ \rho^{(0)} s^{(0)} (s^{(0)})^T \right] \left[ V^{(t - m + 2)} V^{(2)} \dots V^{(t)} \right]
\\
& \ + \ \cdots \cdots
\\
& \ + (V^{(t)})^T \left[ \rho^{(t - 1)} s^{(t - 1)} (s^{(t - 1)})^T \right] V^{(t)}
\\
& \ + \rho^{(t)} s^{(t)} (s^{(t)})^T
\end{aligned}
$

$D^{(t)}$ 的迭代计算很繁琐，但计算 $D^{(t)}$ 的目的是为了得到搜索方向 $d^{(t)} = D^{(t)} g^{(t)}$，因此可以设计快速计算 $D^{(t)} g^{(t)}$ 的方法：

双向循环快速求解搜索方向

算法中还有部分不太明白，暂时只放上图片，待研究透彻后改为 Python 代码形式。另有几篇关于该双向循环快速求解 $D^{(t)} g^{(t)}$ 算法的参考文章如下：

机器学习中牛顿法凸优化的通俗解释

牛顿法与拟牛顿法学习笔记（五）L-BFGS 算法

spark L-BFGS实现

LBFGS.scala

LBFGS方法推导

3. Logistic优化求解算法

Logistic 回归采用 Logistic 损失 / 交叉熵损失：

$
L(y, \mu (X)) = - y \log (\mu (X)) - (1 - y) \log (1 - \mu (X))
$

其中 $y$ 为真值， $\mu (X)$ 为预测值为 1 的概率。

与其他机器学习一样，Logistic 回归的目标函数也包括两项：训练集上的损失和 + 正则项。

$
J(W; \lambda) = \sum^N_{i = 1} L(y_i, \mu (X_i; W)) + \lambda R(W)
$

由于 L1 正则在零点不可导，因此当正则项中含有 L1 正则时，不能直接使用基于梯度、Hessian 矩阵的优化求解算法，而通常使用坐标轴下降法求解，或也可使用次梯度法。

在给定正则参数 $\lambda$ 的情况下，目标函数的最优解为：$\hat{W} = \arg_W \min J(W, \lambda)$，取得最优解的必要条件即一阶导数为零：$\dfrac {\partial J(W, \lambda)} {\partial W} = 0$。

与线性回归模型不同的是，Logistic 回归模型的参数无法用解析法求解，因此可使用迭代法逼近求解。

其中一阶近似有与梯度相关的几个算法：

梯度下降（Logistic 使用梯度下降法收敛速度较慢）
随机梯度下降（SGD）
随机平均梯度法（SAG）
随机平均梯度法改进版（SAGA）
共轭梯度
坐标轴下降

二阶近似有：

牛顿法
拟牛顿法（BFGS、L-BFGS）

观察目标函数的损失和部分：

$
J_1 (W) = \sum^N_{i = 1} \left( - y_i \log (\mu (X_i; W)) - (1 - y_i) \log (1 - \mu (X_i; W)) \right)
$

3.1 梯度

最优化问题的求解离不开梯度的计算，由于目标函数中包含了 $\mu$，即 Sigmoid 变换，记 $\mu_i = \mu (X_i; W) = \sigma (W^T X_i)$，令 $z_i = W^T X_i$，根据复合函数的求导，$\dfrac {\partial \mu_i} {\partial W}$ 的求解如下：

$
\begin{aligned}
\dfrac {\partial \mu_i} {\partial W} &= \dfrac {\partial \sigma (W^T X_i)} {\partial W} = \dfrac {\partial \sigma {z_i}} {\partial W}
\\
&= \dfrac {d \sigma (z_i)} {d z_i} \dfrac {\partial z_i} {\partial W}
\end{aligned}
$

其中 ① 复合函数外层 $\sigma (z_i)$ 求导：

$
\begin{aligned}
\dfrac {d \sigma (z_i)} {d z_i} &= \dfrac {d (\dfrac {1} {1 + e^{- z_i}})} {d z_i}
\\
&= - \dfrac {1} {(1 + e^{- z_i})^2} \times \dfrac {d (1 + e^{- z_i})} {d z_i}
\\
&= - \dfrac {1} {(1 + e^{- z_i})^2} \times (- e^{- z_i})
\\
&= \dfrac {1} {(1 + e^{- z_i})^2} \times (e^{- z_i})
\\
&= \dfrac {1} {1 + e^{- z_i}} \times \dfrac {(1 + e^{- z_i}) - 1} {1 + e^{- z_i}}
\\
&= \sigma (z_i) \times (1 - \sigma (z_i))
\end{aligned}
$

② 复合函数内层 $z_i$ 求导：

$
\dfrac {\partial z_i} {\partial W} = \dfrac {\partial (W^T X_i)} {\partial W} = X_i
$

综合上述 2 式得：

$
\begin{aligned}
\dfrac {\partial \mu_i} {\partial W} &= \dfrac {d \sigma (z_i)} {d z_i} \dfrac {\partial z_i} {\partial W}
\\
&= \sigma (z_i) \times (1 - \sigma (z_i)) \times X_i
\\
&= \mu_i \times (1 - \mu_i) \times X_i
\end{aligned}
$

回到目标函数梯度，将上式代入求导得：

$
\begin{aligned}
J_1 (W) &= \sum^N_{i = 1} \left( - y_i \log (\mu (X_i; W)) - (1 - y_i) \log (1 - \mu (X_i; W)) \right)
\\
&= \sum^N_{i = 1} \left( - y_i \log \mu_i - (1 - y_i) \log (1 - \mu_i) \right)
\\
&= - \sum^N_{i = 1} \left(y_i \log \mu_i + (1 - y_i) \log (1 - \mu_i) \right)
\end{aligned}
$

$
\begin{aligned}
g_1 (W) &= \nabla J_1 (W) = \dfrac {d J_1 (W)} {d W}
\\
&= - d \left( \sum^N_{i = 1} \left( y_i \log \mu_i + (1 - y_i) \log (1 - \mu_i) \right) \right) / d W
\\
&= - \sum^N_{i = 1} \left( y_i \dfrac {1} {\mu_i} \times \dfrac {\partial \mu_i} {\partial W} + (1 - y_i) \dfrac {1} {1 - \mu_i} \times - \dfrac {\partial \mu_i} {\partial W} \right)
\\
&= - \sum^N_{i = 1} \left( y_i \dfrac {1} {\mu_i} - (1 - y_i) \dfrac {1} {1 - \mu_i} \right) \times \mu_i (1 - \mu_i) X_i
\\
&= - \sum^N_{i = 1} \left( y_i (1 - \mu_i) - (1 - y_i) \mu_i \right) \times X_i
\\
&= - \sum^N_{i = 1} (y_i - \mu_i) \times X_i
\\
&= (\mu - y) X
\end{aligned}
$

整理可得 Logistic 回归损失和部分的梯度表达式：

$$
g_1 (W) = X^T (\mu - y)
$$

这里直接解出来的结果是 $(\mu - y) X$，但表达式使用的是 $X^T (\mu - y)$，对此我有些自己的理解方式如下。

首先要提一下矩阵的形式。在机器学习中有一个很重要的工具包 numpy，这个工具包里其中两个很重要的类：numpy.matrix 和 numpy.array也即矩阵和向量，通常手动创建一个矩阵的时候可以用如下方式：

>>> import numpy as np

>>> A = np.matrix([[1, 2], [3, 4]])

>>> print(A)
[[1 2]
 [3 4]]

这就可以视为一个 2 x 2 的矩阵：$\left[ \begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix} \right]$，但是当我们如下创建一个向量时：

>>> B = np.array([1, 2])

>>> print(B)
[1 2]

>>> print(B.shape)
(2,)

可以看到，尽管输出的 B 的形式还是 [1 2]，但这是一个 2 行（也即 1 列）的向量（之所以没有标出 (2, 1) 是因为向量要么只有 1 行要么只有 1 列，防止与矩阵的 n x m 搞混了）。也就是说，默认情况下的向量是列向量，这也符合机器学习中的直觉，例如标签 y，预测值，或单个特征等都是列向量。

回到原问题的梯度表达式中，假设训练数据共有 D 维特征 N 个样本，则 X 是 N x D 维的矩阵，$\mu$ 和 y 均为 N 行的列向量，都知道矩阵的乘法 (A x B) 需要满足 A 的列数 = B 的行数时才有意义，而当作为矩阵运算 $(\mu - y) X$ 时，$(\mu - y)$ 列数为 1，此时无论 $X$ 是否转置，$(\mu - y) X$ 都无意义，不能做乘法运算。因此 $\sum^N_{i = 1} (\mu_i - y_i) \times X_i$ 转换为矩阵表达式时，将 $X$ 提到左乘并转置，不仅 $X^T$ 的列数恰好为 N，与 $\mu$ 和 $y$ 的行数相等，可以做乘法，而且从矩阵的乘法运算规则上符合直观地计算过程（$X$ 转置后，每一列为一个样本，分别与每个 $(\mu_i - y_i)$ 相乘）。

因此，梯度的表达式为：$g_1 (W) = X^T (\mu - y)$。

3.2 Hessian矩阵

Logistic 损失和部分的梯度为 $g_1 (W) = X^T (\mu - y)$，由此求解 Hessian 矩阵：

$
\begin{aligned}
H_1 (W) &= \dfrac {\partial g_1 (W)} {\partial W} = \dfrac {\partial \left( \sum^N_{i = 1} (\mu_i - y_i) X_i \right)} {\partial W} = \dfrac {\partial \left( \sum^N_{i = 1} X^T_i (\mu_i - y_i) \right)} {\partial W}
\\
&= \sum^N_{i = 1} X^T_i \dfrac {\partial (\mu_i)} {\partial W}
\\
&= \sum^N_{i = 1} X^T_i \mu_i (1 - \mu_i) X_i
\end{aligned}
$

说明：$\dfrac {\partial (a^T y)} {\partial y} = a^T$

令矩阵 $S \triangleq \mu_i (1 - \mu_i)$，即 $S$ 为对角阵，对角元素为 $\mu_i (1 - \mu_i)$。可得 Hessian 矩阵表达式：

$$
H_1 (W) = X^T S X
$$

3.3 牛顿法求解Logistic损失函数和极小值：IRLS

当解得梯度 $g_1 (W) = X^T (\mu - y)$ 和 Hessian 矩阵 $H_1 (W) = X^T S X$ 后，即可代入牛顿迭代公式中：

$
\begin{aligned}
W^{(t + 1)} &= W^{(t)} - H(W^{(t)})^{-1} g(W)
\\
&= W^{(t)} - (X^T S^{(t)} X)^{-1} X^T (\mu - y)
\\
&= (X^T S^{(t)} X)^{-1} (X^T S^{(t)} X) \times W^{(t)} - (X^T S^{(t)} X)^{-1} X^T (\mu - y)
\\
&= (X^T S^{(t)} X)^{-1} \times \left( X^T S^{(t)} X W^{(t)} - X^T (\mu - y) \right)
\\
&= (X^T S^{(t)} X)^{-1} \times \left( X^T S^{(t)} X W^{(t)} - X^T S^{(t)} {S^{(t)}}^{-1} (\mu - y) \right)
\\
&= (X^T S^{(t)} X)^{-1} X^T S^{(t)} \times \left( X W^{(t)} - {S^{(t)}}^{-1} (\mu - y) \right)
\end{aligned}
$

令 $z^{(t)} \triangleq X W^{(t)} - {S^{(t)}}^{-1} (\mu - y)$，则得到牛顿法迭代求 Logistic 损失函数和极小值的解：

$$
W^{(t + 1)} = (X^T S^{(t)} X)^{-1} X^T S^{(t)} z^{(t)}
$$

对比线性回归 $X W = y$ 方程中最小二乘的解：$\hat{W}_{OLS} = (X^T X)^{-1} X^T y$，给每个样本加权（每个样本的权重为 $S_i$）即可得加权最小二乘的解：

$
\hat{W}_{OLS_weight} = (X^T S X)^{-1} X^T y
$

牛顿迭代法求得的 Logistic 损失函数和极小值的解与加权最小二乘的解形式类似，因此也称为 迭代加权最小二乘（Iteratively Reweighted Least Squares, IRLS），其中每个样本的权重为 $S_i = \mu_i (1 - \mu_i)$。而 IRLS 又通过共轭梯度（Conjugate Gradient）法求解，因此 Scikit-Learn 中采用牛顿法求解的优化算法为 'newton-cg'。

4. Logistic的优化求解器Solver

Logistic 回归有多种优化求解方法。当使用 L2 正则时，可采用所有优化算法，而由于 L1 正则在零点处不可导，因此次不能使用需要计算梯度 / Hessian 矩阵的方法，此时可以类似 Lasso 求解，采用坐标轴下降法。

Scikit-Learn 中的 Logistic 类已在：《ML入门——Logistic回归简介》文中介绍，此处展开参数 solver 的一些可选项：

'liblinear'：
线性求解器，适用于小数据集，支持 L1 正则和 L2 正则。
内部使用了坐标轴下降法来迭代优化损失函数，如果模型的特征非常多，希望一些不重要的特征系数归零从而让模型系数稀疏的话，可以使用 L1 正则化。在多类 Logistic 回归任务中仅支持 OvR，不支持多项分布损失（MvM），但 MVM 相对精确。

'lbfgs'：
拟牛顿法，适用于较大数据集，仅支持 L2 正则。
支持 OvR 和 MvM 两种多类 Logistic 回归。

'newton-cg'：
牛顿法，适用于较大数据集，仅支持 L2 正则。
每个大迭代中的加权最小二乘回归部分采用共轭梯度算法实现。支持 OvR 和 MvM 两种多类 Logistic 回归。

'sag'：
随机平均梯度下降，适用于很大（如大于 5 万）的数据集，仅支持 L2 正则。
梯度下降法的变种，支持 OvR 和 MvM 两种多类 Logistic 回归。

'saga'：
改进的随机平均梯度下降，适用于非常大的数据集，支持 L1 正则。
当数据量很大，且选择 L1 正则时，只能采用 'saga' 优化求解器。支持 OvR 和 MvM 两种多类 Logistic 回归。

其中，'sag' 和 'saga' 只有在特征尺度大致相等时才能保证收敛，因此需要对数据做缩放（class sklearn.preprocessing 可以实现如：标准化、MinMaxScaler、MaxAbsScaler 等）。在实际任务中，大部分情况下数据预处理时都最好做标准化。实际上，加正则项本身也要求对每维特征做缩放。

另外，对于大数据集的训练任务，可以使用 SGDClassifier，并使用 LogLoss 作为损失函数。

若 SGDClassifier 使用 HingeLoss 作为损失函数，则为实现随机梯度下降的 SVM。在回归任务中还有 SDGRegressor。