线性回归与最小二乘法

线性回归与最小二乘法

最小二乘法是我们经常用到的求解模型的回归方法，是一种优化方法，对于一个线性模型

$$X\theta=y$$ $X$为$m$行 $n$列的矩阵，代表$m$个样本，每个样本有 $n$个变量，即有$m$个等式，未知数 $n$个，$m>n$， $\theta$和 $y$为一维向量，该方程没有解，但可以选择适合的$\theta$使等式尽可能接近，如果以残差平方和评价，代价函数为 $$loss=||X\theta - y||^2$$ 如何使代价最小是一个优化问题，可以通过梯度下降法求解，也可以通过最小二乘求解，最小二乘解为 $$\hat\theta = (X^TX)^{-1}X^Ty$$ 最小二乘解是未知变量的无偏估计，但当X不是列满秩或某些列之间的线性相关性比较大时， $X^TX$的这些行相同或线性相关性也比较大，使行列式接近于0，即 $X^TX$接近于奇异，此时 $X^TX$条件数很大，计算 $(X^TX)^{-1}$时的误差会很大，其解缺乏稳定性和可靠性。X矩阵中不同列之间具有完全或高度相关关系叫做多重共线性，这种情况下用最小二乘求得的解释不可靠的，一言蔽之，模型无偏，方差大，解决办法有：

1. 找到引起多重共线性的变量，将其排除；
2. 使用逐步回归法，逐步添加和排除变量，使最终使用的变量集合最优；(最常用，最有效)
3. 通过加入扰动(正则项)，增加解的稳定性，即岭回归。
4. 通过主成分分析得到新的变量集合。

可以用最小二乘法求解的线性回归模型为

$$y=X\theta + b$$ 对于y的取值为{0,1}或在[0,1]区间的问题，可以用逻辑回归模型建模 $$y=\sigma(X\theta + b)$$ 其中 $\sigma$是sigmoid函数

岭回归

岭回归(ridge regression)，是线性模型一种有偏回归方法，相比最小二乘法增加了偏差，减小了方差，将最小二乘法中的$loss$增加l2正则项（二范数的平方$||\theta||_2^2$，向量内积，向量所有元素的平方和$\sum_{i=1}^n\theta_i^2$），损失函数为

$$loss=||X\theta - y||^2 + \alpha\sum_{i=1}^n\theta_i^2$$ 岭回归的解为 $$\theta(\alpha) = (X^TX + \alpha I)^{-1}X^Ty$$ 其中 $I$为单位矩阵，随着$\alpha$的增大， $\theta$的各项相对于正确值的偏差越来越大， $\alpha$无穷大时， $\theta$趋近于0，可以取多组不同的 $\alpha$，画出每个 $\theta$分量的变化轨迹，也叫岭迹图，借助岭迹图可以确定最佳 $\alpha$。

Lasso回归：(Least Absolute Shrinkage and Selection Operator)，在最小二乘法的损失函数中添加l1正则项（一范数$||\theta||_1$，向量所有元素的绝对值之和$\sum_{i=1}^n|\theta_i|$），损失函数为

$$loss=||X\theta - y||^2 + \alpha\sum_{i=1}^n|\theta_i|$$ Lasso回归的特点，具有特征选择性，对于不重要的特征，对应的 $\theta$容易收缩为0，比岭回归更快速。lasso回归的解法主要有坐标轴下降法和最小角回归法。
坐标轴下降法是一种非梯度优化方法，与梯度下降法沿一个特定方向迭代不同，在每次迭代过程中按一定的次序分别沿着各个坐标轴方向搜索极小值。坐标轴下降法在稀疏矩阵上的计算速度非常快，同时也是lasso回归最快的解法。
最小角回归法是前向选择算法和前向梯度算法的折中版本，详情看 Lasso Regression，利用lasso回归也可以做特征选择。Lasso惩罚项超参数可以基于交叉验证，或者基于信息准则。

ElasticNet回归：同时包含l1和l2正则项的最小二乘法，损失函数为

$$loss=||X\theta - y||^2 + \alpha\sum_{i=1}^n|\theta_i| + \beta\sum_{i=1}^n\theta_i^2$$ 另外，可以将ridge,lasso和elasticNet表示成带约束的优化问题 $$\min||X\theta - y||^2, \qquad st.\sum_{i=1}^n\theta_i^2 \le s$$ $$\min||X\theta - y||^2, \qquad st.\sum_{i=1}^n|\theta_i| \le s$$ $$\min||X\theta - y||^2, \qquad st.\sum_{i=1}^n|\theta_i| \le s1 \space and \space \sum_{i=1}^n\theta_i^2 \le s2$$

