通过输入的协方差快速缩放高斯核答案

【问题标题】：Fast scaling of Gaussian Kernel by the Covariance of the Inputs通过输入的协方差快速缩放高斯核
【发布时间】：2015-10-08 04:56:37
【问题描述】：

我目前正在摆弄多元核密度估计，以使用 Matlab 估计水文数据集的概率密度函数 (PDF)。我最熟悉使用高斯内核进行内核密度估计，如 Sharma（2000 和 2014）中所述（其中内核带宽使用高斯参考规则 (GRR) 设置）。 GRR 编写如下（Sharma，2000）：

其中 lambda_ref = 内核的 GRR 带宽，n 是样本大小，d 是我们用于密度估计的数据集的维度。为了估计我们的数据集 X 的多元密度，我们使用以下公式 (Sharma, 2000)：

其中 lamda 与上面的 lamda_ref 相同，S 是 X 的样本协方差，det() 代表行列式。

我的问题是：我知道有很多“快速”的方法可以计算术语 exp() 表示的高斯核函数，例如这里提出的方法（使用 Matlab）：@987654321 @。由于我将使用样本量很大（1000-10,000）的数据集，因此我正在寻找一个快速代码。有谁知道我如何为考虑到样本协方差矩阵 (S^-1) 的逆的第二个方程编写快速代码？

我非常感谢可以在此问题上提供的任何帮助。谢谢！

备注：

我知道有一个用于计算第二个方程的 Matlab 代码，作为子函数在：http://www.mathworks.com/matlabcentral/fileexchange/29039-mutual-information-2-variablle/content/MutualInfo.m 中找到。然而，这段代码在计算核矩阵方面存在瓶颈。

参考资料：

1 A. Sharma，改进供水管理的季节性到年际降雨概率预测：第 3 部分 — 非参数概率预测模型，水文学杂志，第 239 卷，第 1-4 期，2000 年 12 月 20 日，第 249- 258，ISSN 0022-1694，http://dx.doi.org/10.1016/S0022-1694(00)00348-6。

2 Sharma, A. 和 R. Mehrotra（2014 年），仅使用观测信息对自然系统进行建模的信息论替代方案，水资源。 Res., 50, 650–660, doi:10.1002/2013WR013845.

【问题讨论】：

您是否考虑过将fspecial 与'gaussian' 输入一起使用？
感谢@jucestain 的建议，我之前没有使用过'fspecial' 命令，但我找到了一个可以修改以服务于我的目的的代码，并将其发布为对此的回答问题。

标签： matlab multidimensional-array kernel-density

【解决方案1】：

我找到了一个可以根据我的目的进行修改的代码。原始代码列于以下链接：http://www.kernel-methods.net/matlab/kernels/rbf.m。

代码

function K = rbf(coord,sig)

%function K = rbf(coord,sig)
%
% Computes an rbf kernel matrix from the input coordinates
%
%INPUTS
% coord =  a matrix containing all samples as rows
% sig = sigma, the kernel width; squared distances are divided by
%       squared sig in the exponent
%
%OUTPUTS
% K = the rbf kernel matrix ( = exp(-1/(2*sigma^2)*(coord*coord')^2) )
%
%
% For more info, see www.kernel-methods.net

%
%Author: Tijl De Bie, february 2003. Adapted: october 2004 (for speedup).

n=size(coord,1);
K=coord*coord'/sig^2;
d=diag(K);
K=K-ones(n,1)*d'/2;
K=K-d*ones(1,n)/2;
K=exp(K);

包含样本协方差缩放的修改代码：

xcov = cov(x.'); % sample covariance of the data
invxc = pinv(xcov); % inversion of data sample covariance
coord = x.';
sig = sigma; % kernel bandwidth
n = size(coord,1);
K = coord*invxc*coord'/sig^2;
d = diag(K);
K = K-ones(n,1)*d'/2;
K = K-d*ones(1,n)/2;
K = exp(K); % kernel matrix

我希望这可以帮助其他人研究同样的问题。

【讨论】：