本文介绍如何对非网格、非均匀、含缺失值的二维散点数据（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）纳入线性框架，通过最小二乘获得全局最优参数。掌握此方法，即可稳健处理大量二维响应建模任务。