R线性相关系数：描述两个变量相对于各自均值同增减/浮沉的关系，接近1或-1代表强烈的正/负相关关系，接近0表示基本没有相关关系，或者说是独立的。

$$R=\frac{Cov(X,Y)} {\sqrt{Var(X)}\sqrt{Var(Y)}}$$

R2可决系数(coefficient of determination)：也叫判定系数或拟合优度，如果只有一个自变量，是自变量与因变量的相关系数的平方，总平方和TSS (total sum of squares)

$$SS_{tot}=\sum_i(y_i-\bar y)^2$$ 回归平方和RSS (regression sum of squares) $$SS_{reg}=\sum_i(f_i-\bar y)^2$$ 残差平方和ESS (error sum of squares) $$SS_{res}=\sum_i(y_i-f_i)^2$$ RSS在TSS中所占的比重称为可决系数，可决系数可以作为综合度量回归模型对样本观测值拟合优度的度量指标。可决系数越大，说明在总平方和中由模型作出了解释的部分占的比重越大，模型拟合优度越好。反之可决系数小，说明模型对样本观测值的拟合程度越差。R2可决系数最好的效果是1，也可能是过拟合。 $$R2=\frac {RSS}{TSS} = 1-\frac{ESS}{TSS}$$

主成分分析

主成分分析(Pricipal Component Analysis)是一种数据降维方法，作用是去除噪声和冗余。PCA将原来的特征通过正交线性变换形成新的特征，使得新特征的第一大方差在第一坐标上(第一主成分)，第二大方差在第二坐标上(第二主成分)，以此类推。PCA变换会损失了一部分信息，但主成分上方差大的特性使样本点尽可能分散，以此保留更多的信息。其计算过程如下：

1. 去中心化，减去均值$X \leftarrow X-X_{mean}$，$X$为$m$行$n$列的矩阵
2. 求原始数据协方差矩阵$Cov = \frac1{m-1}(X^TX)$
3. 计算协方差矩阵的特征值$\lambda_i$和特征向量$w_i$
4. 取最大的$N$个特征值，其对应的特征向量形成标准正交基$W=(w_1,w_2,...,w_N)$,$n$行$N$列，每列一个特征向量
5. 进行线性变换$Y=XW$

通过以上过程，$n$维的$X$变成了$N$维$Y$($N<n$)，$N$的取值可以固定，一般使$N$个特征值的和与所有特征值的和之比超过一定阈值$t(t<1.0)$，如$t=0.85$，t越大，保留的信息越多。

PCA原理推导（为什么要用协方差矩阵的特征值和特征向量）

PCA的目标是求一种正交变换$Y=XW$，使变换后各个维度的方差和最大，满足$W^TW=I$，这是一个优化问题

$$\max tr(W^TX^TXW),\qquad st.W^TW=I$$ 因为 $X^TX$是半正定对称矩阵，所以用拉格朗日乘子法求解 $$L = tr(W^TX^TXW) + \lambda(I-W^TW)$$ 取极值的条件是一阶导数等于0，而且该问题是凸函数，极值就是最值。 $$\frac{\partial L}{\partial W} = 2X^TXW-2\lambda W = 0$$ 因此 $\lambda$是 $X^TX$的特征值， $W$是$X^TX$的特征向量，另外，由于 $tr(W^TX^TXW)=tr(W^Tdiag(\lambda_1,\lambda_2,...,\lambda_n) W) = \sum_{i=1}^n\lambda_i$，因此特征值 $\lambda$越大，该特征的方差越大。推导过程用到了迹的求导公式 $$d(tr(AXBX^T))=AXB+A^TXB^T$$ 其中 $A=X^TX,B=I$

矩阵的迹及其求导公式

矩阵的迹trace是矩阵对角线元素的和

