C#实战：如何用最小二乘法将两组数据拟合成二次函数

共计 2449 个字符，预计需要花费 7 分钟才能阅读完成。

数据拟合的应用场景

在工程实践中，数据拟合是常见的基础需求。比如在传感器校准中，我们需要将传感器的原始输出值与实际物理量建立准确的数学关系；在机械臂运动轨迹预测中，通过离散的轨迹点拟合出平滑的运动曲线。二次函数拟合作为最基础的曲线拟合方法，可以很好地描述具有单峰特性的数据关系。

最小二乘法原理与库选择

最小二乘法通过最小化误差平方和来寻找最佳函数匹配。对于二次函数拟合问题，主要有两种解法：

1. 解析法（正规方程）：直接求解线性方程组，计算复杂度 O(n^3)
2. 迭代法：如梯度下降，适合大规模数据但需要调参

我们选择 MathNet.Numerics 库，因为：

• 提供完整的线性代数运算功能
• 优化过的矩阵运算性能
• MIT 许可的开源项目
• 良好的.NET 生态集成

核心实现步骤

1. 数据预处理

// 数据标准化处理
public static (double[] x, double[] y) PreprocessData(List<(double x, double y)> rawData)
{
    // 异常值检测（3σ 原则）var stats = new DescriptiveStatistics(rawData.Select(d => d.y));
    double threshold = 3 * stats.StandardDeviation;

    var filtered = rawData
        .Where(d => Math.Abs(d.y - stats.Mean) <= threshold)
        .ToList();

    return (filtered.Select(d => d.x).ToArray(), 
            filtered.Select(d => d.y).ToArray());
}

2. 矩阵构建与求解

二次函数的一般形式：$y = ax^2 + bx + c$

对应的正规方程：

$$
\begin{bmatrix}
Σx^4 & Σx^3 & Σx^2 \
Σx^3 & Σx^2 & Σx \
Σx^2 & Σx & n \
\end{bmatrix}
\begin{bmatrix}
a \
b \
c \
\end{bmatrix}
=
\begin{bmatrix}
Σx^2y \
Σxy \
Σy \
\end{bmatrix}
$$

实现代码：

public static (double a, double b, double c) FitQuadratic(double[] x, double[] y)
{
    // 构建矩阵元素
    double sx4 = x.Sum(xi => Math.Pow(xi, 4));
    double sx3 = x.Sum(xi => Math.Pow(xi, 3));
    double sx2 = x.Sum(xi => xi * xi);
    double sx = x.Sum();
    double n = x.Length;

    double sx2y = x.Zip(y, (xi, yi) => xi * xi * yi).Sum();
    double sxy = x.Zip(y, (xi, yi) => xi * yi).Sum();
    double sy = y.Sum();

    // 构建矩阵
    var A = Matrix<double>.Build.DenseOfArray(new[,] {{sx4, sx3, sx2},
        {sx3, sx2, sx},
        {sx2, sx, n}
    });

    var b = Vector<double>.Build.Dense(new[] {sx2y, sxy, sy});

    // 解方程
    var coeff = A.Solve(b);

    return (coeff[0], coeff[1], coeff[2]);
}

3. 拟合优度计算

R²计算公式：

$$
R^2 = 1 – \frac{SS_{res}}{SS_{tot}}
$$

实现代码：

public static double CalculateRSquared(double[] x, double[] y, 
    Func<double, double> model)
{double yMean = y.Average();
    double ssTot = y.Sum(yi => Math.Pow(yi - yMean, 2));
    double ssRes = x.Zip(y, (xi, yi) => Math.Pow(yi - model(xi), 2)).Sum();

    return 1 - (ssRes / ssTot);
}

完整项目结构

建议的 Visual Studio 解决方案结构：

QuadraticFitter/
├── FitterCore/           # 核心算法库
│   ├── IFitter.cs        # 拟合接口
│   ├── QuadraticFitter.cs
│   └── Models/
├── DataServices/         # 数据服务
│   ├── IDataReader.cs
│   ├── CsvDataReader.cs
│   └── JsonDataReader.cs
├── Visualization/        # 可视化
│   └── PlotBuilder.cs    # 使用 OxyPlot
└── ConsoleApp/           # 示例应用

性能优化

稀疏矩阵处理

当数据点超过 10,000 时，应考虑稀疏矩阵：

var sparseA = Matrix<double>.Build.Sparse(3, 3);
// 只设置非零元素...

多线程计算

Parallel.For(0, x.Length, i => {// 并行计算矩阵元素...});

常见问题与解决方案

病态矩阵识别

计算条件数：

var evd = A.Evd();
double condNumber = evd.D.Maximum / evd.D.Minimum;
if(condNumber > 1e10) 
    Console.WriteLine("警告：病态矩阵！");

浮点精度处理

1. 使用 decimal 类型处理极端数值
2. 数据归一化（将 x 映射到 [0,1] 区间）
3. 增加正则化项

延伸思考

1. 三次函数拟合：只需在矩阵中增加 x^5, x^4, x^3 项
2. 实时数据流处理：
3. 使用递推最小二乘法（RLS）
4. 维护滑动窗口数据
5. 增量更新矩阵元素

总结

通过本文介绍的方法，我们实现了：

1. 使用正规方程法求解二次函数拟合
2. 完整的数据预处理流程
3. 拟合结果的可视化展示
4. 生产环境中的性能优化方案

完整示例代码已上传 GitHub（示例仓库）。在实际项目中，建议结合具体业务需求调整异常值检测阈值和精度处理策略。