Python中快速核回归实现咨询：解决KernelReg运行缓慢问题

我正在进行一维非参数回归的重复模拟，用于评估多种非参数方法的性能，其中Nadaraya-Watson核估计是重点测试对象之一。但目前遇到了KernelReg运行速度过慢的问题，无法满足模拟需求。

起初我猜测速度慢是因为带宽参数的自适应选择，由于模拟中信号的正则性固定，我尝试指定理论上可达最优收敛速率的带宽，但运行速度反而更慢，原因不明。

想请教：是否有办法提升当前代码的计算速度？或者有没有其他能提供更快核回归估计的Python包？

简化示例代码 #

import numpy as np
from numpy import linalg as LA
from statsmodels.nonparametric.kernel_regression import KernelReg

##define true signal
def f(x):
  return(np.cos(x*(2*np.pi)))

L=[250,500,750,1000,1250,1500,1750,2000] ##list of increasing number of observations 
l_NW=[] ##list of errors of all simulations (structured as len(L) by 100 array)
for m in L:
  S_NW=[] ##list of errors for all simulations for a fixed number of observations
  for j in range(100): ## repeat 100 times for each number of observations
    x = np.zeros(m+1)
    l = np.zeros(m+1)
    for i in range(m+1):
      x[i] = f(i/m) + np.random.randn() ##construct noisy signal
      l[i]=i/m
    ##compute Nadayar-Watson estimates
    kde = KernelReg(x, l,var_type='c') 
    f_NW=kde.fit(l)[0]
    err=LA.norm(np.array(f_NW)-np.array([f(x) for x in l])) ##compute error of estimation
    S_NW.append(err)
  l_NW.append(S_NW)

提速方案建议 #

1. 用向量运算替代循环（最直接的优化） #

原代码中大量逐元素循环可以用numpy向量运算替换，这能大幅降低运行时间：

for m in L:
    S_NW = []
    for j in range(100):
        # 向量化生成自变量与响应变量
        l = np.linspace(0, 1, m+1)
        x = f(l) + np.random.randn(m+1)
        # 向量化计算真实信号值
        f_true = f(l)
        # 核回归拟合
        kde = KernelReg(x, l, var_type='c')
        f_NW = kde.fit(l)[0]
        # 计算误差
        err = LA.norm(f_NW - f_true)
        S_NW.append(err)
    l_NW.append(S_NW)

2. 正确指定固定带宽 #

之前指定带宽后变慢大概率是参数传递错误，KernelReg的bw参数需要传入对应自变量维度的列表。对于一维回归，你可以根据理论最优公式预计算带宽：

# 针对高斯核的最优带宽公式（可根据信号调整常数C）
def optimal_bw(n):
    C = 1.0  # 常数可根据信号光滑性微调
    return C * n ** (-1/5)

for m in L:
    sample_size = m + 1
    h = optimal_bw(sample_size)
    S_NW = []
    for j in range(100):
        l = np.linspace(0, 1, sample_size)
        x = f(l) + np.random.randn(sample_size)
        # 传入固定带宽初始化模型
        kde = KernelReg(x, l, var_type='c', bw=[h])
        f_NW = kde.fit(l)[0]
        err = LA.norm(f_NW - f(l))
        S_NW.append(err)
    l_NW.append(S_NW)

3. 替换为更快的核回归实现包 #

如果statsmodels的KernelReg仍不满足速度需求，可尝试以下替代方案：

• scikit-learn KernelRidge：经过高度优化，支持向量化运算

  from sklearn.kernel_ridge import KernelRidge
  
  # 用RBF核模拟Nadaraya-Watson（alpha=0等价于无正则化的核回归）
  kr = KernelRidge(kernel='rbf', alpha=0, gamma=1/(2*h**2))
  kr.fit(l.reshape(-1,1), x)
  f_NW = kr.predict(l.reshape(-1,1))
  

• 手动实现Nadaraya-Watson：一维场景下手动实现简单且速度极快

  def nadaraya_watson(x, l, h, x_eval):
      # 高斯核实现
      def kernel(u):
          return np.exp(-u**2 / 2) / np.sqrt(2*np.pi)
      # 向量化计算权重并归一化
      weights = kernel((x_eval[:, None] - l[None, :])/h)
      weights = weights / weights.sum(axis=1, keepdims=True)
      # 计算估计值
      return weights @ x
  

4. 并行化模拟循环 #

100次独立模拟可以用并行计算加速：

from joblib import Parallel, delayed

def single_simulation(m, h):
    sample_size = m + 1
    l = np.linspace(0, 1, sample_size)
    x = f(l) + np.random.randn(sample_size)
    f_NW = nadaraya_watson(x, l, h, l)
    return LA.norm(f_NW - f(l))

for m in L:
    h = optimal_bw(m+1)
    # 用所有CPU核心并行执行100次模拟
    S_NW = Parallel(n_jobs=-1)(delayed(single_simulation)(m, h) for _ in range(100))
    l_NW.append(S_NW)

内容的提问来源于stack exchange，提问作者BabaUtah