案例
总是举拟合直线的例子实在太简单了，这里就使用一个更加复杂一点问题模型：双线性变换。具体来说，假设存在两幅地图需要配置，并且找到了各自地图上的同名点，可以使用双线性变换模型来进行快速、初步的校正。也就是说一组同名点满足如下关系：
{
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
其中，
(
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, d0
params_y_true << 0.4, 1.8, -0.05, -1.0; // a1, b1, c1, d1
cout << "真实参数 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 + d0
MatrixXd A_x(N, 4); // N 行，4 列（a0, b0, c0, d0）
VectorXd b_x(N);
// 对于 y 方向：y1 = a1*x0 + b1*y0 + c1*x0*y0 + d1
MatrixXd 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]; // a0
A_x(i, 1) = y0[i]; // b0
A_x(i, 2) = x0[i] * y0[i]; // c0
A_x(i, 3) = 1.0; // d0
b_x(i) = x1[i];
A_y(i, 0) = x0[i]; // a1
A_y(i, 1) = y0[i]; // b1
A_y(i, 2) = x0[i] * y0[i]; // c1
A_y(i, 3) = 1.0; // d1
b_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次，得到：
x
=
[
10.1
,
9.9
,
10.0
,
10.2
,
9.8
]
cm
那么平均值（估计值）是：
¯
x
=
1
5
∑
x
i
=
10.0
cm
这个值就是棍子“真实长度”的最佳估计，那么方差（不确定性）是：
s
2
=
1
n
−
1
∑
(
x
i
−
¯
x
)
2
=
0.025
标准差（精度）则是：
s
=
√
0.025
≈
0.158
cm
那么可以这样评估精度：“长度是
10.0
±
0.16
cm”。
那么现在进行升级，不是估计一个数，而是在估多个参数，例如本例中的双线性变换中，需要估计：
θ
=
[
a
0
,
b
0
,
c
0
,
d
0
]
这就好比需要同时测量四根棍子的长度，并且由于数据的相关性和估计过程的耦合，这些参数并不独立。为了更好的描述多个参数的估计，就引入了协方差矩阵（Covariance Matrix）：当有多个参数
θ
=
[
θ
1
,
θ
2
,
…
,
θ
p
]
，它们的不确定性用一个矩阵表示：
C
o
v
(
θ
)
=
⎡
⎢
⎢
⎣
V
a
r
(
θ
1
)
C
o
v
(
θ
1
,
θ
2
)
⋯
C
o
v
(
θ
2
,
θ
1
)
V
a
r
(
θ
2
)
⋯
⋮
⋮
⋱
⎤
⎥
⎥
⎦
对角线：每个参数的方差 → 表示它自己的不确定性
非对角线：两个参数之间的协方差 → 表示它们的估计是否相互影响
方差我们还好理解，那么协方差是什么的呢？这个概念也有点抽象，假设有两个随机变量
X
和
Y
（比如身高和体重）：
如果
X
大时
Y
也大 → 正协方差
如果
X
大时
Y
小 → 负协方差
如果无关 → 协方差接近 0
也就是说，协方差衡量的是两个变量一起变化的趋势。如果还要说的更加具体一点，那么可以这样定义：对于两个参数
^
a
0
和
^
b
0
，它们的协方差是：
C
o
v
(
^
a
0
,
^
b
0
)
=
1
N
−
1
N
∑
i
=
1
(
^
a
(
i
)
0
−
¯
a
0
)
(
^
b
(
i
)
0
−
¯
b
0
)
其中：
N
：重复实验次数
^
a
(
i
)
0
：第
i
次实验估计的
a
0
¯
a
0
：所有
^
a
(
i
)
0
的平均值
¯
b
0
：所有
^
b
(
i
)
0
的平均值
也即是两个参数的协方差就是两者误差乘积的平均，如果误差经常同号，乘积为正，那么协方差就大于0；如果误差经常异号，乘积为负，那么协方差就小于0。
3.3 计算
那么，如何得到这个协方差矩阵呢？这里进行推导说明有点复杂，先直接给出结论，以后有机会再进行详细说明：
C
o
v
(
^
θ
)
=
σ
2
(
A
T
A
)
−
1
其中
σ
2
是噪声方差的估计：
^
σ
2
=
∥
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;
