Modifier une boucle <<for>>

**Mataka** · 26/05/2007, 01h20

Bonjour,

J'ai un programme qui fait beaucoup de minimisation de fonction, et il est donc assez lent. Je cherche une facon de l'optimiser, et pour cela j'aurais grandement avantage a changer ce genre de boucle :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
for j=2:round(ll(1)/dl)
    sigmaT(:,:,j)=drift'*sigmaT(:,:,j-1)*drift;
end

sous forme purement matriciel, mais je sais pas trop comment m'y prendre. J'avais penser a quelque chose du genre :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

sigmaT(:,:,2:round(ll(1)/dl))=drift'*sigmaT(:,:,2:round(ll(1)/dl)-1)*drift;

Mais bien evidemment ca ne marche pas parce que j'utilise la multiplication matricielle, et non element par element.

Avez vous des idees ?

**ol9245** · 26/05/2007, 16h04

peux-tu nous décrire sigmaT et drift ?
OL

**rostomus** · 26/05/2007, 17h34

Envoyé par ol9245

peux-tu nous décrire sigmaT et drift ?
OL

Il me semple que sigmaT(:,:,i) et drift sont deux matrices carrées de même taille.

**Mataka** · 26/05/2007, 18h38

Drift est une matrice carré 6x6.

SigmaT(:,:,i) est aussi une matrice carrée 6x6. J'ai ajouté un dimension à SigmaT pour stocker l'information des différentes multiplications.

**Jerome Briot** · 26/05/2007, 18h43

As-tu essayé la préallocation de mémoire sur ta matrice SigmaT ?

A lire dans la

: Pré-allocation mémoire des matrices

Mais en y réfléchissant, vu la taille de la matrice (6x6), je ne pense pas que cela aura beaucoup d'effet sur la rapidité du code

**Mataka** · 26/05/2007, 20h24

Oui toutes mes matrices sont déjà préaloué.

Ça fait un grosse diférence parceque même si ma matrice sigmatT(:,:,i) est une matrice 6x6, ma matrice sigmaT(:,:,: ) est une matrice 6x6x100000.

**ol9245** · 27/05/2007, 03h31

bjr,
je ne vois pas de solution radicale a ton prob car la multiplication matricielle est definie uniquement sur des matrices a 2 dimensions. donc de toutes façons tu va etre obligé d'itérer shift' * sigmaT * shift le long de la troisieme dimension de ton grand tableau.

Il y a quand meme une solution éventuelle : tu veux créer une multiplication matricielle A(:,:,: ) * B(:,:,: ) = C(:,:,: )
tu définis ta multiplication ainsi :
A(:,:,i) * B(:,:,i) = C(:,:,i)
cette définition est équivalente à la suivante, qui itère sur les dimensions 1&2 au mieu d'itérer comme tu le fais sur la dimension 3 :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
function C = mult3(A,B)
s = size(A) ;
if s(1) ~= s(2), error('A(:,:,i) doit etre carree') ; end
if s ~= size(B), error('A & B doivent avoir la meme taille') ; end 
C = zeros(s) ;
for i=1:s(1), for j=1:s(2)
     C(i,j,:) = sum(squeeze(A(i,:,:)) .* squeeze(B(:,j,:))) ;
end, end
return

dans cette fonction, je fais la multiplication matricielle 6x6 à la main, comme à l'école, dans deux boucles for imbriquées. Mais au lieu de manipuler des nombres, je manipule tes 10000 lignes d'un seul coup. donc en toute théorie ça devrait aller plus vite.
Pour utiliser cette fonction, tu dois fabriquer (par la fonction repmat) des matrices de même taille que sigmaT à partir de shift et shift' en les répétant autant de fois que nécessaire.

OL

**Mataka** · 27/05/2007, 04h33

C'est une bonne idée, même très bonne, mais j'y vois deux gros problèmes.

