题解 | #多项式岭回归#

题目链接

多项式岭回归

题目描述

在连续 $N=200$ 天的供水量监测中，有 $M$ 处记录丢失（ $Gap\_1, \dots, Gap\_M$ ）。目标是利用每一处丢失记录前后的连续真实数据段，构建二阶多项式岭回归模型来估算缺失值。

区间确定：对于位于第 $pos$ 天的丢失记录：
- 前方区间 $[left\_start, pos-1]$ ： $left\_start$ 是 $pos$ 之前最近的一个丢失记录的后一天，若无则为第 $1$ 天。
- 后方区间 $[pos+1, right\_end]$ ： $right\_end$ 是 $pos$ 之后最近的一个丢失记录的前一天，若无则为第 $N$ 天。
- 训练样本为这两个区间内所有的真实记录。
模型构建：使用二阶多项式 $y = \beta_2 x^2 + \beta_1 x + \beta_0$ 。系数 $\beta = [\beta_2, \beta_1, \beta_0]^T$ 通过公式求解： $\beta = (X^T X + \lambda I)^{-1} X^T y$ 其中 $\lambda = 0.1$ ， $I$ 为 $3 \cdot 3$ 单位矩阵， $X$ 的每一行为 $[x_i^2, x_i, 1]$ 。

解题思路

数据读取与标记：读取 $M$ 和 $N$ ，记录每一天的供水量。如果是丢失记录，记录其对应的 $Gap$ 编号和位置。
针对每个 Gap 进行预测：
- 确定 $pos$ 对应的 left_start 和 right_end。
- 收集训练集：天序号 $x_i \in [left\_start, pos-1] \cup [pos+1, right\_end]$ ，对应的供水量为 $y_i$ 。
- 构建矩阵 $A = X^T X + \lambda I$ 。由于是 $3 \cdot 3$ 矩阵，可以直接计算各元素：
  - $A_{00} = \sum x_i^4 + \lambda$
  - $A_{01} = A_{10} = \sum x_i^3$
  - $A_{02} = A_{20} = \sum x_i^2$
  - $A_{11} = \sum x_i^2 + \lambda$
  - $A_{12} = A_{21} = \sum x_i$
  - $A_{22} = n + \lambda$ （ $n$ 为训练样本数）
- 构建向量 $B = X^T y$ ：
  - $B_0 = \sum x_i^2 y_i$
  - $B_1 = \sum x_i y_i$
  - $B_2 = \sum y_i$
- 求解线性方程组 $A \beta = B$ 。由于是 $3 \cdot 3$ 规模，可以使用克拉默法则（Cramer's Rule）或求逆矩阵。
- 计算预测值 $\hat{y} = \beta_2 \cdot pos^2 + \beta_1 \cdot pos + \beta_0$ 。
输出结果：按 $Gap\_1$ 到 $Gap\_M$ 的顺序，保留两位小数输出。

代码

c++
java
python

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <iomanip>

using namespace std;

