线性回归/16. 带约束条件的线性回归.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在前面的学习中，我们从同一份样本数据中得到了三个模型，而且三个模型的拟合效果都很好，不过很显然，二次模型和三次模型都过于复杂了，属于过拟合。为了判断哪个模型的性能最优，泛化能力最强，我们将数据集划分为训练集和测试集，然后选择在测试集上表现最好的模型。使用这种方法可以有效的评估模型结果，避免过拟合。\n",
    "\n",
    "这一节我们将学习一种避免过拟合的新方法：**正则化**（regularization）方法。\n",
    "\n",
    "在线性回归模型 $ \\hat{y} = w_0 + w_1x_1 + w_2x_2 + ... + w_nx_n = w^Tx $ 中，每一个变量 $x_i$ 就是模型的一个特征，参数 $w_i$ 代表相应的特征在模型中所占的权重。在前面的例子里，本来只有一个特征，我们人为造出了两个特征来（$x^2$ 和 $x^3$），得到了三个模型，这两个特征对模型来说其实根本没用，过拟合就是因为模型过于复杂，把这些没用的特征也学习到了。所以如果能降低这些没用特征的权重，譬如直接降到 0，那就不会出现过拟合了，为了达到这一点，我们可以稍微改进下我们的损失函数：\n",
    "\n",
    "$$\n",
    "loss = \\sum_{i=1}^m (y_i - w^Tx_i)^2 + \\lambda \\| w \\|^2\n",
    "$$\n",
    "\n",
    "这个函数相比于我们之前的平方损失函数多了一个 $\\lambda \\| w \\|^2$，这被称之为 **正则化项** 或 **惩罚项**，其中正则化参数 $\\lambda > 0$，它是一个超参数，$\\| w \\|^2$ 表示 $L_2$ 范数，如下：\n",
    "\n",
    "$$\n",
    "\\| w \\|^2 = \\sum_{i=1}^m w_i^2 = w_1^2 + w_2^2 + \\dots + w_m^2\n",
    "$$\n",
    "\n",
    "很显然，为了使上面的损失函数最小，模型中不相干的参数越小越好，从而降低过拟合的风险，使用这个损失函数的回归方法称为 **岭回归**（ridge regression）。\n",
    "\n",
    "上面的 $L_2$ 范数可以改成 $L_p$ 范数，当 $p = 1$ 时，有损失函数：\n",
    "\n",
    "$$\n",
    "loss = \\sum_{i=1}^m (y_i - w^Tx_i)^2 + \\lambda \\| w \\|\n",
    "$$\n",
    "\n",
    "其中 $L_1$ 范数如下：\n",
    "\n",
    "$$\n",
    "\\| w \\| = \\sum_{i=1}^m |w_i| = |w_1| + |w_2| + \\dots + |w_m|\n",
    "$$\n",
    "\n",
    "这种回归方法称为 **LASSO 回归**（Least Absolute Shrinkage and Selection Operator）。\n",
    "\n",
    "可以看到这里的损失函数仍然是最小二乘，但是在后面增加了 $L_1$ 或 $L_2$ 范数约束，所以岭回归又被称为 **基于 $L_2$ 约束的最小二乘法**，LASSO 回归又被称为 **基于 $L_1$ 约束的最小二乘法**。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $L_2$ 约束的最小二乘\n",
    "\n",
    "可以看到，为了避免过拟合，我们在模型的参数上施加了一定的约束条件。在普通最小二乘回归中，参数空间没有限制，参数可以取任何值，但是在岭回归和 Lasso 回归中，参数被限制在一定的范围内，这被称为 **部分空间约束的最小二乘学习法**。为了保证参数不会偏移到值域范围之外，可以附加一个 $P\\theta = \\theta$ 的约束条件。\n",
    "\n",
    "$$\n",
    "\\mathop{\\min}_{\\theta}J(\\theta), P\\theta = \\theta\n",
    "$$\n",
    "\n",
    "$P$ 是满足 $P^2 = P$ 和 $P^T = P$ 的 $b \\times b$ 维矩阵，表示的是矩阵 $P$ 的值域的正交投影矩阵。$P$ 通常是手工设置的，通过主成分分析法，也可以基于数据进行设置。\n",
    "\n",
    "> 疑惑：$P\\theta = \\theta$，$P$ 不就是单位矩阵吗？\n",
    "\n",
    "不过正交投影矩阵 $P$ 的设置有很大的自由度，所以改用操作相对容易的 $L_2$ 约束条件，譬如：\n",
    "\n",
    "$$\n",
    "\\mathop{\\min}_{\\theta}J(\\theta), \\|\\theta\\|^2 \\leqslant R\n",
    "$$\n",
    "\n",
    "这时的参数空间如下图所示：\n",
    "\n",
    "![](../images/l2-parameter-space.webp)\n",
    "\n",
    "从图中可以看出来，这种约束方法是以参数空间的原点为圆心，在一定半径范围的圆内（一般为**超球**）进行求解的。$R$ 表示的是圆的半径。\n",
    "\n",
    "利用**拉格朗日对偶问题**，该问题可以转换为求解下式的最优解：\n",
    "\n",
    "$$\n",
    "\\mathop{\\max}_{\\lambda} \\mathop{\\min}_{\\theta} [ J(\\theta) + \\frac{\\lambda}{2}(\\|\\theta\\|^2 - R) ], \\lambda \\geqslant 0\n",
    "$$\n",
    "\n",
    "拉格朗日对偶问题的拉格朗日乘子 $\\lambda$ 的解由圆的半径 $R$ 决定，如果不根据 $R$ 来决定 $\\lambda$，而是直接指定的话，$L_2$ 约束的最小二乘解 $\\hat{\\theta}$ 可以通过下式求得：\n",
    "\n",
    "$$\n",
    "\\hat{\\theta} = \\mathop{\\arg\\min}_{\\theta} [ J(\\theta) + \\frac{\\lambda}{2}\\|\\theta\\|^2 ]\n",
    "$$\n",
    "\n",
    "这就是岭回归。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 一般 $L_2$ 约束的最小二乘\n",
    "\n",
    "上面的 $L_2$ 约束可以写成更一般的情况：\n",
    "\n",
    "$$\n",
    "\\mathop{\\min}_{\\theta}J(\\theta), \\theta^TG\\theta \\leqslant R\n",
    "$$\n",
    "\n",
    "这里的 $G$ 是一个 $b \\times b$ 的正则化矩阵。当 $G$ 是对称正定矩阵的时候，$\\theta^TG\\theta \\leqslant R$ 可以把参数限制在一个椭圆区域内。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $L_1$ 约束的最小二乘\n",
    "\n",
    "和 $L_2$ 约束类似，$L_1$ 约束可以写成下面的形式：\n",
    "\n",
    "$$\n",
    "\\mathop{\\min}_{\\theta}J(\\theta), \\|\\theta\\| \\leqslant R\n",
    "$$\n",
    "\n",
    "这时的参数空间如下图所示：\n",
    "\n",
    "![](../images/l1-parameter-space.webp)\n",
    "\n",
    "可见参数被限制在一个菱形范围内，利用拉格朗日对偶问题，$L_1$ 约束的最小二乘解可以通过下式求得：\n",
    "\n",
    "$$\n",
    "\\hat{\\theta} = \\mathop{\\arg\\min}_{\\theta} ( J(\\theta) + \\lambda\\|\\theta\\| )\n",
    "$$\n",
    "\n",
    "这就是 Lasso 回归。如果我们在参数空间中画出损失函数 $J(\\theta)$ 的等高线，可以把 $L_1$ 约束和 $L_2$ 约束作一个对比。由于 $J(\\theta)$ 是平方损失函数，它是一个二次凸函数，所以它的等高线在参数空间呈椭圆状，椭圆的中心即是最小二乘解，如下图：\n",
    "\n",
    "![](../images/l1-vs-l2-parameter-space.webp)\n",
    "\n",
    "而等高线和约束边界的交点即为带约束条件的最小二乘解。可以直观的看出，由于 $L_1$ 约束边界在各个参数的轴有一个尖角，等高线和它的交点更容易出现在参数的轴上，也就意味着其他的参数等于 0，像这样的解被称为 **稀疏解**（Sparse solution）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 岭回归的几何解释\n",
    "\n",
    "如果上面的等高线看上去不是很直观，我们可以在三维场景下再看看岭回归到底是怎么限制参数空间的。在一元线性回归的求解过程中，我们知道，损失函数实际上就是三维空间的一个曲面，损失函数的最小值就是该曲面的最低点。\n",
    "\n",
    "![](../images/ridge-regression-1.webp)\n",
    "\n",
    "岭回归意味着参数的取值被限制在一个圆内，也就是下面三维空间中的圆柱体，这个圆柱体和曲面的交点就是岭回归的解。\n",
    "\n",
    "![](../images/ridge-regression-2.webp)\n",
    "\n",
    "所以，从参数平面理解，即为抛物面等高线在水平面的投影和圆的交点，如下图所示：\n",
    "\n",
    "![](../images/ridge-regression-3.webp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $L_p$ 约束的最小二乘\n",
    "\n",
    "上面介绍了 $L_1$ 约束和 $L_2$ 约束，我们可以将其进一步推广到 $L_p$ 约束：\n",
    "\n",
    "$$\n",
    "\\|\\theta\\|_p = \\lgroup \\sum_{j=1}^b |\\theta_j|^p \\rgroup ^{\\frac{1}{p}} \\leqslant R\n",
    "$$\n",
    "\n",
    "$L_p$ 范数中有两个比较特殊的场景，一种是当 $p = 0$ 时的 $L_0$ 范数：\n",
    "\n",
    "$$\n",
    "\\|\\theta\\|_0 = \\sum_{j=1}^b I(\\theta_j \\neq 0)\n",
    "$$\n",
    "\n",
    "其中 $I(\\theta_j \\neq 0)$ 为指示函数，$\\theta_j \\neq 0$ 时为 1，否则为 0，所以 $L_0$ 范数表示的是非零元素的个数。\n",
    "\n",
    "另一种场景是当 $p = \\infty$ 时的 $L_\\infty$ 范数：\n",
    "\n",
    "$$\n",
    "\\|\\theta\\|_\\infty = \\max \\{ |\\theta_1|, \\dots, |\\theta_b| \\}\n",
    "$$\n",
    "\n",
    "它表示的是元素绝对值中的最大值，因此 $L_\\infty$ 范数也称为 **最大值范数**。\n",
    "\n",
    "我们画出不同的 $p$ 值参数空间的约束边界，如下图：\n",
    "\n",
    "![](../images/lp-parameter-space.webp)\n",
    "\n",
    "当 $p \\leqslant 1$ 时，约束边界在坐标轴上呈有峰值的尖形；当 $p \\geqslant 1$ 时，呈凸形，在坐标轴上呈尖行是存在稀疏解的秘诀。但是，如果约束边界不是凸形的话，可能存在多个解，最优化往往很困难。因此 $p = 1$ 是稀疏解存在的唯一凸形。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $L_1 + L_2$ 约束的最小二乘\n",
    "\n",
    "$L_1$ 约束在实际使用时存在些许限制：\n",
    "\n",
    "1. 参数 b 比训练样本 n 多时，线性模型可选择的最大特征数被局限为 n；\n",
    "2. 如果有多个基函数相似的集合时，$L_1$ 约束会选择一个而忽略其它的，另外，$L_1$ 约束只能在多个相关性较强中的特征中选择一个；\n",
    "3. 参数 b 比训练样本 n 少时, $L_1$ 的通用性比 $L_2$ 稍差；\n",
    "\n",
    "通过结合 $L_1$ 和 $L_2$ 两个约束可以解决上述问题，这就是 $L_1 + L_2$ 约束的最小二乘，也被称为 **弹性网**。\n",
    "\n",
    "$$\n",
    "(1-\\tau)\\|\\theta\\|_1 + \\tau \\|\\theta\\|^2 \\leqslant R\n",
    "$$\n",
    "\n",
    "其中 $0 \\leqslant \\tau \\leqslant 1$，当 $\\tau = 0$ 时就是 $L_1$ 约束，当 $\\tau = 1$ 时就是 $L_2$ 约束，当 $\\tau = 0.5$ 时，约束边界如下图所示：\n",
    "\n",
    "![](../images/l1-l2-parameterr-space.webp)\n",
    "\n",
    "可以看出 $\\tau = 0.5$ 时的 $L_1 + L_2$ 约束和 $L_{1.4}$ 约束几乎一样，但是如果我们把角的部分放大，可以看出 $L_{1.4}$ 约束和 $L_2$ 约束一样平滑，而 $L_1 + L_2$ 约束却是和 $L_1$ 约束一样呈尖行。所以 $L_1 + L_2$ 约束也和 $L_1$ 约束一样容易求得稀疏解。\n",
    "\n",
    "实验证明 $L_1 + L_2$ 约束比 $L_1$ 约束具有更高的精度，但是为了调整 $L_1$ 和 $L_2$ 的平衡，引入了一个新的超参数 $\\tau$，求解要更复杂一点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考\n",
    "\n",
    "* https://www.jianshu.com/p/1677d27e08a7\n",
    "* https://www.jianshu.com/p/9c0b029478e9"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}