Premier problème, mes matrices DRIFT sont des 6x6, pas des 6x6x100000. C'est pas le plus gros problème parce que au pire je peux les extensionné sur la troisème dimension.

Deuxième problème, le gros problème c'est que ma boucle est structuré de sorte que l'élément sigmaT(:,:,j) est fonction de l'élément sigmaT(:,:,j-1) , et donc dans ce cas puis que matlab exécute toute la multiplication en même temps il ne sait pas ce qu'est devenue l'élément précédent.

Bref comme vous dite je pense pas qu'il y ait une solution radicale. Par contre, mes matrices DRIFT possède une propriété intéressantes que j'essaie d'exploité, c'est-à-dire qu'elles sont linéaire dans le sens que

drift(dl)*drift(dl)=drift(2dl)

où dl est un élément crucial dont dépend drift. Et donc, j'aurais aussi pu faire ma boucle comme :

sigmaT(:,:,j) = drift(j*dl)'*sigmaT(:,:,1)*drift(j*dl);

où drift=@(dl)[drift];

Bon, sous cette forme il y a tu une façon d'enlevé la boucle ?

**ol9245** · 27/05/2007, 10h43

Envoyé par Mataka

C'est une bonne idée, même très bonne, mais j'y vois deux gros problèmes.

Premier problème, mes matrices DRIFT sont des 6x6, pas des 6x6x100000. C'est pas le plus gros problème parce que au pire je peux les extensionné sur la troisème dimension.

Oui, c'est bien ce que je comptais que tu fasse.

Deuxième problème, le gros problème c'est que ma boucle est structuré de sorte que l'élément sigmaT(:,:,j) est fonction de l'élément sigmaT(:,:,j-1) , et donc dans ce cas puis que matlab exécute toute la multiplication en même temps il ne sait pas ce qu'est devenue l'élément précédent.

est-ce que ça veut dire que ton problème est entièrement défini par récurrence ?
dans ce cas je ne comprend plus ton algo car si sigmaT(,,3) dépend de sigmatT(,,2) qui dépend de sigmat(,1) alors il vaut mieux dérouler toutes tes multiplications d'un coup parceque là, tu va simplifier grandement ton problème :

T(3) = D * T(2) * D'
T(2) = D * T(1) * D'

donc
T(3) = D² * T(1) * D'²

et T(1000000) = D^100000 * T(1) * D' ^100000

Comme tu as une manière simple de calculer la puissance nième de Drift, ça peut peut-être te servir à utiliser la fonction que je t'ai donnée : tu calcules une matrice D1(6x6x100000) qui contient à la suite les 100000 première puissances de Drift. Même chose pour D2 avec Drift'. Puis tu repmat 100000 fois la matrice SigmaT(1). Tu as ainsi construit 3 matrices (6x6x100000) que tu multiples entre elles avec la fonction que je t'ai donnée.
OL

**Mataka** · 27/05/2007, 17h11

wow merci ! Je l'ai pas essayé encore mais
ça me semble être une excellente idée !

Mais est-ce que le temps de calcul que je sauve sur ma boucle de récurrence je vais le perdre sur ma boucle d'extension dans la troisième dimension ?

Parce que je vois pas d'autre façon d'extensionné que :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
drift2=zeros(6,6,10000);
for j=1:100000
drift2(:,:,j)=drift^j;
end

**Mataka** · 27/05/2007, 19h19

Hum je viens de penser à ça pis mettons que j'ai pas vrament le choix de faire une boucle pour extensionner drift dans la troisième dimension, du moins je pense.

Mais pour extensionner sigmaT(:,:,1) dans la troisème dimension, où toutes les dimensions sont pareils il doit il y avoir une façon de faire ça sans faire de boucle ... non ?

Il y aurait pas une commande qui dupliquerait un matrice sur une autre dimension déjà préfait ?

**Mataka** · 27/05/2007, 20h32