// 3x3 矩阵求逆并与向量相乘
vector<double> solve_beta(double a[3][3], double b[3]) {
    double det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
                 a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
                 a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);

    double inv[3][3];
    inv[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
    inv[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / det;
    inv[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
    inv[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
    inv[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det;
    inv[1][2] = (a[1][0] * a[0][2] - a[0][0] * a[1][2]) / det;
    inv[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
    inv[2][1] = (a[2][0] * a[0][1] - a[0][0] * a[2][1]) / det;
    inv[2][2] = (a[0][0] * a[1][1] - a[1][0] * a[0][1]) / det;

    vector<double> res(3);
    for (int i = 0; i < 3; ++i) {
        res[i] = inv[i][0] * b[0] + inv[i][1] * b[1] + inv[i][2] * b[2];
    }
    return res;
}

int main() {
    int m, n;
    cin >> m >> n;
    vector<double> data(n + 1, -1.0);
    vector<int> gap_pos(m + 1);
    for (int i = 1; i <= n; ++i) {
        string s;
        cin >> s;
        if (s.find("Gap_") != string::npos) {
            int id = stoi(s.substr(4));
            gap_pos[id] = i;
        } else {
            data[i] = stod(s);
        }
    }

    for (int i = 1; i <= m; ++i) {
        int pos = gap_pos[i];
        int left_start = 1;
        for (int p = pos - 1; p >= 1; --p) {
            if (data[p] < 0) {
                left_start = p + 1;
                break;
            }
        }
        int right_end = n;
        for (int p = pos + 1; p <= n; ++p) {
            if (data[p] < 0) {
                right_end = p - 1;
                break;
            }
        }

        double a[3][3] = {0}, b[3] = {0};
        double lambda = 0.1;
        int count = 0;
        for (int p = left_start; p <= right_end; ++p) {
            if (data[p] >= 0) {
                double x = (double)p;
                double y = data[p];
                double x2 = x * x;
                double x3 = x2 * x;
                double x4 = x3 * x;
                a[0][0] += x4; a[0][1] += x3; a[0][2] += x2;
                a[1][1] += x2; a[1][2] += x;
                b[0] += x2 * y; b[1] += x * y; b[2] += y;
                count++;
            }
        }
        a[1][0] = a[0][1]; a[2][0] = a[0][2]; a[2][1] = a[1][2];
        a[2][2] = (double)count;
        a[0][0] += lambda; a[1][1] += lambda; a[2][2] += lambda;

        vector<double> beta = solve_beta(a, b);
        double ans = beta[0] * pos * pos + beta[1] * pos + beta[2];
        cout << "Gap_" << i << ": " << fixed << setprecision(2) << ans << endl;
    }

    return 0;
}

import java.util.*;

public class Main {
    public static void main(String[] args) {
        Scanner sc = new Scanner(System.in).useLocale(Locale.US);
        int m = sc.nextInt();
        int n = sc.nextInt();
        double[] data = new double[n + 1];
        int[] gapPos = new int[m + 1];
        Arrays.fill(data, -1.0);

        for (int i = 1; i <= n; i++) {
            String s = sc.next();
            if (s.startsWith("Gap_")) {
                int id = Integer.parseInt(s.substring(4));
                gapPos[id] = i;
            } else {
                data[i] = Double.parseDouble(s);
            }
        }

        for (int i = 1; i <= m; i++) {
            int pos = gapPos[i];
            int leftStart = 1;
            for (int p = pos - 1; p >= 1; p--) {
                if (data[p] < 0) {
                    leftStart = p + 1;
                    break;
                }
            }
            int rightEnd = n;
            for (int p = pos + 1; p <= n; p++) {
                if (data[p] < 0) {
                    rightEnd = p - 1;
                    break;
                }
            }

            double[][] a = new double[3][3];
            double[] b = new double[3];
            double lambda = 0.1;
            int count = 0;
            for (int p = leftStart; p <= rightEnd; p++) {
                if (data[p] >= 0) {
                    double x = p;
                    double y = data[p];
                    double x2 = x * x;
                    a[0][0] += x2 * x2; a[0][1] += x2 * x; a[0][2] += x2;
                    a[1][1] += x2; a[1][2] += x;
                    b[0] += x2 * y; b[1] += x * y; b[2] += y;
                    count++;
                }
            }
            a[1][0] = a[0][1]; a[2][0] = a[0][2]; a[2][1] = a[1][2];
            a[2][2] = count;
            a[0][0] += lambda; a[1][1] += lambda; a[2][2] += lambda;

            double[] beta = solve(a, b);
            double ans = beta[0] * pos * pos + beta[1] * pos + beta[2];
            System.out.printf("Gap_%d: %.2f\n", i, ans);
        }
    }

    private static double[] solve(double[][] a, double[] b) {
        double det = a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
                     a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
                     a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);

        double[][] inv = new double[3][3];
        inv[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
        inv[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / det;
        inv[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
        inv[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
        inv[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det;
        inv[1][2] = (a[1][0] * a[0][2] - a[0][0] * a[1][2]) / det;
        inv[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
        inv[2][1] = (a[2][0] * a[0][1] - a[0][0] * a[2][1]) / det;
        inv[2][2] = (a[0][0] * a[1][1] - a[1][0] * a[0][1]) / det;

        double[] res = new double[3];
        for (int i = 0; i < 3; i++) {
            res[i] = inv[i][0] * b[0] + inv[i][1] * b[1] + inv[i][2] * b[2];
        }
        return res;
    }
}

def solve():
    line1 = input().split()
    m, n = int(line1[0]), int(line1[1])
    data = [None] * (n + 1)
    gap_pos = {}
    
    for i in range(1, n + 1):
        s = input().strip()
        if s.startswith("Gap_"):
            idx = int(s[4:])
            gap_pos[idx] = i
        else:
            data[i] = float(s)
            
    for i in range(1, m + 1):
        pos = gap_pos[i]
        
        # 确定 left_start
        left_start = 1
        for p in range(pos - 1, 0, -1):
            if data[p] is None:
                left_start = p + 1
                break
        
        # 确定 right_end
        right_end = n
        for p in range(pos + 1, n + 1):
            if data[p] is None:
                right_end = p - 1
                break
        
        # 训练集
        x_list = []
        y_list = []
        for p in range(left_start, pos):
            if data[p] is not None:
                x_list.append(p)
                y_list.append(data[p])
        for p in range(pos + 1, right_end + 1):
            if data[p] is not None:
                x_list.append(p)
                y_list.append(data[p])
        
        # 构建矩阵 A = X^T * X + lambda * I
        lambd = 0.1
        a00, a01, a02 = 0.0, 0.0, 0.0
        a11, a12, a22 = 0.0, 0.0, 0.0
        b0, b1, b2 = 0.0, 0.0, 0.0
        
        for x, y in zip(x_list, y_list):
            x2 = x * x
            x3 = x2 * x
            x4 = x3 * x
            a00 += x4
            a01 += x3
            a02 += x2
            a11 += x2
            a12 += x
            b0 += x2 * y
            b1 += x * y
            b2 += y
            
        num = len(x_list)
        a22 = float(num)
        
        a = [
            [a00 + lambd, a01, a02],
            [a01, a11 + lambd, a12],
            [a02, a12, a22 + lambd]
        ]
        b = [b0, b1, b2]
        
        # 3x3 矩阵求逆
        det = (a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
               a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
               a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]))
        
        inv = [[0.0]*3 for _ in range(3)]
        inv[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det
        inv[0][1] = (a[0][2] * a[2][1] - a[0][1] * a[2][2]) / det
        inv[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det
        inv[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det
        inv[1][1] = (a[0][0] * a[2][2] - a[0][2] * a[2][0]) / det
        inv[1][2] = (a[1][0] * a[0][2] - a[0][0] * a[1][2]) / det
        inv[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det
        inv[2][1] = (a[2][0] * a[0][1] - a[0][0] * a[2][1]) / det
        inv[2][2] = (a[0][0] * a[1][1] - a[1][0] * a[0][1]) / det
        
        beta = [0.0] * 3
        for r in range(3):
            beta[r] = inv[r][0] * b[0] + inv[r][1] * b[1] + inv[r][2] * b[2]
            
        pred = beta[0] * pos * pos + beta[1] * pos + beta[2]
        print(f"Gap_{i}: {pred:.2f}")

solve()

算法及复杂度

算法：二阶多项式岭回归（最小二乘法的正则化版本）。
时间复杂度： $\mathcal{O}(M \cdot N)$ 。其中 $M$ 是缺失值个数， $N$ 是总天数。对于每个缺失值，我们需要线性扫描确定区间，并计算矩阵元素（线性复杂度），矩阵求逆为 $\mathcal{O}(1)$ （固定 $3 \cdot 3$ ）。
空间复杂度： $\mathcal{O}(N)$ 。用于存储 $N$ 天的监测数据。