{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从一元线性回归和二元线性回归的解法中可以看到，线性回归的一般思路是这样的：\n",
    "\n",
    "1. 定义损失函数：$ loss = (y - \\hat{y})^2 $\n",
    "2. 将线性回归函数 $ \\hat{y} = a_0 + a_1x_1 + a_2x_2 + ... + a_nx_n $ 代入损失函数，并求解系数 $a_0$，$a_1$，...，$a_n$ 使得损失函数最小\n",
    "3. 损失函数的最小值可以通过对损失函数求导来实现，分别对每个系数求偏导，并令其等于 0\n",
    "4. 得到一个关于各个系数的线性方程组，求解这个线性方程组得到各个系数的值\n",
    "\n",
    "在求解二元线性回归时，线性方程组是这样的：\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial}{\\partial a_0} F(a_0, a_1, a_2) &=& 2a_0 + 2a_1\\bar{x_1} + 2a_2\\bar{x_2} - 2\\bar{y} &=& 0 \\\\\n",
    "\\frac{\\partial}{\\partial a_1} F(a_0, a_1, a_2) &=& 2a_1\\bar{x_1^2} + 2a_2\\bar{x_1x_2} + 2a_0\\bar{x_1} - 2\\bar{x_1y} &=& 0 \\\\\n",
    "\\frac{\\partial}{\\partial a_2} F(a_0, a_1, a_2) &=& 2a_2\\bar{x_2^2} + 2a_1\\bar{x_1x_2} + 2a_0\\bar{x_2} - 2\\bar{x_2y} &=& 0 \\\\\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "求得方程组的解是：\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "a_0 &=& \\bar{y} - a_1\\bar{x_1} - a_2\\bar{x_2} \\\\\n",
    "a_1 &=& \\frac{(\\bar{x_1y}-\\bar{x_1}\\bar{y})(\\bar{x_2^2}-(\\bar{x_2})^2) - (\\bar{x_2y}-\\bar{x_2}\\bar{y})(\\bar{x_1x_2}-\\bar{x_1}\\bar{x_2})}{(\\bar{x_1^2}-(\\bar{x_1})^2)(\\bar{x_2^2}-(\\bar{x_2})^2) - (\\bar{x_1x_2}-\\bar{x_1}\\bar{x_2})^2} \\\\\n",
    "a_2 &=& \\frac{(\\bar{x_1y}-\\bar{x_1}\\bar{y})(\\bar{x_1x_2}-\\bar{x_1}\\bar{x_2}) - (\\bar{x_2y}-\\bar{x_2}\\bar{y})(\\bar{x_1^2} - (\\bar{x_1})^2)}{(\\bar{x_1x_2}-\\bar{x_1}\\bar{x_2})^2 - (\\bar{x_2^2}-(\\bar{x_2})^2)(\\bar{x_1^2}-(\\bar{x_1})^2)}\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "可以看到，使用这种方法来求解线性回归非常麻烦。当我们扩展到多元线性回归时，这个计算量是极大的。那么有没有一种通用的方法来求解线性回归问题呢？我们不妨把上面的方程组化简：\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "a_0 + a_1\\bar{x_1} + a_2\\bar{x_2} &=& \\bar{y} \\\\\n",
    "a_1\\bar{x_1^2} + a_2\\bar{x_1x_2} + a_0\\bar{x_1} &=& \\bar{x_1y} \\\\\n",
    "a_2\\bar{x_2^2} + a_1\\bar{x_1x_2} + a_0\\bar{x_2} &=& \\bar{x_2y} \\\\\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "转换为矩阵形式：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}{\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "1 & \\bar{x_1} & \\bar{x_2} \\\\\n",
    "\\bar{x_1} & \\bar{x_1^2} & \\bar{x_1x_2} \\\\\n",
    "\\bar{x_2} & \\bar{x_1x_2} & \\bar{x_2^2}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "a_0 \\\\\n",
    "a_1 \\\\\n",
    "a_2\n",
    "\\end{array} \n",
    "\\right ]\n",
    "=\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "\\bar{y} \\\\\n",
    "\\bar{x_1y} \\\\\n",
    "\\bar{x_2y}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "}\\end{aligned}\n",
    "$$\n",
    "\n",
    "所以：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}{\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "a_0 \\\\\n",
    "a_1 \\\\\n",
    "a_2\n",
    "\\end{array} \n",
    "\\right ]\n",
    "=\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "1 & \\bar{x_1} & \\bar{x_2} \\\\\n",
    "\\bar{x_1} & \\bar{x_1^2} & \\bar{x_1x_2} \\\\\n",
    "\\bar{x_2} & \\bar{x_1x_2} & \\bar{x_2^2}\n",
    "\\end{array} \n",
    "\\right ]^{-1}\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "\\bar{y} \\\\\n",
    "\\bar{x_1y} \\\\\n",
    "\\bar{x_2y}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "}\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bar_x1 = 7600.0\n",
      "bar_x2 = 80.0\n",
      "bar_y = 408000.0\n",
      "bar_x1y = 3188000000.0\n",
      "bar_x2y = 33280000.0\n",
      "bar_x1x1 = 60000000.0\n",
      "bar_x1x2 = 620000.0\n",
      "bar_x2x2 = 6600.0\n",
      "a = [[62105.2631579 ]\n",
      " [   32.10526316]\n",
      " [ 1273.68421053]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X1 = np.array([6000, 6000, 8000, 8000, 10000])\n",
    "X2 = np.array([60, 80, 70, 100, 90])\n",
    "Y = np.array([340000, 350000, 400000, 450000, 500000])\n",
    "\n",
    "x1 = np.sum(X1) / X1.size\n",
    "x2 = np.sum(X2) / X2.size\n",
    "y = np.sum(Y) / Y.size\n",
    "print(\"bar_x1 = {0}\".format(x1))\n",
    "print(\"bar_x2 = {0}\".format(x2))\n",
    "print(\"bar_y = {0}\".format(y))\n",
    "\n",
    "x1y = np.sum(np.multiply(X1, Y)) / Y.size\n",
    "x2y = np.sum(np.multiply(X2, Y)) / Y.size\n",
    "x1x1 = np.sum(np.multiply(X1, X1)) / Y.size\n",
    "x1x2 = np.sum(np.multiply(X1, X2)) / Y.size\n",
    "x2x2 = np.sum(np.multiply(X2, X2)) / Y.size\n",
    "print(\"bar_x1y = {0}\".format(x1y))\n",
    "print(\"bar_x2y = {0}\".format(x2y))\n",
    "print(\"bar_x1x1 = {0}\".format(x1x1))\n",
    "print(\"bar_x1x2 = {0}\".format(x1x2))\n",
    "print(\"bar_x2x2 = {0}\".format(x2x2))\n",
    "\n",
    "# 不使用科学计数法\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "mx = np.matrix([1, x1, x2, x1, x1x1, x1x2, x2, x1x2, x2x2]).reshape(3, 3)\n",
    "my = np.matrix([y, x1y, x2y]).reshape(3, 1)\n",
    "#a = np.linalg.inv(mx).dot(my)\n",
    "a = mx.I * my\n",
    "print(\"a = {0}\".format(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里使用了矩阵求逆和乘法运算，在计算量上比之前的方法并没有减少，但是从编程的角度来说，确实简化了很多人工计算，直接使用 `a = mx.I * my` 就完成了多元线性方程组的求解。\n",
    "\n",
    "我们再回过头来看下线性回归函数： \n",
    "\n",
    "$$\n",
    "\\hat{y} = a_0 + a_1x_1 + a_2x_2 + ... + a_nx_n \n",
    "$$\n",
    "\n",
    "既然线性方程组可以转换为矩阵表示，这里的线性函数当然也可以转换为矩阵表示：\n",
    "\n",
    "$$\n",
    "\\hat{y} =\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "1 & x_1 & x_2 & ... & x_n\n",
    "\\end{array} \n",
    "\\right ]\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "a_0 \\\\\n",
    "a_1 \\\\\n",
    "a_2 \\\\\n",
    "... \\\\\n",
    "a_n\n",
    "\\end{array} \n",
    "\\right ]\n",
    "= \\bf{x}^{T}\\bf{a}\n",
    "$$\n",
    "\n",
    "其中，$\\bf{x}$ 和 $\\bf{a}$ 都是一个 $n \\times 1$ 的矩阵，也就是一个列向量，所以用小写加粗字母表示（一般情况，小写字母 $x$ 表示标量，粗体小写字母 $\\bf{x}$ 表示（列）向量，大写字母 $X$ 表示矩阵）。\n",
    "\n",
    "上面只是对某个样本数据的表示，假设有 k 个样本，那么可以写成矩阵形式：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}{\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "\\hat{y_1} \\\\\n",
    "\\hat{y_2} \\\\\n",
    "... \\\\\n",
    "\\hat{y_k}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "=\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "1 & x_{11} & ... & x_{1n} \\\\\n",
    "1 & x_{21} & ... & x_{2n} \\\\\n",
    "... & ... & ... & ... \\\\\n",
    "1 & x_{k1} & ... & x_{kn}\n",
    "\\end{array} \n",
    "\\right ]\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "a_0 \\\\\n",
    "a_1 \\\\\n",
    "... \\\\\n",
    "a_n\n",
    "\\end{array} \n",
    "\\right ]\n",
    "}\\end{aligned}\n",
    "$$\n",
    "\n",
    "分别用 $\\bf{\\hat{y}}$，$\\bf{a}$、$X$ 代表上面的矩阵，上式可以简写成：\n",
    "\n",
    "$$\n",
    "\\bf{\\hat{y}} = \\rm{X}\\bf{a}\n",
    "$$\n",
    "\n",
    "损失函数仍然和之前一样，使用平方损失：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "loss &= (y - \\hat{y})^2 \\\\\n",
    "&= \\sum_{i=0}^{k}(y_i - \\hat{y_i})^2 \\\\\n",
    "&= \n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "y_1-\\hat{y_1} & y_2-\\hat{y_2} & ... & y_k-\\hat{y_k}\n",
    "\\end{array}\n",
    "\\right ]\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "y_1-\\hat{y_1} \\\\\n",
    "y_2-\\hat{y_2} \\\\\n",
    "... \\\\ \n",
    "y_k-\\hat{y_k}\n",
    "\\end{array}\n",
    "\\right ] \\\\\n",
    "&=\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "y_1-\\hat{y_1} \\\\\n",
    "y_2-\\hat{y_2} \\\\\n",
    "... \\\\ \n",
    "y_k-\\hat{y_k}\n",
    "\\end{array}\n",
    "\\right ]^T\n",
    "\\left[ \n",
    "\\begin{array}{ccc}\n",
    "y_1-\\hat{y_1} \\\\\n",
    "y_2-\\hat{y_2} \\\\\n",
    "... \\\\ \n",
    "y_k-\\hat{y_k}\n",
    "\\end{array}\n",
    "\\right ] \\\\\n",
    "&= (\\bf{y}-\\bf{\\hat{y}})^T (\\bf{y}-\\bf{\\hat{y}}) \\\\\n",
    "&= (\\bf{y}-\\rm{X}\\bf{a})^T (\\bf{y}-\\rm{X}\\bf{a})\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "根据矩阵的基本计算法则：$ (A+B)^T = A^T+B^T $ 和 $ (AB)^T = B^TA^T $ 可得：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "loss &= (\\bf{y}-\\rm{X}\\bf{a})^T (\\bf{y}-\\rm{X}\\bf{a}) \\\\\n",
    "&= (\\bf{y}^T-\\bf{a}^T\\rm{X}^T) (\\bf{y}-\\rm{X}\\bf{a}) \\\\\n",
    "&= \\bf{y}^T\\bf{y} - \\bf{y}^T\\rm{X}\\bf{a} - \\bf{a}^T\\rm{X}^T\\bf{y} + \\bf{a}^T\\rm{X}^T\\rm{X}\\bf{a}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "求解线性回归，也就是求解系数 $\\bf{a}$ 使得 loss 函数取得最小值，为了求函数的最小值，这里也要用到之前说的求导数的方法，不过这里我们是对矩阵求导数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**矩阵求导（Matrix Derivative）** 也叫 **矩阵微分（Matrix Differential）**，在机器学习、图像处理、最优化等领域的公式推导中经常用到。关于矩阵求导的概念和原理这里不做过多的介绍，可以参考[这里](https://en.wikipedia.org/wiki/Matrix_calculus#Vector-by-vector_identities)。我们要求损失函数关于向量 $\\bf{a}$ 的导数：\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial loss}{\\partial \\bf{a}} \n",
    "&= \\frac{\\partial \\bf{y}^T\\bf{y}}{\\partial \\bf{a}} \n",
    "- \\frac{\\partial \\bf{y}^T\\rm{X}\\bf{a}}{\\partial \\bf{a}} \n",
    "- \\frac{\\partial \\bf{a}^T\\rm{X}^T\\bf{y}}{\\partial \\bf{a}} \n",
    "+ \\frac{\\partial \\bf{a}^T\\rm{X}^T\\rm{X}\\bf{a}}{\\partial \\bf{a}} \\\\\n",
    "&= 0 - \\rm{X}^T\\bf{y} - \\rm{X}^T\\bf{y} + 2\\rm{X}^T\\rm{X}\\bf{a} \\\\\n",
    "&= 2\\rm{X}^TX\\bf{a} - 2\\rm{X}^T\\bf{y}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "令 $ \\frac{\\partial loss}{\\partial \\bf{a}} = 0 $ 可得：\n",
    "\n",
    "$$\n",
    "\\rm{X}^T\\bf{y} = \\rm{X}^T\\rm{X}\\bf{a} \\Rightarrow \\bf{a} = (\\rm{X}^T\\rm{X})^{-1}\\rm{X}^T\\bf{y}\n",
    "$$\n",
    "\n",
    "其中，$\\rm{X}^T\\rm{X}$ 必须是可逆的，也就是说必须是**满秩矩阵**（又叫 **非奇异矩阵**）或**正定矩阵**。可以看出，通过一些简单的矩阵运算（转置、求逆、乘法）就可以求出损失函数的最小值，上面这个公式在统计学中被称为 **正规方程**（normal equation）。其中 $(\\rm{X}^T\\rm{X})^{-1}\\rm{X}^T$ 称为矩阵 $\\rm{X}$ 的 **伪逆矩阵**（pseudoinverse）或 **广义逆**（[Generalized inverse](https://zh.wikipedia.org/wiki/%E5%B9%BF%E4%B9%89%E9%80%86%E9%98%B5)），通常记为 $\\rm{X}^{\\dagger}$。\n",
    "\n",
    "$$\n",
    "\\bf{a} = \\rm{X}^{\\dagger}\\bf{y}\n",
    "$$\n",
    "\n",
    "通过这种方法求解线性回归，代码实现非常简便："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = [[62105.26315789]\n",
      " [   32.10526316]\n",
      " [ 1273.68421053]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.matrix([\n",
    "    [1, 1, 1, 1, 1], \n",
    "    [6000, 6000, 8000, 8000, 10000], \n",
    "    [60, 80, 70, 100, 90]]).T\n",
    "Y = np.matrix([340000, 350000, 400000, 450000, 500000]).T\n",
    "\n",
    "a = (X.T * X).I * X.T * Y\n",
    "print(\"a = {0}\".format(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "或者直接使用 numpy 自带的 `pinv()` 函数求伪逆矩阵："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a = [[62105.26315789]\n",
      " [   32.10526316]\n",
      " [ 1273.68421053]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "X = np.matrix([\n",
    "    [1, 1, 1, 1, 1], \n",
    "    [6000, 6000, 8000, 8000, 10000], \n",
    "    [60, 80, 70, 100, 90]]).T\n",
    "Y = np.matrix([340000, 350000, 400000, 450000, 500000]).T\n",
    "\n",
    "a = np.linalg.pinv(X).dot(Y)\n",
    "print(\"a = {0}\".format(a))"
   ]
  }
 ],
 "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
}