A vetorização no MatLab

Novembro 2016

Sumário





O software MatLab


MatLab, contração de Matrixe do Laboratory, é um software otimizado para o cálculo de matrizes. Isso resulta em um programa que pode ser feito tanto com um ou mais loops, ou usando um cálculo de matriz, geralmente se executará mais rápido usando uma matriz de cálculo, já que neste caso, MatLab acessa os dados em uma ordem adequada para o armazenamento na memória.
A gestão dos Loops em MatLab melhorou com as últimas versões, mas muitas vezes é vantajoso vetorializar seus cálculos para obter em um código a maior parte de cálculos de matrizes, em vez de loops.

Configuração material para os testes


Processador: 1,6 GHz
RAM: 2 Go
Swap: 4 Go

Como vetorizar um cálculo


Antes de tudo, a maior parte das funções MatLab podem agir em vetores e não somente em simples números ; por exemplo, para cálculos os valores do logaritmo neperianos para os inteiros de 1 à 10^6, estocados em um vetor linha V, basta fazer:
v=log((1:1e6));

Você pode então verificar, com a ajuda de algumas linhas, que esta operação é mais rápida que um loop.
tic;  
v=zeros(1,1e6);  
for p=1:1e6  
    v(p)=log(p);  
end  
toc;
O tempo de execução para este loop é de 3,7s, enquanto que o tempo de execução do mesmo cálculo escrito desta maneira
tic;v=log((1:1e6));toc;
est de 0,3s.

Seja um fator 19 entre os dois tempos de execução.

Em seguida, muitas das funções MatLab permitem de se livrar de cálculos com a ajuda dos loops. Eis aqui uma lista não exaustiva :
  • As operações aritméticas *, /, ^ podem agir em quadros de mesmo tamanho efetuando a operação termo à termo ; basta adicionar um ponto diante da operação para precisar aquilo que se quer realizar nesta operação, quer dizer *, ./, .^. Por exemplo, (1:10).^2 dará os 10 primeiros quadrados perfeitos.
  • repmat - muito simples de usar, digit help repmat ou doc repmat para ter uma ajuda para esta função - permite repetir um quadro (este quadro pode ser um vetor ou uma matriz) para construir uma maior.
  • reshape pode redimensionar uma imagem, usar esta função requer algum conhecimento de MatLab . Como ele "substituir" os elementos ao redimensionar. Novamente, a ajuda dada pelo MatLab ajuda a remodelar doc, você pode ver a dica no [manipulação / faq/sujet-11091-manipulations-elementaires-des-tableaux-sous-matlab básicos das tabelas] para uma explicação feita de forma diferente.
  • Muitas outras coisas !

Experiência sobre um exemplo concreto


No objetivo de ilustrar a vantagem da vetorização, trabalhemos sobre um exemplo concreto. Saber o que faz o programa que nus vamos utilizar não é necessariamente uma obrigação absoluta, ele foi escolhido de maneira a mostrar que vetorizar um cálculo não é um artifício destinado a melhorar a velocidade de execução de simples "caso de escola", mas pode servir nos casos reais de programação. É necessário, então, passar diretamente as sub-sessões IV.2 t IV.3 sem que a compreensão da dica seja alterada.

. Exemplo concreto escolhido


O programa utilizado para os testes calcula o campo magnético gerado no espaço pelo loop de corrente.

O campo magnético gerado por um tal loop é dado pela lei de Biot e Savart:

Programa não vetorização


function B = compute_B_loop(R,I,x,n)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%% Calcule o campo magnético gerado por um loop de corrente 
%% R é o raio do loop 
%% I é a intensidade da corrente 
%% x é a matriz dos pontos que se quer calcular o campo B  
%%  o tamanho de  x é  3 x número_de_pontos 
%%     ( x1 ... ...)  
%% x = ( y1 ... ...)  
%%     ( z1 ... ...)  
%% n é o número de pontos de descreditação para o loop corrente 
%% B é uma  matriz 3 x número_de_pontos  
%%     ( Bx1 ... ...)  
%% B = ( By1 ... ...)  
%%     ( Bz1 ... ...)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
mu0 = 4*pi*10^(-7);  

nx=size(x,2);  

B=zeros(3,nx);  

