计算科学 - 高效的协方差矩阵计算 MATLAB（来自数据的每行组合） - 吾爱随笔录

高效的协方差矩阵计算 MATLAB（来自数据的每行组合）

计算科学 matlab

2021-12-06 07:49:30

我在统计部门的朋友问我如何有效地进行以下计算。

假设我们有数据 $X\in\mathbb{R}^{N\times 2}$ . 他需要做以下计算：

C_{i, j} = σ^{2} \exp (- \frac{| | X_{i} - X_{j} | |_{2}}{θ})

$C_{i,j}=\sigma^2\exp\left(-\frac{||X_i-X_j||_2}{\theta}\right)$ 在哪里

X_{i}

$X_i$ 是一个有两个元素的元组和

σ, θ

$\sigma,\theta$ 是标量。所以

C

$C$ 是一个

N \times N

$N\times N$ 通过从

X

$X$ .

我认为自己是一个 MATLAB 爱好者，但我不知道如何使用内置的 MATLAB 矩阵运算和函数来有效地计算它。我最快的实现使用标准的for循环：

clear all; close all; clc;

n=10;
X=rand(n,2);
sigma=1;
theta=1;

tic
C=zeros(n,n);
for i=1:n
    for j=1:n
        C(i,j)=sigma*sigma*exp(-norm(X(i,:)-X(j,:),2)/theta);
    end
end
toc

有更快的方法吗？

2个回答

试试下面的片段：

clear all; 
close all; 

n=500;
rand('seed',0)
X=rand(n,2);
sigma=1;
theta=1;

% Original method, A
tic
Ca=zeros(n,n);
% DX1a = zeros(n,n); 
% DX2a = zeros(n,n); 
for i=1:n
    for j=1:n
        Ca(i,j)=sigma*sigma*exp(-norm(X(i,:)-X(j,:),2)/theta);
        % DX1a(i,j) = X(i,1) - X(j,1);
        % DX2a(i,j) = X(i,2) - X(j,2);
    end
end
toc

% Alternative method, B
tic
X1 = X(:,1);
X2 = X(:,2);
[MX1, MX2] = meshgrid(X1,X2);
DX1b = MX1'-MX1;
DX2b = MX2-MX2';
Cb = sigma*sigma*exp(-sqrt(DX1b.^2 + DX2b.^2)/theta);
toc

C_error = norm(Ca-Cb,'fro')

% DX1_error = norm(DX1b-DX1a,'fro')
% DX2_error = norm(DX2b-DX2a,'fro')

这里的技巧是使用 meshgrid 将坐标溢出/展开到一对 NxN 数组中（这里 MX1 表示坐标 1，MX2 表示坐标 2），这样您就可以一次计算坐标差矩阵 DX1 和 DX2 . 然后它只是几个元素操作（squares，sqrts，exp's）来形成C。虽然原始方法没有计算DX1或DX2，但我在那里放了一些注释掉的代码，如果你愿意，你可以重新打开更深入地探索。无论如何，方法 B 在我的盒子上要快得多：

Elapsed time is 7.594 seconds.
Elapsed time is 0.0380321 seconds.
C_error =   1.9128e-014

我鼓励你多测试一些。meshgrid() 函数有点特殊，因此请仔细检查您的数据集/更大的程序以确保。

@rchilton1980 方法的一个稍快的变体，它使用单例扩展而不是meshgrid：

Cc = sigma * sigma * exp(-sqrt((X(:,1) - X(:,1)').^2 + (X(:,2) - X(:,2)').^2) / theta);

还有另一个使用的变体vecnorm（在 R2017b 中引入）。它的速度与方法 B 相当，但我想这将是最有效的解决方案 $X$ 是 $n\times k$ 和 $k>2$ .

Cd = sigma * sigma * exp(-vecnorm(reshape(X, [500,1,2]) - reshape(X, [1,500,2]), 2, 3)/theta);

在我的机器上：

Elapsed time is 0.575712 seconds. % Method A
Elapsed time is 0.013994 seconds. % Method B
Elapsed time is 0.006800 seconds. % Method C
Elapsed time is 0.016893 seconds. % Method D
Cb_error =
   1.3146e-14
Cc_error =
   1.3146e-14
Cd_error =
     0

其它你可能感兴趣的问题

上一篇一种快速检查矩阵是否病态并将其转换为条件良好的方法下一篇谈到数值积分器时，什么是 ABA 和 BAB 方案