$$tr(A) = \sum_{i=1}^nA_{ii}$$ 矩阵的迹及其导数有以下性质，以下 $d(tr(expr))$表示对表达式中的 $X$求导
- $a \in R,\space tr(a)=a$
- $tr(A)=tr(A^T)$
- $tr(AB)=tr(BA)$
- $tr(ABC)=tr(CAB)=tr(BCA)$
- $d(tr(X))=I$
- $d(tr(XB))=d(tr(BX))=B^T$
- $d(tr(X^TB))=d(tr(BX^T))=B$
- $d(tr(AXB))=A^TB^T$
- $d(tr(AX^TB))=BA$
- $d(tr(AXBX^T))=AXB+A^TXB^T$
- $d(tr(AXBX))=BXA+AXB=A^TX^TB^T+B^TX^TA^T$
注意：标量对矩阵求导的结果是矩阵，比如雅克比矩阵

如何证明上述定理，最好理解的办法是将迹展开成求和形式，例如，求该式的导数

$$f=tr(AXBXC^T)=\sum_i \sum_j \sum_k \sum_l \sum_m A_{ij} X_{jk} B_{kl} X_{lm} C_{im}$$ 第一种证明，按角标求： $$\frac{\partial f}{\partial X_{jk}} = \sum_i \sum_l \sum_m A_{ij} B_{kl} X_{lm} C_{im} = [BXC^TA]_{kj}$$ $$\frac{\partial f}{\partial X_{lm}} = \sum_i \sum_j \sum_k A_{ij} X_{jk} B_{kl} C_{im} = [C^TAXB]_{ml}$$ $$\frac{\partial tr(AXBXC^T)}{\partial X} = A^TCX^TB^T + B^TX^TA^TC$$ 第二种证明，利用微分乘法法则 + 迹的基本求导公式： $$\frac{\partial tr(AXBXC^T)}{\partial X} = \frac{\partial tr(AXD)}{\partial X} + \frac{\partial tr(EXC^T)}{\partial X} = A^TD^T + E^TC$$ 其中 $D=BXC^T, E=AXB$，该方法十分方便例如求下式关于X的导数，对于复杂的算式，依次将每个 $X$拆分成$AXB$或 $AX^TB$的形式，运用相应的求导公式，将结果相加即可： $$f = tr(AXX^TBCX^TXC)$$

方差的无偏估计

在PCA求协方差时，前面的系数是$\frac1{m-1}$而非$\frac1m$，实际上用的是方差的修正公式，而我们熟知的方差计算公式是有偏的，符号说明：随机变量$X$，$n$个样本，各样本为$x_i$，总体均值$\mu$，总体方差$\sigma^2$，样本均值$\bar X$，样本方差$s^2$，具体推导如下

$$\begin{alignedat}{1} s^2 &=\frac1n\sum_{i=1}^n (x_i - \bar X)^2\\ &=\frac1n\sum_{i=1}^n ((x_i-\mu) + (\mu-\bar X))^2\\ &=\frac1n\sum_{i=1}^n (x_i-\mu)^2 + \frac2n\sum_{i=1}^n(x_i-\mu)(\mu-\bar X) + \frac1n\sum_{i=1}^n(\mu-\bar X)^2\\ &=\frac1n\sum_{i=1}^n (x_i-\mu)^2 - (\bar X - \mu)^2\\ &\le\frac1n\sum_{i=1}^n (x_i-\mu)^2\\ E[s^2] &= E[\frac1n\sum_{i=1}^n (x_i-\mu)^2] - E[(\bar X - \mu)^2]\\ &= D[X] - D[\bar X]\\ &= \sigma^2 - \frac1n\sigma^2\\ &= \frac{n-1}n\sigma^2 \end{alignedat}$$ 化简过程中 $\bar X$和 $\mu$与索引i无关，而且 $$D[\bar X]=D[\frac1n\sum_{i=1}^nx_i]=\frac1{n^2}D[\sum_{i=1}^nx_i]=\frac1{n^2}nD[X]=\frac1nD[X]$$ 从上面可以看出，通过样本均值直接计算得到的样本方差的期望比总体方差小，只要样本均值不等于总体均值，我们总会低估总体的方差，因为样本均值就是使样本方差最小的那个值，为了使方差的估计无偏，就要使用修正后的方差公式 $$s^2 = \frac1{n-1}\sum_{i=1}^n(x_i-\bar X)^2$$

参考链接：

岭回归（ridge regression）
Lasso Regression
【机器学习】汇总详解：矩阵的迹以及迹对矩阵求导
机器学习中对矩阵的迹（trace）求导的一些操作
为什么样本方差里面要除以（n-1）而不是n？
彻底理解样本方差为何除以n-1