for p=1:nx  
    for q=0:n-1  
        normxmy3=((x(1,p)-R*cos(2*pi/n*q))^2+(x(2,p)-R*sin(2*pi/n*q))^2+x(3,p)^2)^(3/2);  
        B(1,p)=B(1,p)+cos(2*pi/n*q)*x(3,p)/normxmy3;  
        B(2,p)=B(2,p)+sin(2*pi/n*q)*x(3,p)/normxmy3;  
        B(3,p)=B(3,p)...  
            +(-sin(2*pi/n*q)*(x(2,p)-R*sin(2*pi/n*q))-...  
            cos(2*pi/n*q)*(x(1,p)-R*cos(2*pi/n*q)))/normxmy3;  
    end  
end  


B = mu0*R*I/(2*n)*B;  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

Programa de vetorização


function B = compute_B(R,I,x,n)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%% Calcule o campo magnético gerado por um loop  corrente 
%% R é o  raio do loop
%% I é a intensidade da corrente 
%% x é a  matriz dos pontos onde se quer calcular o campo B  
%% O tamanho x é  3 x número_de_pontos  
%%     ( x1 ... ...)  
%% x = ( y1 ... ...)  
%%     ( z1 ... ...)  
%% n é o  número de pontos de descretização para o loop corrente  
%% B é uma matriz  3 x número_de_pontos  
%%     ( Bx1 ... ...)  
%% B = ( By1 ... ...)  
%%     ( Bz1 ... ...)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
mu0 = 4*pi*10^(-7);  

nx=size(x,2);  