Ouin je viens de le tester, avec une extension en 3D avec des boucles et c'est pas super. À vrai dire c'est encore plus long qu'avant ...

**Mataka** · 27/05/2007, 22h27

Ok je l'ai essayer avec la fonction ''repmat'' pis je gagne pas vraiment de vitesse, je dirais même que c'est un peu plus lent qu'avant.

Bref c'était une belle approche mais je pense pas que ça augmente vraiment la rapidité du calcul. Je dois avouer que je ne sais pas trop pourquoi par contre ... le calcul devrait être plus vite mais en pratique il ne l'est pas.

Je crois qu'un autre approche devrait êre envisager, mais je sais pas trop laquelle. Je vais essayer de penser à ça et je vous reviens la desssus...

**ol9245** · 28/05/2007, 13h06

Je ne connais pas ton problème, mais quand même, ça m'interpelle. En lisant entre les lignes, on dirat que t représente le temps et que tes 100000 matrices représentent 100000 pas de temps intermédiaires à calculer. J'ai juste ?

si oui, je continue : sachant que pour tout t, sigmaT(t) = Drift^t * sigmaT(1) * drift' ^t, sachant en plus que drift(x)^k = drift(k*x), donc calculable d'un coup d'un seul, pourquoi tu as encore besoin de calculer des pas de temps intermédiaires ?

Je m'explique. Tu dis que tu fais de la minimisation. Je comprend que tu ajustes des valeurs calculées et observées. tes observations sont à des dates t. Donc pour caler ton modèle, il te suffit de calculer ces dates là par la formule directe. Dans la phase de calage, tu n'a aucun usage des autres.

Tout cela n'est que conjectures basées sur ce que mon petit doigt me dit de ton problème...

OL

**Mataka** · 28/05/2007, 16h43

En réalité, t représente l'index de l'axe des Z, c'est un axe d'espace et non de temps. sigmaT représente la matrice qui caractérise un faisceau de particules ; mon programme sert à faire du traitement d'optique de faisceau de particules pour un accélérateur linéaire

Mais tout ça n'est pas vraiment important.

Quand je dit que je fais de la minimisation, c'est que mon programme doit rouler plusieurs fois, disons 1000 fois pour trouver un minimum. Donc ma boucle problématique doit être très rapide.

Le problème c'est que même en utilisant l'algorithme que vous m'avez conseillé, et que je croyais moi même comme étant plus performant, et bien au contraire il est plus lent. La boucle de proche en proche s'avère être plus rapide. Faut-dire que cette dernière manipule des matrice beaucoup moins lourdes, c'est peut-être ça le problème. En fait, j'ai fait des tests et là ou votre algo est plus lent que ma boucle, c'est pour la multiplication matricielle, dans mult3.

**ol9245** · 29/05/2007, 04h34

merci e ces détails très intéressants.

Pour le code, je veux bien regarder ton nouveau code si tu le postes sur le forum.

Autre question qui n'a rien à voir avec Matlab : pourquoi discrétiser la profondeur en autant d'éléments ? Tu as essayé de diviser disons par 10 ton pas d'espace ? ça donne quoi ?

En modélisation, il y a une règle d'or absolue à laquelle je m réfère pour te poser cette question :
A/ un bon modèle est toujours bon. Même avec une discrétisation grossière de l'espace, il doit continuer à donner de bons résultats (mais imprécis).
B/ un mauvais modèle est toujours mauvais. Diminuer les pas d'espace et de temps n'améliore jamais un mauvais modèle.

En application de cette règle, il est toujours souhaitable de mettre au point et de tester les modèls avec des mailles de temps et d'espace grossières. Tant qu'ils ne donne pas satisfaction dans ce cadre là, il n'est pas encore temps de les faire tourner à pas plus fin.
OL

**Mataka** · 29/05/2007, 06h31

