有没有更快的方法来计算 Matlab 中矩阵的元素指数？

Question

我正在分析我的代码，这一行似乎是一个瓶颈。这个想法是基于相同维度的θ向量创建给定维度的n个观测值的空间相关矩阵 - 其中包含要拟合的维度超参数。

在我的整体代码中，在计算似然函数以优化我的 theta 参数时，这个空间相关函数被调用了数千次（每次迭代）。在整个代码的每次迭代中，观察的数量都会增加一个，并且每次迭代都必须重新拟合 theta 参数。因此，您可以看到随着算法的进展，这行代码如何变得至关重要。

我认为最慢的线，即R = exp(sum(bsxfun(@times, -abs(D_X).^corr_model,reshape(theta,[1 1 d])),3)); 花费最多的时间来计算最终的 n × n 矩阵的最终元素指数（n × n 矩阵是在对 3 维求和后产生的）距离矩阵）。但是这条线有很多事情要做，所以我不确定这是否是整体性能最关键的方面。

感谢您提供的任何见解！

我已经用bsxfun 替换了repmat 命令，以将给定的thetas 乘以维度距离矩阵D_X，这显着加快了代码速度。但是，嵌套第二个bsxfun 来执行距离平方，即替换abs(D_X).^corr_model，会使代码运行速度变慢。

n = 500;
dimension = 20;
D_X = repmat(rand(n),[1 1 dimension]);
theta = rand(dimension, 1);
corr_model = 2;

R = corr(corr_model,theta, D_X);

function R = corr(corr_model,theta,D_X)
% calculate the correlation matrix for the modified nugget effect model
% corr_model - the correlation model used for the spatial correlation
% corr_model = 2: gaussian correlation function
% theta - vector of hyperparameters
% D_X - the distance matrix for the input locations

d = size(theta,1);

switch corr_model
    case 2 %Gaussian correlation function  
        R = exp(sum(bsxfun(@times, -abs(D_X).^corr_model,reshape(theta,[1 1 d])),3));
end
end

Answer 1

这是一个更优化的版本：

n = 500;
dimension = 20;
D_X = rand(n);
theta = rand(dimension, 1);
corr_model = 2;

function R = corr(corr_model, theta, D_X)
    switch corr_model
        case 2  
            R = exp(-sum(theta)) .^ ( D_X .^ corr_model );
    end
end

已采取以下步骤：

R = exp(sum(bsxfun(@times, -abs(D_X).^corr_model,reshape(theta,[1 1 d])),3));
R = exp((-abs(D_X) .^ corr_model) .* sum(theta)) ;
R = exp((-D_X .^ corr_model) .* sum(theta)) ;
R = exp((D_X .^ corr_model) .* (-sum(theta)));
R = exp(-sum(theta)) .^ ( D_X .^ corr_model );

• 您可以只设置D_X = rand(n) 而不是repmat，这样您就可以分解出(-abs(D_X) .^ corr_model) 并将其乘以theta 的总和。

• 由于D_X 是一个实矩阵abs(D_X)^2 等价于D_X^2。

• 将(-D_X.^corr_model) .*sum(theta) 转换为(D_X.^corr_model) .*(-sum(theta)) 可能会加快计算速度，因为您只会对一个标量进行求反，而该标量应该比对数组求反更便宜。

• 知道x^(ab) = (x^a)^b 表达式exp(( D_X .^ corr_model ) .* (-sum(theta))) 可以转换为exp(-sum(theta)) .^ ( D_X .^ corr_model )。在第一个表达式中，我们对D_X 进行了三个操作，但在第二个表达式中，对D_X 应用了两个操作。

编辑：

正如D_X矩阵的注释部分所指出的那样，上述方法不适用，因此您可以使用以下方法：

D_X1 = reshape(permute(D_X, [3 2 1]), dimension,[]);
%calling thousands of times
for k = 1 : 10000
    R = corr(corr_model,theta, D_X1);
end
% matrix reshaped to its original size
Result = reshape(R,n,n);

function R = corr(corr_model, theta, D_X)
    switch corr_model
        case 2  
            R = exp(-theta.' * (D_X.^corr_model));
    end
end

这里D_X 被置换并重新整形为[20 x 250000] 矩阵，因此您可以使用矩阵乘法来代替乘积之和。
此外，如果相关矩阵是对称的，您可以使用其一半的元素。现在矩阵大小为[20 x 125250]。

idx = tril(true(n));
D_X1 = reshape(permute(D_X, [3 2 1]), dimension,[]);
D_X1 = D_X1(:,idx);

for k = 1 : 10000
    R = corr(corr_model,theta, D_X1);
end

% matrix reshaped to its original size
Result = zeros(n);
Result(idx) = R;
d = diag(Result);
Result = Result + Result.'
Result(1:n+1:end)=d;