расчет тангенса касательной в матричной форме MATLAB - PullRequest
Новая
       4

расчет тангенса касательной в матричной форме MATLAB

0 голосов
/ 26 ноября 2018

Рассмотрим следующий расчет корреляции тангенса касательной, который выполняется в цикле for 

v1=rand(25,1);
v2=rand(25,1);
n=25;
nSteps=10;
mean_theta = zeros(nSteps,1);
for j=1:nSteps
    theta=[];
    for i=1:(n-j)
        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);
        n1 = norm([v1(i) v2(i)]);
        n2 = norm([v1(i+j) v2(i+j)]);
        theta = [theta acosd(d/n1/n2)];  

    end
    mean_theta(j)=mean(theta);
end
plot(mean_theta)

Как можно использовать вычисления матрицы Matlab для улучшения этой производительности?

Ответы [ 2 ]

0 голосов
/ 26 ноября 2018

Вот полное векторизованное решение: 

i = 1:n-1;
j = (1:nSteps).';
ij= min(i+j,n);
a = cat(3, v1(i).', v2(i).');
b = cat(3, v1(ij), v2(ij));

d = sum(a .* b, 3);
n1 = sum(a .^ 2, 3);
n2 = sum(b .^ 2, 3);
theta = acosd(d./sqrt(n1.*n2));
idx = (1:nSteps).' <= (n-1:-1:1);

mean_theta = sum(theta .* idx ,2) ./ sum(idx,2);

Результат синхронизации октав для моего метода method4 из ответа, предоставленного @CrisLuengo и оригинального метода (n=250): 

Full vectorized    : 0.000864983 seconds
Method4(Vectorize) : 0.002774 seconds
Original(loop)     : 0.340693 seconds
0 голосов
/ 26 ноября 2018

Есть несколько вещей, которые вы можете сделать, чтобы ускорить ваш код.Во-первых, всегда предварительно распределяйте .Это преобразует: 

theta = [];
for i = 1:(n-j)
    %...
    theta = [theta acosd(d/n1/n2)];
end

в: 

theta = zeros(1,n-j);
for i = 1:(n-j)
    %...
    theta(i) = acosd(d/n1/n2);
end

Затем переместите нормализацию из циклов.Нет необходимости снова и снова нормализовать, просто нормализуйте входные данные: 

v = [v1,v2];
v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB
%...
    theta(i) = acosd(dot(v(i,:),v(i+j,:)));

Это очень незначительно меняет выходные данные с точностью до чисел, потому что разный порядок операций приводит к разным округлениям с плавающей запятойошибка.

Наконец, вы можете удалить внутренний цикл, векторизовав это вычисление: 

i = 1:(n-j);
theta = acosd(dot(v(i,:),v(i+j,:),2));

Время (n=25):

  • Оригинал: 0,0019 с
  • Предварительное выделение: 0,0013 с
  • Нормализовать один раз: 0,0011 с
  • Векторизация: 1,4176e-04 с

Время (n=250):

  • Оригинал: 0,0185 с
  • Предварительно выделить: 0,0146 с
  • Нормализовать один раз: 0,0118 с
  • Векторизация: 2,5694e-04 с

Обратите внимание, что векторизованный код - единственный, чье время не растет линейно с n.

Код времени: 

function so
n = 25;
v1 = rand(n,1);
v2 = rand(n,1);
nSteps = 10;
mean_theta1 = method1(v1,v2,nSteps);
mean_theta2 = method2(v1,v2,nSteps);
fprintf('diff method1 vs method2: %g\n',max(abs(mean_theta1(:)-mean_theta2(:))));
mean_theta3 = method3(v1,v2,nSteps);
fprintf('diff method1 vs method3: %g\n',max(abs(mean_theta1(:)-mean_theta3(:))));
mean_theta4 = method4(v1,v2,nSteps);
fprintf('diff method1 vs method4: %g\n',max(abs(mean_theta1(:)-mean_theta4(:))));

timeit(@()method1(v1,v2,nSteps))
timeit(@()method2(v1,v2,nSteps))
timeit(@()method3(v1,v2,nSteps))
timeit(@()method4(v1,v2,nSteps))

function mean_theta = method1(v1,v2,nSteps)
n = numel(v1);
mean_theta = zeros(nSteps,1);
for j = 1:nSteps
    theta=[];
    for i=1:(n-j)
        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);
        n1 = norm([v1(i) v2(i)]);
        n2 = norm([v1(i+j) v2(i+j)]);
        theta = [theta acosd(d/n1/n2)];
    end
    mean_theta(j) = mean(theta);
end

function mean_theta = method2(v1,v2,nSteps)
n = numel(v1);
mean_theta = zeros(nSteps,1);
for j = 1:nSteps
    theta = zeros(1,n-j);
    for i = 1:(n-j)
        d = dot([v1(i) v2(i)],[v1(i+j) v2(i+j)]);
        n1 = norm([v1(i) v2(i)]);
        n2 = norm([v1(i+j) v2(i+j)]);
        theta(i) = acosd(d/n1/n2);
    end
    mean_theta(j) = mean(theta);
end

function mean_theta = method3(v1,v2,nSteps)
v = [v1,v2];
v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB
n = numel(v1);
mean_theta = zeros(nSteps,1);
for j = 1:nSteps
    theta = zeros(1,n-j);
    for i = 1:(n-j)
        theta(i) = acosd(dot(v(i,:),v(i+j,:)));
    end
    mean_theta(j) = mean(theta);
end

function mean_theta = method4(v1,v2,nSteps)
v = [v1,v2];
v = v./sqrt(sum(v.^2,2)); % Can use VECNORM in newest MATLAB
n = numel(v1);
mean_theta = zeros(nSteps,1);
for j = 1:nSteps
    i = 1:(n-j);
    theta = acosd(dot(v(i,:),v(i+j,:),2));
    mean_theta(j) = mean(theta);
end
...