高価なループを高速化するMatlabのベクトル化

Question

ベクトル化を使用して、次のMATLABコードを高速化するにはどうすればよいですか？現在、ループ内の単一の行は、ケースupper = 1e7のために何時間も実行されています。

p = 8; 
lower = 1; 
upper = 1e1; 
n = setdiff(lower:upper,primes(upper)); % contains composite numbers between lower + upper 
x = ones(length(n),p); % Preallocated 2-D array of ones 

% This loop stores the unique prime factors of each composite 
% number from 1 to n, in each row of x. Since the rows will have 
% varying lengths, the rows are padded with ones at the end. 

for i = 1:length(n) 
    x(i,:) = [unique(factor(n(i))) ones(1,p-length(unique(factor(n(i)))))]; 
end 

出力：ここ

サンプル出力とコメントコードである

x = 

1  1  1  1  1  1  1  1 
2  1  1  1  1  1  1  1 
2  3  1  1  1  1  1  1 
2  1  1  1  1  1  1  1 
3  1  1  1  1  1  1  1 
2  5  1  1  1  1  1  1 

例えば、我々はものを無視した場合、最後の行は、10の素因数を含んでいます。私は行列を8列幅にし、1000万までの数の多くの素因数を説明しました。

ありがとうございました！

Answer 1

これはベクトル化ではなく、ループのこのバージョンは、時間の約半分を保存します：

function t = getPrimeTime 
lower = 1; 
upper = 2.^(1:8); 
t = zeros(numel(upper),2); 
for k = 1:numel(upper) 
    n = setdiff(lower:upper(k),primes(upper(k))); % contains composite numbers between lower to upper 
    t(k,1) = timeit(@() getPrime1(n)); 
    t(k,2) = timeit(@() getPrime2(n)); 
    disp(k) 
end 
p = plot(log2(upper),log10(t)); 
p(1).Marker = 'o'; 
p(2).Marker = '*'; 
xlabel('log_2(range of numbers)') 
ylabel('log(time (sec))') 
legend({'getPrime1','getPrime2'}) 
end 

function x = getPrime1(n) % the originel function 
p = 8; 
x = ones(length(n),p); % Preallocated 2-D array of ones 
for k = 1:length(n) 
    x(k,:) = [unique(factor(n(k))) ones(1,p-length(unique(factor(n(k)))))]; 
end 
end 

function x = getPrime2(n) 
p = 8; 
x = ones(numel(n),p); % Preallocated 2-D array of ones 
for k = 1:numel(n) 
    tmp = unique(factor(n(k))); 
    x(k,1:numel(tmp)) = tmp; 
end 
end 

Answer 2

：ここ

for k = 1:numel(n) 
    tmp = unique(factor(n(k))); 
    x(k,1:numel(tmp)) = tmp; 
end 

は、このための迅速なベンチマークです別の方法があります：

p = 8; 
lower = 1; 
upper = 1e1; 
p = 8; 
q = primes(upper); 
n = setdiff(lower:upper, q); 
x = bsxfun(@times, q, ~bsxfun(@mod, n(:), q)); 
x(~x) = inf; 
x = sort(x,2); 
x(isinf(x)) = 1; 
x = [x ones(size(x,1), p-size(x,2))]; 

これは、他の2つのオプションより速いようですが（しかし、より多くのメモリを使用しています）。借入EBH's benchmarking code：あなたは並列コンピューティングツールボックスへのアクセス権を持っている場合

function t = getPrimeTime 
lower = 1; 
upper = 2.^(1:12); 
t = zeros(numel(upper),3); 
for k = 1:numel(upper) 
    n = setdiff(lower:upper(k),primes(upper(k))); 
    t(k,1) = timeit(@() getPrime1(n)); 
    t(k,2) = timeit(@() getPrime2(n)); 
    t(k,3) = timeit(@() getPrime3(n)); 
    disp(k) 
end 
p = plot(log2(upper),log10(t)); 
p(1).Marker = 'o'; 
p(2).Marker = '*'; 
p(3).Marker = '^'; 
xlabel('log_2(range of numbers)') 
ylabel('log(time (sec))') 
legend({'getPrime1','getPrime2','getPrime3'}) 
grid on 
end 

function x = getPrime1(n) % the originel function 
p = 8; 
x = ones(length(n),p); % Preallocated 2-D array of ones 
for k = 1:length(n) 
    x(k,:) = [unique(factor(n(k))) ones(1,p-length(unique(factor(n(k)))))]; 
end 
end 

function x = getPrime2(n) 
p = 8; 
x = ones(numel(n),p); % Preallocated 2-D array of ones 
for k = 1:numel(n) 
    tmp = unique(factor(n(k))); 
    x(k,1:numel(tmp)) = tmp; 
end 
end 

function x = getPrime3(n) % Approach in this answer 
p = 8; 
q = primes(max(n)); 
x = bsxfun(@times, q, ~bsxfun(@mod, n(:), q)); 
x(~x) = inf; 
x = sort(x,2); 
x(isinf(x)) = 1; 
x = [x ones(size(x,1), p-size(x,2))]; 
end