Je dois discrétisé en autant d'élément pour une raison bien simple; en fait mon programme est un clône d'un autre programme bien connu parmis les physiciens d'accélérateur, soit Trace3D. Et pour arriver aux même résultats que Trace3D je dois discrétisé en beaucoup d'éléments (en fait trace3D doit être écrit en C++ je crois et donc c'est normal qu'il soit plus vite). De fait, même avec une résolution de 100 000 j'ai encore une petite erreur comparativement à trace3d, mais à mesure que j'augmente la résolution je diminue l'erreur. Je m'arrête à 100 000 car l'erreur entre trace3d et mon programme est suportable et la vitesse de mon programme est potable, sans être très bonne.

Le problème c'est que trace3d n'est pas en code ouvert, et je n'ai donc pas accès à leurs astuces de programmation

Pour l'instant, voici mon code de ma fonction de calcul d'une des lignes de l'accélérateur (avec nelm=100 000). Vous y reconnaîtrez mes boucles très longue d'exécution :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
function  [sigmaTcal]=calcu(I)
 
global Aq A energy beta gamma nelm numq c mu q z dl lengthQ f ex ey ez Xini ll firsttimecal m zer sigmaTcal
 
if firsttimecal == 0
 
    %  Line construction
    SP=0.01*[65 19 84.5 19 19 69];   % length between the quadpole
    ll=zeros(1,11);
    ll(1)=SP(1);
    for j=2:length(ll)
        if (-1)^j==1
            ll(j)=ll(j-1)+lengthQ(1);
        else
            ll(j)=ll(j-1)+SP((j+1)/2);
        end
    end
    L=ll(length(ll));   % length of the line
    dl=L/nelm; %infinitesimal element
    z=linspace(0,L,nelm); % z axis
 
    % GRADIENT CONVERSION
    m=zeros(1,4);
    m(1) = 0;
    m(2) = (1/(10.8))/2.6;
    m(3) = (-2.699029619469e-5)/2.6;
    m(4) = (-3.81181036994644e-7)/2.6;
 
    zer = zeros(2,2);
    sigmaTcal = zeros(6,6,nelm);
 
    firsttimecal=1;
 
end
 
gradient=zeros(1,numq);
qdp=zeros(6,6,numq);
 
% Calculate the sigma matrix :
ext = ex/1e6;
eyt = ey/1e6;
ezt =(beta*c/f)/(gamma*(gamma+1)*A*energy*1e3*360)*ez;
BB = (beta*c/f)*gamma*(gamma+1)/360*A*energy*1e3*Xini(6);
 
sigmax = ext*[Xini(2) -Xini(1); -Xini(1) ((1+Xini(1)^2)/Xini(2))];
sigmay = eyt*[Xini(4) -Xini(3); -Xini(3) ((1+Xini(3)^2)/Xini(4))];
sigmaz = ezt*[BB Xini(5); Xini(5) (1+Xini(5)^2)/BB];
 
sigma = [ sigmax zer zer; zer sigmay zer; zer zer sigmaz];
 
%%%%%%%%%%%%%%%%   Calculate the matrix for the transformastion of the beam %%%%%%%%%%%%%%%%%%%%%%%
 
%  Free field drift matrix
drift=[1 dl 0 0 0 0; 0 1 0 0 0 0; 0 0 1 dl 0 0; 0 0 0 1 0 0; 0 0 0 0 1 (dl/gamma^2); 0 0 0 0 0 1]';
 
for j=1:numq
    gradient(j) = (m(1) + m(2)*I(j) + m(3)*I(j)^2 + m(4)*I(j)^3)*10;
end
 
% Polarity conversion :
gradient = [-gradient(1) gradient(2) gradient(3) -gradient(4) gradient(5)];
% end of polarity conversion.
 
kq = sqrt(gradient*q/(Aq*c*beta*gamma*mu));
 
%%%%%%%%%%%%%%%%%%%% ZERO PROTECTION %%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:numq
    if abs(kq(j))==0
        kq(j)=eps;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Quadrupole matrixs
