Есть ли более быстрый способ вычисления поэлементных экспонент матрицы в Matlab?

Question

Я профилирую свой код, и эта строка кажется узким местом.Идея состоит в том, чтобы создать пространственную корреляционную матрицу из n наблюдений заданного измерения, основанную на тета-векторе того же измерения, который содержит размерные гиперпараметры для подгонки.

В моем общем коде эта функция пространственной корреляции вызывается тысячи раз (за одну итерацию) при вычислении функции правдоподобия для оптимизации моих тета-параметров.Число наблюдений увеличивается на единицу в каждой итерации всего кода, и тета-параметры должны быть обновлены в каждой итерации.Таким образом, вы можете видеть, как эта строка кода становится критически важной по мере продвижения алгоритма.

Я думаю, что самой медленной строке, т. Е. R = exp(sum(bsxfun(@times, -abs(D_X).^corr_model,reshape(theta,[1 1 d])),3));, требуется больше времени для вычисления конечной поэлементной экспоненты по полученному nпо матрице n (матрица n по n получается после суммирования по 3 измерениям матрицы расстояний).Но в этой линии многое происходит, поэтому я не уверен, является ли это наиболее важным аспектом общей производительности.

Спасибо за ваше понимание!

Я уже заменил команду repmat на bsxfun, чтобы умножить данные тэты на матрицу размерного расстояния 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 );

• Вместо repmat вы можете просто установить D_X = rand(n), чтобы вы могли вычленить (-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;