Python双线性拟合方法及实现步骤

本文介绍如何对非网格、非均匀、含缺失值的二维散点数据（x, y → z）进行双线性模型 z = ax + by + cxy + d 的最小二乘拟合，使用纯 NumPy 构建解析解，输出可解释的系数，并支持误差评估与预测。

双线性拟合是二维线性回归的自然扩展，适用于建模两个连续变量共同影响一个响应量的场景（如性能调优、传感器校准、图像插值建模等）。与规则网格上的双线性插值不同，本文解决的是任意分布的散点数据（xᵢ, yᵢ, zᵢ）的全局最优拟合问题——即寻找参数 a, b, c, d，使残差平方和
$$ \sum_{i=1}^N \left( a x_i + b y_i + c x_i y_i + d - z_i \right)^2 $$
最小化。该问题本质是线性最小二乘（尽管模型含交叉项 xy，但对参数 a,b,c,d 仍是线性的），因此无需迭代优化，可直接通过正规方程（Normal Equation）求得闭式解析解。

核心实现：构建并求解正规方程组

将模型重写为向量形式：
$$ \mathbf{z} \approx \mathbf{X} \boldsymbol{\beta}, \quad \text{其中 } \boldsymbol{\beta} = [a,\, b,\, c,\, d]^T,\quad \mathbf{X} = \begin{bmatrix} x_1 & y_1 & x_1 y_1 & 1 \ x_2 & y_2 & x_2 y_2 & 1 \ \vdots & \vdots & \vdots & \vdots \ x_N & y_N & x_N y_N & 1 \ \end{bmatrix} $$
则最优参数满足正规方程：
$$ (\mathbf{X}^\top \mathbf{X}) \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{z} $$

为提升数值稳定性与计算效率（尤其当 N 很大时），我们不显式构造大型设计矩阵 X，而是直接累加其 Gram 矩阵 $\mathbf{X}^\top \mathbf{X}$ 和向量 $\mathbf{X}^\top \mathbf{z}$ 的各项元素。对应代码如下：

import numpy as np

def bilinear_fit(data):
    """
    对散点数据 (x, y, z) 进行双线性拟合：z = a*x + b*y + c*x*y + d
    data: list of [x, y, z] triplets
    Returns: tuple (a, b, c, d)
    """
    N = len(data)
    # 累加 Gram 矩阵各元素（4×4 对称阵）
    Sx = Sy = Sxy = Sz = Sxz = Syz = Sxyz = 0.0
    Sxx = Syy = Sxxy = Sxyy = Sxxyy = 0.0

    for x, y, z in data:
        Sx   += x
        Sy   += y
        Sxy  += x * y
        Sz   += z
        Sxx  += x * x
        Syy  += y * y
        Sxz  += x * z
        Syz  += y * z
        Sxyz += x * y * z
        Sxxy += x * x * y
        Sxyy += x * y * y
        Sxxyy += x * x * y * y

    # 构造正规方程系数矩阵 A 和右端向量 RHS
    A = np.array([
        [Sxx,  Sxy,  Sxxy, Sx ],
        [Sxy,  Syy,  Sxyy, Sy ],
        [Sxxy, Sxyy, Sxxyy, Sxy],
        [Sx,   Sy,   Sxy,  N  ]
    ])
    RHS = np.array([Sxz, Syz, Sxyz, Sz])

    # 求解线性系统：A @ [a,b,c,d] = RHS
    coeffs = np.linalg.solve(A, RHS)
    return coeffs[0], coeffs[1], coeffs[2], coeffs[3]

# 示例：使用提供的实测数据
D = [
    [1056,   8,   50.89124679], [1056,  16,  61.62827273], # ...（完整数据同题）
    # （此处省略中间数据，实际使用时请填入全部30组）
    [4096, 144, 259.193829]
]

a, b, c, d = bilinear_fit(D)
print(f"a = {a:.12f}\nb = {b:.12f}\nc = {c:.12f}\nd = {d:.12f}")

使用说明与关键注意事项

• ✅ 适用性广：不要求 x/y 构成矩形网格，允许任意采样密度、缺失组合（如某 y 值下无对应 x），天然处理非结构化数据。
• ⚠️ 数值稳定性：当 x 或 y 的量级很大（如示例中 x 达数千），Sxx, Sxxyy 等高阶项可能引起浮点精度损失。强烈建议在拟合前对 x、y 进行标准化（如减均值除标准差），拟合后再将系数逆变换回原始尺度（详见进阶技巧）。
• ? 拟合质量评估：代码末尾提供了预测值与残差输出。建议进一步计算 R² 分数或 RMSE：

  z_pred = np.array([a*x + b*y + c*x*y + d for x,y,_ in D])
  z_true = np.array([z for _,_,z in D])
  r2 = 1 - np.sum((z_true - z_pred)**2) / np.sum((z_true - np.mean(z_true))**2)
  print(f"R² = {r2:.4f}")

• ? 为什么不用 sklearn？ LinearRegression 完全可用——只需手动构造特征矩阵：

  from sklearn.linear_model import LinearRegression
  X = np.array([[x, y, x*y, 1] for x,y,_ in D])
  y = np.array([z for _,_,z in D])
  model = LinearRegression(fit_intercept=False).fit(X, y)
  a, b, c, d = model.coef_

  此方式更简洁、自动处理缩放/正则化，且兼容 Pipeline。

总结

双线性拟合并非黑盒插值，而是可解释的统计建模工具。本文提供的 NumPy 解析解方法透明、高效、零依赖，特别适合嵌入轻量级部署或教学演示；而 sklearn 方案则更适合工程化场景。无论哪种方式，核心思想一致：将非线性特征（xy）纳入线性框架，通过最小二乘获得全局最优参数。掌握此方法，即可稳健处理大量二维响应建模任务。