本文介绍如何对非网格、不规则分布的二维输入（x, y）与输出（z）数据，使用最小二乘法拟合双线性模型 z = a·x + b·y + c·x·y + d，并准确求解系数 a、b、c、d。方法纯 numpy 实现，无需网格化或插值，适用于缺失值、不等采样密度的真实实验数据。
双线性拟合本质上是带交叉项的线性回归：虽然模型在 x 和 y 上是非线性的（含乘积项 x·y），但关于待估参数 a、b、c、d 是线性的。因此，可直接套用多元线性最小二乘法——即最小化残差平方和
$$
S = \sum_{i=1}^N \left( a x_i + b y_i + c x_i y_i + d - z_i \right)^2
$$
对 a、b、c、d 分别求偏导并令其为 0，可导出标准的正规方程组（Normal Equations）：
$$
\mathbf{A} \cdot \begin{bmatrix}a \ b \ c \ d\end{bmatrix} = \mathbf{b}
$$
其中设计矩阵 $\mathbf{A}$ 和右侧向量 $\mathbf{b}$ 的元素由数据的各阶统计量构成（如 $\sum x_i$, $\sum x_i y_i$, $\sum x_i y_i z_i$ 等）。
以下是一个高效、可复用的 bilinear_fit 函数实现，完全基于 NumPy，无外部依赖：
import numpy as np
def bilinear_fit(data):
"""
对散点数据 (x, y, z) 进行双线性最小二乘拟合：z = a*x + b*y + c*x*y + d
Parameters:
-----------
data : list of [x, y, z] or ndarray of shape (N, 3)
输入数据点集合，支持任意数量、非均匀、非网格分布
Returns:
--------
tuple: (a, b, c, d) 拟合系数
"""
data = np.asarray(data)
x, y, z = data[:, 0], data[:, 1], data[:, 2]
N = len(x)
# 计算必要统计量（避免循环，更高效）
Sx = np.sum(x)
Sxx = np.sum(x * x)
Sy = np.sum(y)
Syy = np.sum(y * y)
Sxy = np.sum(x * y)
Sxxy = np.sum(x * x * y)
Sxyy = np.sum(x * y * y)
Sxxyy = np.sum(x * x * y * y)
Sz = np.sum(z)
Sxz = np.sum(x * z)
Syz = np.sum(y * z)
Sxyz = np.sum(x * y * z)
# 构建正规方程组 A @ [a,b,c,d] = 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])
try:
coeffs = np.linalg.solve(A, RHS)
return coeffs[0], coeffs[1], coeffs[2], coeffs[3]
except np.linalg.LinAlgError as e:
raise ValueError("正规方程组矩阵奇异，请检查数据是否共线（如 x 或 y 全为常数）或样本量不足（N < 4）") from e
# 示例：使用用户提供的真实数据
D = [
[1056, 8, 50.89124679], [1056, 16, 61.62827273], [1056, 32, 78.83079982],
[1056, 48, 92.90073197], [1056, 64, 105.103744 ], [1056, 80, 116.0303753],
[1056, 96, 126.0130906], [1056, 112, 135.2610439], [1056, 128, 143.9159512],
[1056, 144, 152.0790946], [2048, 8, 63.71675604], [2048, 16, 77.15971099],
[2048, 32, 98.69757849], [2048, 48, 116.313387 ], [2048, 64, 131.5917779],
[2048, 80, 145.2721136], [2048, 96, 157.7706532], [2048, 112, 169.3492575],
[2048, 128, 180.1853546], [2048, 144, 190.4057615], [4096, 8, 86.7357654 ],
[4096, 16, 105.0352703], [4096, 32, 134.3541477], [4096, 48, 158.334033 ],
[4096, 64, 179.1320602], [4096, 80, 197.7547066], [4096, 96, 214.7686034],
[4096, 112, 230.5302193], [4096, 128, 245.2810877], [4096, 144, 259.193829 ]
]
a, b, c, d = bilinear_fit(D)
print(f"拟合结果：")
print(f"a (x 系数) = {a:.10f}")
print(f"b (y 系数) = {b:.10f}")
print(f"c (xy 交叉项) = {c:.10f}")
print(f"d (截距) = {d:.10f}")
# 验证拟合效果
z_pred = a * np.array(D)[:,0] + b * np.array(D)[:,1] + c * np.array(D)[:,0] * np.array(D)[:,1] + d
residuals = z_pred - np.array(D)[:,2]
print(f"\nR² 决定系数 = {1 - np.sum(residuals**2) / np.sum((np.array(D)[:,2] - np.mean(D, axis=0)[2])**2):.4f}")✅ 关键优势说明：
- ✅ 完全适配散点数据：不要求 x-y 构成矩形网格，允许任意分布、缺失点、重复 x/y 值；
- ✅ 解析解、高精度：通过正规方程直接求解，比迭代法更稳定、更快；
- ✅ 可解释性强：返回明确物理意义的系数 a, b, c, d，便于后续分析或嵌入公式；
- ✅ 轻量可靠：仅依赖 NumPy，无 sklearn/scipy 插值器的黑盒限制（如无法导出系数、边界外推异常等）。
⚠️ 注意事项：
- 数据量应 ≥ 4（参数个数），否则方程组欠定；
- 若 x 或 y 列全为常数（如所有 y=8），会导致设计矩阵秩亏，需预处理或改用正则化（如 np.linalg.lstsq 加 rcond）；
- 对强非线性或存在显著离群点的数据，建议先做标准化（如 (x - x_mean)/x_std）以提升数值稳定性；
- 如需不确定性评估（如系数标准误），可进一步计算协方差矩阵 $(\mathbf{X}^\top \mathbf{X})^{-1} \sigma^2$，其中 $\mathbf{X}$ 为设计矩阵 $[x_i,\, y_i,\, x_i y_i,\, 1]$。
该方法已在多个工程场景（如传感器响应建模、GPU性能预测、图像几何校正）中验证有效，是处理双线性关系最直接、透明且可控的解决方案。
立即学习“Python免费学习笔记（深入）”；
