最小二乘法是一种常用的数学方法，用于通过一组参考点来近似一个多项式函数。该方法的核心思想是找到一个多项式函数，使得它与参考点的误差平方和最小。这种方法在数据分析、信号处理和机器学习等领域有着广泛的应用。本文将详细介绍最小二乘法的原理、实现以及在多项式拟合中的应用。
假设有一组参考点，每个点由一对坐标(x, y)组成。目标是找到一个多项式函数，其阶数低于参考点的数量，以便对这些点进行拟合。例如，如果有7个参考点，并且想要拟合一个3阶的多项式函数，可以建立7个方程来求解3个独立值。这将导致一个矩阵方程Cx = d，其中C是一个7行3列的矩阵，x是要求解的多项式系数向量，d是参考点的y值向量。
正规方程是最小二乘法中的一种求解方法，由卡尔·弗里德里希·高斯提出，用于计算恒星的椭圆轨道。正规方程通过引入一个误差向量r（残差向量）来求解矩阵方程Cx - d = r。通过矩阵的转置和一系列的计算，可以得到一个n阶的矩阵方程，其中n是所求多项式的阶数。这个新的矩阵方程可以通过高斯消元法求解，得到正确的结果向量x。
构建矩阵方程的第一步是填充矩阵C。C的第一列由1组成，因为正在寻找的第一个参数c0是乘以1的。第二列由x1到xn组成，第三列由(x1)^2到(xn)^2组成，以此类推。
// 定义矩阵结构
struct TMatrix {
double[,] c; // 矩阵C
double[] x; // 系数向量x
double[] d; // 参考点y值向量d
}
// 创建矩阵实例
TMatrix m = new TMatrix();
for (int i = 0; i < samples; i++) {
m.d[i] = y_k[i];
for (int k = 0; k < order; k++) {
m.c[i, k] = Math.Pow(x_k[i], k);
}
}
接下来，需要将矩阵方程乘以C的转置。C的转置矩阵是C沿其主对角线镜像的。例如，一个3x7的矩阵在转置后变成了一个7x3的矩阵。然后，需要计算转置(C) * C和转置(C) * d，并将结果放入临时矩阵a中。
for (int i = 0; i < order; i++) {
a.d[i] = 0.0;
for (int k = 0; k < samples; k++) {
a.d[i] = a.d[i] + m.c[k, i] * m.d[k];
}
for (int j = 0; j < order; j++) {
a.c[j, i] = 0.0;
for (int k = 0; k < samples; k++) {
a.c[j, i] = a.c[j, i] + m.c[k, j] * m.c[k, i];
}
}
}
在文献中提到，使用正规方程可能会导致不准确的结果，这主要是由于A * A^T的乘法可能会导致结果矩阵的条件数变差，从而在后续操作中产生舍入误差。为了避免这种情况，可以使用正交变换。正交变换的基本思想是使用Givens旋转矩阵，通过这些矩阵的乘法不会改变向量的长度，因此矩阵方程Cx - d = r可以被这些正交矩阵转换，而不影响残差平方和。
构建上三角矩阵的过程是迭代的，从元素c21开始，构建旋转矩阵以消除C'中的这个元素，然后进行乘法运算，接着构建下一个旋转矩阵以消除下一个元素c31，依此类推，直到得到一个上三角矩阵。
public void RunGivens(int maxOrder) {
for (int i = 0; i < maxOrder; i++) {
for (int k = i + 1; k < maxOrder; k++) {
if (m.c[k, i] != 0.0) {
CalcRMatrix(i, k);
m.MultMatrix(r, i, k);
}
}
}
}
在上述示例中，最终只对列1到m的元素感兴趣，因为正在寻找一个阶数为m的多项式。如果p大于m，旋转矩阵的乘法只会影响大于m的行元素。因此，可以跳过这些乘法，让p只从1运行到m，这可以节省大量的计算时间。
C#代码示例包含一个主窗口。实现了一个高斯消元法的类来解决矩阵方程，并实现了一个函数。在MainWin类顶部定义了public const samples来固定样本数量。这个样本数量可以在4到40之间设置。在窗口中，有10个TextBox用于编辑最多10个样本的x和y。如果样本数量大于10，这些TextBox将变得不可见。在这种情况下，样本值必须在应用程序启动时初始化。在MainWin函数中实现了7个样本点并初始化它们。