Vektorisierung in Matlab zur Beschleunigung teurer Schleife

Question

Wie kann ich den folgenden MATLAB-Code mit Vektorisierung beschleunigen? Momentan benötigt die einzelne Zeile in der Schleife Stunden für den Fall upper = 1e7. Hier

ist der kommentierten Code mit Beispielausgabe:

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 

Ausgang:

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 

Zum Beispiel die letzte Zeile enthält die Primfaktoren von 10, wenn wir diejenigen ignorieren. Ich habe die Matrix 8 Spalten breit gemacht, um die vielen Primfaktoren von Zahlen bis zu 10 Millionen zu berücksichtigen.

Danke für jede Hilfe!

Answer 1

Dies ist nicht Vektorisierung, aber diese Version der Schleife wird etwa die Hälfte der Zeit sparen:

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

Hier ist ein kurzer Maßstab für diesen:

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

Hier ist ein anderer Ansatz:

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))]; 

Dies scheint schneller zu sein als die anderen beiden Optionen (verwendet aber mehr Speicher). Borrowing 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