python - 如何在 Numpy 中实现二次采样 RBF(径向基函数)？

Question

我试图按照 CalTech lecture here 的描述在 Python 和 Numpy 中实现径向基函数.数学对我来说似乎很清楚，所以我觉得它不起作用(或者似乎不起作用)很奇怪。这个想法很简单，为每个高斯分布选择一定数量的二次采样中心形成核矩阵，并尝试找到最佳系数。即求解 Kc = y，其中 K 是具有最小二乘法的高斯核 (gramm) 矩阵。为此我做了:

beta = 0.5*np.power(1.0/stddev,2)
Kern = np.exp(-beta*euclidean_distances(X=X,Y=subsampled_data_points,squared=True))
#(C,_,_,_) = np.linalg.lstsq(K,Y_train)
C = np.dot( np.linalg.pinv(Kern), Y )

但是当我尝试用原始数据绘制插值时，它们看起来完全不同:

有 100 个随机中心(来自数据集)。我还尝试了 10 个中心，它们生成的图形与使用训练集中的每个数据点产生的图形基本相同。我假设使用数据集中的每个数据点应该或多或少完美地复制曲线，但事实并非如此(过度拟合)。它产生:

这似乎不正确。我将提供完整的代码(运行时没有错误):

import numpy as np
from sklearn.metrics.pairwise import euclidean_distances

from scipy.interpolate import Rbf

import matplotlib.pyplot as plt

## Data sets
def get_labels_improved(X,f):
    N_train = X.shape[0]
    Y = np.zeros( (N_train,1) )
    for i in range(N_train):
        Y[i] = f(X[i])
    return Y

def get_kernel_matrix(x,W,S):
    beta = get_beta_np(S)
    #beta = 0.5*tf.pow(tf.div( tf.constant(1.0,dtype=tf.float64),S), 2)
    Z = -beta*euclidean_distances(X=x,Y=W,squared=True)
    K = np.exp(Z)
    return K

N = 5000
low_x =-2*np.pi
high_x=2*np.pi
X = low_x + (high_x - low_x) * np.random.rand(N,1)
# f(x) = 2*(2(cos(x)^2 - 1)^2 -1
f = lambda x: 2*np.power( 2*np.power( np.cos(x) ,2) - 1, 2) - 1
Y = get_labels_improved(X , f)

K = 2 # number of centers for RBF
indices=np.random.choice(a=N,size=K) # choose numbers from 0 to D^(1)
subsampled_data_points=X[indices,:] # M_sub x D
stddev = 100

beta = 0.5*np.power(1.0/stddev,2)
Kern = np.exp(-beta*euclidean_distances(X=X,Y=subsampled_data_points,squared=True))
#(C,_,_,_) = np.linalg.lstsq(K,Y_train)
C = np.dot( np.linalg.pinv(Kern), Y )

Y_pred = np.dot( Kern , C )

plt.plot(X, Y, 'o', label='Original data', markersize=1)
plt.plot(X, Y_pred, 'r', label='Fitted line', markersize=1)
plt.legend()
plt.show()

由于绘图看起来很奇怪，我决定阅读绘图函数的文档，但我找不到任何明显的错误。

Answer 1

插值函数的缩放

主要问题是用于插值的函数标准偏差的不幸选择:

stddev = 100

你的函数的特征(它的驼峰)的大小约为 1。所以，使用

stddev = 1

X值的顺序

之所以会出现红线，是因为 matplotlib 中的 plt 按照给定的顺序连接了连续的数据点。由于您的 X 值是随机排列的，因此会导致困惑的左右移动。使用排序的 X:

X = np.sort(low_x + (high_x - low_x) * np.random.rand(N,1), axis=0)

效率问题

您的 get_labels_improved 方法效率低下，循环遍历 X 的元素。使用 Y = f(X)，将循环留给低级 NumPy 内部。

此外，超定系统的最小二乘解的计算应使用 lstsq 完成。而不是计算伪逆(计算成本高)并乘以它。

这是清理后的代码；使用 30 个中心可以很好地拟合。

import numpy as np
from sklearn.metrics.pairwise import euclidean_distances

import matplotlib.pyplot as plt

N = 5000

low_x =-2*np.pi
high_x=2*np.pi
X = np.sort(low_x + (high_x - low_x) * np.random.rand(N,1), axis=0)
f = lambda x: 2*np.power( 2*np.power( np.cos(x) ,2) - 1, 2) - 1
Y = f(X) 

K = 30 # number of centers for RBF
indices=np.random.choice(a=N,size=K) # choose numbers from 0 to D^(1)
subsampled_data_points=X[indices,:] # M_sub x D
stddev = 1

beta = 0.5*np.power(1.0/stddev,2)
Kern = np.exp(-beta*euclidean_distances(X=X, Y=subsampled_data_points,squared=True))
C = np.linalg.lstsq(Kern, Y)[0]

Y_pred = np.dot(Kern, C)

plt.plot(X, Y, 'o', label='Original data', markersize=1)
plt.plot(X, Y_pred, 'r', label='Fitted line', markersize=1)
plt.legend()
plt.show()