跳到正文
格致开物
MATHWIKI
探索
学科导航
学习路径
搜索
☾
登录
探索
学科导航
学习路径
随机漫游
希腊字母
关于本站
管理员登录
搜索
数学百科
/
知识地图
查看“︁最小二乘法”︁的源代码
←
最小二乘法
因为以下原因，您没有权限编辑该页面：
您请求的操作仅限属于这些用户组的用户执行：
管理员
、aipublisher
您可以查看和复制此页面的源代码。
'''最小二乘法'''（least squares）通过最小化观测值与模型预测值之间的残差平方和，选择模型参数。它既是一种拟合方法，也可以解释为在允许的预测中寻找离观测最近的点。 下面用四个水位数据拟合一条直线，逐步推导要最小化的量、系数方程及几何解释。数据是合成教学数据，与[[数学建模]]条目的复现附件一致。 == 一条直线怎样接近四个数据点 == 时间为 <math>t=(0,1,2,3)</math> 分钟，水位为 <math>h=(1,2,2,4)</math> 厘米。选择直线模型 <math display="block">\hat h=b+at,</math> 其中 <math>b</math> 是初始水位估计，单位为厘米；<math>a</math> 是每分钟的水位变化量，单位为厘米每分钟。 四个点不在同一直线上。例如前两个点的连线为 <math>\hat h=1+t</math>，它在第三个时间 <math>t=2</math> 预测三厘米，而观测为二。因此不能要求一条直线同时精确经过全部点，只能规定一种比较拟合好坏的方法。 对每个点，采用“观测减预测”的残差约定： <math display="block">r_i=h_i-(b+at_i).</math> 直接把残差相加会相互抵消。若取 <math>b=a=1</math>，残差为 <math>(0,0,-1,0)</math>；若取 <math>b=a=0.9</math>，残差为 <math>(0.1,0.2,-0.7,0.4)</math>，后者和为零，却仍没有精确经过各点。 平方后再相加，每个误差都贡献非负数。本例的目标为 <math display="block">S(b,a)=(1-b)^2+(2-b-a)^2+(2-b-2a)^2+(4-b-3a)^2.</math> 第一组参数的平方和为一，第二组为 <math>0.01+0.04+0.49+0.16=0.7</math>，所以按这一标准第二条线更好。接下来求出所有直线中平方和最小的那条。 [[File:Gezhi-teaching-fit-residuals.svg|frame|center|alt=左边四个残差有正有负且和为零，右边平方后全部非负，零点七残差的平方零点四九占主要部分|先比较左图的正负抵消，再看右图的平方贡献。残差和为零，不等于每个残差为零。]] 本例第三个时间点的残差绝对值最大，因此对平方和的贡献也最大。取平方既防止正负抵消，也使大误差受到更重惩罚。 == 从残差推导两个系数方程 == 先固定斜率，只改变截距 <math>b</math>。每个残差对 <math>b</math> 的导数都是负一，因此 <math display="block">\frac{\partial S}{\partial b}=-2\sum_{i=1}^{4}r_i.</math> 再改变斜率 <math>a</math>，第 <math>i</math> 个残差的变化率为 <math>-t_i</math>，所以 <math display="block">\frac{\partial S}{\partial a}=-2\sum_{i=1}^{4}t_i r_i.</math> 在最小点两者为零，得到 <math>\sum r_i=0</math> 和 <math>\sum t_i r_i=0</math>。代入数据： <math display="block">4b+6a=9,\qquad6b+14a=18.</math> 第一式乘 <math>3/2</math>，得到 <math>6b+9a=13.5</math>；第二式减去它，得 <math>5a=4.5</math>，所以 <math>a=0.9</math>。代回第一式，<math>4b=9-5.4=3.6</math>，所以 <math>b=0.9</math>。 因此候选拟合线为 <math>\hat h=0.9+0.9t</math>。它使残差总和为零，且 <math display="block">\sum t_ir_i=0\cdot0.1+1\cdot0.2+2\cdot(-0.7)+3\cdot0.4=0.</math> 这两个条件分别说明，继续平移或转动这条拟合线，都没有一阶降低平方和的方向。 [[File:Gezhi-model-linear-fit-theme.svg|frame|center|alt=四个实心训练点与拟合直线，竖线为各点到相同时间预测的残差，空心点为未参与拟合的两个留出点|先看四个实心点到直线的竖直距离，再看右侧两个空心点。空心点用于检验，没有参与确定斜率和截距。]] 图中采用竖直距离，因为模型预测的是给定时间下的水位。若时间坐标本身也有显著误差，需另行建立误差模型，不能直接把这个目标解释成点到直线的垂直距离。 == 为什么求出的驻点确实最优 == 导数为零一般不足以证明最小值，本例可以直接比较任何另一条直线。把参数改为 <math>b+u,a+v</math>，新预测比旧预测增加 <math>u+vt_i</math>，新残差为 <math>r_i-u-vt_i</math>。展开平方和： <math display="block">S(b+u,a+v)=\sum r_i^2-2u\sum r_i-2v\sum t_ir_i+\sum(u+vt_i)^2.</math> 在刚才求得的参数处，中间两项都为零，于是 <math display="block">S(0.9+u,0.9+v)=0.7+\sum_{i=1}^{4}(u+vt_i)^2\ge0.7.</math> 所以没有任何参数能更好。要让等号成立，每个 <math>u+vt_i</math> 都须为零；取 <math>t_1=0</math> 得 <math>u=0</math>，再取 <math>t_2=1</math> 得 <math>v=0</math>。因此最优参数也唯一。 含截距拟合还有一个方便性质。从 <math>\sum r_i=0</math> 得 <math display="block">\bar h=b+a\bar t,\qquad b=\bar h-a\bar t.</math> 所以拟合线经过样本平均点。本例 <math>\bar t=1.5,\bar h=2.25</math>。把 <math>b</math> 消去后，另一条系数方程给出 <math display="block">a=\frac{\sum (t_i-\bar t)(h_i-\bar h)}{\sum(t_i-\bar t)^2}=\frac{4.5}{5}=0.9.</math> 分母描述时间点的分散程度；如果所有时间都相同，分母为零，就不能用这些数据分别确定斜率和截距。 == 矩阵形式把同样的推导用于更多参数 == 把截距列与时间列排成设计矩阵： <math display="block">X=\begin{pmatrix}1&0\\1&1\\1&2\\1&3\end{pmatrix},\qquad\beta=\begin{pmatrix}b\\a\end{pmatrix},\qquad y=\begin{pmatrix}1\\2\\2\\4\end{pmatrix}.</math> 每行对应一次观测，每列对应一个模型项，所有预测为 <math>X\beta</math>。目标就是 <math display="block">\min_\beta\|y-X\beta\|^2.</math> 前面的两条残差条件统一写成 <math>X^{\mathsf T}(y-X\hat\beta)=0</math>，即'''正规方程''' <math display="block">X^{\mathsf T}X\hat\beta=X^{\mathsf T}y.</math> 本例直接计算得到 <math display="block">X^{\mathsf T}X=\begin{pmatrix}4&6\\6&14\end{pmatrix},\qquad X^{\mathsf T}y=\begin{pmatrix}9\\18\end{pmatrix},</math> 与上一节逐项求出的两个方程完全相同。 一般线性最小二乘的“线性”指参数的出现方式。例如 <math>\hat y=\beta_0+\beta_1t+\beta_2t^2</math> 虽然画出抛物线，对三个参数仍是线性的，只需使用三列 <math>1,t,t^2</math>。模型 <math>ae^{bt}</math> 对参数 <math>b</math> 非线性，则不能直接使用同一组线性正规方程。 == 几何上是在寻找投影 == 所有可能预测 <math>X\beta</math> 构成一个向量子空间，即 <math>X</math> 的列空间。观测向量 <math>y</math> 未必在其中；最小二乘寻找该空间中离 <math>y</math> 最近的预测 <math>p</math>。 [[File:Gezhi-linear-projection-theme.svg|frame|center|alt=平面点二一投影到横轴上的二零，残差垂直于允许的输出方向|以横轴代表允许的预测空间：移动预测点会在原有竖直误差之外，再增加水平误差。]] 图是高维投影的二维示意。水位数据有四个坐标，不能直接把它们画成普通平面点，但“残差垂直于所有可拟合方向”的关系完全相同。 正规方程说，残差 <math>r=y-p</math> 与每一列正交，因此与列空间中任意方向都正交。若换一个预测 <math>p+Xd</math>，勾股关系给出 <math display="block">\|y-(p+Xd)\|^2=\|r-Xd\|^2=\|r\|^2+\|Xd\|^2.</math> 这就是上一节平方展开证明的几何形式。预测点是'''正交投影'''；原始数据可能有四个或更多坐标，但关系与平面上作垂线相同。 预测唯一，参数却未必唯一。例如让两个参数总以和 <math>\beta_1+\beta_2</math> 出现，去拟合两次观测二和四。令 <math>s=\beta_1+\beta_2</math>，目标变为 <math display="block">(2-s)^2+(4-s)^2=2(s-3)^2+2.</math> 最优预测都是三，但参数 <math>(0,3),(1,2),(1.5,1.5)</math> 都符合要求。设计矩阵的两列相同，无法区分两个参数各自的贡献。只有列线性无关时，参数才唯一。 == QR 分解怎样利用这幅几何图景 == 先考虑设计矩阵列线性无关、观测数不少于参数数的情形。为计算投影，可以把设计矩阵的列改写为互相垂直的单位方向。将 <math>X=QR</math>，其中 <math>Q</math> 的列正交归一，<math>R</math> 是上三角矩阵，这称为薄 QR 分解。 仍用水位数据。第一列 <math>(1,1,1,1)</math> 的长度为二，归一化得 <math>q_1=(1,1,1,1)/2</math>。时间列 <math>t=(0,1,2,3)</math> 沿它的投影系数是 <math>q_1^{\mathsf T}t=3</math>。减去投影后，剩余 <math display="block">t-3q_1=(-1.5,-0.5,0.5,1.5),</math> 长度为 <math>\sqrt5</math>，所以 <math display="block">q_2=\frac{(-3,-1,1,3)}{2\sqrt5},\qquad R=\begin{pmatrix}2&3\\0&\sqrt5\end{pmatrix}.</math> 于是第一列是 <math>2q_1</math>，第二列是 <math>3q_1+\sqrt5q_2</math>。观测在这两个单位方向上的坐标为 <math display="block">q_1^{\mathsf T}y=\frac92,\qquad q_2^{\mathsf T}y=\frac{-3-2+2+12}{2\sqrt5}=\frac9{2\sqrt5}.</math> 让预测具有相同的两个坐标，就得到三角方程 <math display="block">2b+3a=\frac92,\qquad\sqrt5a=\frac9{2\sqrt5}.</math> 先解第二式得 <math>a=0.9</math>，再解第一式得 <math>b=0.9</math>。无法投影到这两个方向上的观测部分，正是无法被模型消除的残差。 QR 也是实用的数值计算方法。直接形成 <math>X^{\mathsf T}X</math> 会使满列秩矩阵的二范数条件数平方，列接近相关时可能加重精度损失；实际数值库通常使用 Householder 反射等方式计算 QR。这里的逐列投影用于展示含义，数值实现细节见 [https://fncbook.com/house/ FNC 的 QR 计算章节]。 == 权重与异常值会怎样影响答案 == 平方会放大大残差的影响。只拟合一个常数时，数据 <math>(0,0,0,10)</math> 的平方损失最小点为均值 <math>2.5</math>；把最后一项改为一百，均值变成二十五。若目标换成绝对偏差和，中位数零在两个例子中都最优。选择平方损失，就是选择了一种具体的误差权衡。 不同观测若需要不同权重，可最小化 <math>\sum_iw_ir_i^2</math>，其中 <math>w_i>0</math>。例如用常数 <math>c</math> 拟合零和十，权重分别为九和一，目标为 <math display="block">9c^2+(10-c)^2=10(c-1)^2+90.</math> 答案是 <math>c=1</math>；等权时则为五。较大的权重使答案更靠近对应观测。在独立误差、已知方差的模型下常取方差倒数为权重，权重如何估计及其限制见 [https://itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm NIST 的加权最小二乘说明]。 == 拟合误差与预测检验 == 水位模型的训练平方和为 <math>0.7</math> 平方厘米，训练均方根误差为 <math display="block">\mathrm{RMSE}_{\mathrm{train}}=\sqrt{0.7/4}\approx0.4183\ \mathrm{cm}.</math> 另外保留未参与拟合的两点：时间四、五分钟，水位 <math>4.6,5.8</math> 厘米。模型预测 <math>4.5,5.4</math> 厘米，残差为 <math>0.1,0.4</math>，留出平方和 <math>0.17</math>，均方根误差约为 <math>0.2915</math> 厘米。两点只是说明检验流程，不足以判断一般预测性能。 计算最小二乘解不需要假设正态分布。若要从参数进一步推断真实变化率或给出置信区间，则需指定误差模型。固定、满列秩设计与零均值误差支持无偏性；方差、相关性及分布假设决定更进一步的统计结论，见[[统计推断]]。训练均方误差的除数是观测数，估计噪声方差时常用观测数减去独立参数数，两者用途不同。 == 历史与复现 == Legendre 在 1805 年发表最小二乘法。Gauss 在 1809 年的著作中发表相关方法，并主张自己更早已使用；发表时间与更早使用的主张需要区分。相关记载见 [https://mathshistory.st-andrews.ac.uk/Biographies/Legendre/ MacTutor 的 Legendre 传记]。 本条数据、标准库 Python 程序与预期结果可从[[数学建模#下载与复现|数学建模的下载说明]]取得。修改某个训练水位后，既可以重新求系数，也可以检查残差之和与时间加权残差之和是否仍为零，再观察留出误差如何变化。 == 参考来源与延伸阅读 == * [https://web.stanford.edu/~boyd/vmls/ Boyd 与 Vandenberghe：Introduction to Applied Linear Algebra]，最小二乘、QR、模型拟合。 * [https://fncbook.com/qr/ Driscoll 与 Braun：The QR factorization]、[https://fncbook.com/house/ Computing QR factorizations]：薄 QR 与稳定实现。 * [https://itl.nist.gov/div898/handbook/pmd/section1/pmd143.htm NIST：Weighted Least Squares Regression]：权重估计与统计条件。 * [https://mathshistory.st-andrews.ac.uk/Biographies/Legendre/ MacTutor：Legendre]，最小二乘发表与优先权背景。 * 先修：[[矩阵]]、[[向量空间]]、[[导数]]；相关：[[优化]]、[[统计推断]]、[[数学建模]]。 [[分类:应用与建模]] [[分类:数值分析]]
返回
最小二乘法
。
格致开物MATHWIKI
登录