Векторизация короткой петли так, чтобы Octave соответствовал скорости Matlab - PullRequest
Новая
       6

Векторизация короткой петли так, чтобы Octave соответствовал скорости Matlab

1 голос
/ 14 марта 2012

Мой код требует ~ 0,003 с для запуска в Октаве и ~ 0,0007 в Matlab. Поскольку в Octave нет JIT-компиляции, я думаю, что Matlab выполняет закулисную оптимизацию, которую я должен делать сам. 

numerators = zeros(1, 64);
for c = 1 : C
  numerators(c) = py(c) * prod(diag(px{c}(:, x)));
end

py - это 1xC вектор строки. px - это массив с C элементами, каждый из которых представляет собой DxV матрицу. x - это вектор столбца Dx1, значения которого дискретны в [1-V].

То, что prod(diag(...)) нечетность - это векторизованный способ умножения px{c}(d, x(d)) для всех d: 

p = 1;
for d = 1 : D
  p *= px{c}(d, x(d))
endfor

cellfun может работать, но я увяз в деталях передачи аргументов. (Если это можно сделать, просто скажи так, и я пойму, как себя). Другим вариантом может быть использование 3-D матрицы для px, однако, я попробовал это, и мне достаточно новичка, что я ничего не смог заставить работать.

1 Ответ

1 голос
/ 14 марта 2012

Я сделаю это на основе предоставленной вами информации. У меня нет информации о времени.

  1. Создание px трехмерной матрицы DxVxC (вместо 1xC массива ячеек из DxV 2D матриц)
  2. px_new=px(:,x,:) реорганизует матрицу, помещая интересующие вас значения в основную диагональ первых двух измерений
  3. Создать логический индекс с индексами по диагонали (mask=eye(D,D);); repmat, так что это размер C в третьем измерении.
  4. Индексируйте в px_new и reshape, чтобы в столбцах было C.
  5. prod столбцы, оставляя 1xC вектор

  6. Умножьте это (поэлементно) на py, чтобы получить результат 

    px = nan(3,4,5);    %# create test 3D matrix  
px(:, :, 1) = [1 2 3 4; 4 5 6 4; 7 8 9 4]; 
px(:,:,2)=px(:,:,1)*1.5;  
px(:,:,3)=px(:,:,2)*1.5;  
px(:,:,4)=px(:,:,3)*1.5;      
px(:,:,5)=px(:,:,4)*1.5;  
x = [4 2 3];            %# 1xD vector discrete on 1-V  
px_new=px(:,x,:);       %# reorganize into DxDxC  

idx=logical(repmat(eye(size(pd_new,1))),[1,1,size(pd_new,3)])); %# logical index  
P = prod(reshape(pd_new(idx),[],size(pd_new,3))); %#P is now 1xC vector

Этот код был протестирован на http://www.online -utility.org / math / math_calculator.jsp

Редактировать: Первоначально я делал несколько ненужных шагов. Я обновил его, чтобы сделать его более кратким.

...