当前位置：首页 > news >正文

最小二乘问题详解3：线性最小二乘实例

news 2025/10/6 7:47:35

1. 引言

在上一篇文章《最小二乘问题详解2：线性最小二乘求解》中笔者详细介绍了如何求解线性最小二乘问题，一般使用QR分解或者SVD分解法，这里笔者就实现一个具体的案例来验证一下。

2. 案例

总是举拟合直线的例子实在太简单了，这里就使用一个更加复杂一点问题模型：双线性变换。具体来说，假设存在两幅地图需要配置，并且找到了各自地图上的同名点，可以使用双线性变换模型来进行快速、初步的校正。也就是说一组同名点满足如下关系：

${x1=a0x0+b0y0+c0x0y0+d0y1=a1x0+b1y0+c1x0y0+d1\begin{cases} x_1 = a_0 x_0 + b_0 y_0 + c_0 x_0 y_0 + d_0 \\ y_1 = a_1 x_0 + b_1 y_0 + c_1 x_0 y_0 + d_1 \end{cases}$

其中， $x_0,y_0)$ 和 $x_1,y_1)$ 是对应的同名点，也就是观测值。而 $a_0,b_0,c_0,d_0$ 和 $a_1,b_1,c_1,d_1$ 则是X和Y方向上的双线性变换系数，也就是待定值。对于观测值来说，这个函数是非线性的，但是对于待定值来说这个问题模型则是线性的。这也是笔者在《最小二乘问题详解1：线性最小二乘》中强调的一点：最小二乘问题是线性还是非线性，需要通过待定值来判断。

要实现这个案例，那么就需要准备一组同名点， $x_0,y_0)$ 可以通过随机数生成， $x_1,y_1)$ 则可以通过双线性变换公式加上一点噪声值来得到。另外，也不需要从头实现矩阵运算的，一定规模的矩阵运算库就可以实现线性方程组的求解，比如这里笔者使用的Eigen。具体的案例代码实现如下：

#include <Eigen/Dense>
#include <random>
#include <vector>using namespace std;
using namespace Eigen;int main() {// ========================// 1. 设置真实参数（我们想要求解的目标）// ========================Vector4d params_x_true, params_y_true;params_x_true << 1.5, -0.3, 0.1, 2.0;    // a0, b0, c0, d0params_y_true << 0.4, 1.8, -0.05, -1.0;  // a1, b1, c1, d1cout << "真实参数 x: a0=" << params_x_true(0) << ", b0=" << params_x_true(1)<< ", c0=" << params_x_true(2) << ", d0=" << params_x_true(3) << endl;cout << "真实参数 y: a1=" << params_y_true(0) << ", b1=" << params_y_true(1)<< ", c1=" << params_y_true(2) << ", d1=" << params_y_true(3) << endl;// ========================// 2. 生成 N 个随机点 (x0, y0) 并计算 (x1, y1)// ========================int N = 100;  // 点的数量vector<double> x0(N), y0(N), x1(N), y1(N);// 随机数生成器random_device rd;  //真随机数生成器，生成不可预测的随机种子mt19937 gen(rd());                                   //伪随机数生成器uniform_real_distribution<double> dis(-10.0, 10.0);  // 均匀分布，模拟观测值normal_distribution<double> noise(0.0, 0.1);  // 正态分布，模拟噪声for (int i = 0; i < N; ++i) {x0[i] = dis(gen);y0[i] = dis(gen);// 应用双线性变换double x1_true = params_x_true(0) * x0[i] + params_x_true(1) * y0[i] +params_x_true(2) * x0[i] * y0[i] + params_x_true(3);double y1_true = params_y_true(0) * x0[i] + params_y_true(1) * y0[i] +params_y_true(2) * x0[i] * y0[i] + params_y_true(3);// 添加噪声x1[i] = x1_true + noise(gen);y1[i] = y1_true + noise(gen);}// ========================// 3. 构造设计矩阵 A 和观测向量 b// ========================// 对于 x 方向：x1 = a0*x0 + b0*y0 + c0*x0*y0 + d0MatrixXd A_x(N, 4);  // N 行，4 列（a0, b0, c0, d0）VectorXd b_x(N);// 对于 y 方向：y1 = a1*x0 + b1*y0 + c1*x0*y0 + d1MatrixXd A_y(N, 4);  // N 行，4 列（a1, b1, c1, d1）VectorXd b_y(N);for (int i = 0; i < N; ++i) {A_x(i, 0) = x0[i];          // a0A_x(i, 1) = y0[i];          // b0A_x(i, 2) = x0[i] * y0[i];  // c0A_x(i, 3) = 1.0;            // d0b_x(i) = x1[i];A_y(i, 0) = x0[i];          // a1A_y(i, 1) = y0[i];          // b1A_y(i, 2) = x0[i] * y0[i];  // c1A_y(i, 3) = 1.0;            // d1b_y(i) = y1[i];}// ========================// 4. 使用 Eigen 求解最小二乘// ========================Vector4d theta_x = A_x.colPivHouseholderQr().solve(b_x);Vector4d theta_y = A_y.colPivHouseholderQr().solve(b_y);// ========================// 5. 输出结果// ========================cout << "\n--- 拟合结果 ---" << endl;cout << "估计参数 x: a0=" << theta_x(0) << ", b0=" << theta_x(1)<< ", c0=" << theta_x(2) << ", d0=" << theta_x(3) << endl;cout << "估计参数 y: a1=" << theta_y(0) << ", b1=" << theta_y(1)<< ", c1=" << theta_y(2) << ", d1=" << theta_y(3) << endl;// ========================// 6. 验证：计算残差// ========================VectorXd residual_x = b_x - A_x * theta_x;VectorXd residual_y = b_y - A_y * theta_y;cout << "\n残差平方和 (x): " << residual_x.squaredNorm() << endl;cout << "残差平方和 (y): " << residual_y.squaredNorm() << endl;cout << "总残差平方和: "<< (residual_x.squaredNorm() + residual_y.squaredNorm()) << endl;return 0;
}

