求拟合类正弦数据下界的Linear/Polynomial回归算法（所有点不低于曲线）

Hey there! 不用客气，研究初期遇到这类问题太正常了～你需要的是带不等式约束的回归算法——常规回归是找贴合数据中心的曲线，而我们要强制所有数据点都在拟合曲线的上方，也就是对每个数据点$(x_i, y_i)$，必须满足$\hat{y}(x_i) \leq y_i$（$\hat{y}(x)$是拟合出的下界曲线）。下面给你拆解具体思路和代码示例：

核心思路：约束最小二乘回归 #

常规多项式回归的目标是最小化$\sum_{i=1}^n (y_i - \hat{y}(x_i))^2$（整体拟合误差），现在我们要给这个目标加上硬约束：所有$\hat{y}(x_i)$都不能超过$y_i$。这属于带不等式约束的优化问题，有两种常用实现方式：

1. 线性规划（适合线性回归场景） #

如果选择线性拟合$\hat{y}(x) = ax + b$，约束条件就是$ax_i + b \leq y_i$对所有数据点成立。我们可以把目标设为最小化曲线的整体"偏差"（比如最小化$\sum_{i=1}^n (y_i - (ax + b))$），转化为线性规划问题求解，很多数值库都支持这类求解。

2. 带约束的最小二乘（推荐用于多项式回归） #

如果要拟合更高阶的多项式（比如二次、三次，对应你说的"更佳方案"），目标还是最小化平方误差，但要满足$\hat{y}(x_i) \leq y_i$的约束。这类问题可以用拉格朗日乘数法推导，或者直接借助成熟的数值优化库（比如Eigen、GSL）来快速实现，不用自己造轮子。

C++ 示例代码（二次多项式下界拟合） #

下面用Eigen库（轻量级线性代数库，很适合C++做数值计算）实现二次多项式的下界拟合，代码里包含了数据模拟、约束构建和结果验证：

#include <Eigen/Dense>
#include <Eigen/ConstrainedLeastSquares>
#include <vector>
#include <iostream>
#include <cmath>

using namespace Eigen;

// 构建多项式特征矩阵：给每个x生成[1, x, x², ..., x^degree]
MatrixXd buildPolyFeatures(const std::vector<double>& x, int degree) {
    MatrixXd features(x.size(), degree + 1);
    for (int i = 0; i < x.size(); ++i) {
        double xi = x[i];
        features(i, 0) = 1.0;
        for (int d = 1; d <= degree; ++d) {
            features(i, d) = pow(xi, d);
        }
    }
    return features;
}

int main() {
    // 模拟你的类正弦离散数据（蓝色点）
    std::vector<double> x = {0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0};
    std::vector<double> y = {0.1, 0.6, 0.9, 0.7, 0.3, 0.2, 0.5, 0.8, 0.9};

    // 选择二次多项式拟合（degree=2，你可以改成3、4等更高阶）
    int degree = 2;
    MatrixXd X = buildPolyFeatures(x, degree);
    VectorXd Y(y.data(), y.size());

    // 构建约束：要求拟合值X*w ≤ Y（所有数据点都在曲线上方）
    ConstrainedLeastSquaresProblem<MatrixXd, VectorXd> clsp;
    clsp.setUnconstrainedProblem(X, Y); // 基础最小二乘目标
    clsp.addInequalityConstraint(X, Y); // 添加不等式约束

    // 求解带约束的最小二乘
    VectorXd w = clsp.solve();

    // 输出拟合得到的多项式系数
    std::cout << "拟合的二次多项式系数：" << std::endl;
    std::cout << "常数项 w0 = " << w(0) << "\n一次项 w1 = " << w(1) << "\n二次项 w2 = " << w(2) << std::endl;

    // 验证每个数据点是否满足约束
    std::cout << "\n验证所有数据点是否在拟合曲线上方：" << std::endl;
    for (int i = 0; i < x.size(); ++i) {
        double y_hat = w(0) + w(1)*x[i] + w(2)*pow(x[i], 2);
        // 加1e-6是避免浮点精度误差导致的误判
        bool constraint_satisfied = (y_hat <= y[i] + 1e-6);
        std::cout << "x=" << x[i] << " | 原始y=" << y[i] << " | 拟合下界y_hat=" << y_hat << " | 约束满足：" << (constraint_satisfied ? "是" : "否") << std::endl;
    }

    return 0;
}

额外说明 #

• 如果你不想用Eigen库，也可以手动实现带约束的梯度下降：每次更新多项式参数后，检查所有数据点的约束是否被违反，若违反则调整参数回到约束范围内，反复迭代直到收敛；
• 你提到的"理想解决方案"（绿色曲线），其实就是约束最小二乘的最优结果——在满足所有点在上方的前提下，让曲线尽可能贴近数据的下界，平衡了约束和拟合精度；
• 如果是线性回归场景，约束逻辑更简单，甚至可以用线性规划的单纯形法手动实现，但用库会更高效稳定。

内容的提问来源于stack exchange，提问作者Kowalski Paweł