线性回归里的几何,优化和概率
方程求解
给定 Xw=y 求解 w,如果 X 是可逆的方阵,直接得到唯一的解: \( w=X^{-1}y \) 。
由于 X 可逆,它的列就作为列空间的基,因此 w 是列空间以 X 的列作为基时的坐标系数。
这是只使用代数语言下对求解方程 Xw=y 的核心解释。
列满秩过约束问题
如果 X 形状是 (m,n) ,且 m>n ,即它的行数大于列数(高瘦型)。
由于只有 n 列,但列空间里向量的长度是 m ,那么列空间的维度不可能大于 n, 即矩阵的秩 r<= n。
y 要么恰好在维度为 r 的列空间中,要么不在。
但 y 是属于 \( R^{n} \) 的,由于 r<=n<m, 某个点高维空间点恰好在特定低维平面的概率几乎为 0 ,于是 y 几乎不会在 X 的列空间。
这种情况方程 Xw=y 是无解的。
如果一般的 f(x)=y 方程无解,人们还能期望什么呢?
这往往不是数学本身能回答的,需要知道 f(x)=y 在解决什么具体的问题。
在数学之外,也就是对外部世界进行建模的思考中,每个线性方程被解释成一个约束,而行多于列意味着约束数量比变量数更多,这会导致约束无法被全部满足。
如果 f(x)=1 是一个逻辑的 sat 问题,我们可能要追求的是尽可能多得满足其中的逻辑子表达式。
一群因为风暴而卷入到荒凉海岛上的人,如果食物总量平分无法保证所有人都能活下去,那么他们可能就得决定一小部分人能得到食物。
以上两个例子都是把原本求解的问题用离散的方式放松了,既然无法满足所有约束,那么删掉一些,问题可能就有解了。
这是一种权衡,是对原问题的修改。
对于 Xw=y ,我们当然也能删除行使得约束越来越少最终有解,甚至可以每次随机删除一些方程,然后求解,最后把多次随机删除后得到解求平均,只要你能说服这种做法对于解决问题是合理的,那都可以尝试。
投影后的问题
在线性空间的几何视角下,Xw=y 无解意味着 y 不在 X 的列空间,比如 X 是 x-z 平面,而 y 是 (1,2,3) ,它不在该平面上。
由于这里本身就有几何的解释,那么可以从这种线性几何中找出一种“放松”问题的解释:
既然 y 不在 X 的列空间里(一个超平面),而是悬浮在平面之外,那么我们可以找 X 平面中和 y 最接近的向量。
线性空间里已经有对“距离”的自然定义,即 \( |x-y|^2 \), 它来自于现实中人类对物理空间的某种稳定感知。
即便不使用微积分,直觉上可以推测出点到平面的最短距离是该点到平面投影点的距离,而投影和正交有关
根据正交性质, y-Xv 应该垂直于投影 Xv 本身,即: \( (Xv)^T(y-Xv)=0 \)
于是得到了方程 \( v^TX^TXv = v^TXy \) 。
然而这个方程是约束不足的,我们只讨论了投影向量 Xv 和垂线 y-Xv 垂直,但投影的完整语义是 y-Xv 和整个列空间垂直,而 X 的列空间是由于它的列张成的,因此 y-Xv 只需要垂直于 X 的所有列即可:
\( X^T(y-Xv)=0 \) ,然后得到 \( X^Ty=X^TXv \) ;
这称为正规方程。
从代数形式上看,当我们无法求解 Xw=y 的时候,两边乘以一个 \( X^T \) 就自然得到了一个后备的方程。
向量微分补充
我们也可以用微积分去求解 \( \min |Xw-y|^2 \) 的问题,为此要做一点向量微积分的延伸:
给定方程 \( f(x) = ax \) ,可以从物理上的变化率解释去计算其导数 \( \lim \frac{a(x+\Delta x) - ax}{\Delta x} \), 结果是 a 。
注意这种解释是动态的,我们用极限这种语言将其描述出来,并得到方程。
而如果把它们扩展到向量,写成 \( f(x) = a^Tx \)
这里 a 是 n 个实数 [a1,a2,… ] 组成的列向量,x 是 n 个位置变量 [x1,x2,…] 组成的向量,那么 \( f(x) = a_1x_1+a_2x_2+\dots a_n x_n \)
对 x 中各个变量求偏导时,也可以回到物理解释,比如对 x1 求偏导,相当于固定 x2 到 \( x_n \), 然后给 x1 一点微小扰动,看 f(x) 的变化率,它就等于 a1 。
其他变量的偏导是并行的,只是最后把各个变量的偏导数打包成一个向量,称为梯度 \( \frac{\partial}{\partial{x}}(a^Tx) = a \) ;
通过方向导数的数学性质,可以从几何上解释出,某点的梯度是该点上扰动导致函数 f(x) 变化最大的方向。
类似地,可以对向量范数(一个标量)求梯度: \( \frac{\partial}{\partial{x}}(||x||^2_2) = 2x \) ;
这相当于并行地对 \( x_i^2 \) 求梯度,然后拼成列向量。
继续把 a 扩展到矩阵,有 \( f(x)=Ax \), 此时 f(x) 是一个向量,其中第 i 个元素是 A 的第 i 行和 x 的内积。
而这个元素针对 x 的梯度就是 A 的第 i 行的转置(因为我们约定梯度最终都写成列向量)。
于是我们还是可以定义广义的偏导,即对 f(x) 里每个元素求偏导,得到一个个梯度列矩阵,再按行拼起来: \( \frac{\partial}{\partial{x}}(Ax) = A^T \) ;
注意这里不会把 \( A^T \) 说成是梯度,因为梯度是针对针对标量函数的。
对于 \( f(x) = || Ax ||^2 \), 它是标量,于是可以谈论梯度,应用链式法则,它是:
\( \nabla f(x) = 2A^TAx \)
对于 \( ||Xw-y||^2 \) 对 w 的梯度就是 \( 2X^T(Xw-y) \) ,置 0 后得到正规方程。
另一种封闭解
如果 \( X^TX \) 可逆,那么 \( v = (X^TX)^{-1}X^Ty \) 就是需要求解的目标系数。
\( X^TX \) 可逆意味着不存在 0 特征根,即 \( X^TXv = \lambda v \) 中 \( \lambda \neq 0 \) ;
我们想知道什么情况下它是可逆的,可以两边乘以 \( v^T \) 得到 \( v^TX^TXv=\lambda v^Tv \)
即 \( (Xv)^2 = \lambda v^Tv \)
如果某个特征值 \( \lambda \) 为 0 ,说明向量 \( Xv \) 模长为 0, 而特征向量 v 从定义上要求是非 0 的。
因此如果 \( \lambda \) 为 0, 说明存在非 0 向量 v 使得 Xv=0 。
这意味着 X 的列不相互独立,也就是 X 的 rank 小于宽度 n 。
因此,如果在 X 是列满秩的情况下, \( X^TX \) 就是可逆的, \( v = (X^TX)^{-1}X^Ty \) 是合法的表达式。
将 X 的 SVD 分解 \( U\Sigma V^T \) 代入 \( v = (X^TX)^{-1}X^Ty \)
得到: \( v=(V\Sigma^T U^TU\Sigma V^T)^{-1} (V\Sigma^TU^T)y \)
整理后 \( v=(V \Sigma^T \Sigma V^T)^{-1} (V\Sigma^T U^T) y \)
继续整理为 \( v = V(\Sigma^T\Sigma)^{-1}V^T V \Sigma^T U^T y\)
因为 \( \Sigma^T\Sigma \) = \( \operatorname{diag} (\sigma_1^2,\dots,\sigma_n^2). \) 这是 (n,n) 形状的
它的逆等于 \( \operatorname{diag} (\frac{1}{\sigma_1^2},\dots,\frac{1}{\sigma_n^2}) \), 和 (n,m) 形状的 \( \Sigma^{-1} \) 相乘相当于给后者每行乘以上了对应的对角系数得到:
\( (\Sigma^T\Sigma)^{-1}\Sigma^T \) = \(diag \left( \frac1{\sigma_1}, \dots, \frac1{\sigma_n}, 0,\dots \right )\), 这是一个 nxm 的对角上是奇异值倒数的矩阵,可以写成 \( \Sigma^{-1} \) (被称为伪逆,因为它不是方阵不能求逆,但这是最接近逆的概念的矩阵了)
所以 \( (\Sigma^T\Sigma)^{-1}\Sigma^T = \Sigma^{-1} \)
因此 \( v= V\Sigma^{-1}U^Ty.\)
总结来看,如果 X 的列是相互独立的,那么 \( v = (X^TX)^{-1}X^Ty \) 是 Xw=y 的一个最近似的解。 其中 \( (X^TX)^{-1}X^T \) 经过 SVD 替换可以写成 \( V\Sigma^{-1}U^T \) ,即还有另一种形式的封闭解 \( v= V\Sigma^{-1}U^Ty\) 。
\( (X^TX)^{-1}X^T \) 称为 X 的左逆矩阵,它左乘 X 就是标准的单位矩阵,但右乘 X 则是 \( X(X^TX)^{-1}X^T \), 它是到 X 列空间的投影矩阵 P ,因为 \( Py = X(X^TX)^{-1}X^Ty \), 它等于 \( Xv \) ,也就是我们求的 y 的投影,另外容易验证 \( P^2=P \) 。
因为 \( (X^TX)^{-1}X^T \) 只在左乘 X 的时候得到单位矩阵,所以被称为 X 的左逆。
回归和优化问题
到目前为止还没有引入回归,纯粹是在相对抽象的线性空间几何解释下谈论问题,通过投影把目标映射到 X 的列空间去求解。 如果 X 的列相互独立,那么代数上 \( v = (X^TX)^{-1}X^Ty \) 就可以求解,或者引入 SVD 分解得到另一种封闭解形式: \( v= V\Sigma^{-1}U^Ty.\)
要切换到“回归”视角,需要从行视角入手,Xw=y 的每行是一个线性约束,针对的是单个样本 x ,每一行的 X[i] 和 y[i] 都是样本的一个属性,但 y 是要预测的目标属性,而 x 则是预测的证据。
w 被解释成了 X[i] 和 y[i] 之间的线性规律,通过这个规律,给定 X[i] 就能去预测 y[i] 。
这个规律写成数学则是 \( y_i=w^Tx \) ;这里小写的 x 指代 X 中的行。
对于每个一个行数据都想使得 \( w^Tx \) 预测的结果尽可能接近对应的真实预测值 \( y_i \) 。
而一种衡量“接近”的数学方式是残差的平方和,即最小化 \( \mid Xw-y\mid^2 \) 。
求梯度后置为 0 的方程就是: \( X^T(Xw-y)=0 \)
同样得到了正规方程 \( X^TXw = X^Ty \)
所以回归问题额外出现了误差最小化的概念,这也是几何意义上的,但和向量空间里的几何形式不同,内在逻辑则是一致的。
最小值的证明
前文我们用直觉说 y 距离 X 列空间最近的点是它在 X 的投影,并没有证明,但当我们有了一个最小化的数学目标 \( \mid Xw-y\mid^2 \) 时,可以通过二阶导数来判断 0 临界点的凸性。
而 \( \mid Xw-y\mid^2 \) 对 w 求二阶导数得到的是 \( X^TX \), 通过能量法去判断 \( v^TX^TXv = (Xv)^2 \) ,由于 X 是列满秩的,Xv 对任何非 0 向量 v 都是非 0, 因此始终大于 0, 即优化目标的 Hessain 矩阵 \( X^TX \) 是正定矩阵。于是它是最小值。
也可以不通过线代的高层封装,直接对底层的 \( w_i \) 求二阶导:
\[ \frac{\partial Err}{\partial w_i} = \sum_{i}^{n}2(w_ix_{i}-y_{i})x_{i}\]
二阶偏导数: \[ \frac{\partial^2 Err}{\partial w_i^2} = \sum_{i}^{n}2x_{i}^{2}; \]
对于不同的 \( i,j \) \[ \frac{\partial^2 Err}{\partial w_i \partial w_j} = \sum_{i}^{n} 2x_{i} = \frac{\partial^2 E}{\partial w_j \partial w_i} \]
将以上结果排列成 Hessian 矩阵就是 \( X^TX \) 形式,于是可以用能量法去判断正定性质。
偏置项 b 和特征工程
线性回归模型更常见的形式是 \( y_i=w^Tx_i + b\) ,这有一个偏置 b,此时 Xw=y 中 X 的某一列(一般是最后一列)会是全 1 ,X 是一个被填充后的对象。
然而是否有 b 并不是线性回归问题的核心,而是建模时对模型“好坏”的选择偏好,如果你观察到的数据总体就是能被过原点的直线很好地近似表示,那么可以不需要偏置。
甚至你观察到数据可能要用 \( x^2 \) 而不是 \( x \), 那么 \( y=ax^2+b \) 也是回归,此时 X 中存储的就是各个元素的平方项以及一列纯 1 项的拼接,它仍然是 Xw=y 的问题,对 X 的变换其实可以看作一种特征工程,而这类工程是由于具体问题里数据的语义和预测目标所决定的(特定领域知识)。
特征工程属于传统机器学习的范畴,其核心是已经把线性回归需要的矩阵 X 构造出来,然后就能得到正规方程以及对应的封闭解。
而深度学习是将 X 的构造也参数化,整个学习过程中梯度不会止步于 X, 因此需要继续链式法则,这也意味着无法找到一个封闭的解,只能依赖于梯度下降等数值优化方法。
损失函数的问题
我们为什么选择了 \( |Xw-y|^2 \) 而不是 \( |Xw-y|_1 \) 或者 \( \max_i |(Xw-y)_i| \) 等其他形式?
从数学一侧来看,平方函数光滑且凸,有很好的计算性质,线性回归里直接导出了正规方程及其封闭解。
同时平方误差对应欧氏空间中两个点(向量)之间的距离,当预测目标本身具有空间位置含义时(比如预测路上的车辆),这种度量是自然的。
因此,这种工具被数学家抛出,至于使用者为什么要选择它,那更多是一种经验判断,并且是实践中尝试出来的。
很多时候由于好的数学性质,人们不需要依托强有力的解释,比如预测房价问题中,价格之间的平方损失并不如物体之间距离那么自然,但还是可以用平方损失,因为这样做好算,尽管不是精确的,但至少能拿出一个结果。
甚至我们会因为这个模型本身数学性质好,而尽量把现实问题的数据转成更符合这种模型的形式,比如改变量纲,归一化等。
引入统计后,会有一些新的解释,比如如果观测数据是在真实规律基础上叠加独立高斯噪声产生的,那么最小化平方误差等价于最大化观测数据出现的概率(最大似然估计),后文再详细展开。
行满秩欠约束问题
现在关注欠约束的问题:矩阵 X 的秩 r 小于宽度 n, 这里先假设其高度 m 小于 n ,并且 r=m, 即行满秩,这和前文的列满秩问题完全对偶。
此时线性约束比自由变量少,会有无穷多个解。
在过约束中,我们面临的选择是如何找到一个新的问题,使得它和原始问题最接近。这里遇到的是,如何增加一个强约束,从无数个解中选择一个“最好”的解。
最常见的是最小化 Xw=y 中 w 的模长,即 \( \min_w |w|^2 \), 同时满足 Xw=y 。
为什么选择它呢?
和为什么选择投影来处理过约束的情况类似,数学家本身尝试过很多度量标准,但 \( |w|^2 \) 本身和欧式距离是相容的,并且马上我们会知道,它的计算性质和最小二乘一样好,能得到封闭解。
而具体问题中是否要使用这个标准,是经验问题,需要在实践中去检验。
当然可以引入一些理论上的解释,赋予 Xw=y 一些具体“意义”,比如 w 被解释为模型参数,它用来描述外部世界的某些对象之间的关系,比如房子的面积和所在地区对房价的影响。
然后哲学里的奥卡姆剃刀原则就可以介入,其声称所有解释能力相同的模型中,最简单的那个最好。
但这并没有太多说服力,即便它在实践中表现良好,可以作为一个经验值准则,但 w 的二范数小并不意味着模型简单, (0.2123 0.101) 两个数都很小,但这两个数的复杂度(或者说信息量)比 (5,5) 更高。
所以,在缺乏具体的问题背景时,还是回到数学本身,和最小二乘选择平方损失类似,最小化二范数会使得计算上变得简单,比如它是可微的,凸的有封闭解。
求解 \( \min_w |w|^2, s.t. Xw=y \) 可以用拉格朗日法, \( w^2+\lambda^T(Xw-y)=0 \)
对 w 和 \( \lambda \) 求导,会得到 \( 2w + X^T\lambda = 0 \) 和 \( Xw=y \)
因此有: \( w=-\frac12 X^T\lambda \) ,代入到 Xw=y 有:
\( -\frac{1}{2}XX^T \lambda = y \)
由于假设是行满秩的,那么 \( XX^T \) 可逆,于是有 \( \lambda = -2(XX^T)^{-1}y \) ,代入到 w 的表达式有
\( w = X^T(XX^T)^{-1}y \)
其中 \( X^T(XX^T)^{-1} \) 乘到 X 右边后会得到单位矩阵,乘到 X 左边得到 \( X^T(XX^T)^{-1}X \), 这是 X 行空间的投影矩阵。 所以它称为 X 的右逆。
这个性质和前文最小二乘 \( \min_w |Xw-y|^2 \) 得到的 \( w=(X^TX)^{-1}X^Ty \) 非常对偶。
而把 \( X=U\Sigma V^T \) 代入右逆得到 \( V\Sigma^TU^T (U \Sigma \Sigma^T U^T)^{-1} \) = \( V\Sigma^TU^TU(\Sigma\Sigma^T)^{-1}U^T \) = \( V\Sigma^{-1}U^T \)
这个结果和最小二乘中左逆代入 SVD 后结果完全一样。但它们分别来自两个虽然相似但完全不同的公式的变换,一个要求列满秩,一个要求行满秩。
即便行和列都不满秩,由于任何矩阵 X 都可以进行 SVD 分解,因此 \( V\Sigma^{-1}U^T \) 永远存在,这被称为伪逆。
它给我们一个提示,任何 Xw=y 都可以得到一个 \( w^*=V\Sigma^{-1}U^T \) 形式的解。
在列满秩但列的数量小于行的数量情况下,X 的列无法完全覆盖整个 m 维空间,因此会错过 y, 于是寻找 y 的在 X 上的投影作为替代。
在行列都不满秩的情况下,由于列空间的维度只有 r, 它小于 m, 于是 y 有可能不在 X 的列空间中,于是需要用 y 到 X 列空间的投影来作为 y 的代理。
然而由于秩小于列数 n, 因此 X 的列是线性相关的,即用它作为基是有冗余的, 因此列空间中的元素 y' 可能存在多种不同的表示 Xw1=y' 和 Xw2=y' 。
因此要选出其中 2 范数最小的作为最终解。
因此,对于一般的 X, \( w^*=V\Sigma^{-1}U^T \) 是一个通用的解形式,它天然进行两个阶段的操作:投影并选择最小范数的 w 。
如果 X 是列满秩,那么第二步是唯一的,不需要选 w. 如果是行满秩,第一步投影是一个恒等变换,因此伪逆还编码了某种自动的 if else 判断。
不确定性和岭回归
如果只从数学的角度看, Xw=y 无解就可以去找投影向量,投影后还有多个解,那么找最小范数。
但从建模的角度看, X 和 y 都 是人搜集到的数据,它对于要解决的问题本身来说是有噪声的,比如人们记录的数据有笔误,或者测量结果就是有误差。
如果用伪逆去求解,那么相当于把不准确地 y 精确地投影到了不准确的 X 的列空间,然后再选出最小的 w ,这每一步都可能放大了最初的噪音。
伪逆的解 \( w^* = V \operatorname{diag} \left( \frac1{\sigma_i} \right) U^Ty \) 中,如果某个奇异值 \( \sigma_i \) 很小,那么解的数值会爆炸,而这个奇异值可能在精确的 X 下为 0 (为 0 的奇异值不会被求倒数),是噪声导致它有微小值,这里就把噪声急剧放大了。
因此一种思路是打破伪逆解形式中蕴含的两阶段优化方案,直接把两个要优化的目标融合成一个,这就得到岭回归:
\[ \hat{w} = \arg\min_{w} \left\{ \| \mathbf{y} - \mathbf{X}w \|_2^2 + \lambda \| w \|_2^2 \right\}, \quad \lambda > 0 \]
展开损失: \[ L(w) = (\mathbf{y} - \mathbf{X}w)^T(\mathbf{y} - \mathbf{X}w) + \lambda w^T w \]
w 求梯度并置零: \[ \frac{\partial L}{\partial w} = -2\mathbf{X}^T(\mathbf{y} - \mathbf{X}w) + 2\lambda w = 0 \]
整理得: \[ (\mathbf{X}^T\mathbf{X} + \lambda \mathbf{I}) w = \mathbf{X}^T \mathbf{y} \]
最终封闭解为: \[ \boxed{\hat{w} = (\mathbf{X}^T\mathbf{X} + \lambda \mathbf{I})^{-1}\mathbf{X}^T \mathbf{y}} \]
它和左逆解的区别仅仅在于在 \( X^TX \) 上增加了一个 \( \lambda I \), 将 SVD 分解替换 X:
\[ \hat{w} = (\mathbf{V} \boldsymbol{\Sigma^T \Sigma} \mathbf{V}^T + \lambda \mathbf{I})^{-1} \mathbf{V} \boldsymbol{\Sigma} \mathbf{U}^T \mathbf{y} \]
\( \Sigma^T\Sigma \) 写成 \( \Sigma^2 \) ,并利用正交性 \( \mathbf{I} = \mathbf{V}\mathbf{V}^T \) 和矩阵乘法结合率: \[ \mathbf{V} \boldsymbol{\Sigma}^2 \mathbf{V}^T + \lambda \mathbf{I} = \mathbf{V} (\boldsymbol{\Sigma}^2 + \lambda \mathbf{I}) \mathbf{V}^T \]
其逆为: \[ (\mathbf{V} (\boldsymbol{\Sigma}^2 + \lambda \mathbf{I}) \mathbf{V}^T)^{-1} = \mathbf{V} (\boldsymbol{\Sigma}^2 + \lambda \mathbf{I})^{-1} \mathbf{V}^T \]
得到: \[ \hat{w} = \mathbf{V} (\boldsymbol{\Sigma}^2 + \lambda \mathbf{I})^{-1} \boldsymbol{\Sigma} \mathbf{U}^T \mathbf{y} \]
展开中间的对角矩阵: \[ \boxed{\hat{w} = \mathbf{V} \cdot \text{diag}\left( \frac{\sigma_i}{\sigma_i^2 + \lambda} \right) \cdot \mathbf{U}^T \mathbf{y}} \]
它和普通二阶段优化的线性回归形式的区别仅仅在于把 \( \text{diag}(\frac{1}{\sigma_i}) \) 变成了 \( \frac{\sigma_i}{\sigma_i^2 + \lambda} \)。
此时即便 sigma 是噪音或固有的很小的值,选择适当的 \( \lambda \) 也可以起到平滑的作用,使得噪声不会被倒数放大,更具体的:
- 当 \( \sigma_i \gg \lambda \) 时,回到一般线性回归解。
- 当 \( \sigma_i \ll \lambda \) 时,\( \frac{\sigma_i}{\sigma_i^2+\lambda} \approx \frac{\sigma_i}{\lambda} \to 0 \) ,这意味着,对于很小的 \( \sigma \), 原本被急剧放大,这里则被压制,可以解释为忽视权重过小的特征。
概率统计语言
从一般的 Xw=y 到岭回归,我们开始考虑模型和数据都可能存在误差的问题,不再把方程看作必须严格满足的约束,而是同时去拟合误差和选择最优解。
不过,这一阶段仍然没有真正引入数学意义上的概率。所谓"不确定性",更多是一种建模思想,而不是一个具有严格数学定义的随机过程。
当随机变量及其概率分布被引入后,可以进一步假设观测值并不是固定产生的,而是由某个随机机制生成:
\( y=Xw+\epsilon, \)
其中 \( \epsilon \) 是具有给定分布的随机变量。
这里建模对象发生了变化。我们不再把数据看成一组孤立的几何点,而是某个随机过程的一次观测结果。模型需要描述的不只是数据之间的空间关系,还包括这些数据为什么会以这样的方式出现。
这里概率模型提供了“一次性”的数据生成的视角,我们可以观察到外部世界中某些数据的产生过程,然后把这个产生过程的信息也近似地映射到数学中去,这比直接观察到静态的空间结构然后映射成几何语言更丰富一些。
一个类比是,如果只能看到夜空中某一时刻星星的位置,我们最多只能根据空间分布进行分类,例如划分星座,这是一种静态的几何描述。
当持续观测星体的位置变化时,我们开始能够讨论轨道、速度以及周期等运动规律。但这些仍然只是对现象本身的几何描述,它告诉我们"发生了什么",也可以用这种轨迹去预测,即谈论物体接下来会如何运动,但却没有回答"为什么会发生"。
等到牛顿提出万有引力定律后,人们开始把这些运动理解为某种普遍规律的结果。虽然引力本身为什么存在仍然未知,但至少各种不同的天体运动可以统一地解释为同一种作用规律的表现,人们可以用这个规律去谈论行星为什么运动。
概率统计是在牛顿定律层面吗?
并不是,它是在描述运动轨道的层面构造的一种描述信息生成的语言,而不是描述确定性机制的语言,比如测量行星轨迹的时候,第一天测量了 1.001, 第二天 1.002, 第三天 0.997. 概率能给这种不确定性建立一个函数描述,不同的值对应一个生成的概率,但它不会解释为什么会有一个如此稳定的轨道,也不会解释为什么有噪声,这些噪声可能来自各方面,比如抄写误差,设备误差,但它们都被抽象掉变成数据之间信息的规律。
不过,仅仅这一层抽象就已经极大增强了我们描述现实的能力。很多时候,我们并不需要知道世界为何如此运行,也不需要建立完整的物理机制,只要能够准确描述观测数据之间的统计规律,就足以完成预测、估计和决策。例如房价预测、图像识别等问题,它们真正利用的是信息层面的“运动”规律,而不是底层物理机制。
回归的概率解释
有了这种能描述“信息动态”的语言,我们看它能如何重新去描述最小二乘法。
我们把观测值 y 不再看作确定性的计算结果,而是一个随机变量,给定输入 \(X\) 后,它服从某个条件概率分布 \(p(y|X; w)\)。
而这个分布可以用以下更细的公式决定:
\[ y_i = x_i^T w + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma^2) \]
在真实线性规律 \(x_i^T w\) 的基础上,叠加了一个独立的高斯(正态)噪声 \(\epsilon_i\) ,它不解释噪声来自抄写误差还是设备误差,只是从信息层面统一描述为围绕均值的随机波动。
\(y_i\) 在给定 \(x_i\) 和参数 \(w\) 后,就服从一个以 \(x_i^T w\) 为中心的高斯分布:
\[ p(y_i | x_i; w) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - x_i^T w)^2}{2\sigma^2} \right) \]
对于全部 \(m\) 个独立观测样本,观测数据出现的联合概率(即似然函数 \(L(w)\))为各点概率的乘积:
\[ L(w) = \prod_{i=1}^{m} p(y_i | x_i; w) = \prod_{i=1}^{m} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - x_i^T w)^2}{2\sigma^2} \right) \]
为了计算方便,取对数得到对数似然(log-likelihood):
\[ \ell(w) = \log L(w) = \text{C} - \frac{1}{2\sigma^2} \sum_{i=1}^{m} (y_i - x_i^T w)^2 \]
要最大化观测数据出现的概率,等价于最小化 \(\sum_{i=1}^{m} (y_i - x_i^T w)^2\) ,也就是最小二乘目标函数 \(\|y - Xw\|^2\)。
注意这里 \( \epsilon_i \) 的方差可以是任意的,比如是 \( \lambda \) ,只要保证每个样本上的方差都相同。
岭回归的概率解释
在贝叶斯统计视角下,我们不仅关注数据,还关注参数 \(w\) 自身出现的可能性。我们可以给 \(w\) 赋予一个先验分布(prior)。假设 \(w\) 的每一个分量都独立地服从均值为 0 方差为 1 的高斯分布,w 的联合概率分布就是
\[ p(w) = \prod_{j=1}^{n} \frac{1}{\sqrt{2\pi\tau^2}} \exp\left( -\frac{w_j^2}{2} \right) = \text{C} \cdot \exp\left( -\frac{1}{2} \|w\|^2 \right) \]
且模型为 \( Y=x^TW+\epsilon_{\lambda} \), 这里 \( \epsilon_{\lambda} \) 是方差为 \( \lambda \) ,期望为 0 的高斯分布。
根据贝叶斯定理,看到数据后 w 的后验概率(posterior)为:
\[ p(w | X, y) \propto p(y | X; w) \cdot p(w) \]
取对数后:
\[ \log p(w | X, y) = C - \frac{1}{2\lambda } \|y - Xw\|^2 - \frac{1}{2} \|w\|^2 \]
最大化它就等于最小化 \[ \|y - Xw\|^2 + \lambda \|w\|^2 \]
这就是岭回归的目标函数。
概率解释的灵活性
至此,我们用两个特定的概率假设模型分别覆盖了一般的回归模型和岭回归:
- 假设在线性机制 Xw 上叠加了一个独立的同分布的标准高斯噪声就得到一般回归
- 继续给参数赋予一个先验分布,就得到了岭回归
因此从表达能力上,概率模型可以重新解释前文提到的几何结果,只要稍微一变化,比如你观察到噪声的分布不是平滑的钟型,而有点尖锐,那么可能会选用 Laplace 噪声 \( p(\epsilon) ~ e^{-\mid \epsilon \mid} \), 最大似然概率就变成了最小化 Xw-y 的 L1 范数。