有以下几点需要注意：

随机生成观测值数据的实现

// 随机数生成器
random_device rd;  //真随机数生成器，生成不可预测的随机种子
mt19937 gen(rd());                                   //伪随机数生成器
uniform_real_distribution<double> dis(-10.0, 10.0);  // 均匀分布，模拟观测值
normal_distribution<double> noise(0.0, 0.1);  // 正态分布，模拟噪声

中噪声的随机值不是随意生成的，而要符合正态分布（高斯分布）。这是因为偶然误差（随机误差）通常符合正态分布。

这里分成X方向和Y方向求解了两组线性方程组，因为有两组不相关的待定系数 $a_0,b_0,c_0,d_0)$ 和 $a_1,b_1,c_1,d_1)$ 。

本例使用的QR分解法求解的线性最小二乘问题，如果想使用SVD也很简单，可以将colPivHouseholderQr替换成如下接口：

Vector4d theta_x = A_x.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b_x);
Vector4d theta_y = A_y.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b_y);

最终运行结果如下：

真实参数 x: a0=1.5, b0=-0.3, c0=0.1, d0=2
真实参数 y: a1=0.4, b1=1.8, c1=-0.05, d1=-1--- 拟合结果 ---
估计参数 x: a0=1.5006, b0=-0.299571, c0=0.100173, d0=1.98956
估计参数 y: a1=0.398688, b1=1.80047, c1=-0.0501635, d1=-0.994459残差平方和 (x): 0.745402
残差平方和 (y): 0.945377
总残差平方和: 1.69078

3. 精度

3.1 引出

虽然把最小二乘解求出来了，不过笔者更加关心一个问题，那就是求解的精度是多少？如果我只需要求解一个大致的解，那么随便取四组点求解出来就可以了，反正不能精确求解，得到结果也大差不差——其实这就是最小二乘的意义：我不仅仅求解出来了，还可以明确计算出求解的精度误差，使得观测值与求解的符合度始终在这个误差范围之内，从而实现科学上从定性到定量的质变。

以这个例子来说，由于是先设置的真实参数再进行仿真求解，那么可以直接计算估计值与真实值的差值来确定误差。但是在真实场景中，肯定是不知道真实值的，那么就只能估计精度，具体来说就是计算待定参数的协方差矩阵。

不过要说清楚协方差矩阵不是一件容易的事情，因为这个概念的背景知识是《概率论与数理统计》和《误差理论》和这两门课程的内容。大概论述一下就是，协方差矩阵描述了待定参数的估值的不确定性有多大，它的对角线元素是每个参数的方差，开根号就是标准差，而标准差越小，说明估计越精确，所以标准差就是“精度”的量化。

3.2 协方差

从头来讲，可能我们大学学习的《概率论与数理统计》都忘光了，但是高中数学学习的概率与统计知识应该还是记得的。假设我们使用尺子测量一根棍子的长度，一共测了5次，得到：

$\text{ cm}$

那么平均值（估计值）是：

$xˉ=15∑xi=10.0cm\bar{x} = \frac{1}{5} \sum x_i = 10.0 \text{ cm}$

这个值就是棍子“真实长度”的最佳估计，那么方差（不确定性）是：

$s2=1n−1∑(xi−xˉ)2=0.025s^2 = \frac{1}{n-1} \sum (x_i - \bar{x})^2 = 0.025$