%% x é a matriz dos pontos do espaço, a gente a repete 
x=repmat(reshape(x.',1,nx,3),[n 1 1]);  

%% y é a matriz dos pontos do loop  
y=R*[cos(2*pi/n*(0:n-1));sin(2*pi/n*(0:n-1));zeros(1,n)];  

%% Repete-se a matriz  y  
y=repmat(reshape(y.',n,1,3),[1 nx 1]);  

%% Pode-se então calcular matricialmente todos os vetores SM (i.e. yx)  
xmy=x-y; clear x y;  

%% Calcula-se todos os vetores dl da formula (multiplica-se por R no fim)  
dy=[-sin(2*pi/n*(0:n-1));cos(2*pi/n*(0:n-1));zeros(1,n)];  

%% repete-se a matriz
dy=repmat(reshape(dy.',n,1,3),[1 nx 1]);  

%% Pode-se então calcular matricialmente todos os vetores  dl^SM da formula (i.e. dy^yx)  
cdyxmy=cross(dy,xmy,3); clear dy;  

%%Calcula-se na passagem todo as as normas de SM  (yx)  
norm_xmy_3=repmat(dot(xmy,xmy,3).^(3/2),[1 1 3]); clear xmy;  

%% Só falta agora a juntar B 
B=cdyxmy./norm_xmy_3; clear cdyxmy norm_xmy_3;  
B=sum(B);  
B = mu0*R*I/(2*n)*B;  
B = squeeze(B);  
B = B.';  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

Tempos de execução para cálculos importantes


Executemos agora estes dois programas para os pontos de [-5,5]x[-5,5]x[-5,5] na grade regular 50x50x50 com 50 pontos de descreditação no loop de corrente
A execução do programa não vetorizado
[X Y Z]=meshgrid(linspace(-5,5,50));  
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];  
clear X Y Z;  
tic;B_loop=compute_B_loop(1,1,x,50);toc;  
clear B_loop;
nos dá um tempo de execução de 20,8s.

A execução do programa de vetorização
[X Y Z]=meshgrid(linspace(-5,5,50));  
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];  
clear X Y Z;  
tic;B=compute_B(1,1,x,50);toc;  
clear B;
nos dá um tempo de execução de 14,4s.

O programa não vetorização toma então perto de 1,5 vezes mais de tempo que o cálculo de vetorização .

Imagine então o resultado sobre um programa que roda em algumas horas. É sempre mais agradável ter os resultados depois de somente 3h de espera ao invés de 4h30. Ainda mais se os resultados são falsos e que é preciso relançar depois da correção. :-)

Desvantagens de vetorização

  • Como podemos ver com o compute_B programa mesmo comentado , um cálculo de vetorização é difícilmente legível para alguém que gostaria de usar, modificar, complementar ou traduzi-lo em outro idioma (C, Fortran, etc .. .) para inclusão em um código maior.
  • A segunda desvantagem é, evidentemente , a memória. De fato, um programa vectorizado às vezes necessita armazenar Matrizes maiores do que aquelas necessária para armazenar o resultado final, é em particular verdadeiro para o nosso exemplo. Tente por exemplo lançar o seguinte cálculo (peguemos agora um a grade regular 100x100x100).


[X Y Z]=meshgrid(linspace(-5,5,100));  
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];  
clear X Y Z;  
tic;B=compute_B(1,1,x,50);toc;  
clear B;
O lugar memória é insuficiente para a execução do programa.

Enquanto que o programa não vetoriza
[X Y Z]=meshgrid(linspace(-5,5,100));  
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];  
clear X Y Z;  
tic;B_loop=compute_B_loop(1,1,x,50);toc;  
clear B_loop;
s'exécute sans problème.

Soluções aos problemas precedentes

  • Para o problema de legibilidade, é claro que o comentário ao máximo de seu códigovetorização (mesmo um código vectorizado não é claro) pode ser uma ajuda valiosa para qualquer usuário futuro.
  • Introduzir o segundo problema é um pouco mais complicado. Trata-se, geralmente, de manter a vectorialisação para fazer passar os cálculos "em bloco", mas para chegar ao segmento de blocos de tamanho adequado para não saturar a memória. Aqui você pode ver ver proceder para o nosso exemplo. Aqui o maior número de pontos (aquele que vai colocar o problema) é o que corresponde à descreditação do volume. Nós vamos segmentar o bloco de Pontos enviados para o cálculo do campo magnético.


function B = compute_B(R,I,x,n)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%% Calcule le campo magnético gerado pelo loop de corrente 
%% R est le raio do loop 
%% I é a intensidade da corrente 
%% x é a matriz dos pontos onde queremos calcular o campo B 
%% O tamanho de X é 3 X número _de_pontos  
%%     ( x1 ... ...)  
%% x = ( y1 ... ...)  
%%     ( z1 ... ...)  
%% n é o  número de pontos de descreditação para o loop da corrente  
%% B é uma  matriz 3 x número_de_pontos  
%%     ( Bx1 ... ...)  
%% B = ( By1 ... ...)  
%%     ( Bz1 ... ...)  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
mu0 = 4*pi*10^(-7);  

nx=size(x,2);  
%% O tamanho  l de um bloco será de 32768 pontos  
l = 32768;  
%% Os cálculos que seguem serão idênticos ao código precedente. 
%% salvo que se trabalha por blocos 
seg=(1:l:nx);  
if seg(end)~=nx, seg(end+1)=nx; end  
seg(end)=seg(end)+1;  
nseg=size(seg,2)-1;  
display([num2str(nseg), ' blocs']);  
B = zeros(nx,3);  
for pq=1:nseg  
    display(['Bloc ',num2str(pq)]);  
    npq=seg(pq+1)-seg(pq);  
    xpq=repmat(reshape(x(:,seg(pq):seg(pq+1)-1).',1,npq,3),[n 1 1]);  

    y=R*[cos(2*pi/n*(0:n-1));sin(2*pi/n*(0:n-1));zeros(1,n)];  
    y=repmat(reshape(y.',n,1,3),[1 npq 1]);  
    xpqmy=xpq-y; clear xpq y;  
    dy=[-sin(2*pi/n*(0:n-1));cos(2*pi/n*(0:n-1));zeros(1,n)];  
    dy=repmat(reshape(dy.',n,1,3),[1 npq 1]);  
    cdyxpqmy=cross(dy,xpqmy,3); clear dy;  
    norm_xpqmy_3=repmat(dot(xpqmy,xpqmy,3).^(3/2),[1 1 3]); clear xmy;  

    B(seg(pq):seg(pq+1)-1,:)=squeeze(sum(cdyxpqmy./norm_xpqmy_3)); clear cdyxpqmy norm_xmy_3;  
end  
B = mu0*R*I/(2*n)*B;  
B = B.';  
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Lancemos agora o cálculo precedente que não se executará, pois o espaço memória é insuficiente.
[X Y Z]=meshgrid(linspace(-5,5,100));  
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];  
clear X Y Z;  
tic;B=compute_B(1,1,x,50);toc;  
clear B;
É preciso mais ou menos 124s para se executar.

O programa não vetorizado
[X Y Z]=meshgrid(linspace(-5,5,100));  
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];  
clear X Y Z;  
tic;B_loop=compute_B_loop(1,1,x,50);toc;  
clear B_loop;
s'exécute quant à lui en 169s.

Aqui ainda, é preciso ao programa não vetorizado mais ou menos 1,5 vezes mais de tempo que ao programa « semi-vetorização » para se executar

Tradução feita por Ana Spadari

Veja também :
Este documento, intitulado «  A vetorização no MatLab »a partir de CCM (br.ccm.net) está disponibilizado sob a licença Creative Commons. Você pode copiar, modificar cópias desta página, nas condições estipuladas pela licença, como esta nota aparece claramente.