高维kriging代理模型代码编写先装环境咱们用numpy打底scipy做优化import numpy as np from scipy.linalg import cholesky, solve_triangular from scipy.optimize import minimize重点来了——核函数计算。二维情况下直接算距离就行高维得玩点矩阵魔法def rbf_kernel(X1, X2, theta): 高斯核函数theta包含方差和长度尺度参数 支持(n,d)维数据计算 variance theta[0] length_scales theta[1:] # 维度加权欧氏距离 dist_sq np.sum(((X1[:, None] - X2) / length_scales)**2, axis2) return variance * np.exp(-0.5 * dist_sq)这里有个魔鬼细节length_scales是各向异性的每个维度都有自己的长度尺度。但直接这么算100维数据的内存直接爆炸。实测发现超过20维时得改用内存映射文件处理矩阵。参数优化是重头戏。传统极大似然估计得这么玩def neg_log_likelihood(theta, X, y): K rbf_kernel(X, X, theta) 1e-6*np.eye(len(X)) # 加jitter防奇异 L cholesky(K, lowerTrue) alpha solve_triangular(L.T, solve_triangular(L, y, lowerTrue)) # 对数行列式计算 log_det 2 * np.sum(np.log(np.diag(L))) return 0.5 * (y.T alpha log_det)但高维情况下协方差矩阵的条件数会变得极其糟糕。建议改用谱分解截断处理或者上低秩近似。这里有个小技巧把长度尺度参数限制在合理范围内避免优化器跑飞。高维kriging代理模型代码编写预测阶段更要命def predict(X_train, y_train, X_test, theta): K rbf_kernel(X_train, X_train, theta) L cholesky(K 1e-6*np.eye(len(X_train)), lowerTrue) alpha solve_triangular(L.T, solve_triangular(L, y_train, lowerTrue)) K_s rbf_kernel(X_train, X_test, theta) mu K_s.T alpha # 方差计算部分在高维时建议注释掉计算量太大 # v solve_triangular(L, K_s, lowerTrue) # var rbf_kernel(X_test, X_test, theta) - v.T v return mu #, var实测50维数据时预测方差的计算时间比均值多出两个数量级。工程应用建议直接放弃方差计算除非真需要不确定性量化。最后给个优化器调用示例# 参数初始化要讲究别用全零 init_theta np.concatenate([[np.var(y)], np.std(X, axis0)]) res minimize(neg_log_likelihood, init_theta, args(X_train, y_train), bounds[(1e-3, 1e3)] [(1e-2, 1e2)]*X_train.shape[1], methodL-BFGS-B)注意这里给每个长度尺度单独设定了边界这对高维问题至关重要。曾经有个项目因为没设边界优化器直接把某个维度的长度尺度缩到1e-20整个模型直接废掉。最后说句大实话超过100维不如转战神经网络Kriging的优势在中等维度10-30维且数据量少千级样本的场景。不过搞懂了这些门道应付大多数工程优化问题足够了。