# LogisticRegression 逻辑回归

逻辑回归是一种统计学中的经典分类模型，其本质是假设数据服从某个分布，然后使用极大似然来估计这个分布的参数。逻辑回归可以解决二分类问题也可以解决多分类问题。

逻辑回归 python 实现：[这里](https://github.com/TryFire/Machine-Learning-Algos-theory-and-implementation/blob/master/LogisticRegression/logistic_regression.py)

**二分类逻辑回归定义**假设数据集 $(x_1,y_1),\dots,(x_i,y_i),\dots,(x_N,y_N),x_i \in \mathbb{R^n}, y_i \in \{0,1\}$， $x_i$ 代表第 $i$ 个样本的特征向量，$y_i$ 代表第 $i$ 个样本的标签，二分逻辑回归模型是如下的条件概率分布：
$$
\begin{align}
& P(y=1|x) = \frac{\exp(w^Tx+b)}{1+\exp(w^Tx+b)} = \frac{1}{1+\exp(-(w^Tx+b))} \\
& P(y=0|x) = 1- P(y=1|x)= \frac{1}{1+\exp(w^Tx+b)} \\
\end{align}
$$
其中$(w \in \mathbb{R^n},b \in \mathbb{R})$是参数。为了简化，可以令$w = [b,w_1,\dots,w_n],x=[1,x_1,\dots,x_n]$， 并且令 $\sigma(t) = \frac{1}{1+\exp(-t)}$，那么如上条件概率就可以简化为：
$$
\begin{align}
& P(y=1|x) = \sigma(w^Tx) \\
& P(y=0|x) = 1-\sigma(w^Tx) = \sigma(-w^Tx) \\
\end{align}
$$

逻辑回归模型就比较两个条件概率值的大小，将实例 $x$ 归类到概率比较大的那一类。

其实逻辑回归和感知器有相似之处。感知器中把$w^Tx=0$看做是分类的超平面，$w^Tx>0$和$$w^Tx<0$$ 分别归到两类。在逻辑回归中，如果$w^Tx>0$，那么$P(y=1|x)>0.5，P(y=0|x)<0.5$，如果$$w^Tx<0$$，那么$P(y=1|x)<0.5，P(y=0|x)>0.5$，正好 $w^Tx=0$ 是二分类的界限。逻辑回归针对一个输入，输出一个把输入划分给某个类概率值，很直观。

## 用极大似然估计来学习参数$w$

这里假设每个输入 $x_{i}$ ，其标签服从参数为 $p(x_i)=\sigma(w^Tx)$ 的伯努利分布，也就是说针对 $x_i$，$P(y_i=1|x_i) = p(x_i)$，所以对$x_i$ 来说，其标签为$y_i$ 的概率为$P(y_i|x_i) = p(x_i)^{y_i}*(1-p(x_i))^{1-y_i}$。训练集 $(x_1,y_1),\dots,(x_i,y_i),\dots,(x_N,y_N),x_i \in \mathbb{R^n}, y_i \in \{0,1\}$的似然函数为：
$$
\begin{align*}
L(w) & = P(y_1,\dots,y_i,\dots,y_N | x_i,\dots,x_i,\dots,x_N) \\
& = \prod_{i=1}^{N} P(y_i|x_i) \ (假设每个样本独立) \\
& = \prod_{i=1}^{N} p(x_i)^{y_i}*(1-p(x_i))^{1-y_i}\\
& = \prod_{i=1}^{N} \sigma(w^Tx_i)^{y_i}*(1-\sigma(w^Tx_i))^{1-y_i}
\end{align*}
$$
我们是想找到参数$w$使得$L(w)$ 越大越好，即 $\max L(w)$，这又等价于 $\min -\log(L(w))$，所以我们的损失函数就出来了：

$$
\begin{aligned}
\text{Loss}(w) & = -\log(L(w)) \\
& = -\log(\prod_{i=1}^{N} \sigma(w^Tx_i)^{y_i}*(1-\sigma(w^Tx_i))^{1-y_i}) \\
& = -[ \sum_{i=1}^{N} y_i \log(\sigma(w^Tx_i)) + (1-y_i)\log(\sigma(-w^Tx_i)) ] \\
& = -[\sum_{i=1}^{N} y_i \log(\sigma(w^Tx_i)) - y_i\log(\sigma(-w^Tx_i)) + \log(\sigma(-w^Tx_i)) ] \\
& = -[\sum_{i=1}^{N} y_i \log(\frac{\sigma(w^Tx_i)}{\sigma(-w^Tx_i)}) + \log(\sigma(-w^Tx_i)) ] \\
&= -[\sum_{i=1}^{N} y_i \log(\frac{P(y_i=1|x_i)}{P(y_i=0|x_i)}) + \log(P(y=0|x)) ] \\
&= -[\sum_{i=1}^{N} y_i \log(\frac{\frac{\exp(w^Tx_i)}{1+\exp(w^Tx_i)}}{\frac{1}{1+\exp(w^Tx_i)}}) + \log(\frac{1}{1+\exp(w^Tx_i)}) ] \\
&= \sum_{i=1}^{N} \log(1+\exp(w^Tx_i)) - y_i(w^Tx_i) \\
\end{aligned}
$$

我们的目标就是找到合适的$w$ 来最小化无约束损失函数：
$$
\min \text{Loss}(w) = \min_{w \in \mathbb{R^{n+1}}} \sum_{i=1}^{N} \log(1+\exp(w^Tx_i)) - y_i(w^Tx_i)
$$

### 求解 $\min \text{Loss}(w)$

要求解 $\min \text{Loss}(w)$，我们常用的办法是如果函数 $\text{Loss}(w)$ 是凸函数(convex)，并且其一阶导数等于0不可解（没有**闭式解**），那么我们可以用梯度下降法（SGD）来求解可以最小化 $\text{Loss}(w)$的近似解。

证明 $\text{Loss}(w)$是凸函数一般有这几种方法：

- 如果 $f(\frac{a+b}{2}) \leq \frac{f(a)+f(b)}{2},\forall a,b \in I$ ，那么函数 $f$ 在区间 $I$ 上是convex的。
- 如果 $f(a) \geq f(b) + f^{'}(b)(a-b), \forall a,b \in I$，那么函数 $f$ 在区间 $I$ 上是convex的。
- 如果 $f^{''}(x) \geq 0, \forall x \in I$ ，么函数 $f$ 在区间 $I$ 上是convex的。
- 函数$f$由一些基本的凸函数的组合，比如$g(x)$ 是convex，那么$f(x) = c*g(x)，f(x) = g(ax+b)$ 或者多个凸函数的线性组合。

如果我们可以证明 $g(w) = \log(1+\exp(w^Tx_i)) - y_i(w^Tx_i)$ 是convex，那么 $\text{Loss}(w)$ 也是convex的。

用第一种方法好像有点难，我们先求一下$g(w)$的一阶导数：
$$
\begin{align}
\frac{\partial g(w)}{\partial w^{(m)}} = \frac{\exp(w^Tx_i)x_{i}^{(m)}}{1+\exp(w^Tx_i)} - y_ix_i^{(m)}
\end{align}
$$
其中$w^{(m)}$是$w$的第$(m)$维的值，$x_{i}^{(m)}$ 是 $x_i$ 这个向量的第$(m)$维的值。

看起来用第三种方法会比较简便，就是求解$g(w)$的**Hessian矩阵 $\frac{\partial^2 g(w)}{\partial w^{(m)}\partial w^{(n)}}$** ：
$$
\begin{aligned}
\frac{\partial^2 g(w)}{\partial w^{(m)}\partial w^{(n)}} & = - \frac{\exp(-w^Tx_i)(-x_i^{(n)})}{(1+\exp(-w^Tx_i))^2}x_i^{(m)} \\
& = \frac{\exp(-w^Tx_i)}{1+\exp(-w^Tx_i)}\frac{1}{1+\exp(-w^Tx_i)}x_i^{(m)}x_i^{(n)} \\
& = \sigma(w^Tx_i)(1-\sigma(w^Tx_i))x_i^{(m)}x_i^{(n)}
\end{aligned}
$$

可以得到**Hessian矩阵**：
$$
\nabla^2g(w) = \alpha (x_ix_i^T)
$$
其中 $\alpha =\sigma(w^Tx_i)(1-\sigma(w^Tx_i)) > 0, \forall w,x_i$ 。其中 $x_i x_i^T$ 是熟悉的二次函数，所以 $x_i x_i^T >= 0\forall x_i$ ，由此我们可以得出 $\nabla^2g(w) \geq 0, \forall w$ ，所以 $g(w)$ 是convex的。因为 $\text{Loss}(w)$ 是由多个 $g(w)$ 线性组合而成，所以 $\text{Loss}(w)$ 也是convex的。

### SGD

因为$\text{Loss}(w)$是convex，我们可以用$\nabla_w \text{Loss}(w) = 0$ 来获得最优的 $w$ ：
$$
\begin{align}
\nabla_w \text{Loss}(w) & = \sum_{i=1}^N (\frac{\exp(w^Tx_i)}{1+\exp(w^Tx_i)} - y_i)x_i \\
& = \sum_{i=1}^N (\sigma(w^Tx_i) - y_i)x_i \\
\end{align}
$$
可惜我们无法求出使得$\nabla_w \text{Loss}(w) = 0 $时候的解，所以可以用梯度下降法（SGD）来求解可以最小化 $\text{Loss}(w)$的近似解。

令$Q_i = \sigma(w^Tx_i), Y_i = y_i, X_i = x_i, Q \in \mathbb{R}^{N \times 1},Y \in \mathbb{R}^{N \times 1}, X \in \mathbb{R}^{N \times (n+1)}$，则更新规则：

$$
\begin{align}
w^{\text{new}} & = w^{(\text{old})} - \eta \nabla_w \text{Loss}(w) \\
& = w^{(\text{old})} - \eta X^T(Q-Y)
\end{align}
$$

其中$\eta \in (0,1]$是学习速率。

### 牛顿法

我们发现$\text{Loss}(w)$在任意一处是高阶可导的，所以可以对其在当前$w^{\text{old}}$处进行泰勒展开：
$$
\text{Loss2}(w) = \text{Loss}(w^{\text{old}}) + \text{Loss}^{(1)}(w^{\text{old}})(w - w^{\text{old}}) + \frac{\text{Loss}^{(2)}(w^{\text{old}})(w - w^{\text{old}})^2 }{2}
$$
这样我们得到了一个新的损失函数$\text{Loss2}(w)$的表达式。而且我们发现它是convex的，所以可以试着求解使得$\nabla_w\text{Loss2}(w) = 0$ 的解：
$$
\nabla_w\text{Loss2}(w) = \text{Loss}^{(1)}(w^{\text{old}}) + \text{Loss}^{(2)}(w^{\text{old}})(w - w^{\text{old}}) = 0 \\
w = w^{\text{old}} - h^{-1}g
$$
其中$h = \text{Loss}^{(1)}(w^{\text{old}})$ 是原损失函数$\text{Loss}(w)$的一阶导数，$g=\text{Loss}^{(2)}(w^{\text{old}})$是原损失函数$\text{Loss}(w)$的二阶导数，$h^{-1}$ 是矩阵 $h$ 的逆。

原损失函数$\text{Loss}(w)$的一阶和二阶导数很好求啊。其中一阶导数在随机梯度下降部分已经求过，我们copy一下：
$$
\begin{align}
g = \text{Loss}^{(1)}(w^{\text{old}}) & = \sum_{i=1}^N (\sigma({w^{\text{old}}}^Tx_i) - y_i)x_i \\
\end{align}
$$
二阶导数在证明 $\text{Loss}(w)$ 是convex的时候求解了一部分了，我们只需要加上累加符号：
$$
\begin{aligned}
\frac{\partial \text{Loss}(w)}{\partial w^{(m)}\partial w^{(n)}} & = \sum_{i=1}^N \sigma(w^Tx_i)(1-\sigma(w^Tx_i))x_i^{(m)}x_i^{(n)} \\
h & = \text{Loss}^{(2)}(w^{\text{old}}) = X^T A X
\end{aligned}
$$
其中$X_i = x_i, X \in \mathbb{R}^{N \times (n+1)}$, $A_{i,i} = \sigma(w^Tx_i)(1-\sigma(w^Tx_i)), A \in \mathbb{R}^{N \times N}$ , A 是对角矩阵。这里的$w$ 都用该用当前的$w^{\text{old}}$ 。

有了 $\text{Loss}(w)$ 的一阶和二阶导数，我们就可以利用更新规则对$w$ 进行更新了：
$$
w = w^{\text{old}} - \eta h^{-1}g
$$
其中$\eta \in (0,1]$是学习速率。

#### Iteratively reweighted least squares (IRLS) 法

有空补

## 正则化

当模型复杂度过大时候，模型会产生过拟合现象，所以需要选择复杂度适当的模型。可以使用正则化，在损失函数上加一个正则化项，或者对模型结构复杂度的一个惩罚项，来防止过拟合。常用的正则化有$L1$ 和$L2$ 正则。

### L1正则

我们直接在 $\text{Loss}(w)$上加一个代表模型复杂度的 $w$ 的$L1$ 范式。
$$
\text{Loss}(w) = \sum_{i=1}^{N} \left ( \log(1+\exp(w^Tx_i)) - y_i(w^Tx_i) \right ) + \lambda \sum_{k=1}^{n+1} |w_k|
$$
其中$\lambda > 0$ 用于调整经验风险和惩罚项的关系。那么其梯度为：
$$
\begin{align}
\nabla_w \text{Loss}(w)
& = \sum_{i=1}^N (\sigma(w^Tx_i) - y_i)x_i + I\\
I^{(m)} & = \left \{ \begin{aligned} \lambda,w^{(m)} > 0 \\ -\lambda,w^{(m)} < 0 \\ 0,w^{(m)} = 0 \end{aligned} \right.
\end{align}
$$
我们依然可以用梯度下降来优化损失函数。然而$L1$ 正则对牛顿法没什么用，因为其Hessian矩阵没有变化。

$L1$ 正则使得梯度下降法让$w$ 向着$0$的方向下降，而且$w$ 向量中任一维度的下降对$\text{Loss}(w)$的影响是一样的，$w$ 向量中在0附近的梯度和远离0的梯度是一样的，所以最后得出来的 $w$ 向量中某些维度可能是$0$.

### L2正则

我们直接在 $\text{Loss}(w)$上加一个代表模型复杂度的 $w$ 的$L2$范式。
$$
\text{Loss}(w) = \sum_{i=1}^{N} \left ( \log(1+\exp(w^Tx_i)) - y_i(w^Tx_i) \right ) + \frac{\lambda}{2} w^Tw
$$
那么其梯度为：
$$
\begin{align}
\nabla_w \text{Loss}(w)
& = \sum_{i=1}^N (\sigma(w^Tx_i) - y_i)x_i + \lambda w\\
\end{align}
$$
我们依然可以用梯度下降来优化损失函数。

$L2$ 正则使得梯度下降法让$w$ 向着$0$的方向下降，但是$w$ 向量中 值比较大 的维度下降 比 值比较小 的维度下降对$\text{Loss}(w)$的影响要大，而且$w$某些维度越接近 $0$ ，梯度越小，所以最后得出来的$w$ 中的维度的值很小的概率是$0$，他们会比较平滑得在$0$ 附近。

如果用牛顿法，一阶导数如上，Hessian 矩阵如下
$$
\begin{align}
h & = \text{Loss}^{(2)}(w^{\text{old}}) = X^T A X + E \\
\end{align}
$$
$X, A$ 的解释见上面的牛顿法那儿，$E_{m,m} = \lambda, E \in \mathbb{R}^{(n+1) \times (n+1)}$ 是对角矩阵。

## 逻辑回归多类分类

## 逻辑回归于其他模型的关系

- 最大熵模型
-

## 逻辑回归优缺点，适用场景

- 输出概率值，可解释性强，分给哪个类很直观。
- 处理离散特征还挺好，因为每个离散特征有独立的权重
- 对特征维度大的数据集处理比较快，因为用的是 权重和特征的内积$(w^Tx)$，**梯度下降**或者**牛顿法**(见上面更新规则的表达式)，用的都是向量内积，不是一个一个特征循环，和树模型不一样，特征越多，树越大，越深，计算速度会下降很多。

- 只能处理$x, y$ 线性关系的数据集

- 依赖特征工程，对某些特征要进行离散化处理，比如年龄，不能处理类别特征，需要手动one-hot

- 对连续特征进行归一化，梯度下降收敛更快