for j=1:numq
    qdp(:,:,j)=[cos(kq(j)*dl) 1/kq(j)*sin(kq(j)*dl) 0 0 0 0; -kq(j)*sin(kq(j)*dl) cos(kq(j)*dl) 0 0 0 0; 0 0 cosh(kq(j)*dl) 1/kq(j)*sinh(kq(j)*dl) 0 0; 0 0 kq(j)*sinh(kq(j)*dl) cosh(kq(j)*dl) 0 0; 0 0 0 0 1 dl/(gamma^2); 0 0 0 0 0 1]';
end
 
%%%%%%% Processing to the transformation for extent and for divergence %%%%%%
 
sigmaTcal(:,:,1) = sigma;
 
% DRIFT TRANSFORMATION
 
for j=2:round(ll(1)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%  QUADPOLE TRANSFORMATION
 
for j=round(ll(1)/dl)+1:round(ll(2)/dl)
    sigmaTcal(:,:,j)=qdp(:,:,1)'*sigmaTcal(:,:,j-1)*qdp(:,:,1);
end
 
% DRIFT TRANSFORMATION
 
for j=round(ll(2)/dl)+1:round(ll(3)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%  QUADPOLE TRANSFORMATION
 
for j=round(ll(3)/dl)+1:round(ll(4)/dl)
    sigmaTcal(:,:,j)=qdp(:,:,2)'*sigmaTcal(:,:,j-1)*qdp(:,:,2);
end
 
% DRIFT TRANSFORMATION
 
for j=round(ll(4)/dl)+1:round(ll(5)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%  QUADPOLE TRANSFORMATION
 
for j=round(ll(5)/dl)+1:round(ll(6)/dl)
    sigmaTcal(:,:,j)=qdp(:,:,3)'*sigmaTcal(:,:,j-1)*qdp(:,:,3);
end
 
% DRIFT TRANSFORMATION
 
for j=round(ll(6)/dl)+1:round(ll(7)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%  QUADPOLE TRANSFORMATION
 
for j=round(ll(7)/dl)+1:round(ll(8)/dl)
    sigmaTcal(:,:,j)=qdp(:,:,4)'*sigmaTcal(:,:,j-1)*qdp(:,:,4);
end
 
% DRIFT TRANSFORMATION
 
for j=round(ll(8)/dl)+1:round(ll(9)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%  QUADPOLE TRANSFORMATION
 
for j=round(ll(9)/dl)+1:round(ll(10)/dl)
    sigmaTcal(:,:,j)=qdp(:,:,5)'*sigmaTcal(:,:,j-1)*qdp(:,:,5);
end
 
% DRIFT TRANSFORMATION
 
for j=round(ll(10)/dl)+1:round(ll(11)/dl)
    sigmaTcal(:,:,j)=drift'*sigmaTcal(:,:,j-1)*drift;
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% END OF TRANSFORMATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

**ol9245** · 29/05/2007, 08h40

A/ tu as une solution analytique pour calculer directement drift^n. Tu êux même calculer d'un seul coup avec une seule formule matlab toutes les puissances niemes de sigmaT que tu veux.

B/ est-ce que tu utilises effectivement chacun des pas de temps intermédiaires dans ton optimisation ? j'ai l'impression que non. Seul le résultat finala compte, n'est-ce pas ? Sachant que analytiquement, tu as une solution équivalente à XX pas de temps en une seule opération, moyenant le calcul de drift^n, donc tu peux enchainer la succession e tes quadpole et drift avezc une ligne matlab à chaque fois, et un résultat analytiquement équivalent non ?

OL

**Mataka** · 29/05/2007, 15h55

Hum ... c'est une bonne idée. J'ai normalement besoin du résultat intermédiaire mais pour l'optimisation le résultat final suffit.

Je vais tester ça aujourd'hui, mais je suis assez confiant que ça marche !

Merci !