标准差（精度）则是：

$\sqrt{0.025} \approx 0.158 \text{ cm}$

那么可以这样评估精度：“长度是 $10.0 \pm 0.16$ cm”。

那么现在进行升级，不是估计一个数，而是在估多个参数，例如本例中的双线性变换中，需要估计：

$θ=[a0,b0,c0,d0]\theta = [a_0, b_0, c_0, d_0]$

这就好比需要同时测量四根棍子的长度，并且由于数据的相关性和估计过程的耦合，这些参数并不独立。为了更好的描述多个参数的估计，就引入了协方差矩阵（Covariance Matrix）：当有多个参数 $θ=[θ1,θ2,…,θp]\theta = [\theta_1, \theta_2, \dots, \theta_p]$ ，它们的不确定性用一个矩阵表示：

$Cov(θ)=[Var(θ1)Cov(θ1,θ2)⋯Cov(θ2,θ1)Var(θ2)⋯⋮⋮⋱]\mathrm{Cov}(\theta) = \begin{bmatrix} \mathrm{Var}(\theta_1) & \mathrm{Cov}(\theta_1,\theta_2) & \cdots \\ \mathrm{Cov}(\theta_2,\theta_1) & \mathrm{Var}(\theta_2) & \cdots \\ \vdots & \vdots & \ddots \end{bmatrix}$

对角线：每个参数的方差 → 表示它自己的不确定性
非对角线：两个参数之间的协方差 → 表示它们的估计是否相互影响

方差我们还好理解，那么协方差是什么的呢？这个概念也有点抽象，假设有两个随机变量 $X$ 和 $Y$ （比如身高和体重）：

如果 $X$ 大时 $Y$ 也大 → 正协方差
如果 $X$ 大时 $Y$ 小 → 负协方差
如果无关 → 协方差接近 0

也就是说，协方差衡量的是两个变量一起变化的趋势。如果还要说的更加具体一点，那么可以这样定义：对于两个参数 $a^0\hat{a}_0$ 和 $b^0\hat{b}_0$ ，它们的协方差是：

$Cov(a^0,b^0)=1N−1∑i=1N(a^0(i)−aˉ0)(b^0(i)−bˉ0)\mathrm{Cov}(\hat{a}_0, \hat{b}_0) = \frac{1}{N-1} \sum_{i=1}^N (\hat{a}_0^{(i)} - \bar{a}_0) (\hat{b}_0^{(i)} - \bar{b}_0)$

其中：

$N$ ：重复实验次数
$a^0(i)\hat{a}_0^{(i)}$ ：第 $i$ 次实验估计的 $a_0$
$aˉ0\bar{a}_0$ ：所有 $a^0(i)\hat{a}_0^{(i)}$ 的平均值
$bˉ0\bar{b}_0$ ：所有 $b^0(i)\hat{b}_0^{(i)}$ 的平均值

也即是两个参数的协方差就是两者误差乘积的平均，如果误差经常同号，乘积为正，那么协方差就大于0；如果误差经常异号，乘积为负，那么协方差就小于0。

3.3 计算

那么，如何得到这个协方差矩阵呢？这里进行推导说明有点复杂，先直接给出结论，以后有机会再进行详细说明：

$Cov(θ^)=σ2(ATA)−1\mathrm{Cov}(\hat{\theta}) = \sigma^2 (A^T A)^{-1}$

其中 $σ2\sigma^2$ 是噪声方差的估计：

$σ^2=∥r∥2n−p\hat{\sigma}^2 = \frac{\|r\|^2}{n - p}$

$r$ : 残差
$n$ : 数据点数
$p$ : 参数个数

回到问题，要评价最小二乘解的精度，就可以取协方差矩阵对角线的值，开平方得到每个待定参数的标准差。

具体的代码实现如下：

// 假设求解了 theta_x
VectorXd residual = b_x - A_x * theta_x;
int n = A_x.rows();  // 数据点数
int p = A_x.cols();  // 参数数 = 4
double sigma_squared = residual.squaredNorm() / (n - p);// 计算 (A^T A)^{-1}
MatrixXd AtA_inv = (A_x.transpose() * A_x).inverse();// 协方差矩阵
MatrixXd cov_theta = sigma_squared * AtA_inv;// 每个参数的标准差（即“精度”的量化）
VectorXd std_error = cov_theta.diagonal().cwiseSqrt();cout << "参数标准差 (std error):" << endl;
cout << "a0: ±" << std_error(0) << endl;
cout << "b0: ±" << std_error(1) << endl;
cout << "c0: ±" << std_error(2) << endl;
cout << "d0: ±" << std_error(3) << endl;