Parallelize matlab algorithm

Question

I'm working on a tomographic algorithm, using SART method. I'm using 2D projections in order to get 3D volume.

My main problem is that computation time is really high..

For example, with projections of size 88*75 pixels my program need 5 days to complete.

By using the matlab profiler, we can see that the lines with the multiplication of the radon matrix (which is heavy) are taking so much time.

That's why I want to parallelised my algorithm but I'm not really familiar with it.

This is the part of my algorithme that I want to parallelize :

line which need a lot of time: Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED');


for k = 1: ite % Boucle définie par le nombre d 'itérations choisi par l'
utilisateur
Soustraction = Projection2 - Projection_old;
Erreur(ligne, k) = mean(Soustraction);
Err = Erreur(ligne, k);
for z = 1: L * Nproj
for a = 1: L * L
Somme = Somme + Matrice_Radon(z, a) * Matrice_Radon(z, a);
end
facteur = landa / Somme;
Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED'); % Formule de correction additive de l
  'image

            Correction = Correction*Soustraction;
            Ik1 = Ik + Correction;
            Ik=Ik1;                                                             % Stockage de l'
image corrigée dans l
  'image_itération_précédente

             Projection_old=Matrice_Radon*Ik;                          % Calcul de la projection à partir de l'
image_itération_précédente
end
pos3 = 1;
for incr2 = 1: L
for incr3 = 1: L
Ik2(incr2, incr3) = Ik1(pos3, 1); % Transformation vecteur en matrice
pos3 = pos3 + 1;
end
end
Coupe_SIRT(ligne, : , : ) = Ik2; % Stockage des coupes générées
end
fprintf('Progression : %2d/%d\n', ligne, H); % Affichage de l
  'avancement du traitement
 end

please format your code to a more readable format: indentation and comments in **English**. — Shai, Jul 07 '14 at 14:30
And this is why people (ASTRA,TIGRE,RTK,....) build GPU tomographic reconstruction algorithms! They are massively slow otherwise, no CPU can do the job fast enough. — Ander Biguri, Feb 16 '17 at 17:02

Eli Duenisch · Answer 1 · 2014-07-07T17:19:20.037

I see a lot of nested for loops here. You should use vectorized code instead of for loops whenever it is possible in MATLAB in order to speed up your algorithm. Also it is not possible to use nested parfor.

Parallelization is possible whenever your loops are independent from each other (and of course you need to have the Parallel Computing Toolbox installed). In this case the use it very esay:

  result = zeros(1,n); % preallocate result storing array
  parfor i=1:n % parallel for loop
    %your code here, e.g.:
    result(i) = ...;
  end

A pool of MATALB workers can be set up manually, but it is also possible to just use parfor and let MATLAB handle the pool automatically (at least in r2014a).

Please note that parallelization does not always speed things up.

user3812144 · Answer 2 · 2014-07-08T08:41:59.037

I really don't know how to vectorize this code..

I don't need to vectorize all the code, just this 2 lines :

Correction = mtimesx(facteur, Matrice_Radon_Transposee, 'SPEED'); 

Correction = Correction*Soustraction;

So maybe it would be interesting to parallelize the for loop ? :

for k = 1: ite

Is it correct to write this :

result_Projection_old = zeros(1,ite);
result_Ik = zeros(1,ite);

parfor  k = 1 : ite                

    Somme = 0;      
    Projection_old = result_Projection_old(k);
    Ik = result_Ik(k);

    Soustraction = Projection2 - Projection_old;


    for z=1:LNproj
        for a=1:L_carre

            Somme = Somme + Matrice_Radon(z,a)*Matrice_Radon(z,a);

        end

        facteur = landa/Somme;
        Correction = mtimesx(facteur,Matrice_Radon_Transposee,'SPEED');    
        Correction = Correction*Soustraction;
        Ik1 = Ik + Correction;
        Ik=Ik1; 
        result_Ik(k) = Ik;                                                          

        Projection_old=Matrice_Radon*Ik;
        result_Projection_old(k) = Projection_old;                     

    end
end

With this code, the reconstruction doesn't work, instead of writing result_Ik(k) = Ik I need to write Result_Ik(k+1) = Ik in order to get at the next iteration the value of Ik from the itération k-1 (previous iteration)

Parallelize matlab algorithm

2 Answers2